API Reference

Spheroidal angular wave functions of the first or second kind.

Computes prolate or oblate angular functions and their coordinate derivatives. kind=1 (default) selects Ps; kind=2 selects the Ferrers-branch Qs of DLMF 30.5 and 30.8(ii), on the real interval -1 < η < 1.

The first-kind functions are eigenfunction solutions to the angular wave equation, orthogonal on the interval η ∈ [-1, 1]. They closely resemble associated Legendre polynomials when c is small.

Mathematical Background:

• The spheroidal wave equation separates into radial and angular parts
• Prolate angular equation: d/dη[(1-η²)dS/dη] + [λ - c²η² - m²/(1-η²)]S = 0
• For oblate angular functions, replace -c²η² by +c²η²
• Here λ(c,m,n) is the eigenvalue (separation constant)
• Angular values and derivatives include the Condon–Shortley phase (-1)^m
• Order m: 0 ≤ m ≤ n
• Degree n: defines which eigenvalue/eigenfunction to compute

Parameters:

m::Integer
    Order parameter (must satisfy 0 ≤ m ≤ n)
n::Integer  
    Degree parameter (defines eigenfunction; n ≥ 0)
c::Real or Complex
    Size parameter (prolate: c = kd/2 with k=wavenumber, d=interfocal distance)
    Real c and real coordinates give real angular values.
η::Union{Number,AbstractVector{<:Number}}
    Real points in [-1, 1], or complex points off the real rays outside [-1, 1].
    Complex values continue the real-coordinate normalization.
    Scalar inputs are forwarded as a one-point batch

Keyword Arguments:
    spheroid::Symbol = :prolate
        :prolate or :oblate geometry
        - :prolate: elongated spheroid (ζ > 0 in ξ,η,ϕ coordinates)
        - :oblate: flattened spheroid (ζ < 0)
    
    precision::Symbol = :double
        :double (default) or :quad precision
        - :double: Float64 outputs (about 16 significant decimal digits)
        - :quad: BigFloat outputs, targeting quad precision (slower)
    
    normalize::Bool = false
        Whether to scale by sqrt((2n+1)(n-m)! / (2(n+m)!))
        - false: Meixner-Schafke normalization (norm → ∞ as m → ∞)
        - true: Unity normalization (constant norm, convenient for applications)
        - Only available for kind=1; kind=2 rejects normalize=true

    kind::Integer = 1
        1: regular first-kind Ps, with the existing normalization
        2: second-kind Qs, normalized by DLMF 30.5.4 and 30.11.4
           Qs has parity (-1)^(n-m+1), and equals Ferrers Q_n^m at c=0.
           Endpoint values and derivatives are signed, one-sided infinities.
           Complex c uses the same eigenmode branch as kind=1.

    scaled::Bool = false
        Return value and derivative as (mantissa, exponent) arrays, with
        value = mantissa * 10^exponent, preserving native scaling.
    logderivative::Bool = false
        Append S′/S with one-sided endpoint limits. Other computed zeros give NaN.
    second_derivative::Bool = false
        Append S″ from the ODE, using one-sided limits at singular endpoints.
    derivatives::Integer = 1
        Return coordinate derivatives through order 1–4, with one-sided endpoint limits.

Returns: NamedTuple with fields: - .value::Vector{Real or Complex} Function values Smₙ(η, c) at each η point - .derivative::Vector{Real or Complex} First derivatives dSmₙ/dη at each η point

Accuracy: Accuracy depends on parameter, order and coordinate. The precision setting selects the working/output precision, not a guaranteed number of accurate digits.

Special Cases: - c = 0: Returns Ferrers Pn^m (kind=1) or Qn^m (kind=2) - kind=1, c = 0, η = ±1: Uses one-sided endpoint derivatives; m = 1 gives signed infinities - c → ∞: Functions oscillate rapidly; use appropriate η resolution - Parity in η is (-1)^(n-m) for kind=1 and (-1)^(n-m+1) for kind=2 - Qs is logarithmically singular at ±1 for m=0 and algebraically singular for m>0. Its logarithmic derivative has limits -Inf at -1 and Inf at +1. Singular or unresolved second-kind normalization raises an error.

