]>
Commit | Line | Data |
---|---|---|
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 | // |