API

A spheroidal scatterer.

Attributes:

• a: the length of the semi-major axis.
• c: the length of the semi-minor axis.
• m: the relative complex refractive index.

gaussquad(s::Spheroid{T}, ngauss) where {T}

Evaluate the quadrature points, weights and the corresponding radius and radius derivative (to $\vartheta$) for a spheroid.

Returns: (x, w, r, r′)

• x: the quadrature points.
• w: the quadrature weights.
• r: the radius at each quadrature point.
• r′: the radius derivative at each quadrature point.

A cylindrical scatterer.

Attributes:

• r: the radius of the cylinder base.
• h: the height of the cylinder.
• m: the relative complex refractive index.

gaussquad(c::Cylinder{T}, ngauss) where {T}

Evaluate the quadrature points, weights and the corresponding radius and radius derivative (to $\vartheta$) for a cylinder.

Returns: (x, w, r, r′)

• x: the quadrature points.
• w: the quadrature weights.
• r: the radius at each quadrature point.
• r′: the radius derivative at each quadrature point.

A Chebyshev scatterer defined by

\[r(\theta, \phi)=r_0(1+\varepsilon T_n(\cos\theta))\]

where $T_n(\cos\theta)=\cos n\theta$.

Attributes:

• r₀: the radius of the base sphere.
• ε: the deformation parameter, which satisfies $-1\le\varepsilon<1$.
• n: the degree of the Chebyshev polynomial.
• m: the relative complex refractive index.

gaussquad(c::Chebyshev{T}, ngauss) where {T}

Evaluate the quadrature points, weights and the corresponding radius and radius derivative (to $\vartheta$) for a Chebyshev particle.

Returns: (x, w, r, r′)

• x: the quadrature points.
• w: the quadrature weights.
• r: the radius at each quadrature point.
• r′: the radius derivative at each quadrature point.

Prism{N, T, CT}(a, h, m)

A prism scatterer. The number of base edges is represented by the type parameter N. The prism is assumed to be aligned with the z-axis, and one of The lateral edges is assumed to have φ=0.

Attributes:

• a: the length of the prism base.
• h: the height of the prism.
• m: the relative complex refractive index.

Supertype of fixed-discretization solvers (one build, no iteration).

AbstractSolver

Supertype of every T-matrix solver. A solver is either a fixed-discretization solver (EBCM, IITM, ShMatrix — subtypes of AbstractFixedSolver) or an Iterative wrapper that converges one.

ConvergencePolicy(; threshold, ndgs, maxiter, nₘₐₓ_only, nₛₜₐᵣₜ, Ngₛₜₐᵣₜ, routine_generator)

Convergence settings for an Iterative solver. threshold is the relative Qsca/Qext tolerance, ndgs the quadrature points added per order, nₛₜₐᵣₜ/Ngₛₜₐᵣₜ the starting resolution (0 ⇒ auto from $k \cdot r_{\max}$), nₘₐₓ_only stops once nₘₐₓ converges, and routine_generator builds the stepping routine (default routine_mishchenko).

EBCM(nₘₐₓ, Ng; stable = false)

Fixed-discretization Extended Boundary Condition Method solver: truncation order nₘₐₓ and Ng Gauss–Legendre quadrature points. stable = true assembles the 𝐔 matrix with the cancellation-free F⁺ formulation (spheroids only). Used as transition_matrix(s, λ, EBCM(nₘₐₓ, Ng)).

IITM(nₘₐₓ, Nr, Nϑ; rₘᵢₙ = nothing)
IITM(nₘₐₓ, Nr, Nϑ, Nφ; rₘᵢₙ = nothing)

Fixed-discretization invariant-imbedding (IITM) solver: order nₘₐₓ, Nr radial and Nϑ zenithal quadrature points (plus Nφ azimuthal points for arbitrary / N-fold shapes). rₘᵢₙ defaults to rmin(s) when left nothing.

Iterative(EBCM; stable = false, threshold = 1e-4, ndgs = 4, maxiter = 20, …)

Iterative solver that sweeps the discretization of a fixed-solver method until convergence. Currently the EBCM method is supported, including stable = true:

transition_matrix(s, λ, Iterative(EBCM))                  # classic, auto-converged
transition_matrix(s, λ, Iterative(EBCM; stable = true))   # stabilized + auto-converged

The non-resolution method options (e.g. stable) and the ConvergencePolicy keywords are passed together; the method is named by its fixed-solver type.

ShMatrix(nₘₐₓ, Ng; B = nothing, momtype = nothing, store = nothing)

Fixed-discretization Sh-matrix (moment-separation) solver. B is the number of radial power-series terms; momtype/store control the moment precision (see prepare_sh; nothing selects the defaults). As a single build, transition_matrix(s, λ, ShMatrix(nₘₐₓ, Ng)); for a parameter sweep, transition_matrix(s, λs, mᵣs, ShMatrix(nₘₐₓ, Ng)) prepares once and reuses.

transition_matrix(s::AbstractAxisymmetricShape, λ)

Compute the T-matrix of an axisymmetric scatterer at wavelength λ, automatically converging the truncation order and quadrature with classic EBCM. Equivalent to transition_matrix(s, λ, Iterative(EBCM)). For a stabilized or fixed build, or another method, pass an explicit solver (see AbstractSolver).

amplitude_matrix(axi::AxisymmetricTransitionMatrix{CT, N, V, T}, ϑᵢ, φᵢ, ϑₛ, φₛ;
                          λ = 2π,
                          rot::Union{Nothing, Rotation{3}} = nothing) where {CT, N, V, T}

Calculate the amplitude matrix of an axisymmetric scatterer.

Parameters:

