Hex
1.0
Hydrogen-electron collision solver
|
#include <cmath>
#include <iostream>
#include <gsl/gsl_sf.h>
#include "arrays.h"
#include "special.h"
#include "complex.h"
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... | |
#define ASYEPS 1e-5 |
double ClebschGordan | ( | int | l1, |
int | m1, | ||
int | l2, | ||
int | m2, | ||
int | L, | ||
int | M | ||
) |
void clipang | ( | double & | theta, |
double & | phi | ||
) |
double computef | ( | int | lambda, |
int | l1, | ||
int | l2, | ||
int | l1p, | ||
int | l2p, | ||
int | L | ||
) |
int coul_F | ( | int | l, |
double | k, | ||
double | r, | ||
double & | F, | ||
double & | Fp | ||
) |
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 |
||
) |
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.
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 | ||
) |
l | Angular momentum. |
k | Wavenumber. |
|
inline |
Derivative of Riccati-Bessel function
n | Degree of the function. |
z | Complex argument. |
Vectorized interface for dric_j.
Returns values of all derivatives of Riccati-Bessel functions of order less than or equal to lmax.
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.
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<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.
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 \ . \]
Hydrogen radial function (radius-multiplied) Evaluate hydrogen radial function (radius-multiplied) for complex argument
n | Principal quantum number. |
l | Orbital quantum number. |
z | Radius: real or complex argument. |
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.
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.
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) . \]
n | Degree of the Riccati-Bessel function. |
z | Complex argument. |
Vectorized interface for ric_j.
Returns values of all Riccati-Bessel functions of order less than or equal to lmax.
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.