All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-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::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 'invIncGamma' is described in section 4.
40 
41 \*---------------------------------------------------------------------------*/
42 
43 #include "mathematicalConstants.H"
44 using namespace Foam::constant::mathematical;
45 
46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 
48 namespace Foam
49 {
50 
51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
52 
53 static scalar minimaxs(const scalar P)
54 {
55  // Eqn. 32
56 
57  static const scalar a[] =
58  {
59  3.31125922108741,
60  11.6616720288968,
61  4.28342155967104,
62  0.213623493715853
63  };
64 
65  static const scalar b[] =
66  {
67  6.61053765625462,
68  6.40691597760039,
69  1.27364489782223,
70  0.03611708101884203
71  };
72 
73  const scalar t = P < 0.5 ? sqrt(-2*log(P)) : sqrt(-2*log(1 - P));
74 
75  const scalar s =
76  t
77  - (a[0] + t*(a[1] + t*(a[2] + t*a[3])))
78  /(1 + t*(b[0] + t*(b[1] + t*(b[2] + t*b[3]))));
79 
80  return P < 0.5 ? -s : s;
81 }
82 
83 
84 static scalar Sn(const scalar a, const scalar x)
85 {
86  //- Eqn. 34
87 
88  scalar Sn = 1;
89  scalar Si = 1;
90 
91  for (int i=1; i<100; i++)
92  {
93  Si *= x/(a + i);
94  Sn += Si;
95 
96  if (Si < 1e-4) break;
97  }
98 
99  return Sn;
100 }
101 
102 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
103 
104 } // End namespace Foam
105 
106 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
107 
108 //- Inverse normalised incomplete gamma function
109 Foam::scalar Foam::invIncGamma(const scalar a, const scalar P)
110 {
111  const scalar Q = 1 - P;
112 
113  if (a == 1)
114  {
115  return -log(Q);
116  }
117  else if (a < 1)
118  {
119  const scalar Ga = tgamma(a);
120  const scalar B = Q*Ga;
121 
122  if (B > 0.6 || (B >= 0.45 && a >= 0.3))
123  {
124  // Eqn. 21
125  const scalar u =
126  (B*Q > 1e-8) ? pow(P*Ga*a, 1/a) : exp((-Q/a) - Eu);
127 
128  return u/(1 - (u/(a + 1)));
129  }
130  else if (a < 0.3 && B >= 0.35)
131  {
132  // Eqn. 22
133  const scalar t = exp(-Eu - B);
134  const scalar u = t*exp(t);
135 
136  return t*exp(u);
137  }
138  else if (B > 0.15 || a >= 0.3)
139  {
140  // Eqn. 23
141  const scalar y = -log(B);
142  const scalar u = y - (1 - a)*log(y);
143 
144  return y - (1 - a)*log(u) - log(1 + (1 - a)/(1 + u));
145  }
146  else if (B > 0.1)
147  {
148  // Eqn. 24
149  const scalar y = -log(B);
150  const scalar u = y - (1 - a)*log(y);
151 
152  return y
153  - (1 - a)*log(u)
154  - log
155  (
156  (sqr(u) + 2*(3 - a)*u + (2 - a)*(3 - a))
157  /(sqr(u) + (5 - a)*u + 2)
158  );
159  }
160  else
161  {
162  // Eqn. 25:
163  const scalar y = -log(B);
164  const scalar c1 = (a - 1)*log(y);
165  const scalar c12 = c1*c1;
166  const scalar c13 = c12*c1;
167  const scalar c14 = c12*c12;
168  const scalar a2 = a*a;
169  const scalar a3 = a2*a;
170  const scalar c2 = (a - 1)*(1 + c1);
171  const scalar c3 = (a - 1)*(-(c12/2) + (a - 2)*c1 + (3*a - 5)/2);
172  const scalar c4 =
173  (a - 1)
174  *(
175  (c13/3)
176  - (3*a - 5)*c12/2
177  + (a2 - 6*a + 7)*c1
178  + (11*a2 - 46*a + 47)/6
179  );
180  const scalar c5 =
181  (a - 1)*(-(c14/4)
182  + (11*a - 17)*c13/6
183  + (-3*a2 + 13*a - 13)*c12
184  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
185  + (25*a3 - 195*a2 + 477*a - 379)/12);
186  const scalar y2 = y*y;
187  const scalar y3 = y2*y;
188  const scalar y4 = y2*y2;
189 
190  return y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
191  }
192  }
193  else
194  {
195  // Eqn. 31:
196  scalar s = minimaxs(P);
197 
198  const scalar s2 = sqr(s);
199  const scalar s3 = s*s2;
200  const scalar s4 = s2*s2;
201  const scalar s5 = s*s4;
202  const scalar sqrta = sqrt(a);
203 
204  const scalar w =
205  a + s*sqrta + (s2 - 1)/3
206  + (s3 - 7*s)/(36*sqrta)
207  - (3*s4 + 7*s2 - 16)/(810*a)
208  + (9*s5 + 256*s3 - 433*s)/(38880*a*sqrta);
209 
210  if (a >= 500 && mag(1 - w/a) < 1e-6)
211  {
212  return w;
213  }
214  else if (P > 0.5)
215  {
216  if (w < 3*a)
217  {
218  return w;
219  }
220  else
221  {
222  const scalar D = max(scalar(2), scalar(a*(a - 1)));
223  const scalar lnGa = lgamma(a);
224  const scalar lnB = log(Q) + lnGa;
225 
226  if (lnB < -2.3*D)
227  {
228  // Eqn. 25:
229  const scalar y = -lnB;
230  const scalar c1 = (a - 1)*log(y);
231  const scalar c12 = c1*c1;
232  const scalar c13 = c12*c1;
233  const scalar c14 = c12*c12;
234  const scalar a2 = a*a;
235  const scalar a3 = a2*a;
236 
237  const scalar c2 = (a - 1)*(1 + c1);
238  const scalar c3 =
239  (a - 1)
240  *(
241  - (c12/2)
242  + (a - 2)*c1
243  + (3*a - 5)/2
244  );
245  const scalar c4 =
246  (a - 1)
247  *(
248  (c13/3)
249  - (3*a - 5)*c12/2
250  + (a2 - 6*a + 7)*c1
251  + (11*a2 - 46*a + 47)/6
252  );
253  const scalar c5 =
254  (a - 1)
255  *(
256  - (c14/4)
257  + (11*a - 17)*c13/6
258  + (-3*a2 + 13*a - 13)*c12
259  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
260  + (25*a3 - 195*a2 + 477*a - 379)/12
261  );
262 
263  const scalar y2 = y*y;
264  const scalar y3 = y2*y;
265  const scalar y4 = y2*y2;
266 
267  return y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
268  }
269  else
270  {
271  // Eqn. 33
272  const scalar u =
273  -lnB + (a - 1)*log(w) - log(1 + (1 - a)/(1 + w));
274 
275  return -lnB + (a - 1)*log(u) - log(1 + (1 - a)/(1 + u));
276  }
277  }
278  }
279  else
280  {
281  scalar z = w;
282  const scalar ap1 = a + 1;
283 
284  if (w < 0.15*ap1)
285  {
286  // Eqn. 35
287  const scalar ap2 = a + 2;
288  const scalar v = log(P) + lgamma(ap1);
289  z = exp((v + w)/a);
290  s = log1p(z/ap1*(1 + z/ap2));
291  z = exp((v + z - s)/a);
292  s = log1p(z/ap1*(1 + z/ap2));
293  z = exp((v + z - s)/a);
294  s = log1p(z/ap1*(1 + z/ap2*(1 + z/(a + 3))));
295  z = exp((v + z - s)/a);
296  }
297 
298  if (z <= 0.01*ap1 || z > 0.7*ap1)
299  {
300  return z;
301  }
302  else
303  {
304  // Eqn. 36
305  const scalar lnSn = log(Sn(a, z));
306  const scalar v = log(P) + lgamma(ap1);
307  z = exp((v + z - lnSn)/a);
308 
309  return z*(1 - (a*log(z) - z - v + lnSn)/(a - z));
310  }
311  }
312  }
313 }
314 
315 
316 // ************************************************************************* //
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const scalar Eu(0.57721566490153286060651209)
Euler&#39;s constant.
dimensionedScalar log(const dimensionedScalar &ds)
static scalar Sn(const scalar a, const scalar x)
Definition: invIncGamma.C:81
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
static scalar minimaxs(const scalar P)
Definition: invIncGamma.C:50
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
dimensionedScalar exp(const dimensionedScalar &ds)
mathematical constants.
scalar invIncGamma(const scalar a, const scalar P)
Inverse normalised incomplete gamma function.
Definition: invIncGamma.C:106
dimensionedScalar lgamma(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
Namespace for OpenFOAM.