coverity 15108 fixed
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowBessel.cxx
1 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
2  *
3  * This library is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU Lesser General Public License
5  * as published by the Free Software Foundation; either version 2.1 of
6  * the License, or (at your option) any later version.
7  *
8  * This library is distributed in the hope that it will be useful, but
9  * WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11  * Lesser General Public License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with this library; if not, write to the Free Software
15  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
16  * USA
17  */
18 /** @file 
19     @brief Implementation of Bessel functions */
20 #include "flow/AliFMDFlowBessel.h"
21 #include <cmath>
22 #include <iostream>
23
24 /** Namespace for utility functions */
25 namespace 
26
27   //__________________________________________________________________
28   /** Utility function to compute 
29       @f[ 
30       envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
31       @f]
32       @param n Order 
33       @param x Argument 
34       @return @f$ envj_n(x)@f$ */
35   Double_t Envj(UInt_t n, Double_t x) 
36   { 
37     // Compute 
38     //   1/2 log10(6.28*n) - n log19(1.36 * x / n)
39     return .5 * log10(6.28 * n) - n * log10(1.36 * x / n);
40   }
41   //__________________________________________________________________
42   /** Utility function to do loop in Msta1 and Msta2 */ 
43   UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp) 
44   {
45     // Utility function to do loop in Msta1 and Msta2
46     Double_t f0 = Envj(n0, a0) - mp;
47     Int_t    n1 = n0 + 5;
48     Double_t f1 = Envj(n1, a0) - mp;
49     Int_t    nn = 0;
50     for (UInt_t i = 0; i < 20; i++) { 
51       nn = Int_t(n1 - (n1 - n0) / (1. - f0 / f1));
52       if (fabs(Double_t(nn - n1)) < 1) break;
53       n0 = n1;
54       f0 = f1;
55       n1 = nn;
56       f1 = Envj(nn, a0) - mp;
57     }
58     return nn;
59   }
60   //__________________________________________________________________
61   /** Determine the starting point for backward recurrence such that
62       the magnitude of @f$J_n(x)@f$ at that point is about @f$
63       10^{-mp}@f$  
64       @param x   Argument to @f$ J_n(x)@f$ 
65       @param mp  Value of magnitude @f$ mp@f$ 
66       @return The starting point */
67   UInt_t Msta1(Double_t x, UInt_t mp) 
68   { 
69     // Determine the starting point for backward recurrence such that
70     // the magnitude of J_n(x) at that point is about 10^{-mp}
71     Double_t a0 = fabs(x);
72     Int_t    n0 = Int_t(1.1 * a0) + 1;
73     return MstaLoop(n0, a0, mp);
74   }
75   //__________________________________________________________________
76   /** Determine the starting point for backward recurrence such that
77       all @f$J_n(x)@f$ has @a mp significant digits.
78       @param x   Argument to @f$ J_n(x)@f$ 
79       @param n   Order of @f$ J_n@f$
80       @param mp  Number of significant digits 
81       @reutnr starting point */
82   UInt_t Msta2(Double_t x, Int_t n, Int_t mp) 
83   { 
84     // Determine the starting point for backward recurrence such that 
85     // all J_n(x) has mp significant digits. 
86     Double_t a0  = fabs(x);
87     Double_t hmp = 0.5 * mp;
88     Double_t ejn = Envj(n, a0);
89     Double_t obj = 0;
90     Int_t    n0  = 0;
91     if (ejn <= hmp) {
92       obj = mp;
93       n0  = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
94     }
95     else { 
96       obj = hmp + ejn;
97       n0  = n;
98     }
99     return MstaLoop(n0, a0, obj) + 10;
100   }
101 }
102 //____________________________________________________________________
103 void 
104 AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0, 
105                       Double_t& bi1, Double_t& di1)
106 {
107   // Compute the modified Bessel functions 
108   // 
109   //   I_0(x) = \sum_{k=1}^\infty (x^2/4)^k}/(k!)^2
110   // 
111   // and I_1(x), and their first derivatives  
112   Double_t       x2 =  x * x;
113   if (x == 0) { 
114     bi0 = 1;
115     bi1 = 0;
116     di0 = 0;
117     di1 = .5;
118     return;
119   }
120     
121   if (x <= 18) { 
122     bi0      = 1;
123     Double_t r = 1;
124     for (UInt_t k = 1; k <= 50; k++) { 
125       r   =  .25 * r * x2 / (k * k);
126       bi0 += r;
127       if (fabs(r / bi0) < 1e-15) break;
128     }
129     bi1 = 1;
130     r   = 1;
131     for (UInt_t k = 1; k <= 50; k++) { 
132       r   =  .25 * r * x2 / (k * (k + 1));
133       bi1 += r;
134     }
135     bi1 *= .5 * x;
136   }
137   else { 
138     const Double_t a[] = { 0.125,            0.0703125, 
139                          0.0732421875,     0.11215209960938, 
140                          0.22710800170898, 0.57250142097473, 
141                          1.7277275025845,  6.0740420012735, 
142                          24.380529699556,  110.01714026925,
143                          551.33589612202,  3038.0905109224 };
144     const Double_t b[] = { -0.375,           -0.1171875, 
145                          -0.1025390625,    -0.14419555664063,
146                          -0.2775764465332, -0.67659258842468,
147                          -1.9935317337513, -6.8839142681099,
148                          -27.248827311269, -121.59789187654,
149                          -603.84407670507, -3302.2722944809 };
150     UInt_t k0 = 12;
151     if (x >= 35) k0 = 9;
152     if (x >= 50) k0 = 7;
153       
154     Double_t ca = exp(x) / sqrt(2 * M_PI * x);
155     Double_t xr = 1. / x;
156     bi0       = 1;
157       
158     for (UInt_t k = 1; k <= k0; k++) {
159       bi0 += a[k-1] * pow(xr, Int_t(k));
160       bi1 += b[k-1] * pow(xr, Int_t(k));
161     }
162       
163     bi0 *= ca;
164     bi1 *= ca;
165   }
166   di0 = bi1;
167   di1 = bi0 - bi1 / x;
168 }
169 //____________________________________________________________________
170 UInt_t 
171 AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di) 
172 {
173   // Compute the modified Bessel functions I_{n/2}(x) and their 
174   // derivatives.  
175   typedef Double_t (*fun_t)(Double_t);
176   Int_t   p = (n > 0 ? 1     : -1);
177   fun_t s = (n > 0 ? &sinh : &cosh);
178   fun_t c = (n > 0 ? &cosh : &sinh);
179   
180   // Temporary buffers - only grows in size. 
181   static Double_t* tbi = 0;
182   static Double_t* tdi = 0;
183   static Int_t     tn  = 0;
184   if (!tbi || tn < p*n/2+2) { 
185     if (tbi) delete [] tbi;
186     if (tdi) delete [] tdi;
187     tn = p*n/2+2;
188     tbi = new Double_t[tn];
189     tdi = new Double_t[tn];
190   }
191   // Temporary buffer 
192   // Double_t tbi[p*n/2+2];
193   // Double_t tdi[p*n/2+2];
194   
195   // Calculate I(-1/2) for n>0 or I(1/2) for n<0 
196   Double_t bim = sqrt(Double_t(2)) * (*c)(x) / sqrt(M_PI * x);
197
198   // We calculate one more than n, that is we calculate 
199   //   I(1/2), I(3/2),...,I(n/2+1) 
200   // because we need that for the negative orders 
201   for (Int_t i = 1; i <= p * n + 2; i += 2) { 
202     switch (i) { 
203     case 1: 
204       // Explicit formula for 1/2
205       tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
206       break;
207     case 3:
208       // Explicit formula for 3/2
209       tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
210       break;
211     default:
212       // Recursive formula for n/2 in terms of (n/2-2) and (n/2-1) 
213       tbi[i/2] = tbi[i/2-2] - (i-2) * tbi[i/2-1] / x;
214       break;
215     }
216   }
217   // We calculate the first one dI(1/2) here, so that we can use it in
218   // the loop. 
219   tdi[0] = bim - tbi[0] / (2 * x);
220
221   // Now, calculate dI(3/2), ..., dI(n/2)
222   for (Int_t i = 3; i <= p * n; i += 2) {
223     tdi[i/2] = tbi[i/2-p*1] - p * i * tbi[i/2] / (2 * x);
224   }
225
226   // Finally, store the output 
227   for (Int_t i = 0; i < p * n / 2 + 1; i++) { 
228     UInt_t j = (p < 0 ? p * n / 2 - i : i);
229     bi[i]    = tbi[j];
230     di[i]    = tdi[j];
231   }
232   return n / 2;
233 }
234                 
235 //____________________________________________________________________
236 UInt_t 
237 AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di) 
238 {
239   // Compute the modified Bessel functions I_n(x) and their
240   // derivatives.  Note, that I_{-n}(x) = I_n(x) and 
241   // dI_{-n}(x)/dx = dI_{n}(x)/dx so this function can be used 
242   // to evaluate I_n for all natural numbers  
243   UInt_t mn = in;
244   if (x < 1e-100) { 
245     for (UInt_t k = 0; k <= in; k++) bi[k] = di[k]  = 0;
246     bi[0] = 1;
247     di[1] = .5;
248     return 1;
249   }
250   Double_t bi0, di0, bi1, di1;
251   I01(x, bi0, di0, bi1, di1);
252   bi[0] = bi0; di[0] = di0; bi[1] = bi1; di[1] = di1; 
253
254   if (in <= 1) return 2;
255     
256   if (x > 40 && Int_t(in) < Int_t(.25 * x)) { 
257     Double_t h0 = bi0;
258     Double_t h1 = bi1;
259     for (UInt_t k = 2; k <= in; k++) { 
260       bi[k] = -2 * (k - 1) / x * h1 + h0;
261       h0    = h1;
262       h1    = bi[k];
263     }
264   }
265   else { 
266     UInt_t      m  = Msta1(x, 100);
267     if (m < in) mn = m;
268     else        m  = Msta2(x, in, 15);
269       
270     Double_t f0 = 0;
271     Double_t f1 = 1e-100;
272     Double_t f  = 0;
273     for (Int_t k = m; k >= 0; k--) { 
274       f  = 2 * (k + 1) * f1 / x + f0;
275       if (k <= Int_t(mn)) bi[k] = f;
276       f0 = f1;
277       f1 = f;
278     }
279     Double_t s0 = bi0 / f;
280     for (UInt_t k = 0; k <= mn; k++) 
281       bi[k] *= s0;
282   }
283   for (UInt_t k = 2; k <= mn; k++) 
284     di[k] =  bi[k - 1] - k / x * bi[k];
285   return mn;
286 }
287
288 //____________________________________________________________________
289 namespace 
290 {
291   Double_t* EnsureSize(Double_t*& a, UInt_t& old, const UInt_t size)
292   {
293     if (a && old < size) { 
294       // Will delete, and reallocate. 
295       if (old != 1) delete [] a;
296       a = 0;
297     }
298     if (!a && size > 0) {
299       // const UInt_t n = (size <= 1 ? 2 : size);
300       a              = new Double_t[size];
301       old            = size; // n;
302     }
303     return a;
304   }
305 }
306       
307 //____________________________________________________________________
308 UInt_t 
309 AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x, 
310                       Double_t* bi, Double_t* di)
311 {
312   // Compute the modified Bessel functions I_{\nu}(x) and
313   // their derivatives for \nu = n_1, n_1+1, \ldots, n_2 for
314   // and integer n_1, n_2 or any half-integer n_1, n_2. 
315   UInt_t  in1 = UInt_t(fabs(n1));
316   UInt_t  in2 = UInt_t(fabs(n2));
317   UInt_t  nt  = UInt_t(n2 - n1 + 1);
318   static Double_t* tbi = 0;
319   static Double_t* tdi = 0;
320   static UInt_t    tn  = 0;
321
322   if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument 
323     if (Int_t(2 * fabs(n2)) % 2 != 1) { 
324       std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
325                 << "must also be half integer" << std::endl;
326       return 0;
327     }
328     // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;
329
330     Int_t s1    = (n1 < 0 ? -1 : 1);
331     Int_t s2    = (n2 < 0 ? -1 : 1);
332     Int_t l1    = (s1 < 0 ? (2*in1+1)/2 : 0);
333     Int_t l2    = (s1 > 0 ? in1 : 0);
334
335     // Temporary buffers - only grows in size. 
336     tbi = EnsureSize(tbi, tn, nt+1);
337     tdi = EnsureSize(tdi, tn, nt+1);
338     // Double_t tbi[nt+1], tdi[nt+1];
339
340     if (s1 < 0) { 
341       Ihalf(Int_t(2 * n1), x, tbi, tdi);
342       for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) { 
343         bi[i] = tbi[i];
344         di[i] = tdi[i];
345       }
346     }
347     if (s2 > 0) { 
348       // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
349       Ihalf(Int_t(2 * n2), x, tbi, tdi);
350       for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) { 
351         UInt_t j = i + l2 - l1;
352         bi[i] = tbi[j];
353         di[i] = tdi[j];
354       }
355     }
356     return nt;
357   }
358   if (Int_t(n1) != n1 || Int_t(n2) != n2) { 
359     std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
360               << "n2 (" << n2 << "/" << in2 << ") must be "
361               << "half-integer or integer" << std::endl;
362     return 0; // Not integer!
363   }
364
365   UInt_t  n   = UInt_t(in1 > in2 ? in1 : in2);
366   // Temporary buffers - only grows in size. 
367   tbi = EnsureSize(tbi, tn, n+1);
368   tdi = EnsureSize(tdi, tn, n+1);
369   // Double_t  tbi[n+1];
370   // Double_t  tdi[n+1];
371
372   UInt_t  r   = Iwhole(n, x, tbi, tdi);
373   if (r < n) 
374     std::cerr << "Only got " << r << "/" << n << " values" 
375               << std::endl;
376   for (UInt_t i = 0; i < nt; i++) { 
377     UInt_t j = UInt_t(i + n1 < 0 ? -n1-i : n1+i);
378     bi[i]    = tbi[j];
379     di[i]    = tdi[j];
380   }
381   return nt;
382 }
383
384     
385 //____________________________________________________________________
386 Double_t 
387 AliFMDFlowBessel::I(Double_t n, Double_t x) 
388
389   // Compute the modified Bessel function of the first kind 
390   // 
391   //   I_n(x) = (x/2)^n \sum_{k=0}^\infty (x^2/4)^k / (k!\Gamma(n+k+1))
392   // 
393   // for arbirary integer order n 
394   Double_t i[2], di[2];
395   Inu(n, n, x, i, di);
396   return i[0];
397 }
398       
399 //____________________________________________________________________
400 Double_t 
401 AliFMDFlowBessel::DiffI(Double_t n, Double_t x) 
402
403   // Compute the derivative of the modified Bessel function of the
404   // first kind 
405   // 
406   //    dI_n(x)/dx = I_{n-1}(x) - n/x I_{n}(x)
407   // 
408   // for arbirary integer order n
409   Double_t i[2], di[2];
410   Inu(n, n, x, i, di);
411   return di[0];
412 }
413
414 //____________________________________________________________________
415 //
416 // EOF
417 //