• axi: the T-Matrix of the scatterer.
• ϑᵢ: the zenith angle of the incident wave.
• φᵢ: the azimuth angle of the incident wave.
• ϑₛ: the zenith angle of the scattered wave.
• φₛ: the azimuth angle of the scattered wave.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.
• rot: the rotation of the scatterer.

ebcm_matrices_m(m, s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Ng) where {T, CT}

Calculate the P and U matrices for the m-th EBCM block.

ebcm_matrices_m₀(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Ng) where {T, CT}

Calculate the P and U matrices for the m=0 EBCM block.

transition_matrix(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Ng; stable = false) where {T, CT}

Calculate the T-Matrix for a given scatterer and wavelength.

Parameters:

• s: the axisymmetric scatterer.
• λ: the wavelength.
• nₘₐₓ: the maximum order of the T-Matrix.
• Ng: the number of Gauss-Legendre quadrature points to be used.
• stable: when true, assemble the 𝐔-matrix integrals with the cancellation-free F⁺ formulation of Somerville et al. (2013). The standard integrands lose all precision for spheroids of high aspect ratio (the irregular $\chi_n\psi_k$ products develop huge Laurent terms that should cancel on integration but do not numerically); stable=true removes that cancellation, recovering a relative accuracy of about 1e-9 in Float64 regardless of aspect ratio. It is only valid for Spheroid (the cancellation relies on the spheroid surface) and costs roughly 2–3× the default assembly, so it is opt-in. For it to help, nₘₐₓ must be large enough that $n_{\max}+1 \gtrsim k\cdot c + 15$, where c is the largest semi-axis — comparable to the order needed for convergence at that size anyway. The remaining Float64 round-off (in particular a ~1e-9 floor as the refractive index $s \to 1$) is orthogonal to the cancellation and is lowered by using an extended-precision element type: stable=true with Double64 reaches ~1e-25 at high aspect ratio.

Returns:

• 𝐓: an AxisymmetricTransitionMatrix struct representing the T-Matrix.

transition_matrix_m(m, s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Ng) where {T, CT}

Calculate the m-th block of the T-Matrix for a given axisymmetric scatterer.

transition_matrix_m₀(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Ng) where {T, CT}

Calculate the m=0 block of the T-Matrix for a given axisymmetric scatterer.

ShPreparation

Precomputed shape-only moment tables for the Sh-matrix moment-separation method, returned by prepare_sh. Holds the moments for every azimuthal index m = 0:nmax, which depend only on the particle geometry — not on the wavelength or refractive index — so a single preparation feeds an arbitrarily long sweep of transition_matrix(prep, λ, mᵣ) evaluations cheaply.

The stable field records whether the analytic 𝐔 stabilization is valid: it is true only for spheroids, whose negative-r-power moments vanish (computed at high precision they contribute ≈0, removing the irregular-product cancellation). For other axisymmetric shapes the moment machinery and 𝐏 are still correct, but the 𝐔 reconstruction is not stabilized.

prepare_sh(shape, nmax, Ng; B, momtype, store) -> ShPreparation

Precompute the Sh-matrix shape-only moment tables for shape up to order nmax with Ng Gauss–Legendre points. The result is reused across many transition_matrix(prep, λ, mᵣ) calls (a wavelength / refractive-index sweep) at the cost of only a cheap coefficient×moment sum each.

Keyword arguments:

