]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TUHKMgen/UHKM/HankelFunction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Do not include from subdirectories
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / HankelFunction.cxx
CommitLineData
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
14enum {kNe = 2, kNCoeff = 9};
15
16const 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
30const 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
44const 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
58const 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
72Double_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
109Double_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
147Double_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
179Double_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
212Double_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