Performance: - Real-coordinate first-kind batches use the native Fortran backend. - Complex coordinates use differential-equation continuation. - Second-kind batches reuse coefficient data and Taylor propagation across sorted coordinates. They are slower, particularly near singular endpoints.

Examples:

# Basic usage: prolate Smn with c=1.5
η = [-0.9, -0.5, 0, 0.5, 0.9]
result = smn(1, 2, 1.5, η)
println(result.value)         # Function values
println(result.derivative)    # Derivatives

# Normalized functions (unity norm)
result_norm = smn(0, 3, 10.0, η; normalize=true)

# Oblate spheroid with higher precision
result_quad = smn(2, 4, 5.0, η; spheroid=:oblate, precision=:quad)

# Complex parameter (advanced)
c_complex = 1.5 + 0.1im
result_c = smn(1, 2, c_complex, η)

# Second kind, with quad precision
qs = smn(1, 2, 1.25, 0.3; kind=2, precision=:quad)

Notes: - Order and degree must satisfy 0 ≤ m ≤ n - All η values must be in [-1, 1] - With kind=1 and normalize=false, uses Meixner–Schäfke normalization - Angular values and derivatives include the Condon–Shortley phase (-1)^m - Complex c continues the eigenmode and normalization sign vertically from real(c), matching neighboring native labels of the same parity when needed. Unresolved continuation paths raise an error. See eigenvalue_sweep for eigenvalue continuation along an explicit complex path.

References: - DLMF §30.2: https://dlmf.nist.gov/30.2 - DLMF §30.4(i), phase after Eq. 30.4.1 and exact-zero identity Eq. 30.4.2: https://dlmf.nist.gov/30.4#i - Second-kind definition and normalization: https://dlmf.nist.gov/30.5 and https://dlmf.nist.gov/30.8#ii - Van Buren & Boisvert (2004): Accurate calculation of prolate spheroidal wave functions

smn(m, n::AbstractUnitRange{<:Integer}, c, eta; kwargs...)

Evaluate angular wave functions with n given as a unit range of degrees. The value and derivative matrices have rows corresponding to eta and columns corresponding to n. Normalization, precision, and spheroid keywords match the single-degree method; a scalar eta gives one row.

Use precision=:quad to request quad precision with the same numeric arguments.

Spheroidal radial wave functions and coordinate derivatives.

Computes prolate or oblate spheroidal radial wave functions and their derivatives. Ordinary calls reconstruct the values. scaled=true instead preserves native mantissas and decimal exponents to avoid reconstruction overflow and underflow.

These are eigenfunction solutions to the radial part of the spheroidal wave equation, defined for x >= 1 (prolate) or x >= 0 (oblate), for real or complex parameters. At the prolate endpoint, only the first kind is defined. Four kinds of radial functions are available.

Mathematical Background:

• Prolate radial functions satisfy: d/dξ[(ξ²-1)dRmₙ/dξ] + [c²ξ² - λ - m²/(ξ²-1)]Rmₙ = 0
• λ(c,m,n) is the eigenvalue from the angular function problem
• Functions are related to spherical Bessel and Hankel functions when c→0
• Characteristic-exponent representation: R = R_characteristic × 10^exponent

Parameters:

m, n, c: See smn() documentation

x::Union{Number,AbstractVector{<:Number}}
    Radial evaluation points
    - Real prolate: x >= 1. Real oblate: x >= 0.
    - Complex prolate: excludes the real segment [-1, 1).
    - Complex oblate: excludes imaginary rays z=iy with abs(y) >= 1.
      Continuation through z=0 fixes the oblate branch.

kind::Integer = 1
    Which radial function kind:
    
    1: First kind R₁ (propagating, finite at x→∞)
       - Regular solution; used for most applications
       - Analogous to spherical Bessel j_ℓ(kr)
    
    2: Second kind R₂
       - Complementary solution
       - Analogous to spherical Bessel y_ℓ(kr)
    
    3: Third kind (linear combination) = R₁ + i·R₂
       - Hankel function analog (outgoing wave)
    
    4: Fourth kind (linear combination) = R₁ - i·R₂
       - Alternate Hankel function analog

