Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Macros | Functions
special.cpp File Reference
#include <cmath>
#include <iostream>
#include <gsl/gsl_sf.h>
#include "arrays.h"
#include "special.h"
#include "complex.h"
Include dependency graph for special.cpp:

Macros

#define ASYEPS   1e-5
 

Functions

void dstev_ (char *jobz, int *n, double *d, double *e, double *z, int *ldz, double *work, int *info)
 Lapack routine: Eigenvalues, eigenvectors of a symmetric tridiagonal matrix. More...
 
cArray ric_jv (int lmax, Complex z)
 
Complex ric_j (int l, Complex z)
 
cArray dric_jv (int lmax, Complex z)
 
Complex dric_j (int l, Complex z)
 
Complex hydro_P_table (unsigned n, unsigned l, Complex z)
 
Complex dhydro_P_table (unsigned n, unsigned l, Complex z)
 
Complex associated_laguerre_poly (int k, int s, Complex z)
 
Complex der_associated_laguerre_poly (int k, int s, Complex z)
 
double hydrogen_wfn_normalization (int n, int l)
 
Complex hydro_P (unsigned n, unsigned l, Complex z)
 
Complex dhydro_P (unsigned n, unsigned l, Complex z)
 
Complex sphY (int l, int m, double theta, double phi)
 Spherical harmonic function. More...
 
void clipang (double &theta, double &phi)
 
Complex sphBiY (int l1, int l2, int L, int M, double theta1, double phi1, double theta2, double phi2)
 Bi-polar spherical harmonic function. 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_sigma (int l, double k)
 Coulomb phase shift. More...
 
double coul_F_asy (int l, double k, double r, double sigma)
 Asymptotic form of the regular Coulomb wave. 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 Wigner3j_2 (int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
 Wigner 3j coefficient. More...
 
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...
 
double computef (int lambda, int l1, int l2, int l1p, int l2p, int L)
 
long double dfact (long double x)
 
double ClebschGordan (int __j1, int __m1, int __j2, int __m2, int __J, 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 maxl)
 
double Hyper2F1 (double a, double b, double c, double x)
 The hypergeometric function \( {}_2F_1 \). More...
 
std::vector< std::vector< int > > FdB_partition (int n)
 Faa di Bruno partitioning. More...
 

Macro Definition Documentation

#define ASYEPS   1e-5

Function Documentation

Complex associated_laguerre_poly ( int  k,
int  s,
Complex  z 
)
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.]
void clipang ( double &  theta,
double &  phi 
)
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 der_associated_laguerre_poly ( int  k,
int  s,
Complex  z 
)
long double dfact ( long double  x)
inline
Complex dhydro_P ( unsigned  n,
unsigned  l,
Complex  z 
)

Derivative of Pnl.

Parameters
n
l
z
Complex dhydro_P_table ( unsigned  n,
unsigned  l,
Complex  z 
)
Complex dric_j ( int  n,
Complex  z 
)

Derivative of Riccati-Bessel function

Parameters
nDegree of the function.
zComplex argument.
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.
void dstev_ ( char *  jobz,
int *  n,
double *  d,
double *  e,
double *  z,
int *  ldz,
double *  work,
int *  info 
)

DSTEV computes the eigenvalues and, optionally, the left and/or right eigenvectors of other matrices.

Parameters
jobz'N' for eigenvalues only, 'V' also eigenvectors.
nOrder of the matrix.
dOn entry, the N diagonal elements of the tridiagonal matrix A. On exit, if INFO = 0, the eigenvalues in ascending order.
eOn entry, the (n-1) subdiagonal elements of the tridiagonal matrix A, stored in elements 1 to N-1 of E. On exit, the contents of E are destroyed.
zIf JOBZ = 'V', then if INFO = 0, Z contains the orthonormal eigenvectors of the matrix A, with the i-th column of Z holding the eigenvector associated with D(i). If JOBZ = 'N', then Z is not referenced.
ldzThe leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N).
workIf JOBZ = 'N', WORK is not referenced.
info= 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; > 0: if INFO = i, the algorithm failed to converge; i off-diagonal elements of E did not converge to zero.
std::vector<std::vector<int> > FdB_partition ( int  n)

The Faa di Bruno partitioning is computed by looping over possible n-tuples

\[ (m_1, m_2, \dots, m_n) \]

of integers and by picking only such that satisfy the Faa di Bruno's sum condition

\[ 1 m_1 + 2 m_2 + \dots + n m_n = n \ . \]

The initial trial partitioning is a zero tuple

\[ (0, 0, \dots, 0) \]

and the further tuples are constructed by incrementing a corresponding multidigit number, that has a number system varying with the position. The number system base for the left-most (least significant) position is \( n + 1 \), for the next position it is \( \lceil (n + 1)/2 \rceil \), for the next it is \( \lceil (n + 1)/3 \rceil \), etc. For example, if \( n = 4 \), the increments are

\[ (0, 0, 0, 0), (1, 0, 0, 0), (2, 0, 0, 0), (3, 0, 0, 0), \mathbf{(4, 0, 0, 0)}, \]

\[ (0, 1, 0, 0), (1, 1, 0, 0), (2, 1, 0, 0), \mathbf{(3, 1, 0, 0)}, (4, 1, 0, 0), \]

\[ (0, 2, 0, 0), (1, 2, 0, 0), \mathbf{(2, 2, 0, 0)}, (3, 2, 0, 0), (4, 2, 0, 0), \]

\[ (0, 0, 1, 0), \mathbf{(1, 0, 1, 0)}, (2, 0, 1, 0), (3, 0, 1, 0), (4, 0, 1, 0), \]

\[ (0, 1, 1, 0), (1, 1, 1, 0), (2, 1, 1, 0), (3, 1, 1, 0), (4, 1, 1, 0), \]

\[ (0, 2, 1, 0), (1, 2, 1, 0), (2, 2, 1, 0), (3, 2, 1, 0), (4, 2, 1, 0), \]

\[ \mathbf{(0, 0, 0, 1)}, \ \mathrm{etc.} \]

Here, the number system are 5, 3, 2 and 2. Only those tuples in bold satisfy the sum condition and will be returned.

This function is needed by the generalized chain_rule.

Todo:
Cache results.
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.
Complex hydro_P_table ( unsigned  n,
unsigned  l,
Complex  z 
)
double hydrogen_wfn_normalization ( int  n,
int  l 
)
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_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.
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.
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.

Generated by   doxygen 1.8.5