• B: number of radial power-series terms per Riccati–Bessel function. It sets the largest size parameter the reconstruction resolves (the series, like the F⁺ evaluation, needs more terms as $k\cdot r_{\max}\cdot|m_r|$ grows); it is fixed at preparation time because it determines the moment band. Default max(30, nmax+15).
• momtype: precision used to accumulate the moments. The default (nothing) selects automatically: for a Spheroid the negative-power moments vanish analytically (the F⁺ result), so they are set to zero and the whole precompute runs in the shape's hardware precision — orders of magnitude faster and with no BigFloat dependency, at the same accuracy as stable=true. For other shapes (whose negative-power moments do not vanish) it falls back to BigFloat so the cancellation is resolved numerically. Pass momtype=BigFloat explicitly to force the high-precision path for a spheroid too (slower, and bit-for-bit closer to the full-precision moments in the highest-order, smallest 𝐓 entries).
• store: element type the moments are stored as (default the shape's real type); reconstruction runs in this precision.

The analytic 𝐔 stabilization is enabled only for Spheroid (see ShPreparation).

transition_matrix(prep::ShPreparation, λ, mᵣ) -> AxisymmetricTransitionMatrix
transition_matrix(prep::ShPreparation, λ)      # uses the prepared shape's mᵣ

Assemble the full T-matrix at wavelength λ and relative refractive index mᵣ from a prepare_sh preparation, reconstructing every m-block from the precomputed moments. This is the cheap inner call of a parameter sweep; the expensive geometry quadrature was done once in prepare_sh.

transition_matrix_spectrum(prep::ShPreparation, λs, mᵣs) -> Vector

Reconstruct a T-matrix at each (λ, mᵣ) pair, reusing the prepared moments. If mᵣs is a single number it is held fixed across all λs; otherwise λs and mᵣs are iterated together (and must have equal length). Pass a dispersion table as mᵣs to sweep a material with wavelength-dependent index.

F⁺(n, k, s, x; nterms, prec) -> Complex

Numerically stable $F^+_{nk}(s,x)$ (Somerville et al. (2013), Eqs. (45)–(46)), evaluated in the precision of x. $F^+/x$ is the cancellation-free replacement for $\chi_n(x)\psi_k(sx)$ in the spheroid EBCM $\mathbf{U}$-matrix integrand.

transition_matrix_iitm(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ, Nr, Nϑ; rₘᵢₙ) where {T, CT}

Use IITM to calculate the T-Matrix for a given scatterer and wavelength.

Parameters:

• s: the axisymmetric scatterer.
• λ: the wavelength.
• nₘₐₓ: the maximum order of the T-Matrix.
• Nr: the number of radial quadrature points to be used.
• Nϑ: the number of zenithal quadrature points to be used.

Keyword arguments:

• rₘᵢₙ: the starting point of the radial quadrature. Default to rmin(s), which is the radius of the maximum inscribed sphere.

Returns:

• 𝐓: an AxisymmetricTransitionMatrix struct representing the T-Matrix.

transition_matrix_iitm(s::AbstractShape{T, CT}, λ, nₘₐₓ, Nr, Nϑ, Nφ; rₘᵢₙ) where {T, CT}

Use IITM to calculate the T-Matrix for a given scatterer and wavelength.

Parameters:

• s: the scatterer.
• λ: the wavelength.
• nₘₐₓ: the maximum order of the T-Matrix.
• Nr: the number of radial quadrature points to be used.
• Nϑ: the number of zenithal quadrature points to be used.
• Nφ: the number of azimuthal quadrature points to be used.

Keyword arguments:

• rₘᵢₙ: the starting point of the radial quadrature. Default to rmin(s), which is the radius of the maximum inscribed sphere.

Returns:

• 𝐓: an AxisymmetricTransitionMatrix struct representing the T-Matrix.

transition_matrix_iitm(s::AbstractNFoldShape{N, T, CT}, λ, nₘₐₓ, Nr, Nϑ, Nφ; rₘᵢₙ) where {T, CT}

Use IITM to calculate the T-Matrix for a given scatterer and wavelength.

Parameters:

• s: the scatterer.
• λ: the wavelength.
• nₘₐₓ: the maximum order of the T-Matrix.
• Nr: the number of radial quadrature points to be used.
• Nϑ: the number of zenithal quadrature points to be used.
• Nφ: the number of azimuthal quadrature points to be used.

Keyword arguments:

• rₘᵢₙ: the starting point of the radial quadrature. Default to rmin(s), which is the radius of the maximum inscribed sphere.

Returns:

• 𝐓: an AxisymmetricTransitionMatrix struct representing the T-Matrix.

bhmie([T=Float64,], x, m; nₘₐₓ)

This is a simplified version of Bohren (1983)'s bhmie routine. Only the Mie coefficients are evaluated and returned.

Inputs:

• T: Type used for calculation. Real numbers will be stored as T, while complex numbers will be stored as C = complex(T).
• x: Size parameter of the sphere scatterer. Defined as $\frac{2\pi r}{\lambda}$
• m: Relative refractive index of the scatterer.

Keyword arguments:

• nₘₐₓ: Maximum order of the Mie coefficients. Default to $\max(x + 4\sqrt{3}{x} + 2, |m|x)$.

Outputs:

• a, b: Mie coefficients. Both are Vector{C} with nₘₐₓ elements.

References:

• C. F. Bohren & D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley (1983).

bhcoat([T=Float64,], xᵢₙ, xₒᵤₜ, mᵢₙ, mₒᵤₜ; nₘₐₓ, tolerance = 1e-8)

Inputs:

• T: Type used for calculation. All real numbers will be stored as T, while complex numbers will be stored as C = complex(T).
• xᵢₙ: Size parameter of the inner sphere. Defined as $\frac{2\pi r}{\lambda}$
• xₒᵤₜ: Size parameter of the coated sphere. xₒᵤₜ >= xᵢₙ should hold.
• mᵢₙ: Refractive index of the inner sphere, relative to the host medium.
• mₒᵤₜ: Refractive index of the mantle, relative to the host medium.

Keyword arguments:

• nₘₐₓ: Maximum order of the Mie coefficients. Default to $\max(x_m + 4\sqrt{3}{x_m} + 2, \max(|m_c|, |m_m|)x_m)$.
• tolerance: Error tolerance. Default to 1e-8.

Outputs:

• a, b: Mie coefficients. Both are of type Vector{C}.

References:

• C. F. Bohren & D. R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley (1983).

According to Mishchenko et al. (2002), Eqs. (5.42)–(5.44), the T-Matrix for a Mie particle can be written as:

\[\begin{array}{l} T_{m n m^{\prime} n^{\prime}}^{12}(P) \equiv 0, \quad T_{m n m^{\prime} n^{\prime}}^{21}(P) \equiv 0, \\ T_{m n m^{\prime} n^{\prime}}^{11}(P)=-\delta_{m m^{\prime}} \delta_{n n^{\prime}} b_n, \\ T_{m n m^{\prime} n^{\prime}}^{22}(P)=-\delta_{m m^{\prime}} \delta_{n n^{\prime}} a_n . \end{array}\]

MieTransitionMatrix{CT, N}(x::Real, m::Number)

Generate the T-Matrix from the Mie coefficients of a homogeneous sphere.

MieTransitionMatrix{CT, N}(x_core::Real, x_mantle::Real, m_core::Number, m_mantle::Number)

Generate the T-Matrix from the Mie coefficients of a coated sphere.

This struct provides the T-Matrix API for a Mie particle.

orientation_average(mie::MieTransitionMatrix, _pₒ; _kwargs...)

The T-Matrix of a Mie scatterer is invariant under rotation. Therefore, the original T-Matrix will be returned.

rotate(mie::MieTransitionMatrix, ::Rotation{3})

The T-Matrix of a Mie scatterer is invariant under rotation. Therefore, the original T-Matrix will be returned.

A general T-Matrix $T_{m n m^{\prime} n^{\prime}}^{k l}$ stored in a 6-dimensional array, in the order $(m, n, m^{\prime}, n^{\prime}, k, l)$.

Concrete type for a general T-Matrix.

absorption_cross_section(𝐓::AbstractTransitionMatrix, λ=2π)

Calculate the absorption cross section from the given T-Matrix.

albedo(𝐓::AbstractTransitionMatrix)

Calculate the single scattering albedo from the given T-Matrix.

amplitude_matrix(𝐓::AbstractTransitionMatrix{CT, N}, ϑᵢ, φᵢ, ϑₛ, φₛ; λ=2π)

Calculate the amplitude matrix of the given T-Matrix 𝐓 at the given incidence and scattering angles.

Parameters:

• 𝐓: the T-Matrix of the scatterer.
• ϑᵢ: the incidence zenith angle in radians.
• φᵢ: the incidence azimuthal angle in radians.
• ϑₛ: the scattering zenith angle in radians.
• φₛ: the scattering azimuthal angle in radians.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.

For a general T-Matrix, Mishchenko et al. (2002), Eqs. (5.11)–(5.17), is used as a fallback.

\[\begin{array}{l} S_{11}\left(\hat{\mathbf{n}}^{\text {sca }}, \hat{\mathbf{n}}^{\text {inc }}\right)=\frac{1}{k_1} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=-n}^n \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \alpha_{m n m^{\prime} n^{\prime}}\left[T_{m n m^{\prime} n^{\prime}}^{11} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right. \\ +T_{m n m^{\prime} n^{\prime}}^{21} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right)+T_{m n m^{\prime} n^{\prime}}^{12} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right) \\ \left.+T_{m n m^{\prime} n^2}^{22} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right] \exp \left[\mathrm{i}\left(m \varphi^{\text {sca }}-m^{\prime} \varphi^{\text {inc }}\right)\right] \text {, } \\ S_{12}\left(\hat{\mathbf{n}}^{\mathrm{sca}}, \hat{\mathbf{n}}^{\mathrm{inc}}\right)=\frac{1}{\mathrm{i} k_1} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=-n}^n \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \alpha_{m n m^{\prime} n^{\prime}}\left[T_{m n m^{\prime} n^{\prime}}^{11} \pi_{m n}\left(\vartheta^{\mathrm{sca}}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right)\right. \\ +T_{m n m^{\prime} n^{\prime}}^{21} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right)+T_{m n m^{\prime} n^{\prime}}^{12} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right) \\ \left.+T_{m n m^{\prime} n^{\prime}}^{22} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right] \exp \left[\mathrm{i}\left(m \varphi^{\text {sca }}-m^{\prime} \varphi^{\text {inc }}\right)\right] \text {, } \\ S_{21}\left(\hat{\mathbf{n}}^{\mathrm{sca}}, \hat{\mathbf{n}}^{\mathrm{inc}}\right)=\frac{\mathrm{i}}{k_1} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=-n}^n \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \alpha_{m n m^{\prime} n^{\prime}}\left[T_{m n m^{\prime} n^{\prime}}^{11} \tau_{m n}\left(\vartheta^{\mathrm{sca}}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right)\right. \\ +T_{m n m^{\prime} n^{\prime}}^{21} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\mathrm{inc}}\right)+T_{m n m^{\prime} n^{\prime}}^{12} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right) \\ \left.+T_{m n m^{\prime} n^{\prime}}^{22} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right] \exp \left[\mathrm{i}\left(m \varphi^{\text {sca }}-m^{\prime} \varphi^{\text {inc }}\right)\right] \text {, } \\ S_{22}\left(\hat{\mathbf{n}}^{\text {sca }}, \hat{\mathbf{n}}^{\mathrm{inc}}\right)=\frac{1}{k_1} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=-n}^n \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \alpha_{m n m^{\prime} n^{\prime}}\left[T_{m n n^{\prime} n^{\prime}}^{11} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right. \\ +T_{m n m^{\prime} n^{\prime}}^{21} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \tau_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)+T_{m n m^{\prime} n^{12}}^{12} \tau_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right) \\ \left.+T_{m n m^{\prime} n^{\prime}}^{22} \pi_{m n}\left(\vartheta^{\text {sca }}\right) \pi_{m^{\prime} n^{\prime}}\left(\vartheta^{\text {inc }}\right)\right] \exp \left[\mathrm{i}\left(m \varphi^{\text {sca }}-m^{\prime} \varphi^{\text {inc }}\right)\right], \\ \end{array}\]

