//#####################################################################
// Copyright 2007, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Header Robust_Functions
//#####################################################################
#ifndef __Robust_Functions__
#define __Robust_Functions__

#include <Core/Math_Tools/sqr.h>
#include <Core/Utilities/TYPE_UTILITIES.h>
#include <cmath>
namespace PhysBAM{

using ::std::abs;
using ::std::atan2;
using ::std::sin;

template<class T>
inline T sinc(T x) // sin(x)/x
{
    static_assert(IS_FLOATING_POINT<T>::value);
    if(abs(x)<(T)1e-8) return 1;
    return sin(x)/x;
}

template<class T>
inline T one_minus_cos_x_over_x_squared(T x) // (1-cos(x))/x^2
{return (T).5*sqr(sinc((T).5*x));}

template<class T>
inline T one_minus_cos_x_over_x(T x) // (1-cos(x))/x
{return one_minus_cos_x_over_x_squared(x)*x;}

// (sin(x)-sin(y))/(x-y)
template<class T>
inline T diff_sin_over_diff(T x, T y)
{return cos(x)*sinc(y-x)-sin(x)*one_minus_cos_x_over_x(y-x);}

// (cos(x)-cos(y))/(x-y)
template<class T>
inline T diff_cos_over_diff(T x, T y)
{return -sin(x)*sinc(y-x)-cos(x)*one_minus_cos_x_over_x(y-x);}

template<class T>
inline T atan2_y_x_over_y(T y,T x) // atan2(y,x)/y
{
    static_assert(IS_FLOATING_POINT<T>::value);
    if(abs(y)<1e-8) return 1;
    return atan2(y,x)/y;
}

template<class T>
inline T log1p_x_over_x(T x) // log(x+1)/x
{
    static_assert(IS_FLOATING_POINT<T>::value);
    if(abs(x)<1e-18) return 1;
    return log1p(x)/x;
}

template<class T>
inline T diff_log_over_diff(T x,T y) // (log(x)-log(y))/(x-y)
{if(x>2*y || y>2*x)return (log(x)-log(y))/(x-y); return log1p_x_over_x((y-x)/x)/x;}

template<class T>
inline T exp_x_minus_1_over_x(T x) // (exp(x)-1)/x
{
    static_assert(IS_FLOATING_POINT<T>::value);
    if(!x) return 1;
    return expm1(x)/x;
}

template<class T>
inline T diff_exp_over_diff(T x,T y) // (exp(x)-exp(y))/(x-y)
{return exp(x)*exp_x_minus_1_over_x(y-x);}
}
#endif
