Hex  1.0 Hydrogen-electron collision solver
special.h File Reference
#include <cmath>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include "arrays.h"
#include "complex.h"
#include "misc.h"
Include dependency graph for special.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

## Data Structures

Base class for radial functions. More...

## Namespaces

special

special::constant

special::integral

## Macros

#define fac(n)   gsl_sf_fact(n)
factorial More...

#define sqr(x)   (gsl_sf_pow_int((x),2))
second power More...

#define ALPHA_PLUS(n, l)   (sqrt(((n)-(l))*((n)+(l)+1)))
Shifting coefficients for Sturmian T-operators. More...

#define ALPHA_MINUS(n, l)   (sqrt(((n)-(l)-1)*((n)+(l))))

#define DELTA(a, b)   ((a) == (b) ? 1 : 0)
Kronecker delta. More...

## Functions

template<class T >
special::integral::trapz (double h, NumberArray< T > const &y)
Uniform trapezoidal integration. More...

template<class T >
special::integral::trapz (NumberArray< T > const &x, NumberArray< T > const &y)
Non-uniform trapezoidal integration. More...

template<class T >
special::integral::simpson (double h, NumberArray< T > const &y)
Uniform Simpson integration. More...

template<class T >
special::integral::pow_exp_hyperg1F1 (T a, T b, T c, T u, T v, T x, double epsrel=1e-10, unsigned maxiter=1000)
Compute integral of the confluent hypergeometric function. More...

template<class T >
NumberArray< T > special::integral::romberg (const ArrayView< T > y)
Romberg integration. More...

int special::coulomb_zeros (double eta, int L, int nzeros, double *zeros, double epsrel=1e-8)
Get zeros of the Coulomb wave function $$F_L(-1/k,kr)$$. More...

std::vector< std::vector< int > > special::FdB_partition (int n)
Faa di Bruno partitioning. More...

template<class T , class OuterFunctionDerivative , class InnerFunctionDerivative >
special::chain_rule (OuterFunctionDerivative Df, InnerFunctionDerivative Dg, int n, double x)
Chain rule for n-th derivative. More...

Complex hydro_P (unsigned n, unsigned l, Complex z)

Complex dhydro_P (unsigned n, unsigned l, Complex z)

Complex ric_j (int n, Complex z)

cArray ric_jv (int lmax, Complex z)

Complex dric_j (int n, Complex z)

cArray dric_jv (int lmax, Complex z)

double ric_j (int n, double x)
Ricatti-Bessel function of the first kind, $$\hat{j}_n(x)$$. More...

double ric_n (int n, double x)
Ricatti-Bessel function of the second kind, $$\hat{y}_n(x)$$. More...

double sph_i_scaled (int n, double x)
Scaled modified spherical Bessel function of the first kind, $$\mathrm{e}^{-x} i_n(x)$$. More...

double ric_i_scaled (int n, double x)
Scaled modified Ricatti-Bessel function of the first kind, $$\mathrm{e}^{-x} \hat{i}_n(x)$$. More...

double sph_i (int n, double x)
Modified spherical Bessel function of the first kind, $$i_n(x)$$. More...

double ric_i (int n, double x)
Modified Ricatti-Bessel function of the first kind, $$\hat{i}_n(x)$$. More...

double sph_k_scaled (int n, double x)
Scaled modified spherical Bessel function of the second kind, $$\mathrm{e}^x k_n(x)$$. More...

double ric_k_scaled (int n, double x)
Scaled modified Ricatti-Bessel function of the second kind, $$\mathrm{e}^x \hat{k}_n(x)$$. More...

double sph_k (int n, double x)
Modified spherical Bessel function of the second kind, $$k_n(x)$$. More...

double ric_k (int n, double x)
Modified Ricatti-Bessel function of the second kind, $$\hat{k}_n(x)$$. More...

Complex ric_h_plus (int n, double x)
Ricatti-Hankel function of the first kind, $$\hat{h}_n^{(+)}(x)$$. More...

Complex sphY (int l, int m, double theta, double phi)
Spherical harmonic function. More...

Complex sphBiY (int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2)
Bi-polar spherical harmonic function. More...

double dric_j (int n, double x)
Derivative of Ricatti-Bessel function of the first kind, $$\hat{j}_n'(x)$$. More...

double dric_n (int n, double x)
Derivative of Ricatti-Bessel function of the second kind, $$\hat{y}_n'(x)$$. More...

double dric_i (int n, double x)
Derivative of the modified Ricatti-Bessel function of the second kind, $$\hat{i}_n'(x)$$. More...

