Thermal-FIST  1.3
Package for hadron resonance gas model applications
xMath.cpp
Go to the documentation of this file.
1 /*
2  * Thermal-FIST package
3  *
4  * Copyright (c) 2014-2018 Volodymyr Vovchenko
5  *
6  * GNU General Public License (GPLv3 or later)
7  */
8 #include "HRGBase/xMath.h"
9 
10 #include <cmath>
11 #include <sstream>
12 #include <stdexcept>
13 #include <cfloat>
14 
15 using namespace std;
16 
17 namespace thermalfist {
18 
19  // Implementation of the special functions below adapted from
20  // CERN-ROOT package: https://root.cern.ch/
21 
22  //______________________________________________________________________________
23  double xMath::BesselI0(double x)
24  {
25  // Compute the modified Bessel function I_0(x) for any real x.
26  //
27  //--- NvE 12-mar-2000 UU-SAP Utrecht
28 
29  // Parameters of the polynomial approximation
30  const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424,
31  p4 = 1.2067492, p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
32 
33  const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
34  q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
35  q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
36 
37  const double k1 = 3.75;
38  double ax = fabs(x);
39 
40  double y = 0, result = 0;
41 
42  if (ax < k1) {
43  double xx = x / k1;
44  y = xx * xx;
45  result = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
46  }
47  else {
48  y = k1 / ax;
49  result = (exp(ax) / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
50  }
51  return result;
52  }
53 
54  //______________________________________________________________________________
55  double xMath::BesselK0(double x)
56  {
57  // Compute the modified Bessel function K_0(x) for positive real x.
58  //
59  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
60  // Applied Mathematics Series vol. 55 (1964), Washington.
61  //
62  //--- NvE 12-mar-2000 UU-SAP Utrecht
63 
64  // Parameters of the polynomial approximation
65  const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
66  p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
67 
68  const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
69  q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3, q7 = 5.3208e-4;
70 
71  if (x <= 0) {
72  //Error("xMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
73  return 0;
74  }
75 
76  double y = 0, result = 0;
77 
78  if (x <= 2) {
79  y = x * x / 4;
80  result = (-log(x / 2.)*xMath::BesselI0(x)) + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
81  }
82  else {
83  y = 2 / x;
84  result = (exp(-x) / sqrt(x))*(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  {
92  // Compute the modified Bessel function I_1(x) for any real x.
93  //
94  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
95  // Applied Mathematics Series vol. 55 (1964), Washington.
96  //
97  //--- NvE 12-mar-2000 UU-SAP Utrecht
98 
99  // Parameters of the polynomial approximation
100  const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869,
101  p4 = 0.15084934, p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
102 
103  const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
104  q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
105  q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
106 
107  const double k1 = 3.75;
108  double ax = fabs(x);
109 
110  double y = 0, result = 0;
111 
112  if (ax < k1) {
113  double xx = x / k1;
114  y = xx * xx;
115  result = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
116  }
117  else {
118  y = k1 / ax;
119  result = (exp(ax) / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
120  if (x < 0) result = -result;
121  }
122  return result;
123  }
124 
125  //______________________________________________________________________________
126  double xMath::BesselK1(double x)
127  {
128  // Compute the modified Bessel function K_1(x) for positive real x.
129  //
130  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
131  // Applied Mathematics Series vol. 55 (1964), Washington.
132  //
133  //--- NvE 12-mar-2000 UU-SAP Utrecht
134 
135  // Parameters of the polynomial approximation
136  const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579,
137  p4 = -0.18156897, p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
138 
139  const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
140  q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3, q7 = -6.8245e-4;
141 
142  if (x <= 0) {
143  //Error("xMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
144  return 0;
145  }
146 
147  double y = 0, result = 0;
148 
149  if (x <= 2) {
150  y = x * x / 4;
151  result = (log(x / 2.)*xMath::BesselI1(x)) + (1. / x)*(p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
152  }
153  else {
154  y = 2 / x;
155  result = (exp(-x) / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
156  }
157  return result;
158  }
159 
160  //______________________________________________________________________________
161  double xMath::BesselK(int n, double x)
162  {
163  if (n < 0) n = -n;
164  // Compute the Integer Order Modified Bessel function K_n(x)
165  // for n=0,1,2,... and positive real x.
166  //
167  //--- NvE 12-mar-2000 UU-SAP Utrecht
168 
169  if (x <= 0 || n < 0) {
170  //Error("xMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
171  return 0;
172  }
173 
174  if (n == 0) return xMath::BesselK0(x);
175  if (n == 1) return xMath::BesselK1(x);
176 
177  // Perform upward recurrence for all x
178  double tox = 2 / x;
179  double bkm = xMath::BesselK0(x);
180  double bk = xMath::BesselK1(x);
181  double bkp = 0;
182  for (int j = 1; j < n; j++) {
183  bkp = bkm + double(j)*tox*bk;
184  bkm = bk;
185  bk = bkp;
186  }
187  return bk;
188  }
189 
190  //______________________________________________________________________________
191  double xMath::BesselI(int n, double x)
192  {
193  // Compute the Integer Order Modified Bessel function I_n(x)
194  // for n=0,1,2,... and any real x.
195  //
196  //--- NvE 12-mar-2000 UU-SAP Utrecht
197  if (n < 0) n = -n;
198 
199  int iacc = 40; // Increase to enhance accuracy
200  const double kBigPositive = 1.e10;
201  const double kBigNegative = 1.e-10;
202 
203  if (n < 0) {
204  //Error("xMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
205  return 0;
206  }
207 
208  if (n == 0) return xMath::BesselI0(x);
209  if (n == 1) return xMath::BesselI1(x);
210 
211  if (x == 0) return 0;
212  if (fabs(x) > kBigPositive) return 0;
213 
214  double tox = 2 / fabs(x);
215  double bip = 0, bim = 0;
216  double bi = 1;
217  double result = 0;
218  int m = 2 * ((n + int(sqrt(float(iacc*n)))));
219  for (int j = m; j >= 1; j--) {
220  bim = bip + double(j)*tox*bi;
221  bip = bi;
222  bi = bim;
223  // Renormalise to prevent overflows
224  if (fabs(bi) > kBigPositive) {
225  result *= kBigNegative;
226  bi *= kBigNegative;
227  bip *= kBigNegative;
228  }
229  if (j == n) result = bip;
230  }
231 
232  result *= xMath::BesselI0(x) / bi; // Normalise with BesselI0(x)
233  if ((x < 0) && (n % 2 == 1)) result = -result;
234 
235  return result;
236  }
237 
238  //______________________________________________________________________________
239  double xMath::BesselJ0(double x)
240  {
241  // Returns the Bessel function J0(x) for any real x.
242 
243  double ax, z;
244  double xx, y, result, result1, result2;
245  const double p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
246  const double p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
247  const double p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
248  const double p10 = 59272.64853, p11 = 267.8532712;
249 
250  const double q1 = 0.785398164;
251  const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
252  const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
253  const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
254  const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
255  const double q10 = 0.934935152e-7, q11 = 0.636619772;
256 
257  if ((ax = fabs(x)) < 8) {
258  y = x * x;
259  result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
260  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
261  result = result1 / result2;
262  }
263  else {
264  z = 8 / ax;
265  y = z * z;
266  xx = ax - q1;
267  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
268  result2 = q6 + y * (q7 + y * (q8 + y * (q9 - y * q10)));
269  result = sqrt(q11 / ax)*(cos(xx)*result1 - z * sin(xx)*result2);
270  }
271  return result;
272  }
273 
274  //______________________________________________________________________________
275  double xMath::BesselJ1(double x)
276  {
277  // Returns the Bessel function J1(x) for any real x.
278 
279  double ax, z;
280  double xx, y, result, result1, result2;
281  const double p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
282  const double p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
283  const double p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
284  const double p10 = 99447.43394, p11 = 376.9991397;
285 
286  const double q1 = 2.356194491;
287  const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
288  const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
289  const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
290  const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
291  const double q10 = 0.105787412e-6, q11 = 0.636619772;
292 
293  if ((ax = fabs(x)) < 8) {
294  y = x * x;
295  result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
296  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
297  result = result1 / result2;
298  }
299  else {
300  z = 8 / ax;
301  y = z * z;
302  xx = ax - q1;
303  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
304  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
305  result = sqrt(q11 / ax)*(cos(xx)*result1 - z * sin(xx)*result2);
306  if (x < 0) result = -result;
307  }
308  return result;
309  }
310 
311  //______________________________________________________________________________
312  double xMath::BesselY0(double x)
313  {
314  // Returns the Bessel function Y0(x) for positive x.
315 
316  double z, xx, y, result, result1, result2;
317  const double p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
318  const double p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
319  const double p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
320  const double p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
321 
322  const double q1 = 0.785398164;
323  const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
324  const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
325  const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
326  const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
327  const double q10 = -0.934945152e-7, q11 = 0.636619772;
328 
329  if (x < 8) {
330  y = x * x;
331  result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
332  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
333  result = (result1 / result2) + p12 * xMath::BesselJ0(x)*log(x);
334  }
335  else {
336  z = 8 / x;
337  y = z * z;
338  xx = x - q1;
339  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
340  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
341  result = sqrt(q11 / x)*(sin(xx)*result1 + z * cos(xx)*result2);
342  }
343  return result;
344  }
345 
346  //______________________________________________________________________________
347  double xMath::BesselY1(double x)
348  {
349  // Returns the Bessel function Y1(x) for positive x.
350 
351  double z, xx, y, result, result1, result2;
352  const double p1 = -0.4900604943e13, p2 = 0.1275274390e13;
353  const double p3 = -0.5153438139e11, p4 = 0.7349264551e9;
354  const double p5 = -0.4237922726e7, p6 = 0.8511937935e4;
355  const double p7 = 0.2499580570e14, p8 = 0.4244419664e12;
356  const double p9 = 0.3733650367e10, p10 = 0.2245904002e8;
357  const double p11 = 0.1020426050e6, p12 = 0.3549632885e3;
358  const double p13 = 0.636619772;
359  const double q1 = 2.356194491;
360  const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
361  const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
362  const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
363  const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
364  const double q10 = 0.105787412e-6, q11 = 0.636619772;
365 
366  if (x < 8) {
367  y = x * x;
368  result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
369  result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y * (p12 + y)))));
370  result = (result1 / result2) + p13 * (xMath::BesselJ1(x)*log(x) - 1 / x);
371  }
372  else {
373  z = 8 / x;
374  y = z * z;
375  xx = x - q1;
376  result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
377  result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
378  result = sqrt(q11 / x)*(sin(xx)*result1 + z * cos(xx)*result2);
379  }
380  return result;
381  }
382 
383  //______________________________________________________________________________
384  double xMath::StruveH0(double x)
385  {
386  // Struve Functions of Order 0
387  //
388  // Converted from CERNLIB M342 by Rene Brun.
389 
390  const int n1 = 15;
391  const int n2 = 25;
392  const double c1[16] = { 1.00215845609911981, -1.63969292681309147,
393  1.50236939618292819, -.72485115302121872,
394  .18955327371093136, -.03067052022988,
395  .00337561447375194, -2.6965014312602e-4,
396  1.637461692612e-5, -7.8244408508e-7,
397  3.021593188e-8, -9.6326645e-10,
398  2.579337e-11, -5.8854e-13,
399  1.158e-14, -2e-16 };
400  const double c2[26] = { .99283727576423943, -.00696891281138625,
401  1.8205103787037e-4, -1.063258252844e-5,
402  9.8198294287e-7, -1.2250645445e-7,
403  1.894083312e-8, -3.44358226e-9,
404  7.1119102e-10, -1.6288744e-10,
405  4.065681e-11, -1.091505e-11,
406  3.12005e-12, -9.4202e-13,
407  2.9848e-13, -9.872e-14,
408  3.394e-14, -1.208e-14,
409  4.44e-15, -1.68e-15,
410  6.5e-16, -2.6e-16,
411  1.1e-16, -4e-17,
412  2e-17, -1e-17 };
413 
414  const double c0 = 2 / xMath::Pi();
415 
416  int i;
417  double alfa, h, r, y, b0, b1, b2;
418  double v = fabs(x);
419 
420  v = fabs(x);
421  if (v < 8) {
422  y = v / 8;
423  h = 2 * y*y - 1;
424  alfa = h + h;
425  b0 = 0;
426  b1 = 0;
427  b2 = 0;
428  for (i = n1; i >= 0; --i) {
429  b0 = c1[i] + alfa * b1 - b2;
430  b2 = b1;
431  b1 = b0;
432  }
433  h = y * (b0 - h * b2);
434  }
435  else {
436  r = 1 / v;
437  h = 128 * r*r - 1;
438  alfa = h + h;
439  b0 = 0;
440  b1 = 0;
441  b2 = 0;
442  for (i = n2; i >= 0; --i) {
443  b0 = c2[i] + alfa * b1 - b2;
444  b2 = b1;
445  b1 = b0;
446  }
447  h = xMath::BesselY0(v) + r * c0*(b0 - h * b2);
448  }
449  if (x < 0) h = -h;
450  return h;
451  }
452 
453  //______________________________________________________________________________
454  double xMath::StruveH1(double x)
455  {
456  // Struve Functions of Order 1
457  //
458  // Converted from CERNLIB M342 by Rene Brun.
459 
460  const int n3 = 16;
461  const int n4 = 22;
462  const double c3[17] = { .5578891446481605, -.11188325726569816,
463  -.16337958125200939, .32256932072405902,
464  -.14581632367244242, .03292677399374035,
465  -.00460372142093573, 4.434706163314e-4,
466  -3.142099529341e-5, 1.7123719938e-6,
467  -7.416987005e-8, 2.61837671e-9,
468  -7.685839e-11, 1.9067e-12,
469  -4.052e-14, 7.5e-16,
470  -1e-17 };
471  const double c4[23] = { 1.00757647293865641, .00750316051248257,
472  -7.043933264519e-5, 2.66205393382e-6,
473  -1.8841157753e-7, 1.949014958e-8,
474  -2.6126199e-9, 4.236269e-10,
475  -7.955156e-11, 1.679973e-11,
476  -3.9072e-12, 9.8543e-13,
477  -2.6636e-13, 7.645e-14,
478  -2.313e-14, 7.33e-15,
479  -2.42e-15, 8.3e-16,
480  -3e-16, 1.1e-16,
481  -4e-17, 2e-17,-1e-17 };
482 
483  const double c0 = 2 / xMath::Pi();
484  const double cc = 2 / (3 * xMath::Pi());
485 
486  int i, i1;
487  double alfa, h, r, y, b0, b1, b2;
488  double v = fabs(x);
489 
490  if (v == 0) {
491  h = 0;
492  }
493  else if (v <= 0.3) {
494  y = v * v;
495  r = 1;
496  h = 1;
497  i1 = (int)(-8. / log10(v));
498  for (i = 1; i <= i1; ++i) {
499  h = -h * y / ((2 * i + 1)*(2 * i + 3));
500  r += h;
501  }
502  h = cc * y*r;
503  }
504  else if (v < 8) {
505  h = v * v / 32 - 1;
506  alfa = h + h;
507  b0 = 0;
508  b1 = 0;
509  b2 = 0;
510  for (i = n3; i >= 0; --i) {
511  b0 = c3[i] + alfa * b1 - b2;
512  b2 = b1;
513  b1 = b0;
514  }
515  h = b0 - h * b2;
516  }
517  else {
518  h = 128 / (v*v) - 1;
519  alfa = h + h;
520  b0 = 0;
521  b1 = 0;
522  b2 = 0;
523  for (i = n4; i >= 0; --i) {
524  b0 = c4[i] + alfa * b1 - b2;
525  b2 = b1;
526  b1 = b0;
527  }
528  h = xMath::BesselY1(v) + c0 * (b0 - h * b2);
529  }
530  return h;
531  }
532 
533 
534  //______________________________________________________________________________
535  double xMath::StruveL0(double x)
536  {
537  // Modified Struve Function of Order 0.
538  // By Kirill Filimonov.
539 
540  const double pi = xMath::Pi();
541 
542  double s = 1.0;
543  double r = 1.0;
544 
545  double a0, sl0, a1, bi0;
546 
547  int km, i;
548 
549  if (x <= 20.) {
550  a0 = 2.0*x / pi;
551  for (i = 1; i <= 60; i++) {
552  r *= (x / (2 * i + 1))*(x / (2 * i + 1));
553  s += r;
554  if (fabs(r / s) < 1.e-12) break;
555  }
556  sl0 = a0 * s;
557  }
558  else {
559  km = int(5 * (x + 1.0));
560  if (x >= 50.0)km = 25;
561  for (i = 1; i <= km; i++) {
562  r *= (2 * i - 1)*(2 * i - 1) / x / x;
563  s += r;
564  if (fabs(r / s) < 1.0e-12) break;
565  }
566  a1 = exp(x) / sqrt(2 * pi*x);
567  r = 1.0;
568  bi0 = 1.0;
569  for (i = 1; i <= 16; i++) {
570  r = 0.125*r*(2.0*i - 1.0)*(2.0*i - 1.0) / (i*x);
571  bi0 += r;
572  if (fabs(r / bi0) < 1.0e-12) break;
573  }
574 
575  bi0 = a1 * bi0;
576  sl0 = -2.0 / (pi*x)*s + bi0;
577  }
578  return sl0;
579  }
580 
581  //______________________________________________________________________________
582  double xMath::StruveL1(double x)
583  {
584  // Modified Struve Function of Order 1.
585  // By Kirill Filimonov.
586 
587  const double pi = xMath::Pi();
588  double a1, sl1, bi1, s;
589  double r = 1.0;
590  int km, i;
591 
592  if (x <= 20.) {
593  s = 0.0;
594  for (i = 1; i <= 60; i++) {
595  r *= x * x / (4.0*i*i - 1.0);
596  s += r;
597  if (fabs(r) < fabs(s)*1.e-12) break;
598  }
599  sl1 = 2.0 / pi * s;
600  }
601  else {
602  s = 1.0;
603  km = int(0.5*x);
604  if (x > 50.0)km = 25;
605  for (i = 1; i <= km; i++) {
606  r *= (2 * i + 3)*(2 * i + 1) / x / x;
607  s += r;
608  if (fabs(r / s) < 1.0e-12) break;
609  }
610  sl1 = 2.0 / pi * (-1.0 + 1.0 / (x*x) + 3.0*s / (x*x*x*x));
611  a1 = exp(x) / sqrt(2 * pi*x);
612  r = 1.0;
613  bi1 = 1.0;
614  for (i = 1; i <= 16; i++) {
615  r = -0.125*r*(4.0 - (2.0*i - 1.0)*(2.0*i - 1.0)) / (i*x);
616  bi1 += r;
617  if (fabs(r / bi1) < 1.0e-12) break;
618  }
619  sl1 += a1 * bi1;
620  }
621  return sl1;
622  }
623 
624 
625 
626  //______________________________________________________________________________
627  double xMath::BesselK0exp(double x)
628  {
629  // Compute the modified Bessel function K_0(x) for positive real x.
630  //
631  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
632  // Applied Mathematics Series vol. 55 (1964), Washington.
633  //
634  //--- NvE 12-mar-2000 UU-SAP Utrecht
635 
636  // Parameters of the polynomial approximation
637  const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
638  p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
639 
640  const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
641  q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3, q7 = 5.3208e-4;
642 
643  if (x <= 0) {
644  //Error("xMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
645  return 0;
646  }
647 
648  double y = 0, result = 0;
649 
650  if (x <= 2) {
651  y = x * x / 4;
652  result = exp(x)*((-log(x / 2.)*xMath::BesselI0(x)) + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
653  }
654  else {
655  y = 2 / x;
656  result = (1. / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
657  }
658  return result;
659  }
660 
661  //______________________________________________________________________________
662  double xMath::BesselK1exp(double x)
663  {
664  // Compute the modified Bessel function K_1(x) for positive real x.
665  //
666  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
667  // Applied Mathematics Series vol. 55 (1964), Washington.
668  //
669  //--- NvE 12-mar-2000 UU-SAP Utrecht
670 
671  // Parameters of the polynomial approximation
672  const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579,
673  p4 = -0.18156897, p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
674 
675  const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
676  q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3, q7 = -6.8245e-4;
677 
678  if (x <= 0) {
679  //Error("xMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
680  return 0;
681  }
682 
683  double y = 0, result = 0;
684 
685  if (x <= 2) {
686  y = x * x / 4;
687  result = exp(x) * ((log(x / 2.)*xMath::BesselI1(x)) + (1. / x)*(p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
688  }
689  else {
690  y = 2 / x;
691  result = (1. / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
692  }
693  return result;
694  }
695 
696  //______________________________________________________________________________
697  double xMath::BesselKexp(int n, double x)
698  {
699  if (n < 0) n = -n;
700  // Compute the Integer Order Modified Bessel function K_n(x)
701  // for n=0,1,2,... and positive real x.
702  //
703  //--- NvE 12-mar-2000 UU-SAP Utrecht
704 
705  if (x <= 0 || n < 0) {
706  //Error("xMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
707  return 0;
708  }
709 
710  if (n == 0) return xMath::BesselK0exp(x);
711  if (n == 1) return xMath::BesselK1exp(x);
712 
713  // Perform upward recurrence for all x
714  double tox = 2 / x;
715  double bkm = xMath::BesselK0exp(x);
716  double bk = xMath::BesselK1exp(x);
717  double bkp = 0;
718  for (int j = 1; j < n; j++) {
719  bkp = bkm + double(j)*tox*bk;
720  bkm = bk;
721  bk = bkp;
722  }
723  return bk;
724  }
725 
726  //______________________________________________________________________________
727  double xMath::BesselI0exp(double x)
728  {
729  // Compute the modified Bessel function I_0(x) for any real x.
730  //
731  //--- NvE 12-mar-2000 UU-SAP Utrecht
732 
733  // Parameters of the polynomial approximation
734  const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424,
735  p4 = 1.2067492, p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
736 
737  const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
738  q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
739  q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
740 
741  const double k1 = 3.75;
742  double ax = fabs(x);
743 
744  double y = 0, result = 0;
745 
746  if (ax < k1) {
747  double xx = x / k1;
748  y = xx * xx;
749  result = exp(-ax) * ( p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))) );
750  }
751  else {
752  y = k1 / ax;
753  result = (1. / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
754  }
755  return result;
756  }
757 
758  //______________________________________________________________________________
759  double xMath::BesselI1exp(double x)
760  {
761  // Compute the modified Bessel function I_1(x) for any real x.
762  //
763  // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
764  // Applied Mathematics Series vol. 55 (1964), Washington.
765  //
766  //--- NvE 12-mar-2000 UU-SAP Utrecht
767 
768  // Parameters of the polynomial approximation
769  const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869,
770  p4 = 0.15084934, p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
771 
772  const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
773  q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
774  q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
775 
776  const double k1 = 3.75;
777  double ax = fabs(x);
778 
779  double y = 0, result = 0;
780 
781  if (ax < k1) {
782  double xx = x / k1;
783  y = xx * xx;
784  result = exp(-ax) * (x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
785  }
786  else {
787  y = k1 / ax;
788  result = (1. / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
789  if (x < 0) result = -result;
790  }
791  return result;
792  }
793 
794  //______________________________________________________________________________
795  double xMath::BesselIexp(int n, double x)
796  {
797  // Compute the Integer Order Modified Bessel function I_n(x)
798  // for n=0,1,2,... and any real x.
799  //
800  //--- NvE 12-mar-2000 UU-SAP Utrecht
801  if (n < 0) n = -n;
802 
803  int iacc = 40; // Increase to enhance accuracy
804  const double kBigPositive = 1.e10;
805  const double kBigNegative = 1.e-10;
806 
807  if (n < 0) {
808  //Error("xMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
809  return 0;
810  }
811 
812  if (n == 0) return xMath::BesselI0exp(x);
813  if (n == 1) return xMath::BesselI1exp(x);
814 
815  if (x == 0) return 0;
816  if (fabs(x) > kBigPositive) return 0;
817 
818  double tox = 2 / fabs(x);
819  double bip = 0, bim = 0;
820  double bi = 1;
821  double result = 0;
822  int m = 2 * ((n + int(sqrt(float(iacc*n)))));
823  for (int j = m; j >= 1; j--) {
824  bim = bip + double(j)*tox*bi;
825  bip = bi;
826  bi = bim;
827  // Renormalise to prevent overflows
828  if (fabs(bi) > kBigPositive) {
829  result *= kBigNegative;
830  bi *= kBigNegative;
831  bip *= kBigNegative;
832  }
833  if (j == n) result = bip;
834  }
835 
836  result *= xMath::BesselI0exp(x) / bi; // Normalise with BesselI0exp(x)
837  if ((x < 0) && (n % 2 == 1)) result = -result;
838 
839  return result;
840  }
841 
842 
843  double xMath::Gamma
844  (
845  double x // We require x > 0
846  )
847  {
848  if (x <= 0.0)
849  {
850  std::stringstream os;
851  os << "Invalid input argument " << x << ". Argument must be positive.";
852  throw std::invalid_argument(os.str());
853  }
854 
855  // Split the function domain into three intervals:
856  // (0, 0.001), [0.001, 12), and (12, infinity)
857 
859  // First interval: (0, 0.001)
860  //
861  // For small x, 1/Gamma(x) has power series x + gamma x^2 - ...
862  // So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
863  // The relative error over this interval is less than 6e-7.
864 
865  const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
866 
867  if (x < 0.001)
868  return 1.0 / (x*(1.0 + gamma * x));
869 
871  // Second interval: [0.001, 12)
872 
873  if (x < 12.0)
874  {
875  // The algorithm directly approximates gamma over (1,2) and uses
876  // reduction identities to reduce other arguments to this interval.
877 
878  double y = x;
879  int n = 0;
880  bool arg_was_less_than_one = (y < 1.0);
881 
882  // Add or subtract integers as necessary to bring y into (1,2)
883  // Will correct for this below
884  if (arg_was_less_than_one)
885  {
886  y += 1.0;
887  }
888  else
889  {
890  n = static_cast<int> (floor(y)) - 1; // will use n later
891  y -= n;
892  }
893 
894  // numerator coefficients for approximation over the interval (1,2)
895  static const double p[] =
896  {
897  -1.71618513886549492533811E+0,
898  2.47656508055759199108314E+1,
899  -3.79804256470945635097577E+2,
900  6.29331155312818442661052E+2,
901  8.66966202790413211295064E+2,
902  -3.14512729688483675254357E+4,
903  -3.61444134186911729807069E+4,
904  6.64561438202405440627855E+4
905  };
906 
907  // denominator coefficients for approximation over the interval (1,2)
908  static const double q[] =
909  {
910  -3.08402300119738975254353E+1,
911  3.15350626979604161529144E+2,
912  -1.01515636749021914166146E+3,
913  -3.10777167157231109440444E+3,
914  2.25381184209801510330112E+4,
915  4.75584627752788110767815E+3,
916  -1.34659959864969306392456E+5,
917  -1.15132259675553483497211E+5
918  };
919 
920  double num = 0.0;
921  double den = 1.0;
922  int i;
923 
924  double z = y - 1;
925  for (i = 0; i < 8; i++)
926  {
927  num = (num + p[i])*z;
928  den = den * z + q[i];
929  }
930  double result = num / den + 1.0;
931 
932  // Apply correction if argument was not initially in (1,2)
933  if (arg_was_less_than_one)
934  {
935  // Use identity gamma(z) = gamma(z+1)/z
936  // The variable "result" now holds gamma of the original y + 1
937  // Thus we use y-1 to get back the orginal y.
938  result /= (y - 1.0);
939  }
940  else
941  {
942  // Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
943  for (i = 0; i < n; i++)
944  result *= y++;
945  }
946 
947  return result;
948  }
949 
951  // Third interval: [12, infinity)
952 
953  if (x > 171.624)
954  {
955  // Correct answer too large to display. Force +infinity.
956  double temp = DBL_MAX;
957  return temp * 2.0;
958  }
959 
960  return exp(LogGamma(x));
961  }
962 
963  double xMath::LogGamma
964  (
965  double x // x must be positive
966  )
967  {
968  if (x <= 0.0)
969  {
970  std::stringstream os;
971  os << "Invalid input argument " << x << ". Argument must be positive.";
972  throw std::invalid_argument(os.str());
973  }
974 
975  if (x < 12.0)
976  {
977  return log(fabs(Gamma(x)));
978  }
979 
980  // Abramowitz and Stegun 6.1.41
981  // Asymptotic series should be good to at least 11 or 12 figures
982  // For error analysis, see Whittiker and Watson
983  // A Course in Modern Analysis (1927), page 252
984 
985  static const double c[8] =
986  {
987  1.0 / 12.0,
988  -1.0 / 360.0,
989  1.0 / 1260.0,
990  -1.0 / 1680.0,
991  1.0 / 1188.0,
992  -691.0 / 360360.0,
993  1.0 / 156.0,
994  -3617.0 / 122400.0
995  };
996  double z = 1.0 / (x*x);
997  double sum = c[7];
998  for (int i = 6; i >= 0; i--)
999  {
1000  sum *= z;
1001  sum += c[i];
1002  }
1003  double series = sum / x;
1004 
1005  static const double halfLogTwoPi = 0.91893853320467274178032973640562;
1006  double logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;
1007  return logGamma;
1008  }
1009 
1010 } // namespace thermalfist
thermalfist::xMath::StruveH0
double StruveH0(double x)
Struve functions of order 0.
Definition: xMath.cpp:384
thermalfist::xMath::BesselKexp
double BesselKexp(int n, double x)
Definition: xMath.cpp:697
thermalfist::xMath::Gamma
double Gamma(double)
Definition: xMath.cpp:844
thermalfist::xMath::BesselI1exp
double BesselI1exp(double x)
Definition: xMath.cpp:759
thermalfist::xMath::BesselI0exp
double BesselI0exp(double x)
Definition: xMath.cpp:727
thermalfist::xMath::Pi
double Pi()
Pi constant.
Definition: xMath.h:24
thermalfist::xMath::BesselJ0
double BesselJ0(double x)
Bessel function J0(x) for any real x.
Definition: xMath.cpp:239
thermalfist::xMath::StruveL0
double StruveL0(double x)
Modified Struve functions of order 0.
Definition: xMath.cpp:535
thermalfist::xMath::BesselJ1
double BesselJ1(double x)
Bessel function J1(x) for any real x.
Definition: xMath.cpp:275
thermalfist::xMath::LogGamma
double LogGamma(double)
Definition: xMath.cpp:964
thermalfist::xMath::BesselI1
double BesselI1(double x)
modified Bessel function I_1(x)
Definition: xMath.cpp:90
thermalfist::xMath::BesselIexp
double BesselIexp(int n, double x)
Definition: xMath.cpp:795
xMath.h
Contains some extra mathematical functions used in the code.
thermalfist::xMath::BesselI
double BesselI(int n, double x)
integer order modified Bessel function I_n(x)
Definition: xMath.cpp:191
thermalfist::xMath::BesselK
double BesselK(int n, double x)
integer order modified Bessel function K_n(x)
Definition: xMath.cpp:161
thermalfist::xMath::StruveL1
double StruveL1(double x)
Modified Struve functions of order 1.
Definition: xMath.cpp:582
thermalfist::xMath::BesselK1
double BesselK1(double x)
modified Bessel function K_1(x)
Definition: xMath.cpp:126
thermalfist::xMath::StruveH1
double StruveH1(double x)
Struve functions of order 1.
Definition: xMath.cpp:454
thermalfist::xMath::BesselK0
double BesselK0(double x)
modified Bessel function K_0(x)
Definition: xMath.cpp:55
thermalfist::xMath::BesselY0
double BesselY0(double x)
Bessel function Y0(x) for positive x.
Definition: xMath.cpp:312
thermalfist::xMath::BesselI0
double BesselI0(double x)
modified Bessel function I_0(x)
Definition: xMath.cpp:23
thermalfist::xMath::BesselY1
double BesselY1(double x)
Bessel function Y1(x) for positive x.
Definition: xMath.cpp:347
thermalfist::xMath::BesselK0exp
double BesselK0exp(double x)
Definition: xMath.cpp:627
thermalfist::xMath::BesselK1exp
double BesselK1exp(double x)
Definition: xMath.cpp:662
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition: BilinearSplineFunction.h:13
Generated by   doxygen 1.8.11