η(a::Number, k::Number)
η(ϵ::Number)

Coulomb parameter. For two arguments, returns 1/(a*k). For one argument, returns 1/sqrt(ϵ).

Both a and k must be nonzero. ϵ must be nonzero.

C(ℓ::Number, η::Number)

Coulomb normalization constant.

For complex parameters, the square root of the gamma product is defined by exp((loggamma(ℓ + 1 + im * η) + loggamma(ℓ + 1 - im * η)) / 2). This selects a local analytic branch using the principal log-gamma branches, away from their cuts and poles, and agrees with the positive normalization for real ℓ ≥ 0 and real η. The full normalization is evaluated in logarithmic form to avoid overflow or underflow of its individual factors.

See Eq. (8) of the implementation paper and DLMF 33.2.5.

θ(ℓ::Number, η::Number, ρ::Number)

Coulomb phase function.

F(ℓ::Number, η::Number, ρ::Number)

Regular Coulomb wave function.

References:

• Coulomb wave function
• Implementation paper

D⁺(ℓ::Number, η::Number)

Coulomb D⁺ normalization factor.

D⁻(ℓ::Number, η::Number)

Coulomb D⁻ normalization factor.

H⁺(ℓ::Number, η::Number, ρ::Number)

Outgoing Coulomb wave function.

References:

• Coulomb wave function
• Implementation paper

H⁻(ℓ::Number, η::Number, ρ::Number)

Incoming Coulomb wave function.

References:

• Coulomb wave function
• Implementation paper

F_imag(ℓ::Number, η::Number, ρ::Number)

Imaginary part of the regular Coulomb wave function.

G(ℓ::Number, η::Number, ρ::Number)

Irregular Coulomb wave function.

References:

• Coulomb wave function
• Implementation paper

M_regularized(α::Number, β::Number, γ::Number)

Regularized confluent hypergeometric function.

Φ(ℓ::Number, η::Number, ρ::Number)

Modified Coulomb function Φ.

