#ifndef RTS_SBESSEL_H
#define RTS_SBESSEL_H
#include <math.h>


namespace rts{

#define RTS_BESSEL_CONVERGENCE_FLOAT	0.00045

template <typename T>
CUDA_CALLABLE void sbesselj(int n, rtcomplex<T> x, rtcomplex<T>* j)
{
    //compute the first bessel function
    if(n >= 0)
        j[0] = sin(x) / x;

    //compute the second bessel function
    if(n >= 1)
        j[1] = j[0] / x - cos(x) / x;

    //use the recurrence relation to compute the rest
    for(int i = 2; i <= n; i++)
    {
        j[i] = ( (2 * i - 1) / x ) * j[i-1] - j[i-2];
    }

    //if x = 0, deal with the degenerate case
    /*if( isnan(j[0].r) )
    {
        j[0] = (T)1.0;
        for(int i = 1; i<=n; i++)
            j[i] = (T)0.0;
    }*/
}

template <typename T>
CUDA_CALLABLE void sbessely(int n, rtcomplex<T> x, rtcomplex<T>* y)
{
    //compute the first bessel function
    if(n >= 0)
        y[0] = -cos(x) / x;

    //compute the second bessel function
    if(n >= 1)
        y[1] = y[0] / x - sin(x) / x;

    //use the recurrence relation to compute the rest
    for(int i = 2; i <= n; i++)
    {
        y[i] = ( (2 * i - 1) / x ) * y[i-1] - y[i-2];
    }

}

//spherical Hankel functions of the first kind
template <typename T>
CUDA_CALLABLE void sbesselh1(int n, rtcomplex<T> x, rtcomplex<T>* h)
{
    //compute j_0 and j_1
    rtcomplex<T> j[2];
    sbesselj(1, x, j);

    //compute y_0 and y_1
    rtcomplex<T> y[2];
    sbessely(1, x, y);

    //compute the first-order Hhankel function
    if(n >= 0)
        h[0] = j[0] + y[0].imul();

    //compute the second bessel function
    if(n >= 1)
        h[1] = j[1] + y[1].imul();

    //use the recurrence relation to compute the rest
    for(int i = 2; i <= n; i++)
    {
        h[i] = ( (2 * i - 1) / x ) * h[i-1] - h[i-2];
    }
}

template <typename T>
CUDA_CALLABLE void init_sbesselj(T x, T* j)
{
	/*if(x < EPSILON_FLOAT)
	{
		j[0] = (T)1;
		j[1] = (T)0;
		return;
	}*/
	
	//compute the first 2 bessel functions
	j[0] = std::sin(x) / x;

	j[1] = j[0] / x - std::cos(x) / x;

	//deal with the degenerate case of x = 0
	/*if(isnan(j[0]))
	{
		j[0] = (T)1;
		j[1] = (T)0;
	}*/
	
}

template <typename T>
CUDA_CALLABLE void shift_sbesselj(int n, T x, T* j)//, T stability = 1.4)
{

	/*if(x < EPSILON_FLOAT)
	{
		j[0] = j[1];
		j[1] = (T)0;
		return;
	}*/

	T jnew;

	//compute the next (order n) Bessel function
	jnew = ((2 * n - 1) * j[1])/x - j[0];

	//if(n > stability*x)
	if(n > x)
		if(jnew < RTS_BESSEL_CONVERGENCE_FLOAT)
			jnew = (T)0;

	//deal with the degenerate case of x = 0
	//if(isnan(jnew))
	//	jnew = (T)0;

	//shift and add the new value to the array
	j[0] = j[1];
	j[1] = jnew;
}



}   //end namespace rts



#endif