Output options: - scaled=true: return (mantissa, exponent) arrays for values and derivatives, representing mantissa*10^exponent. - logderivative=true: append R′/R, calculated before output conversion; use one-sided endpoint limits and return NaN at other computed zeros. - second_derivative=true: append R″ from the ODE, with one-sided limits for the regular prolate first kind at x=1; undefined kinds retain NaN. - derivatives=1: return coordinate derivatives through order 1–4. Endpoint derivatives use finite or infinite one-sided limits. - normalization=:static: evaluate c^(-n)R₁ or c^(n+1)R₂, including their limits at c=0. Supports kinds 1 and 2. The default is :standard.

Returns (default options): NamedTuple with fields: - .value::Vector{ComplexF64 or Complex{BigFloat}} Function values Rmₙ(x, c) (always complex, even for real input) Real parameters give real kinds 1/2 and complex combinations for kinds 3/4. Complex parameters can give complex values for every kind.

- .derivative::Vector{ComplexF64 or Complex{BigFloat}}
    First derivatives dRmₙ/dx

Accuracy & Overflow Protection: - Complex kinds 3/4 use direct spherical Hankel expansions, preserving the outgoing/incoming exponential and its coordinate derivatives - Uses characteristic-exponent representation internally - Ordinary output can overflow or underflow during reconstruction - scaled=true transfers mantissas and exponents before reconstruction - Scaling does not improve the native solver's numerical accuracy

Special Cases: - c = 0: Requires normalization=:static - x = 1 (prolate boundary): kinds 2–4 return NaN - Large c or x: Rapid oscillation; may need fine resolution

Examples:

# Basic radial function (kind 1)
x = [1.5, 2.0, 3.0, 5.0]
r1 = rmn(0, 1, 200, x; kind=1)
println(r1.value)

# Second kind for same parameter set
r2 = rmn(0, 1, 200, x; kind=2)

# Hankel-like combinations
r3_hankel = rmn(0, 1, 200, x; kind=3)  # Outgoing wave
r4_hankel = rmn(0, 1, 200, x; kind=4)  # Incoming wave

# Verification: scaled prolate Wronskian should be one (c=200 here)
W = r1.value .* r2.derivative - r1.derivative .* r2.value
println(200 .* (x.^2 .- 1) .* W)  # Approximately one at each point

# Oblate spheroid
r_oblate = rmn(1, 2, 500, x; spheroid=:oblate, kind=1)

Notes: - Return values are complex even when c and x are real - For real c and x, kinds 1 and 2 are real and kinds 3 and 4 are complex - Scalar coordinates return length-one arrays - For real c > 0: W = 1/(c(x²-1)) (prolate), 1/(c(x²+1)) (oblate)

Performance: - Vectorized batch computation in Fortran - Four kinds computed efficiently in single Fortran call - Typical: 1000 points in ~0.001s

References: - DLMF §30.3: https://dlmf.nist.gov/30.3 - Van Buren & Boisvert (2004): Accurate calculation of prolate spheroidal wave functions

rmn(m, n::AbstractUnitRange{<:Integer}, c, x; kwargs...)

Evaluate radial wave functions with n given as a unit range of degrees. The value and derivative matrices have rows corresponding to x and columns corresponding to n. All four radial kinds and the spheroid and precision keywords are supported. A scalar x gives one row. Use precision=:quad to request quad precision with the same numeric arguments.

eigenvalue(m, n, c; spheroid=:prolate, precision=:double,
           operator=:separation, form=:value)

Compute the spheroidal separation constant λₘₙ(c), or an integral-operator eigenvalue.

Returns the eigenvalue associated with order m, degree n, and size parameter c for either prolate or oblate spheroidal wave functions.

operator=:separation is the default. The integral operators require m=0, n>=0, spheroid=:prolate, and finite real c>=0:

• operator=:concentration: eigenvalue Λ of the kernel sin(c*(x-t))/(pi*(x-t)) integrated over [-1,1] (diagonal value c/pi). form=:value returns Λ, :log returns log(Λ), and :complement returns 1-Λ computed before output rounding. For positive c, 0<Λ<1.
• operator=:fourier: eigenvalue μ of the kernel exp(-im*c*x*t) integrated over [-1,1]. Its phase is (-im)^n, and Λ=c*abs2(μ)/(2pi) for c>0. Only form=:value is supported. Returns ComplexF64 or Complex{BigFloat}.

At c=0, concentration eigenvalues are zero; Fourier eigenvalues are 2 for n=0 and zero otherwise. Concentration results are Float64 or BigFloat. Ordinary output may round to zero or one; use the logarithm or complement forms to retain small tails. Working precision and expansion size are refined internally; failed convergence raises an error. Definitions: DLMF 30.15.3, 30.15.5.

For the separation constant with complex c, follow degree n from real(c) vertically to c, matching the same mode used by smn and rmn. Use eigenvalue_sweep with a complex vector for continuation along another path; closed paths can exchange modes.

Args: m, n, c: spheroidal parameters spheroid: :prolate or :oblate precision: :double or :quad

Returns: For the separation constant: real when c is real, complex when c is complex. Integral-operator return types are specified above.

eigenvalue_sweep(m, n, c_grid; spheroid=:prolate, precision=:double,
                branch_lock=true, branch_window=1, use_jacobian_predictor=true)

Track the separation constant over a nonempty, finite real grid c_grid. The grid must be strictly increasing or decreasing; a single point is allowed. m and the starting degree n are integers with 0 <= m <= n.

Keywords

• spheroid=:prolate: geometry, either :prolate or :oblate.
• precision=:double: :double or :quad; selects Float64 or BigFloat.
• branch_lock=true: choose the candidate eigenvalue closest to the predicted value. With false, evaluate the original degree n at every point.
• branch_window=1: nonnegative number of degrees to search on each side of the previously selected degree k: max(m,k-branch_window):k+branch_window. Zero keeps that degree fixed; larger windows evaluate more candidates. Has no effect when branch_lock=false.
• use_jacobian_predictor=true: predict the second grid value using lambda[1] + jacobian_eigen(m,n,c_grid[1]) * (c_grid[2]-c_grid[1]). If disabled or the derivative is unavailable/unreliable, use lambda[1]. Later steps use the slope between the two preceding results. This option affects branch selection only, not the eigenvalue evaluation itself.
• evaluator: optional function (m,n,c) -> real_eigenvalue, defaulting to eigenvalue with the requested geometry and precision. With a custom evaluator, set use_jacobian_predictor=false unless its derivative agrees with jacobian_eigen.

Returns

A named tuple with one entry per grid point:

• c: grid coordinates at the requested precision.
• lambda: selected eigenvalues at the requested precision.
• selected_n: integer degree selected at each point.
• switched_branch: true where the selected degree differs from the preceding point; the first entry is always false.

eigenvalue_sweep(m, n, c_grid::AbstractVector{<:Complex};
                spheroid=:prolate, precision=:double,
                branch_lock=true, branch_window=1)

Continue an eigenvalue along the ordered, piecewise-linear complex path c_grid. Initialize degree n by continuation from real(first(c_grid)) to the first point, as in scalar eigenvalue. Subsequent points continue from the preceding state, with adaptive subdivision and angular-profile matching. Repeated points and closed paths are allowed.

Only native degrees of the same parity can be selected. branch_window counts neighbors in that parity class (one means the current degree and its valid neighbors n-2 and n+2; zero keeps the label fixed). The returned selected_n contains native labels; a label change need not be a discontinuity of the continued eigenvalue. switched_branch marks those label changes. c and lambda are complex vectors at the requested precision.

With branch_lock=false, evaluate the original native label independently at each point. Complex continuation uses profiles, not the real sweep's Jacobian predictor or custom eigenvalue-only evaluator. Unresolved paths raise an error. use_jacobian_predictor is accepted but has no effect for complex grids; evaluator must be nothing. Other keywords and return fields follow the real method, with ComplexF64 or Complex{BigFloat} for c and lambda. Paths around branch points can return a different eigenvalue at their starting coordinate; this function does not impose a globally single-valued branch.

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