double dric_k (int n, double x)
Derivative of the modified Ricatti-Bessel function of the second kind, $$\hat{k}_n'(x)$$. More...

double dric_i_scaled (int n, double x)
Derivative of scaled modified Ricatti-Bessel function of the second kind, $$(\mathrm{e}^{-x} \hat{i}_n(x))'$$. More...

double dric_k_scaled (int n, double x)
Scaled derivative of modified Ricatti-Bessel function of the second kind, $$(\mathrm{e}^x \hat{k}_n(x))$$. More...

Complex dric_h_plus (int n, double x)
Derivative of Ricatti-Hankel function of the first kind, $${\hat{h}_n^{(+)}}'(x)$$. More...

int coul_F_michel (int l, double k, double r, double &F, double &Fp)
Uniform approximation to the Coulomb wave function. More...

int coul_F (int l, double k, double r, double &F, double &Fp)
Evaluate Coulomb wave function (and its derivative). More...

double coul_F_asy (int l, double k, double r, double sigma=special::constant::Nan)
Asymptotic form of the regular Coulomb wave. More...

double coul_F_sigma (int l, double k)
Coulomb phase shift. More...

bool makes_triangle (int two_j1, int two_j2, int two_j3)
Check triangle inequality. More...

double logdelta (int two_j1, int two_j2, int two_j3)
Auxiliary function used in coupling coefficient computations. More...

double computef (int lambda, int l1, int l2, int l1p, int l2p, int L)

double ClebschGordan (int l1, int m1, int l2, int m2, int L, int M)
Clebsch-Gordan coefficient. More...

double Gaunt (int l1, int m1, int l2, int m2, int l, int m)
Gaunt's integral. More...

int triangle_count (int L, int maxell)

double Hyper2F1 (double a, double b, double c, double x)
The hypergeometric function $${}_2F_1$$. More...

## Variables

const double special::constant::e = M_E

const double special::constant::pi = M_PI

const double special::constant::pi_half = M_PI_2

const double special::constant::pi_quart = M_PI_4

const double special::constant::pi_sqrt = M_SQRTPI

const double special::constant::pi_inv = M_1_PI

const double special::constant::pi_two_inv = M_2_PI

const double special::constant::pi_two_invsqrt = M_2_SQRTPI

const double special::constant::sqrt_two = M_SQRT2

const double special::constant::sqrt_three = M_SQRT3

const double special::constant::sqrt_half = M_SQRT1_2

const double special::constant::log2e = M_LOG2E

const double special::constant::log10e = M_LOG10E

const double special::constant::ln2 = M_LN2

const double special::constant::ln10 = M_LN10

const double special::constant::lnpi = M_LNPI

const double special::constant::euler = M_EULER

const double special::constant::Inf = std::numeric_limits<double>::infinity()

const double special::constant::Nan = std::numeric_limits<double>::quiet_NaN()

#define Wigner3j(a, b, c, d, e, f)   Wigner3j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))

