incGamma.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2019-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Global
25  Foam::incGamma
26 
27 Description
28  Calculates the upper and lower incomplete gamma functions as well as their
29  normalised versions.
30 
31  The algorithm is described in detail in DiDonato et al. (1986).
32 
33  \verbatim
34  DiDonato, A. R., & Morris Jr, A. H. (1986).
35  Computation of the incomplete gamma function ratios and their inverse.
36  ACM Transactions on Mathematical Software (TOMS), 12(4), 377-393.
37  \endverbatim
38 
39  All equation numbers in the following code refer to the above paper.
40  The algorithm in function 'incGammaRatio_Q' is described in section 3.
41  The accuracy parameter IND is set to a value of 1.
42 
43 \*---------------------------------------------------------------------------*/
44 
45 #include "mathematicalConstants.H"
46 
47 using namespace Foam::constant::mathematical;
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 namespace Foam
52 {
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 
56 // Eqn. (13)
57 static scalar calcQE11(const scalar a, const scalar x, const int e = 30)
58 {
59  scalar a_2n = 0;
60  scalar b_2n = 1;
61 
62  scalar a_2np1 = 1;
63  scalar b_2np1 = x;
64 
65  int n = 1;
66  for (n = 1; (2*n) <= e; n++)
67  {
68  const scalar a_2nm1 = a_2np1;
69  const scalar b_2nm1 = b_2np1;
70 
71  a_2n = a_2nm1 + (n - a)*a_2n;
72  b_2n = b_2nm1 + (n - a)*b_2n;
73 
74  a_2np1 = x*a_2n + n*a_2nm1;
75  b_2np1 = x*b_2n + n*b_2nm1;
76  }
77 
78  if (2*(n - 1) < e)
79  {
80  return a_2np1/b_2np1;
81  }
82  else
83  {
84  return a_2n/b_2n;
85  }
86 }
87 
88 
89 // Eqn. (15)
90 static scalar calcPE15(const scalar a, const scalar x, const int nmax = 20)
91 {
92  scalar prod = 1;
93  scalar sum = 0;
94 
95  for (int n = 1; n <= nmax; n++)
96  {
97  prod *= (a + n);
98  sum += pow(x, n)/prod;
99  }
100 
101  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
102 
103  return R/a*(1 + sum);
104 }
106 
107 // Eq. (16)
108 static scalar calcQE16(const scalar a, const scalar x, const int N = 20)
109 {
110  scalar an = 1;
111  scalar sum = 0;
112 
113  for (int n = 1; n <= (N - 1); n++)
114  {
115  an *= (a - n);
116  sum += an/pow(x, n);
117  }
118 
119  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
120 
121  return R/x*(1 + sum);
122 }
123 
125 // Eq. (18)
126 static scalar calcTE18
127 (
128  const scalar a,
129  const scalar e0,
130  const scalar x,
131  const scalar lambda,
132  const scalar sigma,
133  const scalar phi
134 )
135 {
136  static const scalar D0[] =
137  {
138  -0.333333333333333E-00,
139  0.833333333333333E-01,
140  -0.148148148148148E-01,
141  0.115740740740741E-02,
142  0.352733686067019E-03,
143  -0.178755144032922E-03,
144  0.391926317852244E-04,
145  -0.218544851067999E-05,
146  -0.185406221071516E-05,
147  0.829671134095309E-06,
148  -0.176659527368261E-06,
149  0.670785354340150E-08,
150  0.102618097842403E-07,
151  -0.438203601845335E-08
152  };
153 
154  static const scalar D1[] =
155  {
156  -0.185185185185185E-02,
157  -0.347222222222222E-02,
158  0.264550264550265E-02,
159  -0.990226337448560E-03,
160  0.205761316872428E-03,
161  -0.401877572016461E-06,
162  -0.180985503344900E-04,
163  0.764916091608111E-05,
164  -0.161209008945634E-05,
165  0.464712780280743E-08,
166  0.137863344691572E-06,
167  -0.575254560351770E-07,
168  0.119516285997781E-07
169  };
170 
171  static const scalar D2[] =
172  {
173  0.413359788359788E-02,
174  -0.268132716049383E-02,
175  0.771604938271605E-03,
176  0.200938786008230E-05,
177  -0.107366532263652E-03,
178  0.529234488291201E-04,
179  -0.127606351886187E-04,
180  0.342357873409614E-07,
181  0.137219573090629E-05,
182  -0.629899213838006E-06,
183  0.142806142060642E-06
184  };
185 
186  const scalar u = 1/a;
187  scalar z = sqrt(2*phi);
188 
189  if (lambda < 1)
190  {
191  z = -z;
192  }
193 
194  if (sigma > (e0/sqrt(a)))
195  {
196  const scalar C0 =
197  D0[6]*pow6(z) + D0[5]*pow5(z) + D0[4]*pow4(z)
198  + D0[3]*pow3(z) + D0[2]*sqr(z) + D0[1]*z + D0[0];
199 
200  const scalar C1 =
201  D1[4]*pow4(z) + D1[3]*pow3(z) + D1[2]*sqr(z) + D1[1]*z + D1[0];
202 
203  const scalar C2 = D2[1]*z + D2[0];
204 
205  return C2*sqr(u) + C1*u + C0;
206  }
207  else
208  {
209  const scalar C0 = D0[2]*sqr(z) + D0[1]*z + D0[0];
210  const scalar C1 = D1[1]*z + D1[0];
211  const scalar C2 = D2[1]*z + D2[0];
212 
213  return C2*sqr(u) + C1*u + C0;
214  }
215 }
216 
217 
218 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
219 
220 } // End namespace Foam
222 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
223 
224 Foam::scalar Foam::incGammaRatio_Q(const scalar a, const scalar x)
225 {
226  const scalar BIG = 14;
227  const scalar x0 = 17;
228  const scalar e0 = 0.025;
229 
230  if (a < 1)
231  {
232  if (a == 0.5)
233  {
234  // Eqn. (8)
235  if (x < 0.25)
236  {
237  return 1 - erf(sqrt(x));
238  }
239  else
240  {
241  return erfc(sqrt(x));
242  }
243  }
244  else if ( x < 1.1)
245  {
246  // Eqn. (12)
247  scalar alpha = x/2.59;
248 
249  if (x < 0.5)
250  {
251  alpha = log(sqrt(0.765))/log(x);
252  }
253 
254  scalar sum = 0;
255 
256  for (int n = 1; n <= 10; n++)
257  {
258  sum += pow((-x), n)/((a + n)*factorial(n));
259  }
260 
261  const scalar J = -a*sum;
262 
263  if (a > alpha || a == alpha)
264  {
265  // Eqn. (9)
266  return 1 - (pow(x, a)*(1 - J))/tgamma(a + 1);
267  }
268  else
269  {
270  // Eqn. (10)
271  const scalar L = exp(a*log(x)) - 1;
272  const scalar H = 1/(tgamma(a + 1)) - 1;
273 
274  return (pow(x, a)*J - L)/tgamma(a + 1) - H;
275  }
276  }
277  else
278  {
279  // Eqn. (11)
280  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
281 
282  return R*calcQE11(a, x);
283  }
284  }
285  else if (a >= BIG)
286  {
287  const scalar sigma = fabs(1 - x/a);
288 
289  if (sigma <= e0/sqrt(a))
290  {
291  // Eqn. (19)
292  const scalar lambda = x/a;
293  const scalar phi = lambda - 1 - log(lambda);
294  const scalar y = a*phi;
295 
296  const scalar E = 0.5 - (1 - y/3)*sqrt(y/pi);
297 
298  if (lambda <= 1)
299  {
300  return
301  1
302  - (
303  E
304  - (1 - y)/sqrt(2*pi*a)
305  *calcTE18(a, e0, x, lambda, sigma, phi)
306  );
307  }
308  else
309  {
310  return
311  E
312  + (1 - y)/sqrt(2*pi*a)
313  *calcTE18(a, e0, x, lambda, sigma, phi);
314  }
315  }
316  else
317  {
318  if (sigma <= 0.4)
319  {
320  // Eqn. (17)
321  const scalar lambda = x/a;
322  const scalar phi = lambda - 1 - log(lambda);
323  const scalar y = a*phi;
324 
325  if (lambda <= 1)
326  {
327  return
328  1
329  - (0.5*erfc(sqrt(y))
330  - exp(-y)/sqrt(2*pi*a)
331  *calcTE18(a, e0, x, lambda, sigma, phi));
332  }
333  else
334  {
335  return
336  0.5*erfc(sqrt(y))
337  + exp(-y)/sqrt(2*pi*a)
338  *calcTE18(a, e0, x, lambda, sigma, phi);
339  }
340  }
341  else
342  {
343  if (x <= max(a, log(10.0)))
344  {
345  // Eqn. (15)
346  return 1 - calcPE15(a, x);
347  }
348  else if (x < x0)
349  {
350  // Eqn. (11)
351  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
352 
353  return R*calcQE11(a, x);
354  }
355  else
356  {
357  // Eqn. (16)
358  return calcQE16(a, x);
359  }
360  }
361  }
362  }
363  else
364  {
365  if (a > x || x >= x0)
366  {
367  if (x <= max(a, log(10.0)))
368  {
369  // Eqn. (15)
370  return 1 - calcPE15(a, x);
371  }
372  else if ( x < x0)
373  {
374  // Eqn. (11)
375  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
376 
377  return R*calcQE11(a, x);
378  }
379  else
380  {
381  // Eqn. (16)
382  return calcQE16(a, x);
383  }
384  }
385  else
386  {
387  if (floor(2*a) == 2*a)
388  {
389  // Eqn. (14)
390  if (floor(a) == a)
391  {
392  scalar sum = 0;
393 
394  for (int n = 0; n <= (a - 1); n++)
395  {
396  sum += pow(x, n)/factorial(n);
397  }
398 
399  return exp(-x)*sum;
400  }
401  else
402  {
403  int i = a - 0.5;
404  scalar prod = 1;
405  scalar sum = 0;
406 
407  for (int n = 1; n <= i; n++)
408  {
409  prod *= (n - 0.5);
410  sum += pow(x, n)/prod;
411  }
412 
413  return erfc(sqrt(x)) + exp(-x)/sqrt(pi*x)*sum;
414  }
415  }
416  else if (x <= max(a, log(10.0)))
417  {
418  // Eqn. (15)
419  return 1 - calcPE15(a, x);
420  }
421  else if ( x < x0)
422  {
423  // Eqn. (11)
424  const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
425 
426  return R*calcQE11(a, x);
427  }
428  else
429  {
430  // Eqn. (16)
431  return calcQE16(a, x);
432  }
433  }
434  }
435 }
436 
437 
438 Foam::scalar Foam::incGammaRatio_P(const scalar a, const scalar x)
439 {
440  return 1 - incGammaRatio_Q(a, x);
441 }
442 
443 
444 Foam::scalar Foam::incGamma_Q(const scalar a, const scalar x)
445 {
446  return incGammaRatio_Q(a, x)*tgamma(a);
447 }
448 
449 
450 Foam::scalar Foam::incGamma_P(const scalar a, const scalar x)
451 {
452  return incGammaRatio_P(a, x)*tgamma(a);
453 }
454 
455 
456 // ************************************************************************* //
scalar incGammaRatio_P(const scalar a, const scalar x)
Normalised lower incomplete gamma function.
Definition: incGamma.C:435
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow5(const dimensionedScalar &ds)
static scalar calcPE15(const scalar a, const scalar x, const int nmax=20)
Definition: incGamma.C:87
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
scalar y
dimensionedScalar exp(const dimensionedScalar &ds)
mathematical constants.
label factorial(label n)
Return n! : 0 < n <= 12.
Definition: label.C:63
phi
Definition: correctPhi.H:3
scalar incGammaRatio_Q(const scalar a, const scalar x)
Normalised upper incomplete gamma function.
Definition: incGamma.C:221
dimensionedScalar erf(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar erfc(const dimensionedScalar &ds)
static scalar calcQE16(const scalar a, const scalar x, const int N=20)
Definition: incGamma.C:105
#define R(A, B, C, D, E, F, K, M)
dimensionedScalar pow4(const dimensionedScalar &ds)
dimensionedScalar pow6(const dimensionedScalar &ds)
label N
Definition: createFields.H:9
label n
scalar incGamma_P(const scalar a, const scalar x)
Lower incomplete gamma function.
Definition: incGamma.C:447
scalar incGamma_Q(const scalar a, const scalar x)
Upper incomplete gamma function.
Definition: incGamma.C:441
static scalar calcTE18(const scalar a, const scalar e0, const scalar x, const scalar lambda, const scalar sigma, const scalar phi)
Definition: incGamma.C:124
Namespace for OpenFOAM.
static scalar calcQE11(const scalar a, const scalar x, const int e=30)
Definition: incGamma.C:54