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.

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.

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.

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.

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.