invIncGamma.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) 2016 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::invIncGamma
29 
30 Description
31  Implementation of the inverse incomplete gamma function.
32 
33 \*---------------------------------------------------------------------------*/
34 
35 #include "MathFunctions.H"
36 #include "mathematicalConstants.H"
37 #include "error.H"
38 
39 using namespace Foam::constant::mathematical;
40 
41 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45 
46 static scalar minimaxs(const scalar P)
47 {
48  // (DM:Eq. 32)
49 
50  constexpr scalar a_0 = 3.31125922108741;
51  constexpr scalar a_1 = 11.6616720288968;
52  constexpr scalar a_2 = 4.28342155967104;
53  constexpr scalar a_3 = 0.213623493715853;
54 
55  constexpr scalar b_0 = 6.61053765625462;
56  constexpr scalar b_1 = 6.40691597760039;
57  constexpr scalar b_2 = 1.27364489782223;
58  constexpr scalar b_3 = 0.03611708101884203;
59 
60  const scalar t = P < 0.5 ? sqrt(-2*log(P)) : sqrt(-2*log(1 - P));
61 
62  const scalar s =
63  t
64  - (a_0 + t*(a_1 + t*(a_2 + t*a_3)))
65  /(1 + t*(b_0 + t*(b_1 + t*(b_2 + t*b_3))));
66 
67  return P < 0.5 ? -s : s;
68 }
69 
70 
71 static scalar Sn(const scalar a, const scalar x)
72 {
73  // (DM:Eq. 34)
74 
75  scalar Sn = 1;
76  scalar Si = 1;
77 
78  for (int i=1; i<100; ++i)
79  {
80  Si *= x/(a + i);
81  Sn += Si;
82 
83  if (Si < 1e-4) break;
84  }
85 
86  return Sn;
87 }
88 
89 } // End namespace Foam
90 
91 
92 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
93 
94 Foam::scalar Foam::Math::invIncGamma(const scalar a, const scalar P)
95 {
96  #ifdef FULLDEBUG
97  if (a <= 0)
98  {
99  WarningInFunction
100  << "The parameter (i.e. a) cannot be negative or zero"
101  << " a = " << a
102  << endl;
103 
104  return std::numeric_limits<scalar>::infinity();
105  }
106 
107  if (P < 0 || P > 1)
108  {
109  WarningInFunction
110  << "The domain of the parameter (i.e. P) should be limited to [0,1]"
111  << " P = " << P
112  << endl;
113 
114  return std::numeric_limits<scalar>::infinity();
115  }
116  #endif
117 
118  const scalar Q = 1 - P;
119 
120  if (a == 1)
121  {
122  return -log(Q);
123  }
124  else if (a < 1)
125  {
126  const scalar Ga = tgamma(a);
127  const scalar B = Q*Ga;
128 
129  if (B > 0.6 || (B >= 0.45 && a >= 0.3))
130  {
131  // (DM:Eq. 21)
132  const scalar u =
133  (B*Q > 1e-8) ? pow(P*Ga*a, 1/a) : exp((-Q/a) - Eu);
134 
135  return u/(1 - (u/(a + 1)));
136  }
137  else if (a < 0.3 && B >= 0.35)
138  {
139  // (DM:Eq. 22)
140  const scalar t = exp(-Eu - B);
141  const scalar u = t*exp(t);
142 
143  return t*exp(u);
144  }
145  else if (B > 0.15 || a >= 0.3)
146  {
147  // (DM:Eq. 23)
148  const scalar y = -log(B);
149  const scalar u = y - (1 - a)*log(y);
150 
151  return y - (1 - a)*log(u) - log(1 + (1 - a)/(1 + u));
152  }
153  else if (B > 0.1)
154  {
155  // (DM:Eq. 24)
156  const scalar y = -log(B);
157  const scalar u = y - (1 - a)*log(y);
158 
159  return y
160  - (1 - a)*log(u)
161  - log
162  (
163  (sqr(u) + 2*(3 - a)*u + (2 - a)*(3 - a))
164  /(sqr(u) + (5 - a)*u + 2)
165  );
166  }
167  else
168  {
169  // (DM:Eq. 25)
170  const scalar y = -log(B);
171  const scalar c1 = (a - 1)*log(y);
172  const scalar c12 = c1*c1;
173  const scalar c13 = c12*c1;
174  const scalar c14 = c12*c12;
175  const scalar a2 = a*a;
176  const scalar a3 = a2*a;
177  const scalar c2 = (a - 1)*(1 + c1);
178  const scalar c3 = (a - 1)*(-(c12/2) + (a - 2)*c1 + (3*a - 5)/2);
179  const scalar c4 =
180  (a - 1)
181  *(
182  (c13/3)
183  - (3*a - 5)*c12/2
184  + (a2 - 6*a + 7)*c1
185  + (11*a2 - 46*a + 47)/6
186  );
187  const scalar c5 =
188  (a - 1)*(-(c14/4)
189  + (11*a - 17)*c13/6
190  + (-3*a2 + 13*a - 13)*c12
191  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
192  + (25*a3 - 195*a2 + 477*a - 379)/12);
193  const scalar y2 = y*y;
194  const scalar y3 = y2*y;
195  const scalar y4 = y2*y2;
196 
197  return y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
198  }
199  }
200  else
201  {
202  // (DM:Eq. 31)
203  scalar s = minimaxs(P);
204 
205  const scalar s2 = sqr(s);
206  const scalar s3 = s*s2;
207  const scalar s4 = s2*s2;
208  const scalar s5 = s*s4;
209  const scalar sqrta = sqrt(a);
210 
211  const scalar w =
212  a + s*sqrta + (s2 - 1)/3
213  + (s3 - 7*s)/(36*sqrta)
214  - (3*s4 + 7*s2 - 16)/(810*a)
215  + (9*s5 + 256*s3 - 433*s)/(38880*a*sqrta);
216 
217  if (a >= 500 && mag(1 - w/a) < 1e-6)
218  {
219  return w;
220  }
221  else if (P > 0.5)
222  {
223  if (w < 3*a)
224  {
225  return w;
226  }
227  else
228  {
229  const scalar D = max(scalar(2), scalar(a*(a - 1)));
230  const scalar lnGa = lgamma(a);
231  const scalar lnB = log(Q) + lnGa;
232 
233  if (lnB < -2.3*D)
234  {
235  // (DM:Eq. 25)
236  const scalar y = -lnB;
237  const scalar c1 = (a - 1)*log(y);
238  const scalar c12 = c1*c1;
239  const scalar c13 = c12*c1;
240  const scalar c14 = c12*c12;
241  const scalar a2 = a*a;
242  const scalar a3 = a2*a;
243 
244  const scalar c2 = (a - 1)*(1 + c1);
245  const scalar c3 =
246  (a - 1)
247  *(
248  - (c12/2)
249  + (a - 2)*c1
250  + (3*a - 5)/2
251  );
252  const scalar c4 =
253  (a - 1)
254  *(
255  (c13/3)
256  - (3*a - 5)*c12/2
257  + (a2 - 6*a + 7)*c1
258  + (11*a2 - 46*a + 47)/6
259  );
260  const scalar c5 =
261  (a - 1)
262  *(
263  - (c14/4)
264  + (11*a - 17)*c13/6
265  + (-3*a2 + 13*a - 13)*c12
266  + (2*a3 - 25*a2 + 72*a - 61)*c1/2
267  + (25*a3 - 195*a2 + 477*a - 379)/12
268  );
269 
270  const scalar y2 = y*y;
271  const scalar y3 = y2*y;
272  const scalar y4 = y2*y2;
273 
274  return y + c1 + (c2/y) + (c3/y2) + (c4/y3) + (c5/y4);
275  }
276  else
277  {
278  // (DM:Eq. 33)
279  const scalar u =
280  -lnB + (a - 1)*log(w) - log(1 + (1 - a)/(1 + w));
281 
282  return -lnB + (a - 1)*log(u) - log(1 + (1 - a)/(1 + u));
283  }
284  }
285  }
286  else
287  {
288  scalar z = w;
289  const scalar ap1 = a + 1;
290 
291  if (w < 0.15*ap1)
292  {
293  // (DM:Eq. 35)
294  const scalar ap2 = a + 2;
295  const scalar v = log(P) + lgamma(ap1);
296  z = exp((v + w)/a);
297  s = log1p(z/ap1*(1 + z/ap2));
298  z = exp((v + z - s)/a);
299  s = log1p(z/ap1*(1 + z/ap2));
300  z = exp((v + z - s)/a);
301  s = log1p(z/ap1*(1 + z/ap2*(1 + z/(a + 3))));
302  z = exp((v + z - s)/a);
303  }
304 
305  if (z <= 0.01*ap1 || z > 0.7*ap1)
306  {
307  return z;
308  }
309  else
310  {
311  // (DM:Eq. 36)
312  const scalar lnSn = log(Sn(a, z));
313  const scalar v = log(P) + lgamma(ap1);
314  z = exp((v + z - lnSn)/a);
315 
316  return z*(1 - (a*log(z) - z - v + lnSn)/(a - z));
317  }
318  }
319  }
320 }
321 
322 
323 // ************************************************************************* //
mathematicalConstants.H
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
B
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
Foam::constant::mathematical::e
constexpr scalar e(M_E)
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::lgamma
dimensionedScalar lgamma(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:278
Foam::constant::mathematical::Eu
constexpr scalar Eu(0.57721566490153286060651209)
Euler's constant.
Foam::constant::physicoChemical::c1
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
MathFunctions.H
error.H
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::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:262
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::constant::physicoChemical::c2
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Foam::Math::invIncGamma
scalar invIncGamma(const scalar a, const scalar P)
Inverse of regularised lower incomplete gamma function.
Definition: invIncGamma.C:91
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical
Mathematical constants.
Foam::Sn
static scalar Sn(const scalar a, const scalar x)
Definition: invIncGamma.C:68
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::minimaxs
static scalar minimaxs(const scalar P)
Definition: invIncGamma.C:43
D
const dimensionedScalar & D
Definition: solveBulkSurfactant.H:4
x
x
Definition: LISASMDCalcMethod2.H:52
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
y
scalar y
Definition: LISASMDCalcMethod1.H:14