f1402849
 dmayerich
renewed commit
|
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
|
#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
|
