f1402849
dmayerich
renewed commit
|
1
2
3
4
5
6
7
|
#ifndef RTS_SBESSEL_H
#define RTS_SBESSEL_H
#include <math.h>
namespace rts{
|
0174d823
dmayerich
bug fixes and exits
|
8
9
|
#define RTS_BESSEL_CONVERGENCE_MIN 0.0145
#define RTS_BESSEL_CONVERGENCE_MAX 0.4
|
a47a23a9
dmayerich
added ENVI functions
|
10
|
#define RTS_BESSEL_MAXIMUM_FLOAT -1e33
|
f1402849
dmayerich
renewed commit
|
11
12
|
template <typename T>
|
7006df5f
David Mayerich
reformat of direc...
|
13
|
CUDA_CALLABLE void sbesselj(int n, complex<T> x, complex<T>* j)
|
f1402849
dmayerich
renewed commit
|
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
|
{
//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>
|
7006df5f
David Mayerich
reformat of direc...
|
39
|
CUDA_CALLABLE void sbessely(int n, complex<T> x, complex<T>* y)
|
f1402849
dmayerich
renewed commit
|
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
{
//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>
|
7006df5f
David Mayerich
reformat of direc...
|
59
|
CUDA_CALLABLE void sbesselh1(int n, complex<T> x, complex<T>* h)
|
f1402849
dmayerich
renewed commit
|
60
61
|
{
//compute j_0 and j_1
|
7006df5f
David Mayerich
reformat of direc...
|
62
|
complex<T> j[2];
|
f1402849
dmayerich
renewed commit
|
63
64
65
|
sbesselj(1, x, j);
//compute y_0 and y_1
|
7006df5f
David Mayerich
reformat of direc...
|
66
|
complex<T> y[2];
|
f1402849
dmayerich
renewed commit
|
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
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)
{
|
f1402849
dmayerich
renewed commit
|
87
|
//compute the first 2 bessel functions
|
a47a23a9
dmayerich
added ENVI functions
|
88
|
j[0] = sin(x) / x;
|
f1402849
dmayerich
renewed commit
|
89
|
|
a47a23a9
dmayerich
added ENVI functions
|
90
|
j[1] = j[0] / x - cos(x) / x;
|
f1402849
dmayerich
renewed commit
|
91
92
93
|
}
template <typename T>
|
a47a23a9
dmayerich
added ENVI functions
|
94
|
CUDA_CALLABLE void init_sbessely(T x, T* y)
|
f1402849
dmayerich
renewed commit
|
95
|
{
|
a47a23a9
dmayerich
added ENVI functions
|
96
97
|
//compute the first 2 bessel functions
y[0] = -cos(x) / x;
|
f1402849
dmayerich
renewed commit
|
98
|
|
a47a23a9
dmayerich
added ENVI functions
|
99
100
|
y[1] = y[0] / x - sin(x) / x;
}
|
f1402849
dmayerich
renewed commit
|
101
|
|
a47a23a9
dmayerich
added ENVI functions
|
102
103
104
105
106
|
template <typename T>
CUDA_CALLABLE void shift_sbesselj(int n, T x, T* b)//, T stability = 1.4)
{
T bnew;
|
f1402849
dmayerich
renewed commit
|
107
108
|
//compute the next (order n) Bessel function
|
a47a23a9
dmayerich
added ENVI functions
|
109
|
bnew = ((2 * n - 1) * b[1])/x - b[0];
|
f1402849
dmayerich
renewed commit
|
110
111
|
//if(n > stability*x)
|
a47a23a9
dmayerich
added ENVI functions
|
112
|
if(n > real(x))
|
0174d823
dmayerich
bug fixes and exits
|
113
|
if(real(bnew) < RTS_BESSEL_CONVERGENCE_MIN || real(bnew) > RTS_BESSEL_CONVERGENCE_MAX)
|
7731cf24
David Mayerich
fixed precision p...
|
114
|
bnew = 0;
|
a47a23a9
dmayerich
added ENVI functions
|
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
//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))
{
|
7731cf24
David Mayerich
fixed precision p...
|
133
134
|
bnew = 0;
b[1] = 0;
|
a47a23a9
dmayerich
added ENVI functions
|
135
|
}
|
f1402849
dmayerich
renewed commit
|
136
|
|
f1402849
dmayerich
renewed commit
|
137
138
|
//shift and add the new value to the array
|
a47a23a9
dmayerich
added ENVI functions
|
139
140
|
b[0] = b[1];
b[1] = bnew;
|
f1402849
dmayerich
renewed commit
|
141
142
143
144
145
146
147
148
149
|
}
} //end namespace rts
#endif
|