invIncGamma.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) 2016-2023 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::invIncGamma
26 
27 Description
28  Function to calculate the inverse of the
29  normalised incomplete gamma function.
30 
31  The algorithm is described in detain in reference:
32  \verbatim
33  DiDonato, A. R., & Morris Jr, A. H. (1986).
34  Computation of the incomplete gamma function ratios and their inverse.
35  ACM Transactions on Mathematical Software (TOMS), 12(4), 377-393.
36  \endverbatim
37 
38  All equation numbers in the following code refer to the above text and
39  the algorithm in the function 'invIncGammaRatio_P' is described in section
40  4.
41 
42 \*---------------------------------------------------------------------------*/
43 
44 #include "mathematicalConstants.H"
45 
46 using namespace Foam::constant::mathematical;
47 
48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
49 
50 namespace Foam
51 {
52 
53 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
54 
55 static scalar minimaxs(const scalar P)
56 {
57  // Eqn. 32
58 
59  static const scalar a[] =
60  {
61  3.31125922108741,
62  11.6616720288968,
63  4.28342155967104,
64  0.213623493715853
65  };
66 
67  static const scalar b[] =
68  {
69  6.61053765625462,
70  6.40691597760039,
71  1.27364489782223,
72  0.03611708101884203
73  };
74 
75  const scalar t = P < 0.5 ? sqrt(-2*log(P)) : sqrt(-2*log(1 - P));
76 
77  const scalar s =
78  t
79  - (a[0] + t*(a[1] + t*(a[2] + t*a[3])))
80  /(1 + t*(b[0] + t*(b[1] + t*(b[2] + t*b[3]))));
81 
82  return P < 0.5 ? -s : s;
83 }
84 
85 
86 static scalar Sn(const scalar a, const scalar x)
87 {
88  //- Eqn. 34
89 
90  scalar Sn = 1;
91  scalar Si = 1;
92 
93  for (int i=1; i<100; i++)
94  {
95  Si *= x/(a + i);
96  Sn += Si;
97 
98  if (Si < 1e-4) break;
99  }
100 
101  return Sn;
102 }
103 
104 
105 static scalar R(const scalar a, const scalar x)
106 {
107  //- Eqn. 3
108  return exp(-x)*pow(x, a)/tgamma(a);
109 }
110 
111 
112 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
113 
114 } // End namespace Foam
115 
116 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
117 
118 Foam::scalar Foam::invIncGammaRatio_P(const scalar a, const scalar P)
119 {
120  const scalar Q = 1 - P;
121 
122  // Determine an initial guess (and potentially do a quick return)
123  scalar x = NaN;
124 
125  if (a == 1)
126  {
127  x = -log(Q);
128  }
129  else if (a < 1)
130  {
131  const scalar Ga = tgamma(a);
132  const scalar B = Q*Ga;
133 
134  if (B > 0.6 || (B >= 0.45 && a >= 0.3))
135  {
136  // Eqn. 21
137  const scalar u =
138  (B*Q > 1e-8) ? pow(P*Ga*a, 1/a) : exp((-Q/a) - Eu);
139 
140  x = u/(1 - (u/(a + 1)));
141  }
142  else if (a < 0.3 && B >= 0.35)
143  {
144  // Eqn. 22
145  const scalar t = exp(-Eu - B);
146  const scalar u = t*exp(t);
147 
148  x = t*exp(u);
149  }
150  else if (B > 0.15 || a >= 0.3)
151  {
152  // Eqn. 23
153  const scalar y = -log(B);
154  const scalar u = y - (1 - a)*log(y);
155 
156  x = y - (1 - a)*log(u) - log(1 + (1 - a)/(1 + u));
157  }
158  else if (B > 0.1)
159  {
160  // Eqn. 24
161  const scalar y = -log(B);
162  const scalar u = y - (1 - a)*log(y);
163 
164  x =
165  y
166  - (1 - a)*log(u)
167  - log
168  (
169  (sqr(u) + 2*(3 - a)*u + (2 - a)*(3 - a))
170  /(sqr(u) + (5 - a)*u + 2)
171  );
172  }
173  else
174  {
175  // Eqn. 25:
176  const scalar y = -log(B);
177  const scalar c1 = (a - 1)*log(y);
178  const scalar c12 = c1*c1;
179  const scalar c13 = c12*c1;
180  const scalar c14 = c12*c12;
181  const scalar a2 = a*a;
182  const scalar a3 = a2*a;
183  const scalar c2 = (a - 1)*(1 + c1);
184  const scalar c3 = (a - 1)*(-(c12/2) + (a - 2)*c1 + (3*a - 5)/2);
185  const scalar c4 =
186  (a - 1)
187  *(
188  (c13/3)
189  - (3*a - 5)*c12/2
190  + (a2 - 6*a + 7)*c1
191  + (11*a2 - 46*a + 47)/6
192  );
193  const scalar c5 =
194  (a - 1)*(-(c14/4)
195  + (11*a - 17)*c13/6
196  + (-3*a2 + 13*a - 13)*c12
197  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
198  + (25*a3 - 195*a2 + 477*a - 379)/12);
199  const scalar y2 = y*y;
200  const scalar y3 = y2*y;
201  const scalar y4 = y2*y2;
202 
203  x = y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
204  }
205 
206  // Quick return condition (a)
207  if (B <= 1e-13)
208  {
209  return x;
210  }
211  }
212  else
213  {
214  // Eqn. 31:
215  scalar s = minimaxs(P);
216 
217  const scalar s2 = sqr(s);
218  const scalar s3 = s*s2;
219  const scalar s4 = s2*s2;
220  const scalar s5 = s*s4;
221  const scalar sqrta = sqrt(a);
222 
223  const scalar w =
224  a + s*sqrta + (s2 - 1)/3
225  + (s3 - 7*s)/(36*sqrta)
226  - (3*s4 + 7*s2 - 16)/(810*a)
227  + (9*s5 + 256*s3 - 433*s)/(38880*a*sqrta);
228 
229  // Quick return condition (b)
230  if (a >= 100 && mag(1 - w/a) < 1e-4)
231  {
232  return w;
233  }
234 
235  if (a >= 500 && mag(1 - w/a) < 1e-6)
236  {
237  x = w;
238  }
239  else if (P > 0.5)
240  {
241  if (w < 3*a)
242  {
243  x = w;
244  }
245  else
246  {
247  const scalar D = max(scalar(2), scalar(a*(a - 1)));
248  const scalar lnGa = lgamma(a);
249  const scalar lnB = log(Q) + lnGa;
250 
251  if (lnB < -2.3*D)
252  {
253  // Eqn. 25:
254  const scalar y = -lnB;
255  const scalar c1 = (a - 1)*log(y);
256  const scalar c12 = c1*c1;
257  const scalar c13 = c12*c1;
258  const scalar c14 = c12*c12;
259  const scalar a2 = a*a;
260  const scalar a3 = a2*a;
261 
262  const scalar c2 = (a - 1)*(1 + c1);
263  const scalar c3 =
264  (a - 1)
265  *(
266  - (c12/2)
267  + (a - 2)*c1
268  + (3*a - 5)/2
269  );
270  const scalar c4 =
271  (a - 1)
272  *(
273  (c13/3)
274  - (3*a - 5)*c12/2
275  + (a2 - 6*a + 7)*c1
276  + (11*a2 - 46*a + 47)/6
277  );
278  const scalar c5 =
279  (a - 1)
280  *(
281  - (c14/4)
282  + (11*a - 17)*c13/6
283  + (-3*a2 + 13*a - 13)*c12
284  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
285  + (25*a3 - 195*a2 + 477*a - 379)/12
286  );
287 
288  const scalar y2 = y*y;
289  const scalar y3 = y2*y;
290  const scalar y4 = y2*y2;
291 
292  x = y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
293  }
294  else
295  {
296  // Eqn. 33
297  const scalar u =
298  -lnB + (a - 1)*log(w) - log(1 + (1 - a)/(1 + w));
299 
300  x = -lnB + (a - 1)*log(u) - log(1 + (1 - a)/(1 + u));
301  }
302  }
303  }
304  else
305  {
306  scalar z = w;
307  const scalar ap1 = a + 1;
308 
309  if (w < 0.15*ap1)
310  {
311  // Eqn. 35
312  const scalar ap2 = a + 2;
313  const scalar v = log(P) + lgamma(ap1);
314  z = exp((v + w)/a);
315  s = log1p(z/ap1*(1 + z/ap2));
316  z = exp((v + z - s)/a);
317  s = log1p(z/ap1*(1 + z/ap2));
318  z = exp((v + z - s)/a);
319  s = log1p(z/ap1*(1 + z/ap2*(1 + z/(a + 3))));
320  z = exp((v + z - s)/a);
321  }
322 
323  // Quick return condition (c)
324  if (z < 0.006*ap1 && (a > 1 || a < 100 || mag(1 - w/a) > 1e-4))
325  {
326  return z;
327  }
328 
329  if (z <= 0.01*ap1 || z > 0.7*ap1)
330  {
331  x = z;
332  }
333  else
334  {
335  // Eqn. 36
336  const scalar lnSn = log(Sn(a, z));
337  const scalar v = log(P) + lgamma(ap1);
338  z = exp((v + z - lnSn)/a);
339 
340  x = z*(1 - (a*log(z) - z - v + lnSn)/(a - z));
341  }
342  }
343  }
344 
345  // Iteration
346  static const label nIter = 20;
347  static const scalar eps = 1e-10;
348  static const scalar tau = 1e-5;
349 
350  for (label iter = 0; iter < nIter; ++ iter)
351  {
352  const scalar dP = incGammaRatio_P(a, max(x, vSmall)) - P;
353  const scalar t = dP/max(R(a, x), vSmall);
354  const scalar w = (a - 1 - x)/2;
355 
356  bool halley = mag(t) < 0.1 && mag(w*t) < 0.1;
357 
358  const scalar h = t + halley*w*sqr(t);
359 
360  x *= 1 - h;
361 
362  if
363  (
364  (halley && mag(w) >= 1 && mag(w*sqr(t)) <= eps)
365  || (mag(h) < eps)
366  || (mag(h) < tau && mag(dP) < tau*max(P, 1 - P))
367  )
368  {
369  return x;
370  }
371  }
372 
373  // Iteration failed to converge. Probably at zero.
374  return 0;
375 }
376 
377 
378 // ************************************************************************* //
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
static const Foam::dimensionedScalar D("D", Foam::dimTemperature, 257.14)
scalar y
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField & b
Definition: createFields.H:25
mathematical constants.
const scalar Eu(0.57721566490153286060651209)
Euler's constant.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
const dimensionedScalar h
Planck constant.
Namespace for OpenFOAM.
dimensionedScalar exp(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedScalar lgamma(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static scalar minimaxs(const scalar P)
Definition: invIncGamma.C:52
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
scalar invIncGammaRatio_P(const scalar a, const scalar P)
Inverse normalised lower incomplete gamma function.
Definition: invIncGamma.C:115
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
scalar incGammaRatio_P(const scalar a, const scalar x)
Normalised lower incomplete gamma function.
Definition: incGamma.C:435
static scalar Sn(const scalar a, const scalar x)
Definition: invIncGamma.C:83