New generator: TUHKMgen
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / HankelFunction.cxx
CommitLineData
b1c2e580 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*/
8
9#include <TMath.h>
10
11#ifndef HANKELFUNCTION_INCLUDED
12#include "HankelFunction.h"
13#endif
14
15//compute Hankel function of zeroth order
16enum {kNe = 2, kNCoeff = 9};
17
18const Double_t i0Coefficient[kNCoeff][kNe] =
19 {
20 //coefficients to compute function I0
21 {1.0, 0.39894228},
22 {3.5156229, 0.01328592},
23 {3.0899424, 0.00225319},
24 {1.2067492, -0.00157565},
25 {0.2659732, 0.00916281},
26 {0.0360768, -0.02057706},
27 {0.0045813, 0.02635537},
28 {0., -0.01647633},
29 {0., 0.00392377}
30 };
31
32const Double_t i1Coefficient[kNCoeff][kNe] =
33 {
34 //coefficients to compute function I1
35 {0.5, 0.39894228},
36 {0.87890594, -0.03988024},
37 {0.51498869, -0.00362018},
38 {0.15084934, 0.00163801},
39 {0.02658733, -0.01031555},
40 {0.00301532, 0.02282967},
41 {0.00032411, -0.02895312},
42 {0., 0.01787654},
43 {0., -0.00420059}
44 };
45
46const Double_t k0Coefficient[kNCoeff][kNe] =
47 {
48 //coefficients to compute modified Hankel function of the zeroth order K0
49 {-0.57721566, 1.25331414},
50 {0.42278420, -0.07832358},
51 {0.23069756, 0.02189568},
52 {0.03488590, -0.01062446},
53 {0.00262698, 0.00587872},
54 {0.00010750, -0.00251540},
55 {0.0000074, 0.00053208},
56 {0., 0. },
57 {0., 0. }
58 };
59
60const Double_t k1Coefficient[kNCoeff][kNe] =
61 {
62 //coefficients to compute modified Hankel function of the first order K1
63 {1.0, 1.25331414},
64 {0.15443144, 0.23498619},
65 {-0.67278579, -0.03655620},
66 {-0.18156897, 0.01504268},
67 {-0.01919402, -0.00780353},
68 {-0.00110404, 0.00325614},
69 {-0.00004686, -0.00068245},
70 { 0., 0. },
71 { 0., 0. }
72 };
73
74Double_t BesselI0(Double_t x) {
75 // (C) Copr. 1986-92 Numerical Recipes Software +7.
76 //compute Bessel function of zeroth order
77
78 const Double_t p1 = i0Coefficient[0][0];
79 const Double_t p2 = i0Coefficient[1][0];
80 const Double_t p3 = i0Coefficient[2][0];
81 const Double_t p4 = i0Coefficient[3][0];
82 const Double_t p5 = i0Coefficient[4][0];
83 const Double_t p6 = i0Coefficient[5][0];
84 const Double_t p7 = i0Coefficient[6][0];
85
86 const Double_t q1 = i0Coefficient[0][1];
87 const Double_t q2 = i0Coefficient[1][1];
88 const Double_t q3 = i0Coefficient[2][1];
89 const Double_t q4 = i0Coefficient[3][1];
90 const Double_t q5 = i0Coefficient[4][1];
91 const Double_t q6 = i0Coefficient[5][1];
92 const Double_t q7 = i0Coefficient[6][1];
93 const Double_t q8 = i0Coefficient[7][1];
94 const Double_t q9 = i0Coefficient[8][1];
95
96 Double_t i0 = 0.;
97
98 if(TMath::Abs(x) < 3.75) {
99 Double_t y = (x / 3.75) * (x / 3.75);
100 i0 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
101 }
102 else {
103 Double_t ax = TMath::Abs(x);
104 Double_t y = 3.75 / ax;
105 i0 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
106 }
107
108 return i0;
109}
110
111Double_t BesselI1(Double_t x) {
112 // (C) Copr. 1986-92 Numerical Recipes Software +7.
113 //compute Bessel function of first order
114
115 const Double_t p1 = i1Coefficient[0][0];
116 const Double_t p2 = i1Coefficient[1][0];
117 const Double_t p3 = i1Coefficient[2][0];
118 const Double_t p4 = i1Coefficient[3][0];
119 const Double_t p5 = i1Coefficient[4][0];
120 const Double_t p6 = i1Coefficient[5][0];
121 const Double_t p7 = i1Coefficient[6][0];
122
123 const Double_t q1 = i1Coefficient[0][1];
124 const Double_t q2 = i1Coefficient[1][1];
125 const Double_t q3 = i1Coefficient[2][1];
126 const Double_t q4 = i1Coefficient[3][1];
127 const Double_t q5 = i1Coefficient[4][1];
128 const Double_t q6 = i1Coefficient[5][1];
129 const Double_t q7 = i1Coefficient[6][1];
130 const Double_t q8 = i1Coefficient[7][1];
131 const Double_t q9 = i1Coefficient[8][1];
132
133 Double_t i1 = 0.;
134
135 if (TMath::Abs(x) < 3.75) {
136 Double_t y = (x / 3.75) * (x / 3.75);
137 i1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
138 }
139 else {
140 Double_t ax = TMath::Abs(x);
141 Double_t y = 3.75/ax;
142 i1 = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))));
143 if(x < 0.) i1 = -i1;
144 }
145
146 return i1;
147}
148
149Double_t HankelK0(Double_t x) {
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