Return a vector of estimated decimal digits for the requested function values. Real coordinates use backend estimates. Complex coordinates compare higher-precision continuation with halved steps and include conditioning and output rounding. These are diagnostics, not rigorous error bounds or statistical confidence intervals, and do not certify the returned coordinate derivatives.

Arguments and keywords

• m, n: integer order and degree with 0 <= m <= n.
• c: finite real or complex spheroidal parameter, following rmn. Radial normalization=:static includes zero. Its real-coordinate calls have no digit estimate.
• arg: nonempty vector of finite coordinates, with one estimate returned per entry in the same order. For angular functions use -1 <= eta <= 1; for real radial coordinates use prolate x>=1 or oblate x>=0. Complex coordinates follow the cuts of smn and rmn. Singular coordinates follow the behavior of the corresponding wave function. A single coordinate must be wrapped in a vector, e.g. [0.3].
• target=:radial: estimate radial values from rmn; :angular selects smn.
• spheroid=:prolate: geometry, :prolate or :oblate.
• precision=:double: :double or :quad; evaluate at that precision before reporting integer digit estimates. This does not request a particular digit count.
• kind=1: angular first kind (1) or second kind (2); radial kinds 1:4. For real coordinates, only radial kind=2 has a radial digit estimator. The default radial kind=1 therefore returns -1 at every coordinate.
• normalize=false: use the default angular normalization. For angular kind=1, true selects unit integral normalization for real parameters and its analytic continuation for complex parameters. Angular kind=2 rejects true. This keyword has no effect on radial estimates.
• diagnostics=false: true returns (; digits, diagnostics). For complex coordinates, diagnostics contains per-point refinement changes, coordinate, parameter and eigenvalue condition numbers, and convergence flags. For real coordinates, diagnostics is nothing.

Interpreting the result

An entry of -1 means no estimate is available. Zero means the diagnostic reports no reliable decimal digits. Positive entries are estimates, not guarantees. Invalid inputs raise the same domain/argument errors as smn or rmn.

For complex coordinates, zeros, endpoint limits, nonfinite outputs and unresolved refinement return -1. Condition numbers above 100 receive a :warning flag. The following backend limitations apply to real coordinates:

• Angular abs(c)<=1e-4: return -1; the spherical-limit and small-parameter expansions have no calibrated digit estimator.
• Angular purely imaginary complex parameters return -1.
• Angular kind=2: evaluate Qs, then return -1; its propagation and coefficient convergence checks are not calibrated decimal-digit estimates.
• Refined angular evaluations for real parameters with abs(c) > n + 1 return -1; native estimates do not describe the coefficient-based result.
• Angular endpoint limits and quad endpoint reconstruction: return -1; the native estimate describes the normalization point rather than the requested coordinate.
• Radial c=0: requires normalization=:static, which returns -1.
• Radial kind=2: report the native second-kind estimate, which may use Wronskian consistency and cancellation diagnostics.
• Radial kind=1, 3, or 4: evaluate the requested function, then return -1. A second-kind/pair estimate does not establish the accuracy of these functions or their potentially cancelling Hankel combinations.
• Nonfinite values and unavailable or out-of-range native estimates return -1.

Quad estimates use the same precision-preserving inputs as the evaluation path. For sensitive calculations, use independent references, identities, or precision comparisons in addition to the backend diagnostic.

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

radial_wronskian(m, n, c, x; spheroid=:prolate, precision=:double, form=:raw, normalization=:standard)

Return W = R1 .* R2′ - R1′ .* R2 at the radial coordinates x. Here R1 and R2 are rmn kinds 1 and 2, and primes denote coordinate derivatives. This function always uses that pair; it has no kind or normalize keyword.

Arguments and keywords

