SpecialFunctions Namespace Reference

Functions
double	bessi1 (double x)
double	bessk1 (double x)
double	beta (double a, double b)
double	beta (double x, double a, double b)
double	betaRegularized (double x, double a, double b)
double	dawson (double x)
double	DebyeD (int n, double x)
double	EllipticE (double k)
double	EllipticE (double x, double k)
double	EllipticF (double x, double k)
double	EllipticK (double k)
double	EllipticPi (double nu, double k)
double	EllipticPi (double x, double nu, double k)
double	functionX (double s)
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)
double	rc (double x, double y)
double	rd (double x, double y, double z)
double	rf (double x, double y, double z)
double	rj (double x, double y, double z, double p)

Detailed Description

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.

Function Documentation

bessi1()

double SpecialFunctions::bessi1

(

double

x

)

This function returns the modified Bessel function $I_1(x)$ of the first kind of integer order 1. The implementation is based on a polynomial approximation for small $x$ and the product of a polynomial approximation and the exponential factor at large $x$. The implementation is adapted from Numerical Recipes.

bessk1()

double SpecialFunctions::bessk1

(

double

x

)

This function returns the modified Bessel function $K_1(x)$ of the second kind of integer order 1. The implementation is based on a polynomial approximation for small $x$ and the product of a polynomial approximation and the exponential factor at large $x$. The implementation is adapted from Numerical Recipes.

beta() [1/2]

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)={\displaystyle \frac{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}{\mathrm{\Gamma }(a+b)}}$$

beta() [2/2]

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.

betaRegularized()

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.

dawson()

double SpecialFunctions::dawson

(

double

x

)

Dawson's integral, defined as

$$D(x)=\exp (-x^2){\int }_0^x\exp (t^2)dt.$$

The implementation is based on the series

$$D(x)=\underset{h\to 0}{lim}\frac{1}{\sqrt{\pi }}\sum _{n\text{ odd}}\frac{e^{-(x-nh{)}^2}}{n},$$

which guarantees an exponential increase in accuracy. The implementation is adapted from Numerical Recipes, where additional tricks are included to make the computation even faster.

DebyeD()

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.

EllipticE() [1/2]

double SpecialFunctions::EllipticE

(

double

k

)

This function returns the complete elliptic integral of the second kind, defined as

$$E(k)={\int }_0^1\frac{\sqrt{1-k^2{t}^2}dt}{\sqrt{1-t^2}},$$

where $0\le k\le 1$. The function is no more than a function call to the incomplete elliptic integral of the second kind $E(x,k)$ with $x=1$. See also the note about argument conventions there.

EllipticE() [2/2]

double SpecialFunctions::EllipticE	(	double	x,
		double	k
	)

This function returns the incomplete elliptic integral of the second kind, defined as

$$E(x,k)={\int }_0^x\frac{\sqrt{1-k^2{t}^2}dt}{\sqrt{1-t^2}}.$$

The arguments must satisfy $0\le x\le 1$ and $0\le k\le 1$. For the evaluation of this function, we express it in terms of Carlson's elliptic integrals of the first and second kind,

$$E(x,k)=x{R}_F(1-x^2,1-k^2{x}^2,1)-\frac{1}{3}{x}^3{k}^2{R}_D(1-x^2,1-k^2{x}^2,1).$$

It is important to note that different incompatible argument conventions are being used for the elliptic integrals. Instead of the coordinate $x$ one often uses the amplitude $\varphi$, where both are related through $x=\sin \varphi$. For the second parameter, one can use the modulus $k$ (used in SKIRT), the parameter $m$ or the modular angle $\alpha$. These are related through $m=k^2={\sin }^2\alpha$. We use the conventions adopted in Maple.

EllipticF()

double SpecialFunctions::EllipticF	(	double	x,
		double	k
	)

This function returns the incomplete elliptic integral of the first kind, defined as

$$F(x,k)={\int }_0^x\frac{dt}{\sqrt{(1-t^2)(1-k^2{t}^2)}}.$$

The arguments must satisfy $0\le x\le 1$ and $0\le k\le 1$. For the evaluation of this function, we express it in terms of Carlson's elliptic integral of the first kind,

