 |
The SKIRT project
advanced radiative transfer for astrophysics
|
|
The SKIRT project
advanced radiative transfer for astrophysics
|
Functions | |
| double | beta (double a, double b) |
| double | beta (double x, double a, double b) |
| double | betaRegularized (double x, double a, double b) |
| double | DebyeD (int n, double x) |
| double | gamma (double a) |
| double | gexp (double p, double x) |
| double | gln (double p, double x) |
| double | gln2 (double p, double x1, double x2) |
| double | incompleteGamma (double a, double x) |
| double | LambertW (double z) |
| double | LambertW1 (double z) |
| double | lngamma (double a) |
| double | lnmean (double x1, double x2) |
| double | lnmean (double x1, double x2, double lnx1, double lnx2) |
This namespace contains a number of special mathematical functions. It includes the usual mathematical special functions such as the gamma function and elliptic integrals, but also a number of specific functions that often occur in the program. Many functions are copied or adapted from the famous Numerical Recipes books. The version used is the second edition of Numerical Recipes in C++ (Press et al. 2002). Other functions are taken from the web or are just implemented ad hoc.
| double SpecialFunctions::beta | ( | double | a, |
| double | b ) |
This function returns the complete Beta function \(B(a,b)\), defined as
\[ B(a,b) = \int_0^1 t^{a-1}\,(1-t)^{b-1}\,dt, \]
with the two parameters \(a\) and \(b\) positive real numbers. The implementation is the trivial translation of the formula
\[ B(a,b) = \dfrac{\Gamma(a)\,\Gamma(b)}{\Gamma(a+b)} \]
| double SpecialFunctions::beta | ( | double | x, |
| double | a, | ||
| double | b ) |
This function returns the incomplete Beta function \(B_x(a,b)\), defined as
\[ B_x(a,b) = \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt, \]
with the parameters \(a\) and \(b\) positive real numbers and the argument \(x\) between 0 and 1.
| double SpecialFunctions::betaRegularized | ( | double | x, |
| double | a, | ||
| double | b ) |
This function returns the regularized incomplete Beta function \(I_x(a,b)\), defined as
\[ I_x(a,b) = \frac{B_x(a,b)}{B(a,b)} = \frac{1}{B(a,b)} \int_0^x t^{a-1}\,(1-t)^{b-1}\,dt, \]
with the parameters \(a\) and \(b\) positive real numbers and the argument \(x\) between 0 and 1. The implementation uses the representation of this function in terms of continued fractions, and is adapted from the Numerical Recipes book.
| double SpecialFunctions::DebyeD | ( | int | n, |
| double | x ) |
This function returns the Debye function of order \(n\), defined as
\[ D_n(x) = \frac{n}{x^n} \int_0^x \frac{t^n\,{\text{d}}t}{{\text{e}}^t-1}. \]
The order \(n\) should be an integer number between 1 and 20; the argument \(x\) is a positive real number. The implementation is adapted from Richard Mathar's website at MPIA Heidelberg.
| double SpecialFunctions::gamma | ( | double | a | ) |
This function returns the Gamma function \(\Gamma(a)\).
| double SpecialFunctions::gexp | ( | double | p, |
| double | x ) |
This function returns a generalized exponential function \({\text{gexp}}(p,x)\), defined as the inverse of the function gln. In formula it is defined as
\[ {\text{gexp}}(p,x) = \begin{cases} \; ((1-p)\,x+1)^{\frac{1}{1-p}} & p\ne1 \\ \; {\text{e}}^p & p=1 \end{cases} \]
When \(p\) approaches 1, we use the expansion
\[ {\text{gexp}}(p,x) \approx {\text{e}}^x\left[ 1 - \frac12\,(1-p)\,x^2 + \frac{1}{24}\, (1-p)^2\, (3x+8)\,x^3 - \frac{1}{48}\, (1-p)^3\, (x^2+8x+12)\, x^4 + \ldots \right] \]
| double SpecialFunctions::gln | ( | double | p, |
| double | x ) |
This function returns a generalized logarithmic function \({\text{gln}}(p,x)\), defined for \(x>0\) and arbitrary real \(p\) as
\[ {\text{gln}}(p,x) = \int_1^x t^{-p}\,{\text{d}}t = \begin{cases} \; \dfrac{x^{1-p}-1}{1-p} & p\ne1 \\ \; \ln x & p=1 \end{cases}. \]
This function is included in this library of special functions because the power law expression diverges as \(p\) approaches 1. In that case we can use the expansion
\[ {\text{gln}}(p,x) \approx \ln x + \frac12\,(1-p)\ln^2x + \frac16\,(1-p)^2 \ln^3x + \frac{1}{24}\,(1-p)^3\ln^4x + \ldots \]
| double SpecialFunctions::gln2 | ( | double | p, |
| double | x1, | ||
| double | x2 ) |
This function returns the difference between two values of the generalized logarithmic function \({\text{gln}}(p,x)\) with the same exponent \(p\). Compared to simply subtracting the two generalized logarithms, this function achieves much better accuracy for large arguments \(x_1\) and \(x_2\) by using the identity
\[ {\text{gln2}}(p,x_1,x_2) = {\text{gln}}(p,x_1) - {\text{gln}}(p,x_2) = (x_2)^{1-p} \, {\text{gln}}(p,\frac{x_1}{x_2}). \]
| double SpecialFunctions::incompleteGamma | ( | double | a, |
| double | x ) |
This function returns the (regularized) incomplete gamma function
\[ P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)} \int_0^x {\text{e}}^{-t}\, t^{a-1}\, {\text{d}}t. \]
The implementation is adapted from the Numerical Recipes book.
| double SpecialFunctions::LambertW | ( | double | z | ) |
This function is the Lambert \(W\) function, also known as the product log function for the principle branch. The Lambert \(W\) function is generally defined as the inverse of the function
\[ w \rightarrow f(w) = w\,e^w. \]
For real \(z\) in the interval \([-1/e,0[\), the equation
\[ z=w\,e^w \]
has two real solutions, one smaller than and one greater than \(-1\). These define the two only real branches of this complex function: the principle branch \(W \equiv W_0\) returns a value \(W(z)>-1\), whereas the other branch, labelled \(W_{-1}\) returns a value \(W_{-1}(z)<-1\). For more details on this function, see the help facility in Maple on the command LambertW or in Mathematica on the command ProductLog. The adopted implementation for the function \(W_0(z)\) is a combination of the C-code available from the homepage of Keith Briggs and the code from the GNU scientific library GSL.
| double SpecialFunctions::LambertW1 | ( | double | z | ) |
This function is the Lambert \(W\) function, also known as the product log function, for the branch \(-1\). The Lambert \(W\) function is generally defined as the inverse of the function
\[ w \rightarrow f(w) = w\,e^w. \]
For real \(z\) in the interval \([-1/e,0[\), the equation
\[ z=w\,e^w \]
has two real solutions, one smaller than and one greater than \(-1\). These define the two only real branches of this complex function: the principle branch \(W \equiv W_0\) returns a value \(W(z)>-1\), whereas the other branch, labelled \(W_{-1}\) returns a value \(W_{-1}(z)<-1\). For more details on this function, see the help facility in Maple on the command LambertW or in Mathematica on the command ProductLog. The adopted implementation for the function \(W_{-1}(z)\) is a combination of the C-code available from the homepage of Keith Briggs and the code from the GNU scientific library GSL. The code returns an error message and the program is halted when the input value does not lie within the interval \([-1/e,0[\). Comparison of the results of this implementation with the results from the Maple function yields identical values. Only for values of \(z\) extremely close to the limiting value zero, we get some minor discrepancies, e.g. the difference between the two results is 0.015% for \(z=-10^{-9}\).
| double SpecialFunctions::lngamma | ( | double | a | ) |
This function returns the logarithm of the Gamma function, i.e. \(\ln\Gamma(a)\). The implementation is taken from the Numerical Recipes book.
| double SpecialFunctions::lnmean | ( | double | x1, |
| double | x2 ) |
This function returns the logarithmic mean \(M(x_1,x_2)\) of two nonnegative values \(x_1\geq 0\) and \(x_2\geq 0\), defined as
\[ M(x_1,x_2) = \begin{cases} 0 & \mathrm{if}\, x_1=0 \,\mathrm{or}\, x_2=0 \\ x_1 & \mathrm{if}\, x_1=x_2 \\ \dfrac{x_2 - x_1}{\ln x_2 - \ln x_1}&\mathrm{otherwise}.\end{cases} \]
The function is invariant for swapping of it arguments, and it can be proven that the logarithmic mean lies between the geometric mean and the arithmetic mean,
\[ \sqrt{(x_1 x_2)} \leq M(x_1,x_2) \leq \frac{x_1+x_2}{2}, \]
where the equality realizes for \(x_1=x_2\).
The function is implemented here because the quotient becomes numerically unstable when \(x_1\approx x_2\). In that case, the function value is calculated using the substitution \(x=\dfrac{x_2}{x_1}-1\) and the following expansion:
\[ M(x_1,x_2) = x_1 \, \frac{x}{\ln(1+x)} = \dfrac{x_1}{1 - \dfrac{x}{2} + \dfrac{x^2}{3} - \dfrac{x^3}{4} + ...} \]
| double SpecialFunctions::lnmean | ( | double | x1, |
| double | x2, | ||
| double | lnx1, | ||
| double | lnx2 ) |
This function returns the logarithmic mean of two values as described for the two-argument lnmean() function in this class, given also the natural logarithm of these two values. In cases where these logarithms are available at the call site anyway, this function is more efficient than its two-argument equivalent.