• m, n: integer order and degree, 0 <= m <= n.
• c: finite real or complex parameter, with the same domain as rmn.
• x: scalar or nonempty vector of radial coordinates, with the same complex cuts as rmn. For real coordinates use x>1 for prolate and x>=0 for oblate geometry. The prolate endpoint x=1 is singular for the second-kind solution, so it cannot give a finite Wronskian check.
• spheroid=:prolate: geometry, :prolate or :oblate.
• precision=:double: :double or :quad, selecting Float64 or BigFloat components in the result.
• form=:raw: return W. :normalized returns c*(x^2-σ)*W and :error returns abs(c*(x^2-σ)*W - 1), with σ=1 for prolate and σ=-1 for oblate.
• normalization=:static: use the static-normalized pair, with W=1/(x²-σ) and normalized output (x²-σ)*W. Includes c=0.

Returns

A vector with one entry per coordinate (length one for scalar x). :raw and :normalized return ComplexF64 or Complex{BigFloat} entries; :error returns real Float64 or BigFloat entries.

For real c > 0, the standard normalization gives W(x) = 1/(c*(x^2-1)) for prolate functions and W(x) = 1/(c*(x^2+1)) for oblate functions. The raw Wronskian therefore varies with position. At c=0, standard normalization throws DomainError, as for rmn.

Use the normalized Wronskian error to check consistency:

errors = radial_wronskian(0, 1, 2.0, [1.5, 2.0, 3.0];
    precision=:quad, form=:error)
maximum(errors)

Normalized forms preserve scaling during intermediate products.

jacobian_eigen(m, n, c; spheroid=:prolate, precision=:double,
               operator=:separation, form=:value, order=1, h=nothing,
               diagnostics=false, adaptive=true, rtol=1e-6, atol=1e-10)

Differentiate eigenvalue with respect to c. Order and degree satisfy 0 <= m <= n. The default uses analytic differentiation; an explicit h selects finite differences.

Keywords

• spheroid=:prolate: :prolate or :oblate.
• precision=:double: :double or :quad, with Float64 or BigFloat components in the result.
• operator=:separation: differentiate the separation constant. :concentration and :fourier require real c >= 0, m=0, and :prolate.
• order=1: use order=2 for the second parameter derivative. This differences analytic first derivatives with fourth-order stencils and step refinement. Complex results have d2_dcreal2, d2_dcreal_dcimag, and d2_dcimag2 fields. diagnostics=true appends a shared metadata field for order 2.
• form=:value: differentiate the selected eigenvalue. Concentration also accepts :log and :complement, giving derivatives of log(Λ) and 1-Λ.
• h=nothing: analytic first derivatives or an automatic step for order 2. A positive finite step selects centered differences. Integral operators require h<c for order 1 or 2h<c for order 2 when c>0. Integral stencils use forward differences at zero.
• diagnostics=false: include diagnostic fields when true; see below.
• adaptive=true: allow smaller finite-difference steps when consistency is poor. Has no effect on the default analytic calculation.
• rtol=1e-6, atol=1e-10: positive finite relative and absolute tolerances for finite-difference step-halving checks; not requested accuracy bounds.

Returns

For real c, return a scalar derivative, or (derivative, metadata) when diagnostics=true. Concentration results are real; Fourier results are complex.

For order 1 and complex c=a+ib, return (d_dcreal, d_dcimag), containing ∂λ/∂a and ∂λ/∂b. With metadata, append metadata_dcreal and metadata_dcimag. On a local analytic branch, d_dcimag = im*d_dcreal.

Diagnostics

Metadata includes method, finite_flag, conditioning_flag, and suggested_action. Coefficient differentiation reports expansion refinement, tail size, eigenvalue conditioning, and eigenproblem/sensitivity residuals; step_used and relative_change_when_halving_step are nothing. Finite differences populate those two step fields, with conditioning :good, :warning, or :poor and action :accept, :retry_smaller_h, or :use_quad. Unresolved curvature returns action :unresolved in quad precision. Integral limits use :right_limit.

jacobian_eigen(1, 2, 1.25; precision=:quad)
jacobian_eigen(0, 2, 1.0; operator=:concentration, diagnostics=true)

Numerical Jacobian of smn outputs with respect to c.

