53 static scalar
minimaxs(
const scalar P)
57 static const scalar a[] =
65 static const scalar
b[] =
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]))));
80 return P < 0.5 ? -
s :
s;
84 static scalar
Sn(
const scalar a,
const scalar
x)
91 for (
int i=1; i<100; i++)
111 const scalar Q = 1 - P;
119 const scalar Ga = tgamma(a);
120 const scalar B = Q*Ga;
122 if (B > 0.6 || (B >= 0.45 && a >= 0.3))
126 (B*Q > 1
e-8) ?
pow(P*Ga*a, 1/a) :
exp((-Q/a) -
Eu);
128 return u/(1 - (u/(a + 1)));
130 else if (a < 0.3 && B >= 0.35)
133 const scalar t =
exp(-
Eu - B);
134 const scalar u = t*
exp(t);
138 else if (B > 0.15 || a >= 0.3)
141 const scalar
y = -
log(B);
142 const scalar u = y - (1 - a)*
log(y);
144 return y - (1 - a)*
log(u) -
log(1 + (1 - a)/(1 + u));
149 const scalar
y = -
log(B);
150 const scalar u = y - (1 - a)*
log(y);
156 (
sqr(u) + 2*(3 - a)*u + (2 - a)*(3 - a))
157 /(
sqr(u) + (5 - a)*u + 2)
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);
178 + (11*a2 - 46*a + 47)/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;
190 return y + c1 + (c2/
y) + (c3/y2) + (c4/y3) + (c5/y4);
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);
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);
210 if (a >= 500 &&
mag(1 - w/a) < 1
e-6)
222 const scalar D =
max(scalar(2), scalar(a*(a - 1)));
223 const scalar lnGa =
lgamma(a);
224 const scalar lnB =
log(Q) + lnGa;
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;
237 const scalar
c2 = (a - 1)*(1 + c1);
251 + (11*a2 - 46*a + 47)/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
263 const scalar y2 = y*
y;
264 const scalar y3 = y2*
y;
265 const scalar y4 = y2*y2;
267 return y + c1 + (c2/
y) + (c3/y2) + (c4/y3) + (c5/y4);
273 -lnB + (a - 1)*
log(w) -
log(1 + (1 - a)/(1 + w));
275 return -lnB + (a - 1)*
log(u) -
log(1 + (1 - a)/(1 + u));
282 const scalar ap1 = a + 1;
287 const scalar ap2 = a + 2;
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);
298 if (z <= 0.01*ap1 || z > 0.7*ap1)
305 const scalar lnSn =
log(
Sn(a, z));
307 z =
exp((v + z - lnSn)/a);
309 return z*(1 - (a*
log(z) - z - v + lnSn)/(a - z));
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const scalar Eu(0.57721566490153286060651209)
Euler's constant.
dimensionedScalar log(const dimensionedScalar &ds)
static scalar Sn(const scalar a, const scalar x)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
static scalar minimaxs(const scalar P)
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
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)
scalar invIncGamma(const scalar a, const scalar P)
Inverse normalised incomplete gamma function.
dimensionedScalar lgamma(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)