Mathematical tools

All examples use ordinary numeric arguments. Add precision=:quad for quad precision.

Expansion coefficients and connection factors

d = dmn(1, 2, 1.25)   # angular expansion coefficients
K = kmn(1, 2, 1.25)   # joining factor
A = amn(1, 2, 1.25)   # radial factor, also accepts negative order

dmn returns degrees, coefficients, the coefficient derivatives, and convergence metadata. The coefficients expand smn in the Ferrers basis:

\[S_{mn}(c,\eta)=\sum_l d_l\,\mathsf P_l^m(\eta).\]

The basis includes the Condon-Shortley phase. normalize matches smn. In DLMF 30.8.1, $d_l=(-1)^k a_{n,k}^m$ for $l=n+2k$.

amn returns $A_n^m(\sigma c^2)$ from DLMF 30.11.4, accepting $-n\le m\le n$. Here $\sigma=1$ for prolate and $\sigma=-1$ for oblate. It uses Meixner-Schäfke normalization. Positive order includes the finite coefficient extension below degree $m$.

kmn supplies the exterior angular-radial connection DLMF 30.11.8:

\[S_n^{m(1)}(z,\gamma)=K_n^m(\gamma)\,\mathrm{Ps}_n^m(z,\gamma^2).\]

Use $\gamma=c$, $z=x$ for prolate and $\gamma=ic$, $z=-ix$ for oblate. The exterior Legendre branch uses $[z\sqrt{1-z^{-2}}]^m$ without the Ferrers phase. Oblate joining factors are complex. normalize=true gives the factor for a unit-normalized angular expansion.

All three support real or complex parameters, both geometries and precisions. By default at c=0, amn is one and kmn is zero except for $(m,n)=(0,0)$, where it is one. Radial calls at c=0 require normalization=:static. rtol and max_terms control coefficient convergence, not relative accuracy of every small coefficient.

Eigenvalue operators

eigenvalue(m, n, c; operator=:separation) selects:

Separation constant

The angular and radial equations share λₘₙ(c), with λₘₙ(0) = n(n+1). Both geometries and real or complex parameters are supported.

eigenvalue(1, 2, 1.25; precision=:quad)
eigenvalue(1, 2, 1.25; spheroid=:oblate, precision=:quad)

Concentration and finite Fourier eigenvalues

For the first-kind angular mode Sₙ = S₀ₙ on [-1,1]:

\[\int_{-1}^{1}\frac{\sin(c(x-t))}{\pi(x-t)}S_n(t)\,dt =\Lambda_n(c)S_n(x), \qquad \int_{-1}^{1}e^{-icxt}S_n(t)\,dt=\mu_n(c)S_n(x).\]

The sinc kernel is c/π at x=t. For c>0,

\[0<\Lambda_n<1,\qquad \mu_n=(-i)^n\sqrt{\frac{2\pi\Lambda_n}{c}}.\]

These conventions follow DLMF 30.15 with the interval scaled to [-1,1].

eigenvalue(0, 2, 1.0; operator=:concentration, precision=:quad)
eigenvalue(0, 2, 1.0; operator=:fourier, precision=:quad)
eigenvalue(0, 0, 50.0; operator=:concentration,
    form=:complement, precision=:quad)

Parameter derivatives

The Jacobians differentiate with respect to c; coordinate derivatives are returned by smn and rmn.

jacobian_smn(1, 2, 1.25, [0.3]; kind=2, precision=:quad)
jacobian_rmn(1, 2, 1.25, [2.0]; kind=2, precision=:quad)
jacobian_eigen(0, 2, 1.0; operator=:concentration, precision=:quad)
jacobian_eigen(0, 2, 1.0; operator=:fourier, precision=:quad)
jacobian_smn(1, 2, 1.25, [0.3]; order=2, diagnostics=true)

order=2 selects second parameter derivatives. Wave results use d2value_dc2 and d2derivative_dc2. Complex results use suffixes dcreal2, dcreal_dcimag, and dcimag2. These calculations difference analytic first derivatives with fourth-order stencils. h controls the step and diagnostics=true adds shared step-refinement diagnostics. Standard radial and positive-bandwidth integral stencils require 2h<abs(c).

jacobian_eigen accepts the same operators and forms as eigenvalue. For a real mode ψₙ with unit norm on [-1,1] and c>0:

\[\frac{d\Lambda_n}{dc}=\frac{|\mu_n|^2}{\pi}\psi_n(1)^2, \qquad \frac{d\log\Lambda_n}{dc}=\frac{2\psi_n(1)^2}{c}, \qquad \frac{d(1-\Lambda_n)}{dc}=-\frac{d\Lambda_n}{dc}.\]

For order 1 and complex c=a+ib, results contain real- and imaginary-direction partials (d_dcreal/d_dcimag for eigenvalues; dvalue_dcreal/dvalue_dcimag and dderivative_dcreal/dderivative_dcimag for functions). On a local analytic branch, ∂f/∂b = i ∂f/∂a.

Inverse bandwidth

find_c_for_eigenvalue solves for real c within a supplied bracket:

\[\lambda_{mn}(c)=\lambda_{\mathrm{target}} \quad\text{or}\quad \Lambda_n(c)=\Lambda_{\mathrm{target}}.\]

target = eigenvalue(1, 2, 1.25; precision=:quad)
root = find_c_for_eigenvalue(1, 2, target;
    bracket=(0.5, 2.0), precision=:quad)

root = find_c_for_eigenvalue(0, 0, 0.999;
    operator=:concentration, bracket=(0.0, 10.0), precision=:quad)
