spherical_bessel.h 2.94 KB
``````#ifndef RTS_SBESSEL_H
#define RTS_SBESSEL_H
#include <math.h>

namespace rts{

#define RTS_BESSEL_CONVERGENCE_MIN		0.0145f
#define RTS_BESSEL_CONVERGENCE_MAX		0.4f
#define RTS_BESSEL_MAXIMUM_FLOAT		-1e33f

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

//compute the second bessel function
if(n >= 1)
j[1] = j[0] / x - (T)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, complex<T> x, complex<T>* y)
{
//compute the first bessel function
if(n >= 0)
y[0] = -(T)cos(x) / x;

//compute the second bessel function
if(n >= 1)
y[1] = y[0] / x - (T)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, complex<T> x, complex<T>* h)
{
//compute j_0 and j_1
complex<T> j[2];
sbesselj(1, x, j);

//compute y_0 and y_1
complex<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)
{
//compute the first 2 bessel functions
j[0] = (T)sin(x) / x;

j[1] = j[0] / x - (T)cos(x) / x;
}

template <typename T>
CUDA_CALLABLE void init_sbessely(T x, T* y)
{
//compute the first 2 bessel functions
y[0] = -(T)cos(x) / x;

y[1] = y[0] / x - (T)sin(x) / x;
}

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

T bnew;

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

//if(n > stability*x)
if(n > real(x))
if(real(bnew) < RTS_BESSEL_CONVERGENCE_MIN || real(bnew) > RTS_BESSEL_CONVERGENCE_MAX)
bnew = 0;

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

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

T bnew;

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

if(bnew < RTS_BESSEL_MAXIMUM_FLOAT ||
(n > x && bnew > 0))
{
bnew = 0;
b[1] = 0;
}

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

}   //end namespace rts

#endif
``````