Spherical Bessel Functions

bessels(x::AbstractVector, nℓ; kwargs...)

Convenience wrapper around bessels! that preallocates output matrices of the appropriate dimensions.

bessels!(j, j′, y, y′, x)

Compute all spherical Bessel (regular: $j_n$) and Neumann (irregular: $y_n$) functions and their derivatives at x and store the results in the vectors j, j′, y, y′. If only the Bessel functions or the Neumann functions are of interest, the other pair of arrays can be substituted by nothing. However, it is not possible to compute only the functions but not the derivatives, since they are generated using the following recurrence relations:

\[\begin{aligned} g_{n-1} &= \frac{n+1}{x} g_n + g_n'\\ g'_{n-1} &= \frac{n-1}{x} g_{n-1} - g_n. \end{aligned}\]

These recurrence relations are employed in a downward fashion for the Bessel functions and an upward fashion for the Neumann functions.

It is assumed that all passed arrays are of the same lengths (not checked).

bessels!(j, j′, y, y′, x::AbstractVector; kwargs...)

Loop through all values of x and compute all Bessel and Neumann functions, storing the results in the preallocated matrices j, j′, y, y′.

bessel_fraction(x, n)

Evaluate the continued fraction for the spherical Bessel function

\[\frac{j'_n(x)}{j_n(x)} = \frac{n}{x} - \frac{1}{T_{n+1}(x)-}\frac{1}{T_{n+2}(x)-}...\frac{1}{T_k(x)-...},\]

where

\[T_k(x) = \frac{2k+1}{x}.\]

bessel_downward_recurrence!(j, j′, x⁻¹, sinc, cosc, nmax, cf1, s; large)

Given the logarithmic derivative cf1=j′[end]/j[end] (computed using bessel_fraction), and the sign s, fill in all lower orders using the downward recurrence

\[\begin{aligned} g_{n-1} &= S_{n+1} g_n + g_n', & g_{n-1}' &= S_{n-1} g_{n-1} - g_n, & S_n &= \frac{n}{x}, \end{aligned}\]

and normalize the whole series using the analytically known expressions

\[\tag{DLMF10.49.3} \begin{aligned} j_0(z) &= \frac{\sin z}{z}, & j_0'(z) &= \frac{\cos z}{z} - \frac{\sin z}{z^2}. \end{aligned}\]

If at any point during the recurrence, j[i]>large, renormalize j′[i:end] ./= j[i] and j[i:end] ./= j[i] to avoid overflow, and the continue the recurrence. This is very useful for small $|x|$ and large $n$.

nmax allows starting the recurrence for a higher $n$ than for which there is space allocated in j and j′ (this can help convergence for the terms which are actually of interest).

neumann_upward_recurrence(y, y′, x⁻¹, sinc, cosc)

Generate the irregular Neumann functions using (stable) upward recurrence

\[\begin{aligned} g_{n+1} &= S_n g_n - g_n', & g_{n+1}' &= g_n - S_{n+2} g_{n+1}, \end{aligned}\]

starting from the analytically known expressions

\[\tag{DLMF10.49.5} \begin{aligned} y_0(z) &= -\frac{\cos z}{z}, & y_0'(z) &= \frac{\sin z}{z} + \frac{\cos z}{z^2}. \end{aligned}\]

reflect!(g, g′, offset=0)

Affect the reflection symmetries

\[\tag{DLMF10.47.14} \begin{aligned} j_n(-z) &= (-)^n j_n(z), & y_n(-z) &= (-)^{n+1}y_n(z), \\ j_n'(-z) &= (-)^{n-1} j_n'(z), & y_n'(-z) &= (-)^{n} y_n'(z), \end{aligned}\]

where the $+1$ can be added using offset.

powneg1(m)

Returns an integer power of negative unity, i.e. $(-)^n$.