root.c
root.converged
root.residual

bracket=(c_lo,c_hi) is required: finite endpoints with c_lo<c_hi whose eigenvalues enclose the target. Concentration inversion also requires c_lo>=0. An endpoint matching within tolerance is returned immediately. maxiter=160 limits iterations; use_jacobian=true permits Newton steps. With false, separation inversion uses secant/bisection and concentration uses bisection.

For concentration, form defines the target:

See find_c_for_eigenvalue for tolerance defaults and stopping criteria.

Eigenvalue sweeps

sweep = eigenvalue_sweep(1, 2, [0.5, 1.0, 1.5]; precision=:quad)
sweep.lambda

path = ComplexF64[2, 2+3im, 3im, 0, 2]
continued = eigenvalue_sweep(0, 0, path; precision=:quad)
continued.selected_n

Sweeps track separation constants only. Real grids must be monotone; complex grids specify an ordered piecewise-linear path. selected_n records degree labels and switched_branch flags label changes. branch_lock=false evaluates the native label independently at each point.

branch_window must be nonnegative; zero fixes the degree label. It has no effect with branch_lock=false. use_jacobian_predictor has no effect for complex paths, which use angular-profile matching. Both geometries and precision=:double/:quad are supported.

Multiple degrees

smn(1, 1:24, 1.25, [-0.7, 0.2, 0.8];
    spheroid=:oblate, precision=:quad, normalize=true)
rmn(1, 1:24, 1.25, [0.4, 1.0, 3.0];
    spheroid=:oblate, precision=:quad, kind=3, scaled=true)

Rows correspond to coordinates; columns correspond to degrees. A scalar coordinate gives one row. Single-degree calls return vectors, including length-one vectors for scalar coordinates.

Scaled values and logarithmic derivatives

\[f=M\,10^E,\qquad \ell=\frac{f'}{f}.\]

s = smn(200, 200, 0.0, 0.3; scaled=true, logderivative=true)
s.value.mantissa
s.value.exponent
s.logderivative

r = rmn(1, 2, 1.25, 2.0;
    kind=2, precision=:quad, scaled=true, logderivative=true)

With scaled=true, value and derivative each contain mantissa and integer exponent arrays. For finite nonzero values, 1≤abs(mantissa)<10. Complex mantissas use a common exponent for both components; zeros and nonfinite values use exponent zero.

Higher coordinate derivatives

For λ=λₘₙ(c) and σ=1 (prolate) or σ=-1 (oblate):

\[S''=\frac{2\eta S'- \left(\lambda-\sigma c^2\eta^2-\frac{m^2}{1-\eta^2}\right)S}{1-\eta^2},\]

\[R''=\frac{-2xR'- \left(c^2x^2-\lambda-\frac{\sigma m^2}{x^2-\sigma}\right)R}{x^2-\sigma}.\]

s = smn(1, 2, 1.25, 0.3; second_derivative=true, precision=:quad)
r = rmn(1, 2, 1.25, 2.0; second_derivative=true, precision=:quad)
s.second_derivative
r.second_derivative

Both options also work with scaled outputs and degree ranges.

Use derivatives=3 or 4 to append third_derivative and fourth_derivative. Endpoint derivatives use one-sided limits. Oblate x=0 is supported. derivatives=2 is equivalent to second_derivative=true.

Radial Wronskians

For real c>0 and nonsingular radial coordinates:

\[W=R_1R_2'-R_1'R_2,\qquad \widehat W=c(x^2-\sigma)W=1,\qquad \varepsilon_W=|\widehat W-1|.\]

radial_wronskian(m,n,c,x; spheroid=:prolate, precision=:double, form=:raw) uses radial kinds 1 and 2 and their coordinate derivatives. m,n are integer order and degree with 0≤m≤n. Standard normalization requires finite c≠0. x is a real scalar or nonempty vector: use prolate x>1 or oblate x≥0. spheroid=:oblate selects σ=-1; precision=:quad selects quad precision. normalization=:static includes c=0, with $W=1/(x^2-\sigma)$ and normalized output $(x^2-\sigma)W$.

Every form returns a vector, including for scalar x. Raw and normalized entries are complex; error entries are real. Normalized forms preserve scaling in intermediate products. Complex c is accepted under the same conventions as rmn; the unit-normalization check stated above is for real c>0.

radial_wronskian(1, 2, 1.25, [1.5, 2.0, 4.0];
    precision=:quad, form=:error)

The prolate identity is DLMF 30.11.7; the oblate factor follows its radial equation and normalization.

Accuracy estimates

accuracy(m,n,c,arg; target=:radial, spheroid=:prolate, precision=:double, kind=1, normalize=false, normalization=:standard, diagnostics=false) returns one integer per coordinate.

Positive entries estimate decimal digits; 0 reports no reliable digits and -1 means no estimate is available. Invalid inputs raise the corresponding wave-function errors; nonfinite values receive -1 when evaluation returns them.

accuracy(0, 1, 2.0, [1.5]; target=:radial, kind=2)
accuracy(0, 2, 1.0, [0.3]; target=:angular, normalize=true, precision=:quad)
accuracy(0, 2, 0.0, [0.3]; target=:angular) # [-1]

Complex coordinates use higher precision and halved continuation steps to estimate convergence. Add diagnostics=true for per-point refinement changes, condition numbers, and convergence flags. Zeros, endpoints and unresolved estimates return -1.

a = accuracy(1, 2, 1.25, [0.3 + 0.2im]; target=:angular, diagnostics=true)
a.digits
a.diagnostics[1].conditioning_flag