spherical_bessel.h
2.94 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#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