order=2 returns d2value_dc2 and d2derivative_dc2 for real c. For complex c, use suffixes dcreal2, dcreal_dcimag, and dcimag2. It differences analytic first derivatives using fourth-order stencils. h selects the step and diagnostics=true appends shared step diagnostics.

By default, differentiate the normalized Legendre coefficient eigenvector and evaluate its expansion at the requested coordinates. Eigenvalue and angular sensitivities reuse the same internal coefficient calculation. For complex parameters this uses the analytic bilinear normalization and the local branch selected by smn. Imaginary-direction partials follow analyticity on that branch. For order 1, supplying h selects centered finite differences with local phase transport. At first-kind angular endpoints, the mixed coordinate/parameter derivative uses its one-sided limit, which can be infinite.

kind=2 differentiates Qs, including its DLMF normalization and differential equation, on the same local branch as smn(...; kind=2). It rejects normalize=true. At the singular coordinates eta=±1, both parameter derivatives return NaN; nearby interior coordinates remain supported. Real first-kind derivatives with abs(c) > n + 1 use extra working precision to protect exponentially small interior values from cancellation.

For order 1 and real c, returns:

• dvalue_dc
• dderivative_dc

For order 1 and complex c = a + ib, returns:

• dvalue_dcreal, dvalue_dcimag
• dderivative_dcreal, dderivative_dcimag

Numerical Jacobian of rmn outputs with respect to c.

order=2 returns d2value_dc2 and d2derivative_dc2 for real c. For complex c, use suffixes dcreal2, dcreal_dcimag, and dcimag2. It differences analytic first derivatives using fourth-order stencils. h selects the step and diagnostics=true appends shared step diagnostics.

By default, differentiate the normalized spherical-Bessel expansions and their coefficient eigenproblem. Near boundaries where the second-kind expansion converges slowly, propagate the differentiated radial equation from x=2. Both geometries, all four kinds, and real/complex parameters are supported. For order 1, supplying h selects centered finite differences with local mode transport. The analytic path reports no finite-difference step. Its coefficient residuals and series-tail checks are diagnostics, not rigorous error bounds.

At prolate x=1, first-kind derivatives use one-sided limits for real or complex parameters. Other kinds return NaN. normalization=:static includes c=0 and differentiates the static-normalized functions, including their powers of c.

Order 1 with explicit h reuses function values and coordinate derivatives. Order 2 with standard normalization requires 2h<abs(c). Step diagnostics are not accuracy bounds.

For order 1 and real c, returns:

• dvalue_dc
• dderivative_dc

For order 1 and complex c = a + ib, returns:

• dvalue_dcreal, dvalue_dcimag
• dderivative_dcreal, dderivative_dcimag

dmn(m, n, c; spheroid=:prolate, precision=:double,
                       normalize=false, rtol=nothing, max_terms=512)

Return the Ferrers expansion of smn(m,n,c,x; kind=1, normalize, ...): S(x) = sum(coefficients[i] * P_degrees[i]^m(x)), with Condon–Shortley phase included in the Ferrers basis. Only degrees of the parity of n-m occur. In DLMF 30.8.1 notation, the returned coefficient at degree l=n+2k is (-1)^k * a[n,k]^m, not a[n,k]^m itself.

The named tuple contains degrees, coefficients, dcoefficients_dc, eigenvalue, deigenvalue_dc, terms, tail_norm, relative_change, and converged. Derivatives hold the chosen normalization fixed. Returned arrays are independent copies. normalize=true uses the same unit bilinear norm as smn. Complex coefficients are analytic continuations, not Hermitian-normalized.

Supports integers 0 ≤ m ≤ n, real or complex c, both geometries, and c=0. precision=:double returns Float64 values, :quad returns BigFloat values, with complex counterparts for complex c. rtol controls coefficient-vector convergence, not relative accuracy of every tiny coefficient. An unresolved mode or exhausted max_terms raises an error.

kmn(m, n, c; spheroid=:prolate, precision=:double,
               normalize=false, rtol=nothing, max_terms=512)

