bug fix
[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowBessel.cxx
CommitLineData
97e94238 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 */
39eefe19 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 */
25namespace
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$ */
97e94238 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)
39eefe19 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 */
97e94238 43 UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp)
44 {
45 // Utility function to do loop in Msta1 and Msta2
39eefe19 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++) {
6ce810fc 51 nn = Int_t(n1 - (n1 - n0) / (1. - f0 / f1));
2a4540e6 52 if (fabs(Double_t(nn - n1)) < 1) break;
39eefe19 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 */
97e94238 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}
39eefe19 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
97e94238 77 all @f$J_n(x)@f$ has @a mp significant digits.
39eefe19 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 */
97e94238 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.
39eefe19 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//____________________________________________________________________
103void
97e94238 104AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0,
105 Double_t& bi1, Double_t& di1)
39eefe19 106{
97e94238 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
39eefe19 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//____________________________________________________________________
170UInt_t
171AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di)
172{
97e94238 173 // Compute the modified Bessel functions I_{n/2}(x) and their
174 // derivatives.
39eefe19 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
2219d590 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 }
39eefe19 191 // Temporary buffer
2219d590 192 // Double_t tbi[p*n/2+2];
193 // Double_t tdi[p*n/2+2];
39eefe19 194
195 // Calculate I(-1/2) for n>0 or I(1/2) for n<0
2a4540e6 196 Double_t bim = sqrt(Double_t(2)) * (*c)(x) / sqrt(M_PI * x);
39eefe19 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//____________________________________________________________________
236UInt_t
237AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di)
238{
97e94238 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
39eefe19 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
97e94238 256 if (x > 40 && Int_t(in) < Int_t(.25 * x)) {
39eefe19 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//____________________________________________________________________
bbcd8619 289namespace
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//____________________________________________________________________
39eefe19 308UInt_t
97e94238 309AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
310 Double_t* bi, Double_t* di)
39eefe19 311{
97e94238 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.
39eefe19 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);
bbcd8619 318 static Double_t* tbi = 0;
319 static Double_t* tdi = 0;
320 static UInt_t tn = 0;
321
39eefe19 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);
2219d590 334
335 // Temporary buffers - only grows in size.
bbcd8619 336 tbi = EnsureSize(tbi, tn, nt+1);
337 tdi = EnsureSize(tdi, tn, nt+1);
2219d590 338 // Double_t tbi[nt+1], tdi[nt+1];
339
39eefe19 340 if (s1 < 0) {
6ce810fc 341 Ihalf(Int_t(2 * n1), x, tbi, tdi);
39eefe19 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 << ",...)";
6ce810fc 349 Ihalf(Int_t(2 * n2), x, tbi, tdi);
39eefe19 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);
bbcd8619 366 // Temporary buffers - only grows in size.
367 tbi = EnsureSize(tbi, tn, n+1);
368 tdi = EnsureSize(tdi, tn, n+1);
2219d590 369 // Double_t tbi[n+1];
370 // Double_t tdi[n+1];
bbcd8619 371
39eefe19 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++) {
6ce810fc 377 UInt_t j = UInt_t(i + n1 < 0 ? -n1-i : n1+i);
39eefe19 378 bi[i] = tbi[j];
379 di[i] = tdi[j];
380 }
381 return nt;
382}
383
384
385//____________________________________________________________________
386Double_t
387AliFMDFlowBessel::I(Double_t n, Double_t x)
388{
97e94238 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
39eefe19 394 Double_t i[2], di[2];
395 Inu(n, n, x, i, di);
396 return i[0];
397}
398
399//____________________________________________________________________
400Double_t
401AliFMDFlowBessel::DiffI(Double_t n, Double_t x)
402{
97e94238 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
39eefe19 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//