New generator: TUHKMgen
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / HankelFunction.cxx
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
16 enum {kNe = 2, kNCoeff = 9};
17       
18 const 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
32 const 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
46 const 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
60 const 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
74 Double_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
111 Double_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
149 Double_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
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