Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Data Structures | Namespaces | Macros | Functions | Variables
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

class  RadialFunction< T >
 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()
 
class RadialFunction exception
 
#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,
 
)    (sqrt(((n)-(l)-1)*((n)+(l))))
#define ALPHA_PLUS (   n,
 
)    (sqrt(((n)-(l))*((n)+(l)+1)))
#define DELTA (   a,
 
)    ((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,
 
)    Wigner3j_2(2*(a),2*(b),2*(c),2*(d),2*(e),2*(f))
#define Wigner6j (   a,
  b,
  c,
  d,
  e,
 
)    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
lAngular momentum.
kWavenumber.
rRadial coordinate.
FOutput reference for resulting value.
FpOutput reference for resulting derivative.
double coul_F_asy ( int  l,
double  k,
double  r,
double  sigma = special::constant::Nan 
)
Parameters
lAngular momentum.
kWavenumber.
rRadial coordinate.
sigmaOptionally, 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
lAngular momentum.
kWavenumber.
rRadial coordinate.
FOutput reference for resulting value.
FpOutput reference for resulting derivative.
double coul_F_sigma ( int  l,
double  k 
)
Parameters
lAngular momentum.
kWavenumber.
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
nDegree of the function.
zComplex 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
lmaxAngular momentum limit.
zComplex 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 
)

Hydrogen radial function (radius-multiplied) Evaluate hydrogen radial function (radius-multiplied) for complex argument

Parameters
nPrincipal quantum number.
lOrbital quantum number.
zRadius: 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
nDegree of the Riccati-Bessel function.
zComplex 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
lmaxAngular momentum limit.
zComplex 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.

Variable Documentation