]>
Commit | Line | Data |
---|---|---|
03896fc4 | 1 | ////////////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna // | |
4 | // amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru // | |
5 | // November. 2, 2005 // | |
6 | // // | |
7 | ////////////////////////////////////////////////////////////////////////////////// | |
b1c2e580 | 8 | |
9 | #include <TMath.h> | |
10 | ||
b1c2e580 | 11 | #include "HankelFunction.h" |
b1c2e580 | 12 | |
13 | //compute Hankel function of zeroth order | |
14 | enum {kNe = 2, kNCoeff = 9}; | |
15 | ||
16 | const Double_t i0Coefficient[kNCoeff][kNe] = | |
17 | { | |
18 | //coefficients to compute function I0 | |
19 | {1.0, 0.39894228}, | |
20 | {3.5156229, 0.01328592}, | |
21 | {3.0899424, 0.00225319}, | |
22 | {1.2067492, -0.00157565}, | |
23 | {0.2659732, 0.00916281}, | |
24 | {0.0360768, -0.02057706}, | |
25 | {0.0045813, 0.02635537}, | |
26 | {0., -0.01647633}, | |
27 | {0., 0.00392377} | |
28 | }; | |
29 | ||
30 | const Double_t i1Coefficient[kNCoeff][kNe] = | |
31 | { | |
32 | //coefficients to compute function I1 | |
33 | {0.5, 0.39894228}, | |
34 | {0.87890594, -0.03988024}, | |
35 | {0.51498869, -0.00362018}, | |
36 | {0.15084934, 0.00163801}, | |
37 | {0.02658733, -0.01031555}, | |
38 | {0.00301532, 0.02282967}, | |
39 | {0.00032411, -0.02895312}, | |
40 | {0., 0.01787654}, | |
41 | {0., -0.00420059} | |
42 | }; | |
43 | ||
44 | const Double_t k0Coefficient[kNCoeff][kNe] = | |
45 | { | |
46 | //coefficients to compute modified Hankel function of the zeroth order K0 | |
47 | {-0.57721566, 1.25331414}, | |
48 | {0.42278420, -0.07832358}, | |
49 | {0.23069756, 0.02189568}, | |
50 | {0.03488590, -0.01062446}, | |
51 | {0.00262698, 0.00587872}, | |
52 | {0.00010750, -0.00251540}, | |
53 | {0.0000074, 0.00053208}, | |
54 | {0., 0. }, | |
55 | {0., 0. } | |
56 | }; | |
57 | ||
58 | const Double_t k1Coefficient[kNCoeff][kNe] = | |
59 | { | |
60 | //coefficients to compute modified Hankel function of the first order K1 | |
61 | {1.0, 1.25331414}, | |
62 | {0.15443144, 0.23498619}, | |
63 | {-0.67278579, -0.03655620}, | |
64 | {-0.18156897, 0.01504268}, | |
65 | {-0.01919402, -0.00780353}, | |
66 | {-0.00110404, 0.00325614}, | |
67 | {-0.00004686, -0.00068245}, | |
68 | { 0., 0. }, | |
69 | { 0., 0. } | |
70 | }; | |
71 | ||
72 | Double_t BesselI0(Double_t x) { | |
73 | // (C) Copr. 1986-92 Numerical Recipes Software +7. | |
74 | //compute Bessel function of zeroth order | |
75 | ||
76 | const Double_t p1 = i0Coefficient[0][0]; | |
77 | const Double_t p2 = i0Coefficient[1][0]; | |
78 | const Double_t p3 = i0Coefficient[2][0]; | |
79 | const Double_t p4 = i0Coefficient[3][0]; | |
80 | const Double_t p5 = i0Coefficient[4][0]; | |
81 | const Double_t p6 = i0Coefficient[5][0]; | |
82 | const Double_t p7 = i0Coefficient[6][0]; | |
83 | ||
84 | const Double_t q1 = i0Coefficient[0][1]; | |
85 | const Double_t q2 = i0Coefficient[1][1]; | |
86 | const Double_t q3 = i0Coefficient[2][1]; | |
87 | const Double_t q4 = i0Coefficient[3][1]; | |
88 | const Double_t q5 = i0Coefficient[4][1]; | |
89 | const Double_t q6 = i0Coefficient[5][1]; | |
90 | const Double_t q7 = i0Coefficient[6][1]; | |
91 | const Double_t q8 = i0Coefficient[7][1]; | |
92 | const Double_t q9 = i0Coefficient[8][1]; | |
93 | ||
94 | Double_t i0 = 0.; | |
95 | ||
96 | if(TMath::Abs(x) < 3.75) { | |
97 | Double_t y = (x / 3.75) * (x / 3.75); | |
98 | i0 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))); | |
99 | } | |
100 | else { | |
101 | Double_t ax = TMath::Abs(x); | |
102 | Double_t y = 3.75 / ax; | |
103 | i0 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))); | |
104 | } | |
105 | ||
106 | return i0; | |
107 | } | |
108 | ||
109 | Double_t BesselI1(Double_t x) { | |
110 | // (C) Copr. 1986-92 Numerical Recipes Software +7. | |
111 | //compute Bessel function of first order | |
112 | ||
113 | const Double_t p1 = i1Coefficient[0][0]; | |
114 | const Double_t p2 = i1Coefficient[1][0]; | |
115 | const Double_t p3 = i1Coefficient[2][0]; | |
116 | const Double_t p4 = i1Coefficient[3][0]; | |
117 | const Double_t p5 = i1Coefficient[4][0]; | |
118 | const Double_t p6 = i1Coefficient[5][0]; | |
119 | const Double_t p7 = i1Coefficient[6][0]; | |
120 | ||
121 | const Double_t q1 = i1Coefficient[0][1]; | |
122 | const Double_t q2 = i1Coefficient[1][1]; | |
123 | const Double_t q3 = i1Coefficient[2][1]; | |
124 | const Double_t q4 = i1Coefficient[3][1]; | |
125 | const Double_t q5 = i1Coefficient[4][1]; | |
126 | const Double_t q6 = i1Coefficient[5][1]; | |
127 | const Double_t q7 = i1Coefficient[6][1]; | |
128 | const Double_t q8 = i1Coefficient[7][1]; | |
129 | const Double_t q9 = i1Coefficient[8][1]; | |
130 | ||
131 | Double_t i1 = 0.; | |
132 | ||
133 | if (TMath::Abs(x) < 3.75) { | |
134 | Double_t y = (x / 3.75) * (x / 3.75); | |
135 | i1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))); | |
136 | } | |
137 | else { | |
138 | Double_t ax = TMath::Abs(x); | |
139 | Double_t y = 3.75/ax; | |
140 | i1 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))); | |
141 | if(x < 0.) i1 = -i1; | |
142 | } | |
143 | ||
144 | return i1; | |
145 | } | |
146 | ||
147 | Double_t HankelK0(Double_t x) { | |
03896fc4 | 148 | // (C) Copr. 1986-92 Numerical Recipes Software +7. |
149 | // compute modified Hankel function of the first order | |
b1c2e580 | 150 | const Double_t p1 = k0Coefficient[0][0]; |
151 | const Double_t p2 = k0Coefficient[1][0]; | |
152 | const Double_t p3 = k0Coefficient[2][0]; | |
153 | const Double_t p4 = k0Coefficient[3][0]; | |
154 | const Double_t p5 = k0Coefficient[4][0]; | |
155 | const Double_t p6 = k0Coefficient[5][0]; | |
156 | const Double_t p7 = k0Coefficient[6][0]; | |
157 | ||
158 | const Double_t q1 = k0Coefficient[0][1]; | |
159 | const Double_t q2 = k0Coefficient[1][1]; | |
160 | const Double_t q3 = k0Coefficient[2][1]; | |
161 | const Double_t q4 = k0Coefficient[3][1]; | |
162 | const Double_t q5 = k0Coefficient[4][1]; | |
163 | const Double_t q6 = k0Coefficient[5][1]; | |
164 | const Double_t q7 = k0Coefficient[6][1]; | |
165 | ||
166 | Double_t k0 = 0.; | |
167 | if(x <= 2.0) { | |
168 | Double_t y = x * x / 4.0; | |
169 | k0 = (-TMath::Log(x/2.0)*BesselI0(x)) + (p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))); | |
170 | } | |
171 | else { | |
172 | Double_t y = (2.0 / x); | |
173 | k0 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7)))))); | |
174 | } | |
175 | ||
176 | return k0; | |
177 | } | |
178 | ||
179 | Double_t HankelK1(Double_t x) { | |
180 | // (C) Copr. 1986-92 Numerical Recipes Software +7. | |
181 | // compute modified Hankel function of the first order | |
182 | const Double_t p1 = k1Coefficient[0][0]; | |
183 | const Double_t p2 = k1Coefficient[1][0]; | |
184 | const Double_t p3 = k1Coefficient[2][0]; | |
185 | const Double_t p4 = k1Coefficient[3][0]; | |
186 | const Double_t p5 = k1Coefficient[4][0]; | |
187 | const Double_t p6 = k1Coefficient[5][0]; | |
188 | const Double_t p7 = k1Coefficient[6][0]; | |
189 | ||
190 | const Double_t q1 = k1Coefficient[0][1]; | |
191 | const Double_t q2 = k1Coefficient[1][1]; | |
192 | const Double_t q3 = k1Coefficient[2][1]; | |
193 | const Double_t q4 = k1Coefficient[3][1]; | |
194 | const Double_t q5 = k1Coefficient[4][1]; | |
195 | const Double_t q6 = k1Coefficient[5][1]; | |
196 | const Double_t q7 = k1Coefficient[6][1]; | |
197 | ||
198 | Double_t k1 = 0.; | |
199 | ||
200 | if(x <= 2.0) { | |
201 | Double_t y = x * x / 4.0; | |
202 | k1 = (TMath::Log(x/2.0)*BesselI1(x)) + (1.0/x)*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))); | |
203 | } | |
204 | else { | |
205 | Double_t y = 2.0 / x; | |
206 | k1 = (TMath::Exp(-x)/TMath::Sqrt(x))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*q7)))))); | |
207 | } | |
208 | ||
209 | return k1; | |
210 | } | |
211 | ||
212 | Double_t HankelKn(Int_t n, Double_t x) { | |
213 | // (C) Copr. 1986-92 Numerical Recipes Software +7. | |
214 | // compute modified Hankel function of the first order | |
215 | if(n < 2) throw "Bad argument n in Kn"; | |
216 | ||
217 | Double_t tox = 2.0 / x; | |
218 | Double_t km = HankelK0(x); | |
219 | Double_t k = HankelK1(x); | |
220 | Double_t kp = 0.; | |
221 | ||
222 | for(Int_t c = 1; c <= n-1; ++c) { | |
223 | kp = km + c * tox * k; | |
224 | km = k; | |
225 | k = kp; | |
226 | } | |
227 | ||
228 | return k; | |
229 | } | |
230 |