1 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
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.
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.
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
19 @brief Implementation of Bessel functions */
20 #include "flow/AliFMDFlowBessel.h"
24 /** Namespace for utility functions */
27 //__________________________________________________________________
28 /** Utility function to compute
30 envj_n(x) = \frac12 \log_{10}(6.28 n) - n \log_{10}(1.36 x / n)
34 @return @f$ envj_n(x)@f$ */
35 Double_t Envj(UInt_t n, Double_t x)
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);
41 //__________________________________________________________________
42 /** Utility function to do loop in Msta1 and Msta2 */
43 UInt_t MstaLoop(Int_t n0, Double_t a0, Float_t mp)
45 // Utility function to do loop in Msta1 and Msta2
46 Double_t f0 = Envj(n0, a0) - mp;
48 Double_t f1 = Envj(n1, a0) - mp;
50 for (UInt_t i = 0; i < 20; i++) {
51 nn = n1 - (n1 - n0) / (1. - f0 / f1);
52 if (fabs(nn - n1) < 1) break;
56 f1 = Envj(nn, a0) - mp;
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$
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)
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);
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)
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);
93 n0 = Int_t(1.1 * a0) + 1; //Bug for x<0.1-vl, 2-8.2002
99 return MstaLoop(n0, a0, obj) + 10;
102 //____________________________________________________________________
104 AliFMDFlowBessel::I01(Double_t x, Double_t& bi0, Double_t& di0,
105 Double_t& bi1, Double_t& di1)
107 // Compute the modified Bessel functions
109 // I_0(x) = \sum_{k=1}^\infty (x^2/4)^k}/(k!)^2
111 // and I_1(x), and their first derivatives
124 for (UInt_t k = 1; k <= 50; k++) {
125 r = .25 * r * x2 / (k * k);
127 if (fabs(r / bi0) < 1e-15) break;
131 for (UInt_t k = 1; k <= 50; k++) {
132 r = .25 * r * x2 / (k * (k + 1));
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 };
154 Double_t ca = exp(x) / sqrt(2 * M_PI * x);
155 Double_t xr = 1. / x;
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));
169 //____________________________________________________________________
171 AliFMDFlowBessel::Ihalf(Int_t n, Double_t x, Double_t* bi, Double_t* di)
173 // Compute the modified Bessel functions I_{n/2}(x) and their
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);
180 // Temporary buffers - only grows in size.
181 static Double_t* tbi = 0;
182 static Double_t* tdi = 0;
184 if (!tbi || tn < p*n/2+2) {
185 if (tbi) delete [] tbi;
186 if (tdi) delete [] tdi;
188 tbi = new Double_t[tn];
189 tdi = new Double_t[tn];
192 // Double_t tbi[p*n/2+2];
193 // Double_t tdi[p*n/2+2];
195 // Calculate I(-1/2) for n>0 or I(1/2) for n<0
196 Double_t bim = sqrt(2) * (*c)(x) / sqrt(M_PI * x);
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) {
204 // Explicit formula for 1/2
205 tbi[0] = sqrt(2 * x / M_PI) * (*s)(x) / x;
208 // Explicit formula for 3/2
209 tbi[1] = sqrt(2 * x / M_PI) * ((*c)(x) / x - (*s)(x) / (x * x));
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;
217 // We calculate the first one dI(1/2) here, so that we can use it in
219 tdi[0] = bim - tbi[0] / (2 * x);
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);
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);
235 //____________________________________________________________________
237 AliFMDFlowBessel::Iwhole(UInt_t in, Double_t x, Double_t* bi, Double_t* di)
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
245 for (UInt_t k = 0; k <= in; k++) bi[k] = di[k] = 0;
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;
254 if (in <= 1) return 2;
256 if (x > 40 && Int_t(in) < Int_t(.25 * x)) {
259 for (UInt_t k = 2; k <= in; k++) {
260 bi[k] = -2 * (k - 1) / x * h1 + h0;
266 UInt_t m = Msta1(x, 100);
268 else m = Msta2(x, in, 15);
271 Double_t f1 = 1e-100;
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;
279 Double_t s0 = bi0 / f;
280 for (UInt_t k = 0; k <= mn; k++)
283 for (UInt_t k = 2; k <= mn; k++)
284 di[k] = bi[k - 1] - k / x * bi[k];
288 //____________________________________________________________________
290 AliFMDFlowBessel::Inu(Double_t n1, Double_t n2, Double_t x,
291 Double_t* bi, Double_t* di)
293 // Compute the modified Bessel functions I_{\nu}(x) and
294 // their derivatives for \nu = n_1, n_1+1, \ldots, n_2 for
295 // and integer n_1, n_2 or any half-integer n_1, n_2.
296 UInt_t in1 = UInt_t(fabs(n1));
297 UInt_t in2 = UInt_t(fabs(n2));
298 UInt_t nt = UInt_t(n2 - n1 + 1);
299 if (Int_t(2 * fabs(n1)) % 2 == 1) { // Half-integer argument
300 if (Int_t(2 * fabs(n2)) % 2 != 1) {
301 std::cout << "If n1 is half-integer (" << n1 << ") then n2 "
302 << "must also be half integer" << std::endl;
305 // std::cout << "Half-integers " << n1 << "," << n2 << std::endl;
307 Int_t s1 = (n1 < 0 ? -1 : 1);
308 Int_t s2 = (n2 < 0 ? -1 : 1);
309 Int_t l1 = (s1 < 0 ? (2*in1+1)/2 : 0);
310 Int_t l2 = (s1 > 0 ? in1 : 0);
312 // Temporary buffers - only grows in size.
313 static Double_t* tbi = 0;
314 static Double_t* tdi = 0;
315 static UInt_t tn = 0;
316 if (!tbi || tn < nt+1) {
317 if (tbi) delete [] tbi;
318 if (tdi) delete [] tdi;
320 tbi = new Double_t[tn];
321 tdi = new Double_t[tn];
323 // Double_t tbi[nt+1], tdi[nt+1];
326 Ihalf(2 * n1, x, tbi, tdi);
327 for (Int_t i = 0; i < l1 && i < Int_t(nt); i++) {
333 // std::cout << "Evaluating Ihalf(" << 2 * n2 << "," << x << ",...)";
334 Ihalf(2 * n2, x, tbi, tdi);
335 for (Int_t i = l1; i <= 2 * Int_t(in2) && i < Int_t(nt); i++) {
336 UInt_t j = i + l2 - l1;
343 if (Int_t(n1) != n1 || Int_t(n2) != n2) {
344 std::cerr << "n1 (" << n1 << "/" << in1 << ") and "
345 << "n2 (" << n2 << "/" << in2 << ") must be "
346 << "half-integer or integer" << std::endl;
347 return 0; // Not integer!
350 UInt_t n = UInt_t(in1 > in2 ? in1 : in2);
351 static Double_t* tbi = 0;
352 static Double_t* tdi = 0;
353 static UInt_t tn = 0;
354 if (!tbi || tn < n+1) {
355 if (tbi) delete [] tbi;
356 if (tdi) delete [] tdi;
358 tbi = new Double_t[tn];
359 tdi = new Double_t[tn];
361 // Double_t tbi[n+1];
362 // Double_t tdi[n+1];
363 UInt_t r = Iwhole(n, x, tbi, tdi);
365 std::cerr << "Only got " << r << "/" << n << " values"
367 for (UInt_t i = 0; i < nt; i++) {
368 UInt_t j = (i + n1 < 0 ? -n1-i : n1+i);
376 //____________________________________________________________________
378 AliFMDFlowBessel::I(Double_t n, Double_t x)
380 // Compute the modified Bessel function of the first kind
382 // I_n(x) = (x/2)^n \sum_{k=0}^\infty (x^2/4)^k / (k!\Gamma(n+k+1))
384 // for arbirary integer order n
385 Double_t i[2], di[2];
390 //____________________________________________________________________
392 AliFMDFlowBessel::DiffI(Double_t n, Double_t x)
394 // Compute the derivative of the modified Bessel function of the
397 // dI_n(x)/dx = I_{n-1}(x) - n/x I_{n}(x)
399 // for arbirary integer order n
400 Double_t i[2], di[2];
405 //____________________________________________________________________