CbmRoot
xMath.cxx
Go to the documentation of this file.
1 #include "xMath.h"
2 #include <cmath>
3 
4 using namespace std;
5 
6 //______________________________________________________________________________
7 double xMath::BesselI0(double x) {
8  // Compute the modified Bessel function I_0(x) for any real x.
9  //
10  //--- NvE 12-mar-2000 UU-SAP Utrecht
11 
12  // Parameters of the polynomial approximation
13  const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424, p4 = 1.2067492,
14  p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
15 
16  const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
17  q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
18  q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
19 
20  const double k1 = 3.75;
21  double ax = fabs(x); //fabs(x);
22 
23  double y = 0, result = 0;
24 
25  if (ax < k1) {
26  double xx = x / k1;
27  y = xx * xx;
28  result = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
29  } else {
30  y = k1 / ax;
31  result =
32  (exp(ax) / sqrt(ax))
33  * (q1
34  + y
35  * (q2
36  + y
37  * (q3
38  + y
39  * (q4
40  + y
41  * (q5
42  + y
43  * (q6
44  + y
45  * (q7
46  + y * (q8 + y * q9))))))));
47  }
48  return result;
49 }
50 
51 //______________________________________________________________________________
52 double xMath::BesselK0(double x) {
53  // Compute the modified Bessel function K_0(x) for positive real x.
54  //
55  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
56  // Applied Mathematics Series vol. 55 (1964), Washington.
57  //
58  //--- NvE 12-mar-2000 UU-SAP Utrecht
59 
60  // Parameters of the polynomial approximation
61  const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
62  p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
63 
64  const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
65  q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3,
66  q7 = 5.3208e-4;
67 
68  if (x <= 0) {
69  //Error("xMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
70  return 0;
71  }
72 
73  double y = 0, result = 0;
74 
75  if (x <= 2) {
76  y = x * x / 4;
77  result =
78  (-log(x / 2.) * xMath::BesselI0(x))
79  + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
80  } else {
81  y = 2 / x;
82  result =
83  (exp(-x) / sqrt(x))
84  * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
85  }
86  return result;
87 }
88 
89 //______________________________________________________________________________
90 double xMath::BesselI1(double x) {
91  // Compute the modified Bessel function I_1(x) for any real x.
92  //
93  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
94  // Applied Mathematics Series vol. 55 (1964), Washington.
95  //
96  //--- NvE 12-mar-2000 UU-SAP Utrecht
97 
98  // Parameters of the polynomial approximation
99  const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869, p4 = 0.15084934,
100  p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
101 
102  const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
103  q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
104  q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
105 
106  const double k1 = 3.75;
107  double ax = fabs(x); //fabs(x);
108 
109  double y = 0, result = 0;
110 
111  if (ax < k1) {
112  double xx = x / k1;
113  y = xx * xx;
114  result =
115  x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
116  } else {
117  y = k1 / ax;
118  result =
119  (exp(ax) / sqrt(ax))
120  * (q1
121  + y
122  * (q2
123  + y
124  * (q3
125  + y
126  * (q4
127  + y
128  * (q5
129  + y
130  * (q6
131  + y
132  * (q7
133  + y * (q8 + y * q9))))))));
134  if (x < 0) result = -result;
135  }
136  return result;
137 }
138 
139 //______________________________________________________________________________
140 double xMath::BesselK1(double x) {
141  // Compute the modified Bessel function K_1(x) for positive real x.
142  //
143  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
144  // Applied Mathematics Series vol. 55 (1964), Washington.
145  //
146  //--- NvE 12-mar-2000 UU-SAP Utrecht
147 
148  // Parameters of the polynomial approximation
149  const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579, p4 = -0.18156897,
150  p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
151 
152  const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
153  q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3,
154  q7 = -6.8245e-4;
155 
156  if (x <= 0) {
157  //Error("xMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
158  return 0;
159  }
160 
161  double y = 0, result = 0;
162 
163  if (x <= 2) {
164  y = x * x / 4;
165  result =
166  (log(x / 2.) * xMath::BesselI1(x))
167  + (1. / x)
168  * (p1
169  + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
170  } else {
171  y = 2 / x;
172  result =
173  (exp(-x) / sqrt(x))
174  * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
175  }
176  return result;
177 }
178 
179 //______________________________________________________________________________
180 double xMath::BesselK(int n, double x) {
181  // Compute the Integer Order Modified Bessel function K_n(x)
182  // for n=0,1,2,... and positive real x.
183  //
184  //--- NvE 12-mar-2000 UU-SAP Utrecht
185 
186  if (x <= 0 || n < 0) {
187  //Error("xMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
188  return 0;
189  }
190 
191  if (n == 0) return xMath::BesselK0(x);
192  if (n == 1) return xMath::BesselK1(x);
193 
194  // Perform upward recurrence for all x
195  double tox = 2 / x;
196  double bkm = xMath::BesselK0(x);
197  double bk = xMath::BesselK1(x);
198  double bkp = 0;
199  for (int j = 1; j < n; j++) {
200  bkp = bkm + double(j) * tox * bk;
201  bkm = bk;
202  bk = bkp;
203  }
204  return bk;
205 }
206 
207 //______________________________________________________________________________
208 double xMath::BesselI(int n, double x) {
209  // Compute the Integer Order Modified Bessel function I_n(x)
210  // for n=0,1,2,... and any real x.
211  //
212  //--- NvE 12-mar-2000 UU-SAP Utrecht
213 
214  int iacc = 40; // Increase to enhance accuracy
215  const double kBigPositive = 1.e10;
216  const double kBigNegative = 1.e-10;
217 
218  if (n < 0) {
219  //Error("xMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
220  return 0;
221  }
222 
223  if (n == 0) return xMath::BesselI0(x);
224  if (n == 1) return xMath::BesselI1(x);
225 
226  if (x == 0) return 0;
227  if (fabs(x) > kBigPositive) return 0;
228 
229  double tox = 2 / fabs(x);
230  double bip = 0, bim = 0;
231  double bi = 1;
232  double result = 0;
233  int m = 2 * ((n + int(sqrt(float(iacc * n)))));
234  for (int j = m; j >= 1; j--) {
235  bim = bip + double(j) * tox * bi;
236  bip = bi;
237  bi = bim;
238  // Renormalise to prevent overflows
239  if (fabs(bi) > kBigPositive) {
240  result *= kBigNegative;
241  bi *= kBigNegative;
242  bip *= kBigNegative;
243  }
244  if (j == n) result = bip;
245  }
246 
247  result *= xMath::BesselI0(x) / bi; // Normalise with BesselI0(x)
248  if ((x < 0) && (n % 2 == 1)) result = -result;
249 
250  return result;
251 }
252 
253 //______________________________________________________________________________
254 double xMath::BesselJ0(double x) {
255  // Returns the Bessel function J0(x) for any real x.
256 
257  double ax, z;
258  double xx, y, result, result1, result2;
259  const double p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
260  const double p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
261  const double p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
262  const double p10 = 59272.64853, p11 = 267.8532712;
263 
264  const double q1 = 0.785398164;
265  const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
266  const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
267  const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
268  const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
269  const double q10 = 0.934935152e-7, q11 = 0.636619772;
270 
271  if ((ax = fabs(x)) < 8) {
272  y = x * x;
273  result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
274  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
275  result = result1 / result2;
276  } else {
277  z = 8 / ax;
278  y = z * z;
279  xx = ax - q1;
280  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
281  result2 = q6 + y * (q7 + y * (q8 + y * (q9 - y * q10)));
282  result = sqrt(q11 / ax) * (cos(xx) * result1 - z * sin(xx) * result2);
283  }
284  return result;
285 }
286 
287 //______________________________________________________________________________
288 double xMath::BesselJ1(double x) {
289  // Returns the Bessel function J1(x) for any real x.
290 
291  double ax, z;
292  double xx, y, result, result1, result2;
293  const double p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
294  const double p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
295  const double p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
296  const double p10 = 99447.43394, p11 = 376.9991397;
297 
298  const double q1 = 2.356194491;
299  const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
300  const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
301  const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
302  const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
303  const double q10 = 0.105787412e-6, q11 = 0.636619772;
304 
305  if ((ax = fabs(x)) < 8) {
306  y = x * x;
307  result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
308  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
309  result = result1 / result2;
310  } else {
311  z = 8 / ax;
312  y = z * z;
313  xx = ax - q1;
314  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
315  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
316  result = sqrt(q11 / ax) * (cos(xx) * result1 - z * sin(xx) * result2);
317  if (x < 0) result = -result;
318  }
319  return result;
320 }
321 
322 //______________________________________________________________________________
323 double xMath::BesselY0(double x) {
324  // Returns the Bessel function Y0(x) for positive x.
325 
326  double z, xx, y, result, result1, result2;
327  const double p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
328  const double p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
329  const double p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
330  const double p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
331 
332  const double q1 = 0.785398164;
333  const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
334  const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
335  const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
336  const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
337  const double q10 = -0.934945152e-7, q11 = 0.636619772;
338 
339  if (x < 8) {
340  y = x * x;
341  result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
342  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
343  result = (result1 / result2) + p12 * xMath::BesselJ0(x) * log(x);
344  } else {
345  z = 8 / x;
346  y = z * z;
347  xx = x - q1;
348  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
349  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
350  result = sqrt(q11 / x) * (sin(xx) * result1 + z * cos(xx) * result2);
351  }
352  return result;
353 }
354 
355 //______________________________________________________________________________
356 double xMath::BesselY1(double x) {
357  // Returns the Bessel function Y1(x) for positive x.
358 
359  double z, xx, y, result, result1, result2;
360  const double p1 = -0.4900604943e13, p2 = 0.1275274390e13;
361  const double p3 = -0.5153438139e11, p4 = 0.7349264551e9;
362  const double p5 = -0.4237922726e7, p6 = 0.8511937935e4;
363  const double p7 = 0.2499580570e14, p8 = 0.4244419664e12;
364  const double p9 = 0.3733650367e10, p10 = 0.2245904002e8;
365  const double p11 = 0.1020426050e6, p12 = 0.3549632885e3;
366  const double p13 = 0.636619772;
367  const double q1 = 2.356194491;
368  const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
369  const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
370  const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
371  const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
372  const double q10 = 0.105787412e-6, q11 = 0.636619772;
373 
374  if (x < 8) {
375  y = x * x;
376  result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
377  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y * (p12 + y)))));
378  result = (result1 / result2) + p13 * (xMath::BesselJ1(x) * log(x) - 1 / x);
379  } else {
380  z = 8 / x;
381  y = z * z;
382  xx = x - q1;
383  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
384  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
385  result = sqrt(q11 / x) * (sin(xx) * result1 + z * cos(xx) * result2);
386  }
387  return result;
388 }
389 
390 //______________________________________________________________________________
391 double xMath::StruveH0(double x) {
392  // Struve Functions of Order 0
393  //
394  // Converted from CERNLIB M342 by Rene Brun.
395 
396  const int n1 = 15;
397  const int n2 = 25;
398  const double c1[16] = {1.00215845609911981,
399  -1.63969292681309147,
400  1.50236939618292819,
401  -.72485115302121872,
402  .18955327371093136,
403  -.03067052022988,
404  .00337561447375194,
405  -2.6965014312602e-4,
406  1.637461692612e-5,
407  -7.8244408508e-7,
408  3.021593188e-8,
409  -9.6326645e-10,
410  2.579337e-11,
411  -5.8854e-13,
412  1.158e-14,
413  -2e-16};
414  const double c2[26] = {.99283727576423943,
415  -.00696891281138625,
416  1.8205103787037e-4,
417  -1.063258252844e-5,
418  9.8198294287e-7,
419  -1.2250645445e-7,
420  1.894083312e-8,
421  -3.44358226e-9,
422  7.1119102e-10,
423  -1.6288744e-10,
424  4.065681e-11,
425  -1.091505e-11,
426  3.12005e-12,
427  -9.4202e-13,
428  2.9848e-13,
429  -9.872e-14,
430  3.394e-14,
431  -1.208e-14,
432  4.44e-15,
433  -1.68e-15,
434  6.5e-16,
435  -2.6e-16,
436  1.1e-16,
437  -4e-17,
438  2e-17,
439  -1e-17};
440 
441  const double c0 = 2 / xMath::Pi();
442 
443  int i;
444  double alfa, h, r, y, b0, b1, b2;
445  double v = fabs(x);
446 
447  v = fabs(x);
448  if (v < 8) {
449  y = v / 8;
450  h = 2 * y * y - 1;
451  alfa = h + h;
452  b0 = 0;
453  b1 = 0;
454  b2 = 0;
455  for (i = n1; i >= 0; --i) {
456  b0 = c1[i] + alfa * b1 - b2;
457  b2 = b1;
458  b1 = b0;
459  }
460  h = y * (b0 - h * b2);
461  } else {
462  r = 1 / v;
463  h = 128 * r * r - 1;
464  alfa = h + h;
465  b0 = 0;
466  b1 = 0;
467  b2 = 0;
468  for (i = n2; i >= 0; --i) {
469  b0 = c2[i] + alfa * b1 - b2;
470  b2 = b1;
471  b1 = b0;
472  }
473  h = xMath::BesselY0(v) + r * c0 * (b0 - h * b2);
474  }
475  if (x < 0) h = -h;
476  return h;
477 }
478 
479 //______________________________________________________________________________
480 double xMath::StruveH1(double x) {
481  // Struve Functions of Order 1
482  //
483  // Converted from CERNLIB M342 by Rene Brun.
484 
485  const int n3 = 16;
486  const int n4 = 22;
487  const double c3[17] = {.5578891446481605,
488  -.11188325726569816,
489  -.16337958125200939,
490  .32256932072405902,
491  -.14581632367244242,
492  .03292677399374035,
493  -.00460372142093573,
494  4.434706163314e-4,
495  -3.142099529341e-5,
496  1.7123719938e-6,
497  -7.416987005e-8,
498  2.61837671e-9,
499  -7.685839e-11,
500  1.9067e-12,
501  -4.052e-14,
502  7.5e-16,
503  -1e-17};
504  const double c4[23] = {1.00757647293865641,
505  .00750316051248257,
506  -7.043933264519e-5,
507  2.66205393382e-6,
508  -1.8841157753e-7,
509  1.949014958e-8,
510  -2.6126199e-9,
511  4.236269e-10,
512  -7.955156e-11,
513  1.679973e-11,
514  -3.9072e-12,
515  9.8543e-13,
516  -2.6636e-13,
517  7.645e-14,
518  -2.313e-14,
519  7.33e-15,
520  -2.42e-15,
521  8.3e-16,
522  -3e-16,
523  1.1e-16,
524  -4e-17,
525  2e-17,
526  -1e-17};
527 
528  const double c0 = 2 / xMath::Pi();
529  const double cc = 2 / (3 * xMath::Pi());
530 
531  int i, i1;
532  double alfa, h, r, y, b0, b1, b2;
533  double v = fabs(x);
534 
535  if (v == 0) {
536  h = 0;
537  } else if (v <= 0.3) {
538  y = v * v;
539  r = 1;
540  h = 1;
541  i1 = static_cast<int>((-8. / log10(v)));
542  for (i = 1; i <= i1; ++i) {
543  h = -h * y / ((2 * i + 1) * (2 * i + 3));
544  r += h;
545  }
546  h = cc * y * r;
547  } else if (v < 8) {
548  h = v * v / 32 - 1;
549  alfa = h + h;
550  b0 = 0;
551  b1 = 0;
552  b2 = 0;
553  for (i = n3; i >= 0; --i) {
554  b0 = c3[i] + alfa * b1 - b2;
555  b2 = b1;
556  b1 = b0;
557  }
558  h = b0 - h * b2;
559  } else {
560  h = 128 / (v * v) - 1;
561  alfa = h + h;
562  b0 = 0;
563  b1 = 0;
564  b2 = 0;
565  for (i = n4; i >= 0; --i) {
566  b0 = c4[i] + alfa * b1 - b2;
567  b2 = b1;
568  b1 = b0;
569  }
570  h = xMath::BesselY1(v) + c0 * (b0 - h * b2);
571  }
572  return h;
573 }
574 
575 
576 //______________________________________________________________________________
577 double xMath::StruveL0(double x) {
578  // Modified Struve Function of Order 0.
579  // By Kirill Filimonov.
580 
581  const double pi = xMath::Pi();
582 
583  double s = 1.0;
584  double r = 1.0;
585 
586  double a0, sl0, a1, bi0;
587 
588  int km, i;
589 
590  if (x <= 20.) {
591  a0 = 2.0 * x / pi;
592  for (i = 1; i <= 60; i++) {
593  r *= (x / (2 * i + 1)) * (x / (2 * i + 1));
594  s += r;
595  if (fabs(r / s) < 1.e-12) break;
596  }
597  sl0 = a0 * s;
598  } else {
599  km = int(5 * (x + 1.0));
600  if (x >= 50.0) km = 25;
601  for (i = 1; i <= km; i++) {
602  r *= (2 * i - 1) * (2 * i - 1) / x / x;
603  s += r;
604  if (fabs(r / s) < 1.0e-12) break;
605  }
606  a1 = exp(x) / sqrt(2 * pi * x);
607  r = 1.0;
608  bi0 = 1.0;
609  for (i = 1; i <= 16; i++) {
610  r = 0.125 * r * (2.0 * i - 1.0) * (2.0 * i - 1.0) / (i * x);
611  bi0 += r;
612  if (fabs(r / bi0) < 1.0e-12) break;
613  }
614 
615  bi0 = a1 * bi0;
616  sl0 = -2.0 / (pi * x) * s + bi0;
617  }
618  return sl0;
619 }
620 
621 //______________________________________________________________________________
622 double xMath::StruveL1(double x) {
623  // Modified Struve Function of Order 1.
624  // By Kirill Filimonov.
625 
626  const double pi = xMath::Pi();
627  double a1, sl1, bi1, s;
628  double r = 1.0;
629  int km, i;
630 
631  if (x <= 20.) {
632  s = 0.0;
633  for (i = 1; i <= 60; i++) {
634  r *= x * x / (4.0 * i * i - 1.0);
635  s += r;
636  if (fabs(r) < fabs(s) * 1.e-12) break;
637  }
638  sl1 = 2.0 / pi * s;
639  } else {
640  s = 1.0;
641  km = int(0.5 * x);
642  if (x > 50.0) km = 25;
643  for (i = 1; i <= km; i++) {
644  r *= (2 * i + 3) * (2 * i + 1) / x / x;
645  s += r;
646  if (fabs(r / s) < 1.0e-12) break;
647  }
648  sl1 = 2.0 / pi * (-1.0 + 1.0 / (x * x) + 3.0 * s / (x * x * x * x));
649  a1 = exp(x) / sqrt(2 * pi * x);
650  r = 1.0;
651  bi1 = 1.0;
652  for (i = 1; i <= 16; i++) {
653  r = -0.125 * r * (4.0 - (2.0 * i - 1.0) * (2.0 * i - 1.0)) / (i * x);
654  bi1 += r;
655  if (fabs(r / bi1) < 1.0e-12) break;
656  }
657  sl1 += a1 * bi1;
658  }
659  return sl1;
660 }
661 
662 
663 //______________________________________________________________________________
664 double xMath::BesselK0exp(double x) {
665  // Compute the modified Bessel function K_0(x) for positive real x.
666  //
667  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
668  // Applied Mathematics Series vol. 55 (1964), Washington.
669  //
670  //--- NvE 12-mar-2000 UU-SAP Utrecht
671 
672  // Parameters of the polynomial approximation
673  const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
674  p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
675 
676  const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
677  q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3,
678  q7 = 5.3208e-4;
679 
680  if (x <= 0) {
681  //Error("xMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
682  return 0;
683  }
684 
685  double y = 0, result = 0;
686 
687  if (x <= 2) {
688  y = x * x / 4;
689  result =
690  exp(x)
691  * ((-log(x / 2.) * xMath::BesselI0(x))
692  + (p1
693  + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
694  } else {
695  y = 2 / x;
696  result =
697  (1. / sqrt(x))
698  * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
699  }
700  return result;
701 }
702 
703 //______________________________________________________________________________
704 double xMath::BesselK1exp(double x) {
705  // Compute the modified Bessel function K_1(x) for positive real x.
706  //
707  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
708  // Applied Mathematics Series vol. 55 (1964), Washington.
709  //
710  //--- NvE 12-mar-2000 UU-SAP Utrecht
711 
712  // Parameters of the polynomial approximation
713  const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579, p4 = -0.18156897,
714  p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
715 
716  const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
717  q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3,
718  q7 = -6.8245e-4;
719 
720  if (x <= 0) {
721  //Error("xMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
722  return 0;
723  }
724 
725  double y = 0, result = 0;
726 
727  if (x <= 2) {
728  y = x * x / 4;
729  result =
730  exp(x)
731  * ((log(x / 2.) * xMath::BesselI1(x))
732  + (1. / x)
733  * (p1
734  + y
735  * (p2
736  + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
737  } else {
738  y = 2 / x;
739  result =
740  (1. / sqrt(x))
741  * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
742  }
743  return result;
744 }
745 
746 //______________________________________________________________________________
747 double xMath::BesselKexp(int n, double x) {
748  // Compute the Integer Order Modified Bessel function K_n(x)
749  // for n=0,1,2,... and positive real x.
750  //
751  //--- NvE 12-mar-2000 UU-SAP Utrecht
752 
753  if (x <= 0 || n < 0) {
754  //Error("xMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
755  return 0;
756  }
757 
758  if (n == 0) return xMath::BesselK0exp(x);
759  if (n == 1) return xMath::BesselK1exp(x);
760 
761  // Perform upward recurrence for all x
762  double tox = 2 / x;
763  double bkm = xMath::BesselK0exp(x);
764  double bk = xMath::BesselK1exp(x);
765  double bkp = 0;
766  for (int j = 1; j < n; j++) {
767  bkp = bkm + double(j) * tox * bk;
768  bkm = bk;
769  bk = bkp;
770  }
771  return bk;
772 }
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:136
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
xMath::StruveH1
double StruveH1(double x)
Definition: xMath.cxx:480
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
xMath::StruveL1
double StruveL1(double x)
Definition: xMath.cxx:622
xMath::BesselY0
double BesselY0(double x)
Definition: xMath.cxx:323
xMath::BesselI1
double BesselI1(double x)
Definition: xMath.cxx:90
xMath::BesselI
double BesselI(int n, double x)
Definition: xMath.cxx:208
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
xMath::BesselK0
double BesselK0(double x)
Definition: xMath.cxx:52
xMath::BesselJ0
double BesselJ0(double x)
Definition: xMath.cxx:254
xMath::Pi
double Pi()
Definition: xMath.h:5
xMath::BesselJ1
double BesselJ1(double x)
Definition: xMath.cxx:288
xMath::StruveH0
double StruveH0(double x)
Definition: xMath.cxx:391
xMath::BesselY1
double BesselY1(double x)
Definition: xMath.cxx:356
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
xMath::BesselK1
double BesselK1(double x)
Definition: xMath.cxx:140
xMath::BesselK1exp
double BesselK1exp(double x)
Definition: xMath.cxx:704
xMath::BesselI0
double BesselI0(double x)
Definition: xMath.cxx:7
xMath::BesselK
double BesselK(int n, double x)
Definition: xMath.cxx:180
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
xMath::BesselKexp
double BesselKexp(int n, double x)
Definition: xMath.cxx:747
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
xMath.h
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:137
xMath::BesselK0exp
double BesselK0exp(double x)
Definition: xMath.cxx:664
xMath::StruveL0
double StruveL0(double x)
Definition: xMath.cxx:577