erfInv.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) 2021 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "MathFunctions.H"
29 #include "mathematicalConstants.H"
30 #include "error.H"
31 
32 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
33 
34 Foam::scalar Foam::Math::erfInv(const scalar y)
35 {
36  #ifdef FULLDEBUG
37  if (mag(y) >= scalar(1))
38  {
40  << "The domain of inverse error function argument "
41  << "(i.e. y) should be limited to (-1, 1):" << nl
42  << " y = " << y
43  << endl;
44 
45  return std::numeric_limits<scalar>::infinity();
46  }
47  #endif
48 
49  // (W:p. 2) to reduce the max relative error to O(1e-4)
50  constexpr scalar a = 0.147;
51 
52  const scalar k =
53  scalar(2)/(a*constant::mathematical::pi) + 0.5*log(scalar(1) - sqr(y));
54 
55  const scalar h = log(scalar(1) - sqr(y))/a;
56 
57  // (W:Eq. 7)
58  const scalar x = sqrt(-k + sqrt(sqr(k) - h));
59 
60  if (y < 0)
61  {
62  return -x;
63  }
64 
65  return x;
66 }
67 
68 
69 // ************************************************************************* //
mathematicalConstants.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
MathFunctions.H
error.H
Foam::constant::universal::h
const dimensionedScalar h
Planck constant.
Definition: setRegionSolidFields.H:33
Foam::Math::erfInv
scalar erfInv(const scalar y)
Inverse error function of a real-number argument.
Definition: erfInv.C:34
Foam::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:262
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
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