w(ℓ::Integer, η::Number)
w(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

w_plus(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

w_minus(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

h_plus(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

h_minus(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

g(ℓ::Number, η::Number)

Auxiliary function for Coulomb wave functions.

Φ_dot(ℓ::Number, η::Number, ρ::Number; h=nothing)

Numerical derivative of Φ with respect to ℓ using a central finite difference. h is the step size; defaults to cbrt(eps(T)) where T is the float type of ℓ, which is near-optimal for double-sided finite differences.

F_dot(ℓ::Number, η::Number, ρ::Number; h=nothing)

Numerical derivative of F with respect to ℓ using a central finite difference. h is the step size; defaults to cbrt(eps(T)) where T is the float type of ℓ.

Ψ(ℓ::Number, η::Number, ρ::Number; h=nothing)

Auxiliary function for Coulomb wave functions. h is the finite-difference step size; defaults to cbrt(eps(T)) (see Φ_dot).

I(ℓ::Number, η::Number, ρ::Number; h=nothing)

Auxiliary function for Coulomb wave functions. h is the finite-difference step size; defaults to cbrt(eps(T)) (see Φ_dot).

debye_function(n::T, β::T, x::T; tol=1e-35, max_terms=2000) where {T <: AbstractFloat}

Compute n/x^n * ∫₀ˣ t^n/(exp(t)-1)^β dt for finite n > 0, 0 < β < n+1, and x ≥ 0. At x = 0, returns the limiting value: 0 for β < 1, 1 for β = 1, and Inf for β > 1. Returns 0 at x = Inf. Supports any AbstractFloat type (e.g., Float32, Float64, BigFloat). Array broadcasting is supported: any one argument may be an AbstractArray.

Throws ArgumentError for invalid parameters, tolerances, or work limits. Throws ErrorException if quadrature cannot meet the requested tolerance within the work limit.

• tol: positive finite relative quadrature tolerance (default 1e-35), floored at 8eps(T)
• max_terms: positive integer limit on quadrature subintervals (default 2000)

References:

• Debye function
• Paper

fresnel(z::Number)

Compute the Fresnel cosine and sine integrals and their combination:

S(z) = ∫₀ᶻ sin(π / 2 * t²) dt
C(z) = ∫₀ᶻ cos(π / 2 * t²) dt

Returns (C, S, C + im * S). For complex inputs, C and S are the analytic complex Fresnel integrals. For ComplexF16, ComplexF32, and ComplexF64, the third value is evaluated independently to avoid cancellation. Other complex float types retain the series/asymptotic implementation, whose combination can lose precision when C and S nearly cancel.

FresnelC(z::Number)

Compute the Fresnel cosine integral C(z).

FresnelS(z::Number)

Compute the Fresnel sine integral S(z).

FresnelE(z::Number)

Compute the Fresnel integral combination C(z) + im * S(z).

Clausen(n::Int, θ::Real; N::Int=10, m::Int=20)

Compute the Clausen function Cl_n(θ) of order n at angle θ. Supports any Real input type for θ.

n must be in 1..6; N must be 10 or 20 (throws ArgumentError otherwise).

Note on singularity: Clausen(1, θ) = -log|2sin(θ/2)| has a logarithmic singularity at θ = 0, 2π, 4π, …, where it returns Inf.

• N: number of Gauss-Legendre quadrature nodes used in the Euler-Maclaurin tail sum; 10 gives double precision, 20 gives extended precision.
• m: number of terms in the direct summation before switching to the tail formula.

References:

• Clausen function
• Implementation paper

Ci_complex(z::Complex{T}) where {T <: AbstractFloat}

Complex cosine integral function used in Clausen function calculations. Supports any AbstractFloat type.

f_n(n::Int, k::Int, θ::Real)

Compute the Clausen series summand fₙ(k, θ): sin(kθ)/kⁿ for even n, cos(kθ)/kⁿ for odd n. Supports any Real input type for θ.

F_clausen(n::Int, z::Complex{T}, θ::T) where {T <: AbstractFloat}

Dispatch to the correct primitive function Fₙ(z, θ) for n = 1..6. Supports any AbstractFloat type.

FermiDiracIntegral(j, x)

The Fermi-Dirac integral

Returns the value $F_j(x)$

Supports any AbstractFloat type (e.g., Float32, Float64, BigFloat).

Resources: [1] D. Bednarczyk and J. Bednarczyk, Phys. Lett. A, 64, 409 (1978) [2] J. S. Blakemore, Solid-St. Electron, 25, 1067 (1982) [3] X. Aymerich-Humet, F. Serra-Mestres, and J. Millan, Solid-St. Electron, 24, 981 (1981) [4] X. Aymerich-Humet, F. Serra-Mestres, and J. Millan, J. Appl. Phys., 54, 2850 (1983) [5] H. M. Antia, Rational Function Approximations for Fermi-Dirac Integrals (1993)

https://arxiv.org/abs/0811.0116 https://en.wikipedia.org/wiki/CompleteFermi%E2%80%93Diracintegral https://dlmf.nist.gov/25.12#iii

FermiDiracIntegralNorm(j,x)

The Fermi-Dirac integral

\[ F_j(x) = \frac{1}{\Gamma(j+1)}\int_0^\infty \frac{t^j}{\exp(t-x)+1} \, dt\]

Returns the value $F_j(x)$

Resources: [1] D. Bednarczyk and J. Bednarczyk, Phys. Lett. A, 64, 409 (1978) [2] J. S. Blakemore, Solid-St. Electron, 25, 1067 (1982) [3] X. Aymerich-Humet, F. Serra-Mestres, and J. Millan, Solid-St. Electron, 24, 981 (1981) [4] X. Aymerich-Humet, F. Serra-Mestres, and J. Millan, J. Appl. Phys., 54, 2850 (1983) [5] H. M. Antia, Rational Function Approximations for Fermi-Dirac Integrals (1993)

https://arxiv.org/abs/0811.0116 https://de.wikipedia.org/wiki/Fermi-Dirac-Integral https://dlmf.nist.gov/25.12#iii

BoseEinsteinIntegral(k, η)

Compute the unnormalized Bose–Einstein integral

\[\int_0^\infty \frac{t^k}{\exp(t-\eta)-1}\,dt = \Gamma(k+1) B_k(\eta).\]

Accepts integer and half-integer orders k > -1 and real η ≤ 0. At η = 0, returns gamma(k+1)*zeta(k+1) for k > 0, and Inf for k = -1/2 or k = 0. Returns zero at η = -Inf. Invalid inputs throw DomainError. Supports Float32, Float64, BigFloat, dot broadcasting, and ForwardDiff differentiation with respect to η.

The naming follows FermiDiracIntegral. For normalized values and continuation to negative orders, use BoseEinsteinIntegralNorm.

BoseEinsteinIntegralNorm(k, η)

Compute the normalized Bose–Einstein integral

\[B_k(\eta) = \frac{1}{\Gamma(k+1)}\int_0^\infty \frac{t^k}{\exp(t-\eta)-1}\,dt = \operatorname{Li}_{k+1}(e^\eta).\]

Accepts integer and half-integer orders k ≥ -9/2 and real η ≤ 0. For k ≤ -1, the polylogarithm defines the continuation by differentiation; the displayed integral itself requires k > -1. At η = 0, returns zeta(k+1) for k > 0 and Inf otherwise. Returns zero at η = -Inf. Invalid orders, positive η, and NaN inputs throw DomainError.

Uses Fukushima's minimax rational approximations for half-integer orders -9/2:1:39/2 and integer orders 1:19, elementary formulas for nonpositive integers, and a convergent series for higher orders. Float32 and Float64 use the double precision coefficients; BigFloat uses convergent series at the working precision. The ForwardDiff extension supports differentiation with respect to η, but not the discrete order k.

Examples

julia> round(BoseEinsteinIntegralNorm(0.5, -1.0); digits=12)
0.4284407346

Reference

T. Fukushima (2020 preprint), Analytical computation of Bose–Einstein integral of half integer orders, −9/2, −7/2, ⋯, and 39/2, and integer orders, 1, 2, ⋯, and 19, by minimax rational function approximations, Eqs. (11)–(24), Tables A.3–A.48, doi:10.13140/RG.2.2.21720.65283.

See also BoseEinsteinIntegral, FermiDiracIntegralNorm.

MarcumQ(μ::Real, a::Real, b::Real)

Compute the generalized Marcum Q-function Q_μ(a, b) of order μ with non-centrality parameter a ≥ 0 and threshold b ≥ 0. Returns a value in [0, 1].

Requires μ ≥ 0.5, a ≥ 0, and b ≥ 0; throws ArgumentError otherwise. Supports any Real input type; arguments are promoted to a common floating-point type. Array broadcasting is supported: any one argument may be an AbstractArray.

Reference: [1] https://arxiv.org/pdf/1311.0681v1

dQdb(M, a, b)

Derivative ∂Q_M(a,b)/∂b of the (standard) Marcum Q-function of order M. Requires M integer ≥1 and a>0.

voigt(x::Real, y::Real)

Evaluate the real Voigt function

\[K(x, y) = \frac{y}{\pi} \int_{-\infty}^{\infty} \frac{e^{-t^2}}{(x - t)^2 + y^2} \, \mathrm{d}t, \qquad y \ge 0.\]

For y == 0, K(x, 0) = exp(-x^2). The implementation follows the Fourier-expansion method of Abrarov, Quine, and Jagpal.

U(a::T, x::T) where {T <: AbstractFloat}

Compute the real parabolic cylinder function U(a,x). Supports Float32, Float64, and BigFloat. Uses a precision-guarded convergent series when the order and argument do not permit an accurate asymptotic expansion.

Reference: DLMF, Chapter 12.

V(a::T, x::T) where {T <: AbstractFloat}

Compute the real parabolic cylinder function V(a,x) using a convergent series with extra working precision. Supports Float32, Float64, and BigFloat.

W(a::T, x::T) where {T <: AbstractFloat}

Compute the real parabolic cylinder function W(a,x). Supports Float32, Float64, and BigFloat. Uses a precision-guarded convergent series or an asymptotic expansion whose decreasing terms reach the requested precision.

Reference: DLMF 12.14.

dU(a::T, x::T) where {T <: AbstractFloat}

Compute ∂U(a,x)/∂x from the same expansion as U. Supports Float32, Float64, and BigFloat.

dV(a::T, x::T) where {T <: AbstractFloat}

Compute ∂V(a,x)/∂x from the same convergent series as V. Supports Float32, Float64, and BigFloat.

dW(a::T, x::T) where {T <: AbstractFloat}

Compute ∂W(a,x)/∂x from the same expansion as W. Supports Float32, Float64, and BigFloat.

ParabolicCylinderD(ν::Real, x::Real)

Compute Dν(x) = U(-ν-1/2, x) for finite real order and argument. Supports Float32, Float64 and BigFloat, with mixed inputs promoted. See DLMF 12.2.5.

dParabolicCylinderD(ν::Real, x::Real)

Compute the argument derivative dDν(x)/dx = dU(-ν-1/2, x). The domain and type behavior match ParabolicCylinderD.

U_scaled(a::Real, x::Real)

Compute F(a,x) * U(a,x) for finite real a and x ≥ 0, evaluating the scaling before rounding to the output type. Here F = exp(L), with L = a*log(x/2 + sqrt(x²/4+a)) + x*sqrt(x²/4+a)/2 - a/2 in the nonoscillatory region, L = a*(log(abs(a))-1)/2 in the oscillatory region, and L = x²/4 when a = 0.

Supports Float32, Float64 and BigFloat. Invalid arguments raise DomainError. The series fallback uses extra precision and can be expensive for large orders. This scaling is from Gil, Segura & Temme (2006), equations (7), (11)–(13): paper. It removes growth/decay in both order and argument; it is not simply multiplication by exp(x²/4).

V_scaled(a::Real, x::Real)

Compute V(a,x) / F(a,x) for finite real a and x ≥ 0, with the scaling factor defined in U_scaled. Scaling is applied before rounding, so values can remain finite when V overflows. Supports Float32, Float64 and BigFloat; invalid arguments raise DomainError.

ParabolicCylinderD_scaled(ν::Real, x::Real)

Compute U_scaled(-ν-1/2, x) for finite real ν and x ≥ 0. Uses the order-dependent scaling of U_scaled, and supports Float32, Float64 and BigFloat. Invalid arguments raise DomainError.

WhittakerM(κ::Number, μ::Number, z::Number)

Compute the Whittaker function Mκμ(z) = exp(-z/2) * z^(μ+1/2) * ₁F₁(μ-κ+1/2, 1+2μ, z) on the principal branch.

Inputs must be finite and z ≠ 0; negative real arguments require an explicit complex input. Negative integer 2μ is a parameter pole. Invalid inputs raise DomainError. Real inputs with z > 0 return real values; complex inputs return complex values. Float32, Float64 and BigFloat are supported, with mixed inputs promoted. Use broadcasting for arrays.

Evaluation uses a precision-guarded Kummer series. Extra precision controls cancellation and intermediate overflow, at a cost in speed for large inputs. See DLMF 13.14 and Thompson & Barnett (1986), doi:10.1016/0021-9991(86)90046-X.

WhittakerW(κ::Number, μ::Number, z::Number)

Compute Wκμ(z) = exp(-z/2) * z^(μ+1/2) * U(μ-κ+1/2, 1+2μ, z), where U here is Tricomi's confluent hypergeometric function.

The domain, principal-branch convention and type behavior match WhittakerM, but W also accepts negative integer 2μ. Uses an asymptotic expansion when it reaches working precision, otherwise a guarded connection formula, including limits at integer 2μ. These are series methods described by Thompson & Barnett (1986); this is not a port of their full COULCC algorithm.

dWhittakerM(κ::Number, μ::Number, z::Number)

Compute the argument derivative of WhittakerM using the analytic derivative of its Kummer representation. The domain and type behavior match WhittakerM.

dWhittakerW(κ::Number, μ::Number, z::Number)

Compute the argument derivative of WhittakerW using its asymptotic expansion or the analytic derivative of its Tricomi representation. The domain and type behavior match WhittakerW.