Return the dimensionless connection coefficient K_n^m(γ) in DLMF 30.11.8: S_n^{m(1)}(z,γ) = K_n^m(γ) * Ps_n^m(z,γ^2). Here Ps is the exterior Legendre-branch continuation, not the Ferrers value at a real radial coordinate. For prolate γ=c. For oblate γ=im*c, and the corresponding radial coordinate is z=-im*x. The oblate result is complex even for real c.

normalize=true gives the multiplier for the unit-normalized angular function. Supports 0 ≤ m ≤ n, real or complex c, and both precisions. At c=0, the factor is 1 for (m,n)=(0,0) and 0 for other modes before optional normalization. These limits do not define rmn at zero parameter. Formulas 30.11.10–11 use origin values or derivatives according to parity and avoid division by angular zeros at an arbitrary matching coordinate. Unresolved sums raise an error.

amn(m, n, c; spheroid=:prolate, precision=:double,
              rtol=nothing, max_terms=512)

Return the dimensionless DLMF 30.11.4 radial factor A_n^m(σ*c^2) in Meixner–Schäfke normalization, where σ=1 for prolate and σ=-1 for oblate. This function accepts signed integer order -n ≤ m ≤ n. For negative order -p, it returns sum(d_l * (l+p)!/(l-p)!) * (n-p)!/(n+p)!, where d_l are the unnormalized angular coefficients of order p. Positive order includes the finite extension below degree m, so it is generally not sum(d_l).

A_n^m(0)=1. Both real and complex parameters are supported. The result is Float64 or BigFloat according to precision, or its complex counterpart. An unresolved coefficient expansion or singular finite extension raises an error. See dmn for truncation controls.

find_c_for_eigenvalue(m, n, lambda_target; bracket,
                      spheroid=:prolate, precision=:double,
                      operator=:separation, form=:value,
                      atol=nothing, rtol=nothing, maxiter=160, use_jacobian=true)

Solve eigenvalue(m, n, c) = lambda_target for real c. m and n are integer order and degree with 0 <= m <= n; lambda_target is the real target in the selected operator and form.

With operator=:concentration, solve for real prolate bandwidth with m=0. form=:value interprets the target as Λ in [0,1); :complement as 1-Λ in (0,1]; :log as log(Λ) in [-Inf,0). Exact unit concentration has no finite bandwidth. A zero-concentration target requires a bracket containing zero. Use complement or log targets to preserve extreme tails. The returned residual is the computed value minus the target in the selected form (zero at an exact zero-concentration match, including form=:log). Fourier inversion is unsupported.

Uses a bracketed hybrid strategy with guaranteed bisection fallback and optional Jacobian-guided acceleration through jacobian_eigen when derivative quality is acceptable.

Keywords

• bracket: required finite interval (c_lo, c_hi), with c_lo < c_hi. Its endpoint eigenvalues must enclose the target; an endpoint already within tolerance is returned immediately. Concentration requires c_lo >= 0.
• operator=:separation: invert the separation constant; :concentration inverts the concentration eigenvalue. :fourier is unsupported.
• form=:value: target convention; concentration also accepts :log and :complement.
• spheroid=:prolate: :prolate or :oblate; concentration requires :prolate.
• precision=:double: :double or :quad.
• atol, rtol: positive finite tolerances. Separation defaults are 1e-10 and 1e-8, used for residual and bracket-width checks. Concentration defaults are 1e-12/1e-10 (double), 1e-30/1e-28 (quad); the coordinate tolerance is atol+rtol*abs(c). Concentration additionally requires the residual in log(Λ/(1-Λ)) to be at most rtol, avoiding false convergence in either tail.
• maxiter=160: positive maximum number of iterations.
• use_jacobian=true: permit Newton steps using parameter derivatives. With false, separation inversion uses secant/bisection steps and concentration inversion uses bisection.

Returns

A named tuple with fields:

method=:endpoint means an input endpoint met tolerance; :newton means a derivative-based step; :secant means interpolation between bracket endpoints; :bisection means interval halving. :maxiter means the iteration limit was reached and converged=false. The method may change during the solve; this field records only the final step or termination reason.

T is Float64 for double precision and BigFloat for quad precision.