Where

\[\begin{array}{l} \alpha_{m n m^{\prime} n^{\prime}}=\mathrm{i}^{n^{\prime}-n-1}(-1)^{m+m^{\prime}}\left[\frac{(2 n+1)\left(2 n^{\prime}+1\right)}{n(n+1) n^{\prime}\left(n^{\prime}+1\right)}\right]^{1 / 2}, \\ \pi_{m n}(\vartheta)=\frac{m d_{0 m}^n(\vartheta)}{\sin \vartheta}, \quad \pi_{-m n}(\vartheta)=(-1)^{m+1} \pi_{m n}(\vartheta), \\ \tau_{m n}(\vartheta)=\frac{\mathrm{d} d_{0 m}^n(\vartheta)}{\mathrm{d} \vartheta}, \quad \tau_{-m n}(\vartheta)=(-1)^m \tau_{m n}(\vartheta) \end{array}\]

asymmetry_parameter(𝐓, λ)

Calculate the asymmetry parameter from the given transition matrix, using Mishchenko et al. (2002), Eq. (4.92):

\[\langle\cos\Theta\rangle=\frac{\alpha_1^1}{3}\]

expansion_coefficients(𝐓::AbstractTransitionMatrix{CT, N}, λ) where {CT, N}

Calculate the expansion coefficients from an arbitrary T-Matrix, using Bi & Yang (2014), Eqs. (24)–(74).

Parameters:

• 𝐓: The precalculated T-Matrix of a scatterer.
• λ: The wavelength.

Keyword arguments:

• full: Whether to return the full expansion coefficients (β₃ to β₆). Default to false.

extinction_cross_section(𝐓::AbstractTransitionMatrix{CT, N}, λ=2π) where {CT, N}

Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.102).

\[\left\langle C_{\mathrm{ext}}\right\rangle=-\frac{2 \pi}{k_1^2} \operatorname{Re} \sum_{n=1}^{\infty} \sum_{m=-n}^n\left[T_{m n n n}^{11}(P)+T_{m n m n}^{22}(P)\right]\]

Parameters:

• 𝐓: the T-Matrix of the scatterer.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.

Get the maximum order of a T-Matrix.

orientation_average(𝐓::AbstractTransitionMatrix{CT, N}, pₒ; Nα = 10, Nβ = 10, Nγ = 10) where {CT, N}

Calculate the orientation average of a transition matrix using numerical integration, given the orientation distribution function $p_o(\alpha,\beta,\gamma)$.

\[\langle T_{m n m^{\prime} n^{\prime}}^{p p^{\prime}}(L)\rangle = \int_0^{2\pi}\mathrm{d}\alpha\int_0^{\pi}\mathrm{d}\beta\sin\beta\int_0^{2\pi}\mathrm{d}\gamma p_o(\alpha,\beta,\gamma) T_{m n m^{\prime} n^{\prime}}^{p p^{\prime}}(L; \alpha,\beta,\gamma)\]

Parameters:

• 𝐓: the T-Matrix to be orientation averaged.
• pₒ: the orientation distribution function. Note that the $\sin\beta$ part is already included.
• Nα: the number of points used in the numerical integration of $\alpha$. Default to 10.
• Nβ: the number of points used in the numerical integration of $\beta$. Default to 10.
• Nγ: the number of points used in the numerical integration of $\gamma$. Default to 10.

phase_matrix(𝐒::AbstractMatrix)

Calculate the phase matrix 𝐙 from the amplitude matrix 𝐒, according to Mishchenko et al. (2002), Eqs. (2.106)–(2.121).

rotate(𝐓::AbstractTransitionMatrix{CT, N}, rot::Rotation{3})

Rotate the given T-Matrix 𝐓 by the Euler angle rot and generate a new T-Matrix.

For a general T-Matrix, Mishchenko et al. (2002), Eq. (5.29), is used as a fallback. A TransitionMatrix will be returned, which is the most general yet concrete type.

\[T_{m n m^{\prime} n^{\prime}}^{p p′}(L ; \alpha, \beta, \gamma)=\sum_{m_1=-n}^n \sum_{m_2=-n^{\prime}}^{n^{\prime}} D_{m m_1}^n(\alpha, \beta, \gamma) T_{m_1 n m_2 n^{\prime}}^{p p′}(P) D_{m_2 m^{\prime}}^{n^{\prime}}(-\gamma,-\beta,-\alpha)\quad p,p′=1,2\]

scattering_cross_section(𝐓::AbstractTransitionMatrix{CT, N}, λ=2π) where {CT, N}

Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.140).

\[\left\langle C_{\mathrm{sca}}\right\rangle=\frac{2 \pi}{k_1^2} \sum_{n=1}^{\infty} \sum_{m=-n}^n \sum_{n^{\prime}=1}^{\infty} \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \sum_{k=1}^2 \sum_{l=1}^2\left|T_{m n m^{\prime} n^{\prime}}^{k l}(P)\right|^2\]

Parameters:

• 𝐓: the T-Matrix of the scatterer.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.

scattering_matrix(𝐓, λ, θs)

Calculate expansion coefficients first and then calculate scatterering matrix elements.

Parameters:

• 𝐓: The transition matrix.
• λ: The wavelength.
• θs: The scattering angles to be evaluated in degrees.

scattering_matrix(α₁, α₂, α₃, α₄, β₁, β₂, θs)

Calculate the scatterering matrix elements from the given expansion coefficients.

Parameters:

• α₁, α₂, α₃, α₄, β₁, β₂: The precalculated expansion coefficients.
• θs: The scattering angles to be evaluated in degrees.

scattering_matrix(α₁, α₂, α₃, α₄, β₁, β₂, β₃, β₄, β₅, β₆, θs)

Calculate all 10 independent scatterering matrix elements (F₁₁, F₁₂, F₁₃, F₁₄, F₂₂, F₂₃, F₂₄, F₃₃, F₃₄, F₄₄) from the given expansion coefficients.

The scattering matrix can be expressed as:

\[\mathbf{F}=\left[\begin{array}{cccc} F_{11} & F_{12} & F_{13} & F_{14} \\ F_{12} & F_{22} & F_{23} & F_{24} \\ -F_{13} & -F_{23} & F_{33} & F_{34} \\ F_{14} & F_{24} & -F_{34} & F_{44} \end{array}\right]\]

Parameters:

• α₁, α₂, α₃, α₄, β₁, β₂, β₃, β₄, β₅, β₆: The precalculated expansion coefficients.
• θs: The scattering angles to be evaluated in degrees.

Transition matrix of an axisymmetric scatterer (spheroid, cylinder, or Chebyshev particle). Axial symmetry makes the matrix block-diagonal in the azimuthal order m, so only the per-m blocks are stored (in the field 𝐓) — a far more compact representation than the general TransitionMatrix. Indexing as T[m, n, m′, n′, p, p′] transparently returns the corresponding element of the dense layout.

expansion_coefficients(𝐓::AxisymmetricTransitionMatrix{CT, N, V, T}, λ) where {CT, N, V, T}

Calculate the expansion coefficients from an axisymmetric T-Matrix. Translated from Mishchenko et al.'s Fortran code.

Parameters:

• 𝐓: The precalculated T-Matrix of a scatterer.
• λ: The wavelength.

