Hex  1.0 Hydrogen-electron collision solver
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
 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 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
 n Degree of the function. z Complex 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
 lmax Angular momentum limit. z Complex 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. n Order of the matrix. d On entry, the N diagonal elements of the tridiagonal matrix A. On exit, if INFO = 0, the eigenvalues in ascending order. e On 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. z If 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. ldz The leading dimension of the array Z. LDZ >= 1, and if JOBZ = 'V', LDZ >= max(1,N). work If 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 > 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 )

Parameters
 n Principal quantum number. l Orbital quantum number. z Radius: 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
 n Degree of the Riccati-Bessel function. z Complex 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
 lmax Angular momentum limit. z Complex 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.