double Wigner3j_2 (int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
Wigner 3j coefficient. More...

#define Wigner6j(a, b, c, d, e, f)   Wigner6j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))

double Wigner6j_2 (int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6)
Wigner 6j coefficient. More...

## Macro Definition Documentation

 #define ALPHA_MINUS ( n, l ) (sqrt(((n)-(l)-1)*((n)+(l))))
 #define ALPHA_PLUS ( n, l ) (sqrt(((n)-(l))*((n)+(l)+1)))
 #define DELTA ( a, b ) ((a) == (b) ? 1 : 0)
 #define fac ( n ) gsl_sf_fact(n)
 #define sqr ( x ) (gsl_sf_pow_int((x),2))
 #define Wigner3j ( a, b, c, d, e, f ) Wigner3j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))
 #define Wigner6j ( a, b, c, d, e, f ) Wigner6j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))

## Function Documentation

 double ClebschGordan ( int l1, int m1, int l2, int m2, int L, int M )
Note
In present implementation valid only for integer (not half-integer) angular momenta. [Otherwise one needs to correct the signs.]
 double computef ( int lambda, int l1, int l2, int l1p, int l2p, int L )
Returns
Value of $$f(\lambda,l_1,l_2,l_1',l_2',L)$$.
 int coul_F ( int l, double k, double r, double & F, double & Fp )
Parameters
 l Angular momentum. k Wavenumber. r Radial coordinate. F Output reference for resulting value. Fp Output reference for resulting derivative.
 double coul_F_asy ( int l, double k, double r, double sigma = special::constant::Nan )
Parameters
 l Angular momentum. k Wavenumber. r Radial coordinate. sigma Optionally, the precomputed Coulomb phase shift.
 int coul_F_michel ( int l, double k, double r, double & F, double & Fp )

This routine uses algorithm from the following article: Michel N, Uniform WKB approximation of Coulomb wave functions for arbitrary partial wave, EPL, 83 (2008) 10002.

The method is asymptotically valid for high energies and partial waves.

Parameters
 l Angular momentum. k Wavenumber. r Radial coordinate. F Output reference for resulting value. Fp Output reference for resulting derivative.
 double coul_F_sigma ( int l, double k )
Parameters
 l Angular momentum. k Wavenumber.
 Complex dhydro_P ( unsigned n, unsigned l, Complex z )

Derivative of Pnl.

Parameters
 n l z
 Complex dric_h_plus ( int n, double x )
inline
 double dric_i ( int n, double x )
inline
 double dric_i_scaled ( int n, double x )
inline
 Complex dric_j ( int n, Complex z )

Derivative of Riccati-Bessel function

Parameters
 n Degree of the function. z Complex argument.
 double dric_j ( int n, double x )
inline
 cArray dric_jv ( int lmax, Complex z )

Vectorized interface for dric_j.

Returns values of all derivatives of Riccati-Bessel functions of order less than or equal to lmax.

Parameters
 lmax Angular momentum limit. z Complex argument.
 double dric_k ( int n, double x )
inline
 double dric_k_scaled ( int n, double x )
inline
 double dric_n ( int n, double x )
inline
 double Gaunt ( int l1, int m1, int l2, int m2, int l, int m )

Computes the integral of three spherical harmonic functions.

$\int_{4\pi} Y_{l_1m_1} Y_{l_2m_2} Y^{\ast}_{lm} \mathrm{d}\Omega \ .$

 Complex hydro_P ( unsigned n, unsigned l, Complex z )

Parameters
 n Principal quantum number. l Orbital quantum number. z Radius: real or complex argument.
 double Hyper2F1 ( double a, double b, double c, double x )

Evaluates the Gauss hypergeometric function

${}_2F_1(a,b;c;x)$

for arbitrary real argument $$x$$. It uses the routine gsl_sf_hyperg_2F1 from GSL library, which is defined for |x| < 1. Using the transformations Abramowitz & Stegun 15.3.3-9 the hypergeometric function is transformed so that is can be evaluated by the library function.

Todo:
Document transformations.
 double logdelta ( int two_j1, int two_j2, int two_j3 )

The triangle function, defined as

$\Delta(j_1,j_2,j_3) = \sqrt{ \frac{(j_1+j_2-j_3)!(j_1-j_2+j_3)!(-j_1+j_2+j_3)!}{(j_1+j_2+j_3+1)!} } \ .$

 bool makes_triangle ( int two_j1, int two_j2, int two_j3 )

Verify that the three angular momenta satisfy triangle inequalities, i.e.

$|j_1 - j_2 | \le j_3 \le j_1 + j_2$

and analogously for all (cyclical) permutations of indices.

 Complex ric_h_plus ( int n, double x )
inline

$$\hat{h}_0^{(+)}(x) = -\mathrm{i} \exp (\mathrm{i}x)$$.

 double ric_i ( int n, double x )
inline

$$\hat{i}_0(x) = \sinh x$$

 double ric_i_scaled ( int n, double x )
inline

$$\mathrm{e}^{-x} \hat{i}_0(x) = \mathrm{e}^{-x} \sinh x$$

 Complex ric_j ( int n, Complex z )

Riccati-Bessel function

Evaluate Riccati-Bessel function for complex argument. Function is not suitable for large degrees, it uses the most naïve (and least stable) evaluation method. Starting from the expressions for zeroth and first Riccati-Bessel function

$j_0(z) = \sin z, \qquad j_1(z) = \frac{\sin z}{z} - \cos z$

the function employs the forward(!) recurrence relation

$j_{n+1}(z) = \frac{2n+1}{z} j_n(z) - j_{n-1}(z) .$

Note
Forward recurrence is unstable for small arguments. In the present implementation those are expected to be real. For real arguments the GSL routine is called, which doesn't suffer from this instability.
Parameters
 n Degree of the Riccati-Bessel function. z Complex argument.
 double ric_j ( int n, double x )
inline

$$\hat{j}_0(x) = \sin x$$

 cArray ric_jv ( int lmax, Complex z )

Vectorized interface for ric_j.

Returns values of all Riccati-Bessel functions of order less than or equal to lmax.

Parameters
 lmax Angular momentum limit. z Complex argument.
 double ric_k ( int n, double x )
inline

$$\hat{k}_0(x) = \mathrm{e}^{-x}$$.

 double ric_k_scaled ( int n, double x )
inline

$$\mathrm{e}^{x} \hat{k}_0(x) = 1$$.

 double ric_n ( int n, double x )
inline

$$\hat{n}_0(x) = \cos x$$

 double sph_i ( int n, double x )
inline

$$i_0(x) = \frac{1}{x} \sinh x$$

 double sph_i_scaled ( int n, double x )
inline

$$\mathrm{e}^{-x} i_0(x) = \mathrm{e}^{-x} \frac{1}{x} \sinh x$$

 double sph_k ( int n, double x )
inline

$$k_0(x) = \frac{1}{x} \mathrm{e}^{-x}$$.

 double sph_k_scaled ( int n, double x )
inline

$$\mathrm{e}^{x} k_0(x) = \frac{1}{x}$$.

 Complex sphBiY ( int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2 )
 Complex sphY ( int l, int m, double theta, double phi )
 int triangle_count ( int L, int maxell )

Compute number of angular momenta pairs $$\ell_1$$ and $$\ell_2$$ that are less than or equal to "maxell" and compose the total angular momentum $$L$$.

 double Wigner3j_2 ( int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3 )

Compute the Wigner 3j coefficient. Even though there is a routine "gsl_sf_coupling_3j" in GSL library, it has to be implemented anew, because of factorial overflows. This routine computes all factorials only in logarithms using the standard function "lgamma". The formula is taken from Edmonds, A. R.: Angular momentum in quantum mechanics, Princeton 1968.

$\left( \matrix{j_1 7 j_2 & j_3 \cr m_1 & m_2 & m_3} \right) = \epsilon(j_1,j_2,j_3)\Delta(j_1,j_2,j_3) \delta_{m_1+m_2+m_3}^0 (-1)^{j_1-j_2-m_3} \sqrt{(j_1+m_1)! (j_1-m_1)! (j_2+m_2)! (j_2-m_2)! (j_3+m_3)! (j_3-m_3)!} \sum_k \frac{(-1)^2}{k! (j_1+j_2-j_3-k)! (j_1-m_1-k)! (j_2+m_2-k)! (j_3-j_2+m_1+k)! (j_3-j_1-m_2+k)!} \ ,$

$\epsilon(j_1,j_2,j_3) = \cases{1 & triangle inequality satisfied \cr 0 & otherwise} \ .$

See logdelta for definition of the triangle function $$\Delta(a,b,c)$$. Note that the arguments $$j_1, j_2, j_3, m_1, m_2, m_3$$ need to be supplied doubled, as $$2j_1, 2j_2, 2j_3, 2m_1, 2m_2, 2m_3$$ so that the parameters can be considered integral even though the half-integral angular momentum is allowed.

 double Wigner6j_2 ( int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6 )

Compute the Wigner 6j coefficient. Even though there is a routine "gsl_sf_coupling_6j" in GSL library, it has to be implemented anew, because of factorial overflows. This routine computes all factorials only in logarithms using the standard function "lgamma". The formula is taken from Edmonds, A. R.: Angular momentum in quantum mechanics, Princeton 1968.

$\left\{ \matrix{j_1 & j_2 & j_3 \cr j_4 & j_5 & j_6} \right\} = \Delta(j_1,j_2,j_3) \Delta(j_3,j_4,j_5) \Delta(j_1,j_5,j_6) \Delta(j_2,j_4,j_6) \sum_k \frac{(-1)^k (k+1)!}{(k-j_1-j_2-j_3)! (k-j_1-j_5-j_6)! (k-j_2-j_4-j_6)! (k-j_3-j_4-j_5)! (j_1+j_2+j_4+j_5-k)! (j_1+j_3+j_4+j_6-k)! (j_2+j_3+j_5+j_6-k)!} \ .$

See logdelta for definition of the triangle function $$\Delta(a,b,c)$$. Note that the arguments $$j_1, j_2, j_3, j_4, j_5, j_6$$ need to be supplied doubled, as $$2j_1, 2j_2, 2j_3, 2j_4, 2j_5, 2j_6$$ so that the parameters can be considered integral even though the half-integral angular momentum is allowed.