extinction_cross_section(axi::AxisymmetricTransitionMatrix{CT, N}, λ=2π) where {CT, N}

Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.107).

\[\left\langle C_{\text {ext }}\right\rangle=-\frac{2 \pi}{k_1^2} \operatorname{Re} \sum_{n=1}^{\infty} \sum_{m=0}^n\left(2-\delta_{m 0}\right)\left[T_{m n m n}^{11}(P)+T_{m n m n}^{22}(P)\right]\]

Parameters:

• 𝐓: the T-Matrix of the scatterer.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.

scattering_cross_section(axi::AxisymmetricTransitionMatrix{CT, N}, λ=2π) where {CT, N}

Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.141).

\[\left\langle C_{\text {sca }}\right\rangle=\frac{2 \pi}{k_1^2} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=0}^{\min \left(n, n^{\prime}\right)} \sum_{k=1}^2 \sum_{l=1}^2\left(2-\delta_{m 0}\right)\left|T_{m n m n^{\prime}}^{k l}(P)\right|^2\]

Parameters:

• 𝐓: the T-Matrix of the scatterer.
• λ: the wavelength of the incident wave in the host medium. Default to 2π.

RandomOrientationTransitionMatrix(𝐓::AbstractTransitionMatrix{CT, N}) where {CT, N}

Calculate the random-orientation T-Matrix according to Mishchenko et al. (2002), Eq. (5.96).

\[\left\langle T_{m n m^{\prime} n^{\prime}}^{k l}(L)\right\rangle=\frac{\delta_{n n^{\prime}} \delta_{m m^{\prime}}}{2 n+1} \sum_{m_1=-n}^n T_{m_1 n m_1 n}^{k l}(P), \quad k, l=1,2\]

orientation_average(ro::RandomOrientationTransitionMatrix, _pₒ; _kwargs...)

The random-orientation T-Matrix is invariant under rotation. Therefore, the original T-Matrix will be returned.

rotate(ro::RandomOrientationTransitionMatrix{CT, N, V}, ::Rotation{3}) where {CT, N, V}

The random-orientation T-Matrix is invariant under rotation. Therefore, the original T-Matrix will be returned.

Iterator for the order-degree pairs of the given maximum order nₘₐₓ.

Example of nₘₐₓ=2:

julia> collect(OrderDegreeIterator(2))
8-element Vector{Tuple{Int64, Int64}}:
 (1, -1)
 (1, 0)
 (1, 1)
 (2, -2)
 (2, -1)
 (2, 0)
 (2, 1)
 (2, 2)

vswf(kind::Symbol, m::Integer, n::Integer, k, r, ϑ, φ)

Evaluate the vector spherical wave functions $\mathbf{M}_{mn}$ and $\mathbf{N}_{mn}$ of degree m and order n at the point (r, ϑ, φ) for size parameter k (so the radial argument is x = k·r), returning the pair (𝐌, 𝐍) as SVector{3} of Cartesian field components.

kind selects the radial dependence:

• :regular — spherical Bessel $j_n$ (the regular RgM/RgN, finite at the origin; used for incident and internal fields).
• :outgoing — spherical Hankel $h_n^{(1)}$ (the radiating M/N used for the scattered field; valid outside the circumscribing sphere).

The convention is Mishchenko, Travis & Lacis (2002), Eqs. (5.16)–(5.17), consistent with amplitude_matrix. 𝐌 and 𝐍 satisfy $\nabla\times\mathbf{M}_{mn} = k\,\mathbf{N}_{mn}$ and $\nabla\times\mathbf{N}_{mn} = k\,\mathbf{M}_{mn}$.

vswf_cartesian(kind::Symbol, m, n, k, xyz::AbstractVector)

Convenience wrapper of vswf taking a Cartesian position xyz = (x, y, z) and returning the (𝐌, 𝐍) Cartesian field vectors. The polar axis is +z.

incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3}

The incident plane wave $(E_\theta\hat{\boldsymbol\vartheta} + E_\varphi\hat{\boldsymbol\varphi})\exp(\mathrm{i}k\,\hat{\mathbf n}\cdot\mathbf r)$ evaluated at the Cartesian point r⃗, propagating along n̂ = (ϑ_inc, φ_inc) with k = 2π/λ. Polarization (Eθ, Eφ) is given in the spherical basis at n̂.

internal_coefficients(x, mᵣ, ϑ_inc, φ_inc, Eθ, Eφ; nmax) -> (c, d)

Internal-field expansion coefficients (cₘₙ, dₘₙ) for a homogeneous sphere of size parameter x = k a and relative refractive index mᵣ, under an incident plane wave ((ϑ_inc, φ_inc) propagation, Jones polarization (Eθ, Eφ) in the spherical basis at the incidence direction). Returns two OffsetArrays indexed [m, n]; reuse them across many field points with internal_field.

internal_field(c, d, mᵣ, λ, r⃗) -> SVector{3}
internal_field(x, mᵣ, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗; nmax) -> SVector{3}

Internal electric field (inside a homogeneous sphere) at the Cartesian point r⃗, reconstructed from regular VSWFs at the internal wavenumber $k_\text{int} = m_r k$. The first form reuses precomputed coefficients (c, d) from internal_coefficients (efficient for a dense grid); the second computes them on the fly for size parameter x = k a and refractive index mᵣ.

r⃗ must lie INSIDE the sphere (the regular expansion represents the interior field). The external field is total_field; the two are continuous (tangentially) across the surface.

scattered_field(p, q, λ, r⃗) -> SVector{3}
scattered_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3}

Scattered electric field at the Cartesian point r⃗, reconstructed from the radiating VSWFs. The first form reuses precomputed coefficients (p, q) from scattering_coefficients (efficient for evaluating many points); the second computes them on the fly for the T-matrix 𝐓 and incidence.

