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