$$F(x,k)=x{R}_F(1-x^2,1-k^2{x}^2,1).$$

It is important to note that different incompatible argument conventions are being used for the elliptic integrals. Instead of the coordinate $x$ one often uses the amplitude $\varphi$, where both are related through $x=\sin \varphi$. For the second parameter, one can use the modulus $k$ (used in SKIRT), the parameter $m$ or the modular angle $\alpha$. These are related through $m=k^2={\sin }^2\alpha$. We use the conventions adopted in Maple.

EllipticK()

double SpecialFunctions::EllipticK

(

double

k

)

This function returns the complete elliptic integral of the first kind, defined as

$$K(k)={\int }_0^1\frac{dt}{\sqrt{(1-t^2)(1-k^2{t}^2)}},$$

where $0\le k\le 1$. The function is no more than a function call to the incomplete elliptic integral of the first kind $F(x,k)$ with $x=1$. See also the note about argument conventions there.

EllipticPi() [1/2]

double SpecialFunctions::EllipticPi	(	double	nu,
		double	k
	)

This function returns the complete elliptic integral of the third kind, defined as

$$\mathrm{\Pi }(\nu ,k)={\int }_0^1\frac{dt}{(1-\nu t^2)\sqrt{(1-t^2)(1-k^2{t}^2)}},$$

where $0\le k\le 1$. The function is no more than a function call to the incomplete elliptic integral of the third kind $\mathrm{\Pi }(x,\nu ,k)$ with $x=1$. See also the note about argument conventions there.

EllipticPi() [2/2]

double SpecialFunctions::EllipticPi	(	double	x,
		double	nu,
		double	k
	)

This function returns the incomplete elliptic integral of the third kind, defined as

$$\mathrm{\Pi }(x,\nu ,k)={\int }_0^x\frac{dt}{(1-\nu t^2)\sqrt{(1-t^2)(1-k^2{t}^2)}}.$$

The arguments must satisfy $0\le x\le 1$ and $0\le k\le 1$. For the evaluation of this function, we express it in terms of Carlson's elliptic integrals of the first and third kind,

$$\mathrm{\Pi }(x,\nu ,k)=x{R}_F(1-x^2,1-k^2{x}^2,1)+\frac{1}{3}\nu x^3{R}_J(1-x^2,1-k^2{x}^2,1,1-\nu x^2).$$

It is important to note that different incompatible argument conventions are being used for the elliptic integrals. Instead of the coordinate $x$ one often uses the amplitude $\varphi$, where both are related through $x=\sin \varphi$. For the second parameter, one can use the modulus $k$ (used in SKIRT), the parameter $m$ or the modular angle $\alpha$. These are related through $m=k^2={\sin }^2\alpha$. Finally, there are different sign conventions for the characteristic $\nu$, which is sometimes defined with a different signature (e.g. Abramowitz & Stegun 1965). We use the conventions adopted in Maple.

functionX()

double SpecialFunctions::functionX

(

double

s

)

This function returns the value of the function defined as

$$X(s)=\{\begin{array}{ll}{\displaystyle \frac{1}{\sqrt{1-s^2}}}\text{arccosh}\left({\displaystyle \frac{1}{s}}\right)& \text{if }0<s<1\text{,}\\ {\displaystyle \frac{1}{\sqrt{s^2-1}}}\text{arccos}\left({\displaystyle \frac{1}{s}}\right)& \text{if }s>1\text{.}\end{array}$$

This function is completely smooth (it is the same function if one considers it in the complex plane).

gamma()

double SpecialFunctions::gamma

(

double

a

)

This function returns the Gamma function $\mathrm{\Gamma }(a)$.

gexp()

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{array}{ll}((1-p)x+1{)}^{\frac{1}{1-p}}& p\ne 1\\ {\text{e}}^p& p=1\end{array}$$

When $p$ approaches 1, we use the expansion

$$\text{gexp}(p,x)\approx {\text{e}}^x[1-\frac{1}{2}(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+\dots ]$$

gln()

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{array}{ll}{\displaystyle \frac{x^{1-p}-1}{1-p}}& p\ne 1\\ \ln x& p=1\end{array}.$$

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+\frac{1}{2}(1-p){\ln }^{2}x+\frac{1}{6}(1-p{)}^2{\ln }^{3}x+\frac{1}{24}(1-p{)}^3{\ln }^{4}x+\dots$$

gln2()

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}).$$

