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