The point r⃗ must lie OUTSIDE the sphere circumscribing the particle (the radiating series diverges inside it).

scattering_coefficients(𝐓, ϑ_inc, φ_inc, Eθ, Eφ) -> (p, q)

Scattered-field expansion coefficients (pₘₙ, qₘₙ) for the T-matrix 𝐓 under an incident plane wave propagating along (ϑ_inc, φ_inc) with Jones polarization (Eθ, Eφ) in the spherical basis at the incidence direction. Returns two OffsetArrays indexed [m, n] (m ∈ -N:N, n ∈ 1:N), where (p, q) = 𝐓 (a, b) and (a, b) are the incident plane-wave coefficients. Reuse these across many field points (see scattered_field).

total_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3}

Total external field 𝐄_inc(r⃗) + 𝐄_sca(r⃗) at the Cartesian point r⃗ for the T-matrix 𝐓 under the incident plane wave (see incident_field and scattered_field).

internal_coefficients(shape::AbstractAxisymmetricShape, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ) -> (c, d)

Internal-field expansion coefficients (cₘₙ, dₘₙ) for a general axisymmetric particle, obtained from the EBCM matrices: per azimuthal block $\mathbf{c} = \tfrac{1}{2}\mathbf{Q}^{-1}\mathbf{a}$ with $\mathbf{Q} = \mathbf{P} + \mathrm{i}\mathbf{U}$, under an incident plane wave ((ϑ_inc, φ_inc) propagation, Jones polarization (Eθ, Eφ)). Returns OffsetArrays indexed [m, n]; reuse them across field points with internal_field.

The interior expansion is guaranteed within the inscribed sphere (radius rmin(shape)) and is empirically stable throughout the interior; the practical limit is EBCM conditioning at high aspect ratio with high nmax (keep nmax/aspect where the T-matrix itself is reliable). See internal_field.

internal_field(shape::AbstractAxisymmetricShape, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3}

Internal electric field inside a general axisymmetric particle at the Cartesian point r⃗, via the EBCM internal coefficients (see internal_coefficients). The relative refractive index is taken from shape, and r⃗ should lie inside the particle. Convergence is guaranteed within the inscribed sphere (radius rmin(shape)) and is empirically stable throughout the interior; the binding constraint is EBCM conditioning at high aspect × nmax, not the inscribed-sphere radius.

AbstractLinearizationBackend

Abstract supertype for analytical or linearized differentiation backends. Concrete backends declare which solvers can provide a Jacobian instead of falling back to numerical differentiation.

EBCMLinearization()

Backend marker for linearized EBCM transition matrices. The current supported slice uses analytical P/U derivatives for fixed-order, fixed-quadrature Spheroid, Chebyshev, and Cylinder problems, then propagates those derivatives through the analytical EBCM matrix identity.

IITMLinearization([variant])

Backend marker for future analytical derivatives of IITM transition matrices. variant distinguishes :axisymmetric, :nfold, :arbitrary, or :auto.

LinearizationProblem(rebuild, x; variables)
LinearizationProblem(x; variables, rebuild)

Describe a real parameter vector and the function that rebuilds the physical scattering input from that vector. variables must name each entry of x.

This mirrors the existing numerical differentiation pattern where users write x -> shape, wavelength, config, while giving analytical backends a stable problem object to dispatch on.

LinearizationResult(value, jacobian, variables; metadata=(;))

Container for a computed value and its Jacobian with respect to a LinearizationProblem parameter vector.

LinearizationSupport(supported, reason)

Capability record returned by supports_linearization.

MieLinearization()

Backend marker for future analytical derivatives of Mie transition matrices.

UnsupportedLinearization

Exception thrown when the unified linearization API is asked for an analytical Jacobian that has not been implemented for the chosen problem/backend/output.

linearize_observable(observable, problem, backend; config=nothing)

Compute a scattering observable and its Jacobian with respect to problem. The default method only defines the unsupported boundary; analytical backends should add methods that reuse linearize_transition_matrix.

linearize_transition_matrix(problem, backend; config=nothing)

Compute a transition matrix and its Jacobian with respect to problem. Backends that do not implement analytical derivatives throw UnsupportedLinearization.

supports_linearization(problem, backend; output=:transition_matrix, config=nothing)

Return a LinearizationSupport record describing whether an analytical linearization implementation is available for the requested output.

ricattibesselj!(ψₙ, ψₙ′, z, nₘₐₓ, nₑₓₜᵣₐ, x)

Parameters:

• ψ: output array for $\psi_n(x)$.
• ψ′: output array for $\psi^{\prime}_n(x)$.
• z: auxiliary array for downward recursion.
• nₘₐₓ: maximum order of $\psi_n(x)$.
• nₑₓₜᵣₐ: extra terms for downward recursion of $z$.
• x: argument.

In-place version of ricattibesselj.

ricattibesselj(nₘₐₓ::Integer, nₑₓₜᵣₐ::Integer, x)

Calculate Ricatti-Bessel function of the first kind, $\psi_n(x)$, and its derivative $\psi^{\prime}_n(x)$ for $1\le n\le n_{\max}$.

$\psi_n(x)$ is defined as

\[\psi_n(x) = xj_n(x)\]

where $j_n(x)$ is the Bessel function of the first kind.

Parameters:

• nₘₐₓ: maximum order of $\psi_n(x)$.
• nₑₓₜᵣₐ: extra terms for downward recursion of $z$.
• x: argument.

Returns: (ψ, ψ′)

• ψ: array of $\psi_n(x)$.
• ψ′: array of $\psi^{\prime}_n(x)$.

ricattibessely!(χ, χ′, nₘₐₓ, x)

Parameters:

• χ: output array for $\chi_n(x)$.
• χ′: output array for $\chi^{\prime}_n(x)$.
• nₘₐₓ: maximum order of $\chi_n(x)$.
• x: argument.

In-place version of ricattibessely.

ricattibessely(nₘₐₓ::Integer, x)

Calculate Ricatti-Bessel function of the second kind, $\chi_n(x)$, and its derivative $\chi^{\prime}_n(x)$ for $1\le n\le n_{\max}$.

$\chi_n(x)$ is defined as

\[\chi_n(x) = xy_n(x)\]

where $y_n(x)$ is the Bessel function of the second kind.

Parameters:

• nₘₐₓ: maximum order of $\chi_n(x)$.
• x: argument.

Returns: (χ, χ′)

• χ: Array of $\chi_n(x)$.
• χ′: Array of $\chi^{\prime}_n(x)$.

factorial([T=Float64,], n)

Calculate factorials $n!$. For $n\le20$, the standard Base.factorial is used. For $n>20$, Arblib.gamma! is used instead.

gausslegendre([T=Float64,], n)

Calculate the n-point Gauss-Legendre quadrature nodes and weights.

• For Float64, we use the $\mathcal{O}(n)$ implementation from FastGaussQuadrature.jl.
• For other types, we use the arbitrary precision implementation from Arblib.jl, then convert the results to the desired type.

pi_func([T=Float64,], m::Integer, n::Integer, ϑ::Number; d=nothing)

Calculate

\[\pi_{m n}(\vartheta)=\frac{m d_{0 m}^n(\vartheta)}{\sin \vartheta}\]

• If d is given, it is used as the value of $d_{0 m}^n(\vartheta)$.

tau_func([T=Float64,], m::Integer, n::Integer, ϑ::Number)

Calculate

\[\tau_{m n}(\vartheta)=\frac{\mathrm{d} d_{0 m}^n(\vartheta)}{\mathrm{d} \vartheta}\]

wigner_D(::Type{T}, m::Integer, m′::Integer, n::Integer, α::Number, β::Number, γ::Number) where {T}

Calculate the Wigner D-function $D^j_{mn}(\theta)$, which is defined as:

\[D_{m m^{\prime}}^{n}(\alpha, \beta, \gamma)=\mathrm{e}^{-\mathrm{i} m \alpha} d_{m m^{\prime}}^{n}(\beta) \mathrm{e}^{-\mathrm{i} m^{\prime} \gamma}\]

where

\[0 \leq \alpha<2 \pi, \quad 0 \leq \beta \leq \pi, \quad 0 \leq \gamma<2 \pi\]

wigner_D_recursion!(d::AbstractVector{CT}, m::Integer, m′::Integer, nmax::Integer, α::Number, β::Number, γ::Number)

Calculate the Wigner D-function recursively, in place.

wigner_D_recursion([T=Float64,], m::Integer, m′::Integer, nmax::Integer, α::Number, β::Number, γ::Number)

Calculate the Wigner D-function recursively (use wigner_d_recursion).

wigner_d([T=Float64,], m::Integer, n::Integer, s::Integer, ϑ::Number) where {T}

Calculate Wigner (small) d-function $d_{mn}^s(\theta)$ for a single $(m, n, s)$ combination, using Mishchenko et al. (2002), Eq. (B.1).

\[\begin{aligned} d_{m n}^{s}(\vartheta)=& \sqrt{(s+m) !(s-m) !(s+n) !(s-n) !} \\ & \times \sum_{k=\max(0,m-n)}^{\min(s + m, s - n)}(-1)^{k} \frac{\left(\cos \frac{1}{2} \vartheta\right)^{2 s-2 k+m-n}\left(\sin \frac{1}{2} \vartheta\right)^{2 k-m+n}}{k !(s+m-k) !(s-n-k) !(n-m+k) !} \end{aligned}\]

wigner_d_recursion!(d::AbstractVector{T}, m::Integer, n::Integer, sₘₐₓ::Integer, ϑ::Number; deriv=nothing)

Calculate the Wigner d-function recursively, in place.

wigner_d_recursion([T=Float64,], m::Integer, n::Integer, sₘₐₓ::Integer, ϑ::Number; deriv::Bool = false) where {T}

Calculate Wigner (small) d-function $d_{mn}^s(\theta)$ for $s\in[s_{\min}=\max(|m|, |n|),s_{\max}]$ (alternatively, its derivative as well) via upward recursion, using Mishchenko et al. (2002), Eq. (B.22).

\[\begin{aligned} d_{m n}^{s+1}(\vartheta)=& \frac{1}{s \sqrt{(s+1)^{2}-m^{2}} \sqrt{(s+1)^{2}-n^{2}}}\left\{(2 s+1)[s(s+1) x-m n] d_{m n}^{s}(\vartheta)\right.\\ &\left.-(s+1) \sqrt{s^{2}-m^{2}} \sqrt{s^{2}-n^{2}} d_{m n}^{s-1}(\vartheta)\right\}, \quad s \geq s_{\min } \end{aligned}\]

The initial terms are given by Mishchenko et al. (2002), Eqs. (B.23) and (B.24).

\[\begin{array}{l} d_{m n}^{s_{\min }-1}(\vartheta)=0 \\ d_{m n}^{s_{\min }}(\vartheta)=\xi_{m n} 2^{-s_{\min }}\left[\frac{\left(2 s_{\min }\right) !}{(|m-n|) !(|m+n|) !}\right]^{1 / 2}(1-x)^{|m-n| / 2}(1+x)^{|m+n| / 2} \end{array}\]

where

\[\xi_{m n}=\left\{\begin{array}{ll} 1 & \text { for } n \geq m \\ (-1)^{m-n} & \text { for } n<m \end{array}\right.\]