incompleteGamma()

double SpecialFunctions::incompleteGamma	(	double	a,
		double	x
	)

This function returns the (regularized) incomplete gamma function

$$P(a,x)=\frac{\gamma (a,x)}{\mathrm{\Gamma }(a)}=\frac{1}{\mathrm{\Gamma }(a)}{\int }_0^x{\text{e}}^{-t}{t}^{a-1}\text{d}t.$$

The implementation is adapted from the Numerical Recipes book.

LambertW()

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\to 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.

LambertW1()

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\to 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}$.

lngamma()

double SpecialFunctions::lngamma

(

double

a

)

This function returns the logarithm of the Gamma function, i.e. $\ln \mathrm{\Gamma }(a)$. The implementation is taken from the Numerical Recipes book.

lnmean() [1/2]

double SpecialFunctions::lnmean	(	double	x1,
		double	x2
	)

This function returns the logarithmic mean $M(x_1,x_2)$ of two nonnegative values $x_1\ge 0$ and $x_2\ge 0$, defined as

$$M(x_1,x_2)=\{\begin{array}{ll}0& \mathrm{i}\mathrm{f}{x}_1=0\mathrm{o}\mathrm{r}{x}_2=0\\ x_1& \mathrm{i}\mathrm{f}{x}_1=x_2\\ {\displaystyle \frac{x_2-x_1}{\ln x_2-\ln x_1}}& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}$$

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)}\le M(x_1,x_2)\le \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={\displaystyle \frac{x_2}{x_1}}-1$ and the following expansion:

$$M(x_1,x_2)=x_1\frac{x}{\ln (1+x)}={\displaystyle \frac{x_1}{1-{\displaystyle \frac{x}{2}}+{\displaystyle \frac{x^2}{3}}-{\displaystyle \frac{x^3}{4}}+...}}$$

lnmean() [2/2]

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.

rc()

double SpecialFunctions::rc	(	double	x,
		double	y
	)

This function returns Carlson's degenerate elliptic integral $R_C$, defined as

$$R_C(x,y)=\frac{1}{2}{\int }_0^{\mathrm{\infty }}\frac{dt}{(t+y)\sqrt{t+x}}.$$

Here $x$ must be nonnegative and $y$ must be nonzero. The implementation, based on the duplication theorem, is taken from Numerical Recipes.

rd()

double SpecialFunctions::rd	(	double	x,
		double	y,
		double	z
	)

This function returns Carlson's elliptic integral $R_D$, defined as

$$R_D(x,y,z)=\frac{3}{2}{\int }_0^{\mathrm{\infty }}\frac{dt}{\sqrt{(t+x)(t+y)(t+z{)}^3}}.$$

The first two arguments $x$ and $y$ must be nonnegative and at most one of them is zero, whereas the last argument $z$ must be positive. The implementation, based on the duplication theorem, is taken from Numerical Recipes.

rf()

double SpecialFunctions::rf	(	double	x,
		double	y,
		double	z
	)

This function returns Carlson's elliptic integral $R_F$, defined as

$$R_F(x,y,z)=\frac{1}{2}{\int }_0^{\mathrm{\infty }}\frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}.$$

The three arguments $x$, $y$ and $z$ must be nonnegative and at most one of them is zero. The implementation of this function is based on the duplication theorem and is taken from Numerical Recipes.

rj()

double SpecialFunctions::rj	(	double	x,
		double	y,
		double	z,
		double	p
	)

This function returns Carlson's elliptic integral $R_J$, defined as

$$R_J(x,y,z,p)=\frac{3}{2}{\int }_0^{\mathrm{\infty }}\frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}.$$

The first three arguments $x$, $y$ and $z$ must be nonnegative and at most one of them is zero, whereas the last argument $p$ must be nonzero. The implementation, based on the duplication theorem, is taken from Numerical Recipes.