mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
460 lines
12 KiB
C
460 lines
12 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2019 OpenFOAM Foundation
|
|
Copyright (C) 2021 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Global
|
|
Foam::Math::incGamma
|
|
|
|
Description
|
|
Implementation of the incomplete gamma functions.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "MathFunctions.H"
|
|
#include "mathematicalConstants.H"
|
|
#include "error.H"
|
|
#include <cmath>
|
|
#include <limits>
|
|
|
|
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
|
|
// (DM:Eq. 13)
|
|
static scalar calcQE11(const scalar a, const scalar x, const int e = 30)
|
|
{
|
|
scalar a_2n = 0;
|
|
scalar b_2n = 1;
|
|
|
|
scalar a_2np1 = 1;
|
|
scalar b_2np1 = x;
|
|
|
|
int n = 1;
|
|
for (n = 1; (2*n) <= e; n++)
|
|
{
|
|
const scalar a_2nm1 = a_2np1;
|
|
const scalar b_2nm1 = b_2np1;
|
|
|
|
a_2n = a_2nm1 + (n - a)*a_2n;
|
|
b_2n = b_2nm1 + (n - a)*b_2n;
|
|
|
|
a_2np1 = x*a_2n + n*a_2nm1;
|
|
b_2np1 = x*b_2n + n*b_2nm1;
|
|
}
|
|
|
|
if (2*(n - 1) < e)
|
|
{
|
|
return a_2np1/b_2np1;
|
|
}
|
|
else
|
|
{
|
|
return a_2n/b_2n;
|
|
}
|
|
}
|
|
|
|
|
|
// (DM:Eq. 15)
|
|
static scalar calcPE15(const scalar a, const scalar x, const int nmax = 20)
|
|
{
|
|
scalar prod = 1;
|
|
scalar sum = 0;
|
|
|
|
for (int n = 1; n <= nmax; n++)
|
|
{
|
|
prod *= (a + n);
|
|
sum += pow(x, n)/prod;
|
|
}
|
|
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R/a*(1 + sum);
|
|
}
|
|
|
|
|
|
// (DM:Eq. 16)
|
|
static scalar calcQE16(const scalar a, const scalar x, const int N = 20)
|
|
{
|
|
scalar an = 1;
|
|
scalar sum = 0;
|
|
|
|
for (int n = 1; n <= (N - 1); n++)
|
|
{
|
|
an *= (a - n);
|
|
sum += an/pow(x, n);
|
|
}
|
|
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R/x*(1 + sum);
|
|
}
|
|
|
|
|
|
// (DM:Eq. 18)
|
|
static scalar calcTE18
|
|
(
|
|
const scalar a,
|
|
const scalar e0,
|
|
const scalar x,
|
|
const scalar lambda,
|
|
const scalar sigma,
|
|
const scalar phi
|
|
)
|
|
{
|
|
constexpr scalar D0_0 = -0.333333333333333E-00;
|
|
constexpr scalar D0_1 = 0.833333333333333E-01;
|
|
constexpr scalar D0_2 = -0.148148148148148E-01;
|
|
constexpr scalar D0_3 = 0.115740740740741E-02;
|
|
constexpr scalar D0_4 = 0.352733686067019E-03;
|
|
constexpr scalar D0_5 = -0.178755144032922E-03;
|
|
constexpr scalar D0_6 = 0.391926317852244E-04;
|
|
// unused: constexpr scalar D0_7 = -0.218544851067999E-05;
|
|
// unused: constexpr scalar D0_8 = -0.185406221071516E-05;
|
|
// unused: constexpr scalar D0_9 = 0.829671134095309E-06;
|
|
// unused: constexpr scalar D0_10 = -0.176659527368261E-06;
|
|
// unused: constexpr scalar D0_11 = 0.670785354340150E-08;
|
|
// unused: constexpr scalar D0_12 = 0.102618097842403E-07;
|
|
// unused: constexpr scalar D0_13 = -0.438203601845335E-08;
|
|
|
|
constexpr scalar D1_0 = -0.185185185185185E-02;
|
|
constexpr scalar D1_1 = -0.347222222222222E-02;
|
|
constexpr scalar D1_2 = 0.264550264550265E-02;
|
|
constexpr scalar D1_3 = -0.990226337448560E-03;
|
|
constexpr scalar D1_4 = 0.205761316872428E-03;
|
|
// unused: constexpr scalar D1_5 = -0.401877572016461E-06;
|
|
// unused: constexpr scalar D1_6 = -0.180985503344900E-04;
|
|
// unused: constexpr scalar D1_7 = 0.764916091608111E-05;
|
|
// unused: constexpr scalar D1_8 = -0.161209008945634E-05;
|
|
// unused: constexpr scalar D1_9 = 0.464712780280743E-08;
|
|
// unused: constexpr scalar D1_10 = 0.137863344691572E-06;
|
|
// unused: constexpr scalar D1_11 = -0.575254560351770E-07;
|
|
// unused: constexpr scalar D1_12 = 0.119516285997781E-07;
|
|
|
|
constexpr scalar D2_0 = 0.413359788359788E-02;
|
|
constexpr scalar D2_1 = -0.268132716049383E-02;
|
|
// unused: constexpr scalar D2_2 = 0.771604938271605E-03;
|
|
// unused: constexpr scalar D2_3 = 0.200938786008230E-05;
|
|
// unused: constexpr scalar D2_4 = -0.107366532263652E-03;
|
|
// unused: constexpr scalar D2_5 = 0.529234488291201E-04;
|
|
// unused: constexpr scalar D2_6 = -0.127606351886187E-04;
|
|
// unused: constexpr scalar D2_7 = 0.342357873409614E-07;
|
|
// unused: constexpr scalar D2_8 = 0.137219573090629E-05;
|
|
// unused: constexpr scalar D2_9 = -0.629899213838006E-06;
|
|
// unused: constexpr scalar D2_10 = 0.142806142060642E-06;
|
|
|
|
const scalar u = 1/a;
|
|
scalar z = sqrt(2*phi);
|
|
|
|
if (lambda < 1)
|
|
{
|
|
z = -z;
|
|
}
|
|
|
|
if (sigma > (e0/sqrt(a)))
|
|
{
|
|
const scalar C0 =
|
|
D0_6*pow6(z) + D0_5*pow5(z) + D0_4*pow4(z)
|
|
+ D0_3*pow3(z) + D0_2*sqr(z) + D0_1*z + D0_0;
|
|
|
|
const scalar C1 =
|
|
D1_4*pow4(z) + D1_3*pow3(z) + D1_2*sqr(z) + D1_1*z + D1_0;
|
|
|
|
const scalar C2 = D2_1*z + D2_0;
|
|
|
|
return C2*sqr(u) + C1*u + C0;
|
|
}
|
|
else
|
|
{
|
|
const scalar C0 = D0_2*sqr(z) + D0_1*z + D0_0;
|
|
const scalar C1 = D1_1*z + D1_0;
|
|
const scalar C2 = D2_1*z + D2_0;
|
|
|
|
return C2*sqr(u) + C1*u + C0;
|
|
}
|
|
}
|
|
|
|
} // End namespace Foam
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
Foam::scalar Foam::Math::incGammaRatio_Q(const scalar a, const scalar x)
|
|
{
|
|
using namespace Foam::constant::mathematical;
|
|
|
|
#ifdef FULLDEBUG
|
|
if (a <= 0)
|
|
{
|
|
WarningInFunction
|
|
<< "The parameter (i.e. a) cannot be negative or zero"
|
|
<< " a = " << a
|
|
<< endl;
|
|
|
|
return std::numeric_limits<scalar>::infinity();
|
|
}
|
|
|
|
if (x < 0)
|
|
{
|
|
WarningInFunction
|
|
<< "The parameter (i.e. x) cannot be negative"
|
|
<< " x = " << x
|
|
<< endl;
|
|
|
|
return std::numeric_limits<scalar>::infinity();
|
|
}
|
|
#endif
|
|
|
|
constexpr scalar BIG = 14;
|
|
constexpr scalar x0 = 17;
|
|
constexpr scalar e0 = 0.025;
|
|
|
|
if (a < 1)
|
|
{
|
|
if (a == 0.5)
|
|
{
|
|
// (DM:Eq. 8)
|
|
if (x < 0.25)
|
|
{
|
|
return 1 - erf(sqrt(x));
|
|
}
|
|
else
|
|
{
|
|
return erfc(sqrt(x));
|
|
}
|
|
}
|
|
else if ( x < 1.1)
|
|
{
|
|
// (DM:Eq. 12)
|
|
scalar alpha = x/2.59;
|
|
|
|
if (x < 0.5)
|
|
{
|
|
alpha = log(sqrt(0.765))/log(x);
|
|
}
|
|
|
|
scalar sum = 0;
|
|
|
|
for (label n = 1; n <= 10; ++n)
|
|
{
|
|
sum += pow((-x), n)/((a + n)*factorial(n));
|
|
}
|
|
|
|
const scalar J = -a*sum;
|
|
|
|
if (a > alpha || a == alpha)
|
|
{
|
|
// (DM:Eq. 9)
|
|
return 1 - (pow(x, a)*(1 - J))/tgamma(a + 1);
|
|
}
|
|
else
|
|
{
|
|
// (DM:Eq. 10)
|
|
const scalar L = exp(a*log(x)) - 1;
|
|
const scalar H = 1/(tgamma(a + 1)) - 1;
|
|
|
|
return (pow(x, a)*J - L)/tgamma(a + 1) - H;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// (DM:Eq. 11)
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R*calcQE11(a, x);
|
|
}
|
|
}
|
|
else if (a >= BIG)
|
|
{
|
|
const scalar sigma = fabs(1 - x/a);
|
|
|
|
if (sigma <= e0/sqrt(a))
|
|
{
|
|
// (DM:Eq. 19)
|
|
const scalar lambda = x/a;
|
|
const scalar phi = lambda - 1 - log(lambda);
|
|
const scalar y = a*phi;
|
|
|
|
const scalar E = 0.5 - (1 - y/3)*sqrt(y/pi);
|
|
|
|
if (lambda <= 1)
|
|
{
|
|
return
|
|
1
|
|
- (
|
|
E
|
|
- (1 - y)/sqrt(2*pi*a)
|
|
*calcTE18(a, e0, x, lambda, sigma, phi)
|
|
);
|
|
}
|
|
else
|
|
{
|
|
return
|
|
E
|
|
+ (1 - y)/sqrt(2*pi*a)
|
|
*calcTE18(a, e0, x, lambda, sigma, phi);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (sigma <= 0.4)
|
|
{
|
|
// (DM:Eq. 17)
|
|
const scalar lambda = x/a;
|
|
const scalar phi = lambda - 1 - log(lambda);
|
|
const scalar y = a*phi;
|
|
|
|
if (lambda <= 1)
|
|
{
|
|
return
|
|
1
|
|
- (0.5*erfc(sqrt(y))
|
|
- exp(-y)/sqrt(2*pi*a)
|
|
*calcTE18(a, e0, x, lambda, sigma, phi));
|
|
}
|
|
else
|
|
{
|
|
return
|
|
0.5*erfc(sqrt(y))
|
|
+ exp(-y)/sqrt(2*pi*a)
|
|
*calcTE18(a, e0, x, lambda, sigma, phi);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (x <= max(a, log(10.0)))
|
|
{
|
|
// (DM:Eq. 15)
|
|
return 1 - calcPE15(a, x);
|
|
}
|
|
else if (x < x0)
|
|
{
|
|
// (DM:Eq. 11)
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R*calcQE11(a, x);
|
|
}
|
|
else
|
|
{
|
|
// (DM:Eq. 16)
|
|
return calcQE16(a, x);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (a > x || x >= x0)
|
|
{
|
|
if (x <= max(a, log(10.0)))
|
|
{
|
|
// (DM:Eq. 15)
|
|
return 1 - calcPE15(a, x);
|
|
}
|
|
else if ( x < x0)
|
|
{
|
|
// (DM:Eq. 11)
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R*calcQE11(a, x);
|
|
}
|
|
else
|
|
{
|
|
// (DM:Eq. 16)
|
|
return calcQE16(a, x);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (floor(2*a) == 2*a)
|
|
{
|
|
// (DM:Eq. 14)
|
|
if (floor(a) == a)
|
|
{
|
|
scalar sum = 0;
|
|
|
|
for (label n = 0; n <= (a - 1); ++n)
|
|
{
|
|
sum += pow(x, n)/factorial(n);
|
|
}
|
|
|
|
return exp(-x)*sum;
|
|
}
|
|
else
|
|
{
|
|
int i = a - 0.5;
|
|
scalar prod = 1;
|
|
scalar sum = 0;
|
|
|
|
for (int n = 1; n <= i; n++)
|
|
{
|
|
prod *= (n - 0.5);
|
|
sum += pow(x, n)/prod;
|
|
}
|
|
|
|
return erfc(sqrt(x)) + exp(-x)/sqrt(pi*x)*sum;
|
|
}
|
|
}
|
|
else if (x <= max(a, log(10.0)))
|
|
{
|
|
// (DM:Eq. 15)
|
|
return 1 - calcPE15(a, x);
|
|
}
|
|
else if (x < x0)
|
|
{
|
|
// (DM:Eq. 11)
|
|
const scalar R = (exp(-x)*pow(x, a))/tgamma(a);
|
|
|
|
return R*calcQE11(a, x);
|
|
}
|
|
else
|
|
{
|
|
// (DM:Eq. 16)
|
|
return calcQE16(a, x);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::Math::incGammaRatio_P(const scalar a, const scalar x)
|
|
{
|
|
return 1 - incGammaRatio_Q(a, x);
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::Math::incGamma_Q(const scalar a, const scalar x)
|
|
{
|
|
return incGammaRatio_Q(a, x)*tgamma(a);
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::Math::incGamma_P(const scalar a, const scalar x)
|
|
{
|
|
return incGammaRatio_P(a, x)*tgamma(a);
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|