Go to the documentation of this file.
55 static scalar
minimaxs(
const scalar P)
59 static const scalar a[] =
67 static const scalar
b[] =
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]))));
82 return P < 0.5 ? -
s :
s;
86 static scalar
Sn(
const scalar a,
const scalar
x)
93 for (
int i=1; i<100; i++)
113 const scalar Q = 1 - P;
121 const scalar Ga = tgamma(a);
122 const scalar
B = Q*Ga;
124 if (
B > 0.6 || (
B >= 0.45 && a >= 0.3))
128 (
B*Q > 1
e-8) ?
pow(P*Ga*a, 1/a) :
exp((-Q/a) -
Eu);
130 return u/(1 - (u/(a + 1)));
132 else if (a < 0.3 && B >= 0.35)
135 const scalar t =
exp(-
Eu -
B);
136 const scalar u = t*
exp(t);
140 else if (
B > 0.15 || a >= 0.3)
143 const scalar
y = -
log(
B);
144 const scalar u =
y - (1 - a)*
log(
y);
146 return y - (1 - a)*
log(u) -
log(1 + (1 - a)/(1 + u));
151 const scalar
y = -
log(
B);
152 const scalar u =
y - (1 - a)*
log(
y);
158 (
sqr(u) + 2*(3 - a)*u + (2 - a)*(3 - a))
159 /(
sqr(u) + (5 - a)*u + 2)
165 const scalar
y = -
log(
B);
166 const scalar
c1 = (a - 1)*
log(
y);
167 const scalar c12 =
c1*
c1;
168 const scalar c13 = c12*
c1;
169 const scalar c14 = c12*c12;
170 const scalar a2 = a*a;
171 const scalar a3 = a2*a;
172 const scalar
c2 = (a - 1)*(1 +
c1);
173 const scalar c3 = (a - 1)*(-(c12/2) + (a - 2)*
c1 + (3*a - 5)/2);
180 + (11*a2 - 46*a + 47)/6
185 + (-3*a2 + 13*a - 13)*c12
186 + (2*a3 - 25*a2 + 72*a - 61)*
c1/2
187 + (25*a3 - 195*a2 + 477*a - 379)/12);
188 const scalar y2 =
y*
y;
189 const scalar y3 = y2*
y;
190 const scalar y4 = y2*y2;
192 return y +
c1 + (
c2/
y) + (c3/y2) + (c4/y3) + (c5/y4);
200 const scalar s2 =
sqr(
s);
201 const scalar s3 =
s*s2;
202 const scalar s4 = s2*s2;
203 const scalar s5 =
s*s4;
204 const scalar sqrta =
sqrt(a);
207 a +
s*sqrta + (s2 - 1)/3
208 + (s3 - 7*
s)/(36*sqrta)
209 - (3*s4 + 7*s2 - 16)/(810*a)
210 + (9*s5 + 256*s3 - 433*
s)/(38880*a*sqrta);
212 if (a >= 500 &&
mag(1 - w/a) < 1
e-6)
224 const scalar
D =
max(scalar(2), scalar(a*(a - 1)));
225 const scalar lnGa =
lgamma(a);
226 const scalar lnB =
log(Q) + lnGa;
231 const scalar
y = -lnB;
232 const scalar
c1 = (a - 1)*
log(
y);
233 const scalar c12 =
c1*
c1;
234 const scalar c13 = c12*
c1;
235 const scalar c14 = c12*c12;
236 const scalar a2 = a*a;
237 const scalar a3 = a2*a;
239 const scalar
c2 = (a - 1)*(1 +
c1);
253 + (11*a2 - 46*a + 47)/6
260 + (-3*a2 + 13*a - 13)*c12
261 + (2*a3 - 25*a2 + 72*a - 61)*
c1/2
262 + (25*a3 - 195*a2 + 477*a - 379)/12
265 const scalar y2 =
y*
y;
266 const scalar y3 = y2*
y;
267 const scalar y4 = y2*y2;
269 return y +
c1 + (
c2/
y) + (c3/y2) + (c4/y3) + (c5/y4);
275 -lnB + (a - 1)*
log(w) -
log(1 + (1 - a)/(1 + w));
277 return -lnB + (a - 1)*
log(u) -
log(1 + (1 - a)/(1 + u));
284 const scalar ap1 = a + 1;
289 const scalar ap2 = a + 2;
292 s = log1p(z/ap1*(1 + z/ap2));
293 z =
exp((v + z -
s)/a);
294 s = log1p(z/ap1*(1 + z/ap2));
295 z =
exp((v + z -
s)/a);
296 s = log1p(z/ap1*(1 + z/ap2*(1 + z/(a + 3))));
297 z =
exp((v + z -
s)/a);
300 if (z <= 0.01*ap1 || z > 0.7*ap1)
307 const scalar lnSn =
log(
Sn(a, z));
309 z =
exp((v + z - lnSn)/a);
311 return z*(1 - (a*
log(z) - z - v + lnSn)/(a - z));
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))
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
dimensionedScalar exp(const dimensionedScalar &ds)
scalar invIncGamma(const scalar a, const scalar P)
Inverse normalized incomplete gamma function.
dimensionedScalar lgamma(const dimensionedScalar &ds)
constexpr scalar Eu(0.57721566490153286060651209)
Euler's constant.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedScalar log(const dimensionedScalar &ds)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static scalar Sn(const scalar a, const scalar x)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
static scalar minimaxs(const scalar P)
const dimensionedScalar & D