vSmartMOM

Appearance

Scattering API ​

Scattering computes aerosol optical properties, Fourier/Greek coefficients, and phase-matrix reconstructions used by the RT core.

Mie and Aerosol Optics ​

vSmartMOM.Scattering.make_mie_model Function
julia
make_mie_model(::NAI2, aerosol::Aerosol, λ, polarization, truncation_type, r_max, nquad_radius;
               r_min = DEFAULT_R_MIN, architecture = Architectures.CPU(),
               precision_policy = nothing) -> MieModel

Construct a MieModel configured for the Siewert NAI-2 workflow.

Arguments

  • aerosol: aerosol size-distribution and refractive-index specification.

  • λ: wavelength (must use the same length units as r_max and the aerosol radius scale).

  • polarization: one of Stokes_I, Stokes_IQ, Stokes_IQU, Stokes_IQUV.

  • truncation_type: typically δBGE.

  • r_max: upper radius used in size quadrature.

  • nquad_radius: number of radius quadrature points.

Keyword arguments

  • r_min: lower radius of the size quadrature (default DEFAULT_R_MIN, 1e-3 µm). The radius nodes are fixed by r_min, r_max and nquad_radius and do not move with the size distribution; see radius_quadrature. Requires 0 < r_min < r_max.

  • architecture: Architectures.CPU() (default) or Architectures.GPU(). Selects the CPU vs GPU Mie path at compute_aerosol_optical_properties time.

  • precision_policy: GPU Dₙ-recursion precision (NativeFloat64/DSEmulated). nothing (default) auto-selects on the GPU and is ignored on the CPU.

Returns

source
julia
make_mie_model(::PCW, aerosol::Aerosol, λ, polarization, truncation_type, r_max, nquad_radius,
               wigner_filepath::String; r_min = DEFAULT_R_MIN, ...) -> MieModel

Construct a MieModel configured for the Domke PCW workflow, loading Wigner tables from wigner_filepath via load_wigner_values. The keyword arguments are those of the NAI-2 method.

source
julia
make_mie_model(::PCW, aerosol::Aerosol, λ, polarization, truncation_type, r_max, nquad_radius,
               wigner_A, wigner_B; r_min = DEFAULT_R_MIN, ...) -> MieModel

Construct a MieModel configured for the Domke PCW workflow, using precomputed Wigner tables wigner_A and wigner_B. The keyword arguments are those of the NAI-2 method.

Use this overload when Wigner tensors are already in memory.

source
vSmartMOM.Scattering.DEFAULT_R_MIN Constant
julia
DEFAULT_R_MIN

Default lower bound of the radius quadrature, 1e-3 µm (a radius of 1 nm, the smallest atmospheric aerosol particles).

The bound is fixed: the radius nodes depend only on r_min, r_max and nquad_radius, never on the size-distribution parameters, so the analytic size-parameter tangents are exact derivatives of the discretized optics. For the size distributions in the shipped configurations, at most 1e-14 of the particle number lies below 1e-3 µm; the extinction cross-section falls at least as fast as r² towards small radii, so their optical contribution is smaller still.

source
vSmartMOM.Scattering.radius_quadrature Function
julia
radius_quadrature(model::MieModel, ::Type{T}=typeof(model.r_max); norm=false) -> (r, wᵣ)

Radius nodes and weights of model's size-distribution quadrature: gauleg_log with model.nquad_radius nodes on [model.r_min, model.r_max]. The bounds are converted to T; the nodes and weights come back in promote_type(T, Float64) (the Gauss–Legendre abscissae are Float64).

The nodes depend only on r_min, r_max and nquad_radius, not on the size distribution, so a change of the distribution parameters changes only the weights (compute_wₓ); the linearized NAI-2 path differentiates exactly those weights. Every log-normal Mie path (NAI-2 CPU and GPU, linearized NAI-2, PCW) builds its nodes here.

source
vSmartMOM.Scattering.compute_aerosol_optical_properties Function
julia
compute_aerosol_optical_properties(model::MieModel{<:NAI2,FT}, FT_out::Type) -> AerosolOptics

Compute bulk aerosol optical properties with the Siewert NAI-2 formulation on the CPU.

This method is the internal two-argument form called by the architecture/method router (_dispatch_aerosol_optics) in phase_function_autodiff.jl. The router supplies FT_out = _mie_output_type(model) (the model's own FT type parameter) so that a Float32 model produces Float32 greek coefficients. For architecture = GPU() the GPU KernelAbstractions pipeline is dispatched instead (see compute_aerosol_optical_properties_gpu).

The second argument FT_out sets the output float type of the returned AerosolOptics (greek coefficients, ω̃, k). It is supplied by the router and must be passed explicitly — there is no default.

Internal precision: the Mie cross-section sums (S₁, S₂, Cₑₓₜ, Cₛ arrays) and the Greek-coefficient angular quadrature always run in Float64 on the CPU path, regardless of FT. This is intentional — the Dₙ downward continued-fraction recursion and the size-integral sums are catastrophically unstable in Float32 for large or absorbing particles (|y| ≳ 50). The final convert.(FT_out, …) call narrows back to the requested output precision.

Reference:

  • S. Sanghavi, Revisiting the Fourier expansion of Mie scattering matrices in generalized spherical functions, JQSRT 136 (2014), 16-27. https://doi.org/10.1016/j.jqsrt.2013.12.015

  • Cross-section expressions are the standard Mie-series forms (paper Eq. (1)); Greek-coefficient Fourier expansion follows the Siewert/Domke framework discussed around paper Eq. (17).

For each radius quadrature point, Mie coefficients are used to compute extinction and scattering cross-sections:

with  . Bulk values are quadrature-weighted means over the size distribution. The returned optical scalars are:

Greek coefficients (α, β, γ, δ, ϵ, ζ) are then derived from the bulk phase matrix moments (Sanghavi, 2014).

Returns

Notes

  • Uses the convention nᵢ >= 0.

  • Radius and wavelength units must be consistent.

source
julia
compute_aerosol_optical_properties(lin::LinMode, model::MieModel{FDT,FT}, FT2::Type) where {FDT<:NAI2, FT}

Reference: Suniti Sanghavi 2014, https://doi.org/10.1016/j.jqsrt.2013.12.015

Compute the aerosol optical properties and their derivatives (w.r.t. nᵣ, nᵢ and the native LogNormal log-radius location and width) using the Siewert-NAI2 linearized method.

The radius nodes (radius_quadrature) are fixed by model.r_min, model.r_max and model.nquad_radius, so the size directions act only on the normalized quadrature weights. The returned tangents are therefore the exact derivatives of the discretized forward optics of the same MieModel at any nquad_radius; their accuracy as derivatives of the continuous size integral follows the convergence of the forward quadrature.

The third argument FT2 sets the output float type for both the forward Greek coefficients and the derivative (lin) greek coefficients. It must be passed explicitly — there is no default. The internal computation always runs in Float64 for numerical stability of the Dₙ recursion (same rationale as the forward path).

Input: LinMode, MieModel Output: (AerosolOptics, linAerosolOptics)

source
julia
compute_aerosol_optical_properties(model::MieModel{<:PCW,FT}, FT2::Type) -> AerosolOptics

Compute bulk aerosol optical properties with the Domke PCW formulation (Sanghavi, 2014), using precomputed Wigner tables.

The second argument FT2 sets the output float type of the returned AerosolOptics. It is supplied by the router (_dispatch_aerosol_optics in phase_function_autodiff.jl) and must be passed explicitly — there is no default.

Reference:

The method evaluates terms (Eq. 22 in Sanghavi, 2014), maps them to Greek coefficients (Eq. 24), and returns AerosolOptics with:

Use make_mie_model with PCW() and Wigner inputs before calling this function.

There is no GPU PCW kernel; a GPU()-architecture PCW model falls back to this CPU implementation (with a one-time @warn) via the single-verb compute_aerosol_optical_properties(model) router.

source
julia
compute_aerosol_optical_properties(model::MieModel; autodiff=false) -> AerosolOptics

Single-verb entry point for aerosol optics. Everything (Fourier-decomposition method, compute architecture, GPU precision policy) is read off the model.

  • autodiff=false (default): dispatches on model.computation_type and model.architecture:

    • NAI2 + CPU() → the analytic CPU implementation (compute_aerosol_optical_properties 2-arg method).

    • NAI2 + a GPU architecture with a Mie pipeline (has_gpu_mie true, CUDA GPU() or Float32 MetalGPU()) → the KernelAbstractions GPU pipeline (compute_aerosol_optical_properties_gpu).

    • NAI2/PCW + a non-CPU architecture without a GPU Mie pipeline → one-time @warn then CPU fallback (Mie on CPU; RT arrays still run on that architecture).

    • PCW + GPU → one-time @warn then CPU fallback (no GPU PCW kernel).

    • PCW + CPU() → the analytic PCW implementation.

  • autodiff=true is retired and raises an error. For aerosol Jacobians use compute_aerosol_optical_properties(LinMode(), model, FT), which returns forward and hand-linearized optics. Its native size coordinates are the lognormal distribution parameters, not median radius/geometric width.

source
vSmartMOM.Scattering.compute_aerosol_optical_properties_nodes Function
julia
compute_aerosol_optical_properties_nodes(radii, weights, n_real, n_imag,
    wavelength, polarization, truncation;
    architecture = Architectures.CPU(),
    l_max = nothing,
    precision_policy = nothing) -> AerosolOptics

Bulk Mie aerosol optics for a CALLER-DEFINED set of size-distribution nodes, using the Siewert NAI-2 formulation (same physics/conventions as compute_aerosol_optical_properties for a log-normal MieModel, and the same AerosolOptics/GreekCoefs output shape).

This is the intended seam for sectional or tabulated size distributions that are not naturally expressed as a Distributions.jl continuous distribution — e.g. GCHPIO's TOMAS two-moment sectional bins. No log-normal (or any other analytic) size distribution is assumed anywhere in this path.

Arguments

  • radii::AbstractVector: wet particle radii, one per node, [μm].

  • weights::AbstractVector: non-negative node weights (number-mixing-ratio or number-concentration weights). Need not be normalized — they are divided by their sum internally, so the returned bulk quantities are number-mean properties and are invariant to a uniform rescaling of weights (e.g. weights .* 7.3 gives the identical AerosolOptics).

  • n_real, n_imag: the ONE complex refractive index shared by the whole node ensemble, using vSmartMOM's convention n = n_real - i·n_imag (n_imag ≥ 0; positive n_imag is absorption — same convention as Aerosol and asserted identically).

  • wavelength: wavelength, same units as radii (μm internally).

  • polarization::AbstractPolarizationType: accepted for API parity with make_mie_model / MieModel; as in the existing NAI2 path, all six Greek coefficient arrays are always computed regardless of polarization type (the type does not currently gate the computation — matching compute_aerosol_optical_properties's existing behavior).

  • truncation::AbstractTruncationType: applied to the raw bulk Greek coefficients via truncate_phase, exactly mirroring how the RT pipeline (model_from_parameters/update_model) applies a MieModel's truncation_type downstream of the raw Mie call: a δBGE is fit only if it would shorten the natural (node-driven) Greek series; otherwise (including NoTruncation) the untruncated series passes through with fᵗ = 0.

Keyword arguments

  • architecture: Architectures.CPU() (default) or a GPU architecture with a registered Mie pipeline (Architectures.has_gpu_mie, e.g. CUDA GPU()) — same architecture routing as compute_aerosol_optical_properties(model). Non-CPU architectures without a GPU Mie pipeline warn once and fall back to CPU (mirroring the existing MieModel router).

  • l_max: optional cap on the number of OUTPUT Greek terms, applied by plain array slicing (min(l_max, series_length)) AFTER truncation is applied — a δBGE fit always sees the full natural series, and the cap then shortens whatever the truncation produced. This is an efficiency knob, independent of truncation: because angular resolution here scales with the largest supplied node (no r_max cap), a very wide node ensemble can produce a long natural Greek series; l_max lets a caller who does not need the full series avoid paying for it downstream, without engaging δBGE's forward-peak re-fit. nothing (default) keeps the full length.

  • precision_policy: GPU-only, see MiePrecisionPolicy (NativeFloat64/DSEmulated/NativeFloat32); nothing (default) auto-selects via Architectures.default_mie_precision_policy exactly as the MieModel GPU path does (NativeFloat32 on MetalGPU(), since Metal has no Float64 hardware at all). Ignored on the CPU path.

Angular resolution

n_max is get_n_max(maximum(2π .* radii ./ wavelength)) — i.e. derived from the size parameter of the ACTUAL largest supplied radius, with no artificial cap. n_mu = 2n_max - 1, exactly as the existing NAI2 path computes it from r_max when the size distribution comes from a log-normal MieModel.

Returns

AerosolOptics with greek_coefs, number-mean ω̃, number-mean k, and fᵗ from the requested truncation.

source
vSmartMOM.Scattering.compute_aerosol_extinction_nodes Function
julia
compute_aerosol_extinction_nodes(radii, weights, n_real, n_imag, wavelength;
    architecture = Architectures.CPU(),
    precision_policy = nothing) -> k

EXTINCTION-ONLY fast path for the caller-node Mie seam: returns just the number-mean bulk extinction cross-section k [μm²] — the SAME quantity, same convention, and the SAME arithmetic as .k from compute_aerosol_optical_properties_nodes — without computing any amplitude matrices, phase-matrix elements, or the six Greek coefficient series. This equality is EXACT (not merely close) on the CPU path, on CUDA (real hardware and the KA-CPU KernelAbstractions backend), and on the KA-CPU backend under every precision policy — see "Numerics" below for exactly what "exact" means on each path and how it is tested.

Motivation

A downstream hybrid consumer (GCHPIO's :parametric_exact_tau mode) needs only the number-mean extinction cross-section from ONE node set while taking shapes (Greek coefficients) from another. Today it would have to call the full compute_aerosol_optical_properties_nodes API and discard the amplitude matrices and all six Greek series, which dominate the cost at large angular orders. This function computes only what k actually depends on: the per-node Mie coefficients (aₙ, bₙ) and a weighted sum over nodes — no Gauss-Legendre angular quadrature, π/τ functions, amplitude functions, or Greek-coefficient projection anywhere in this path.

Arguments

Same node-ensemble arguments as compute_aerosol_optical_properties_nodes, EXCEPT there is no polarization/truncation/l_max — none of these are relevant to a bare extinction cross-section (there is no phase matrix or Greek series to polarize, truncate, or length-cap).

  • radii::AbstractVector: wet particle radii, one per node, [μm].

  • weights::AbstractVector: non-negative node weights (number-mixing-ratio or number-concentration weights). Need not be normalized — divided by their sum internally, so k is invariant to a uniform rescaling of weights.

  • n_real, n_imag: the ONE complex refractive index shared by the whole node ensemble, convention n = n_real - i·n_imag (n_imag ≥ 0).

  • wavelength: same units as radii (μm internally).

Keyword arguments

  • architecture: Architectures.CPU() (default) or a GPU architecture with a registered Mie pipeline (Architectures.has_gpu_mie) — identical routing to compute_aerosol_optical_properties_nodes (CPU / GPU / warn-once CPU fallback for an architecture without a GPU Mie pipeline).

  • precision_policy: GPU-only, see MiePrecisionPolicy (NativeFloat64/DSEmulated/NativeFloat32); ignored on the CPU path.

Returns

The number-mean bulk extinction cross-section k [μm²], as a scalar of the promoted float type of the inputs.

Validation

Identical to compute_aerosol_optical_properties_nodes: _validate_node_inputs plus sum(weights) > 0.

Numerics

CPU path: internal computation runs in IC = Float64 (same IC/FT narrowing convention as the rest of this module). Per node, this computes exactly the SAME C_ext expression the shared NAI-2 core uses — 2π/k² * dot(n_v, real.(an_v .+ bn_v)) (the internal _nai2_bulk_optics core's C_ext formula, via the shared _mie_node_ab_and_C_ext! helper) — in the SAME arithmetic order, weighted by the SAME normalized weights (same normalization step), so k is BIT-IDENTICAL (===) to compute_aerosol_optical_properties_nodes(...).k given identical (radii, weights, n_real, n_imag, wavelength).

GPU path (compute_aerosol_extinction_nodes_gpu): shares its preamble (validation, FT/RA, per-node arrays, Kernel 1) with compute_aerosol_optical_properties_nodes_gpu via the internal _prepare_node_mie_gpu helper, so both entry points run cross_sections_kernel! and the single-thread (ndrange=1) weighted_sum_kernel! reduction on IDENTICAL per-node/weight device arrays. k is confirmed EXACT-EQUALITY (==, not merely isapprox) tested for all three precision policies (NativeFloat64, DSEmulated, NativeFloat32) on both the standard and wide-range TOMAS-like node sets — CPU path (===), the KA-CPU KernelAbstractions backend (==), and real CUDA hardware (==), all in test/test_mie_nodes.jl. On MetalGPU() (NativeFloat32 only, see MiePrecisionPolicy) the same exact-equality test is wired in and expected to hold by the same reasoning, but is unverified on real Metal hardware (this package's CI has none); the test is written to activate automatically the first time it runs on a Mac with a functional Metal device.

source
vSmartMOM.Scattering.truncate_phase Function
julia
truncate_phase(::NoTruncation, aero::AerosolOptics; kwargs...) -> AerosolOptics

Identity passthrough. Returns an AerosolOptics with the same Greek coefficients, ω̃ and k, and fᵗ = 0 (the f_tr → 0 limit of Sanghavi & Stephens 2015 Eq. 8 — the truncation-modified τ*, ω*, Z* collapse to the originals).

Note: raw Mie outputs initialise fᵗ = 1 as a "untruncated yet" sentinel — passing them through unchanged would let downstream delta_m_forward interpret the 1 as "everything is in the forward peak" and silently zero out aerosol scattering. We return fᵗ = 0 so the rescaling is a true no-op.

source
julia
truncate_phase(mod::δBGE, aero::AerosolOptics; reportFit=false) -> AerosolOptics

Apply δ-BGE truncation to aerosol Greek coefficients.

The method removes/approximates the forward peak using a least-squares fit over all quadrature angles, then renormalizes with retained scattering fraction :

All six families share the renormalisation (SS2015 Eq. 8:     ). Only the diagonal families additionally have the forward-peak amount   removed — the Dirac spike is diagonal, so    (S2014 Eqs. A.5–A.10; cf. SS2015 Eqs. 27c,d for the δ-m analogue). This keeps the degree of linear polarisation invariant under truncation outside the exclusion cone.

Returns a new AerosolOptics with truncated coefficients and updated fᵗ. The common 1/c₀ normalization is required because all three fits target the unnormalized original matrix outside the forward cone; downstream δ-scaling multiplies the smooth remainder by c₀ = 1-fᵗ.

source
julia
truncate_phase(mod::δBGE, aero::AerosolOptics, lin_aero::linAerosolOptics; reportFit=false)

Apply δ-BGE truncation to aerosol optics and their linearized derivatives (Jacobians).

Performs the same least-squares fit as truncate_phase(mod::δBGE, aero) for β, γ, ϵ, then propagates derivatives through the truncation to obtain truncated Greek-coefficient Jacobians. Reference: Sanghavi & Stephens (2015).

All six truncated families carry the 1/c₀ renormalisation of SS2015 Eq. 8 (c₀ = 1 - fᵗ), so each derivative picks up the matching -xᵗ·ċ₀ chain-rule term with ċ₀ = ẋβ[:,1].

Arguments

Returns

  • (AerosolOptics, linAerosolOptics): truncated optics and their derivatives
source
julia
truncate_phase(::NoTruncation, aero::AerosolOptics, lin_aero::linAerosolOptics; kwargs...)
    -> (AerosolOptics, linAerosolOptics)

Identity passthrough for the linearised (Jacobian) path. Resets the raw Mie sentinel fᵗ = 1 (and the corresponding ḟᵗ) to zero, matching the two-argument truncate_phase(::NoTruncation, aero) so delta_m_truncation_lin doesn't silently zero the aerosol SSA Jacobian when params.truncation = NoTruncation().

source

Batched Caller-Node Mie Seam ​

Exact-nmax-grouped GPU batching + reusable single-owner workspace for many-ensemble columns (e.g. GCHPIO's per-layer TOMAS inventory). See proposals/batched_mie_nodes_seam.md for the full design rationale.

vSmartMOM.Scattering.prepare_mie_node_geometry Function
julia
prepare_mie_node_geometry(radii, weights, offsets;
    architecture = Architectures.CPU(), backend = nothing) -> MieNodeGeometry

Build the λ-independent, immutable MieNodeGeometry backing the batched caller-node Mie API. Concatenates E ensembles' nodes (E = length(offsets) - 1) into flat radii/weights arrays, normalizes each ensemble's weights to sum to 1 (Float64, then narrowed — the same "P2 lesson" convention as compute_aerosol_optical_properties_nodes_gpu: normalizing in Float64 before narrowing avoids a Float32 sum of large-but-finite weights overflowing to Inf), and (for GPU architectures) uploads the result once.

Arguments

  • radii::AbstractVector: concatenated wet particle radii for ALL ensembles, one per node, [μm], ensemble-contiguous (ensemble e's nodes occupy radii[offsets[e]+1:offsets[e+1]]).

  • weights::AbstractVector: concatenated non-negative node weights, same layout as radii. Need not be normalized per ensemble — normalized internally, exactly as the single-ensemble API's weights are.

  • offsets::AbstractVector{<:Integer}: length E+1, offsets[1] == 0, non-decreasing, offsets[end] == length(radii). Every ensemble needs at least one node (offsets[e+1] > offsets[e]).

Keyword arguments

  • architecture: Architectures.CPU() (default: GFT = Float64 always, matching the CPU-threaded path's IC = Float64 convention exactly, no device copy made) or a GPU architecture with a registered Mie pipeline (Architectures.has_gpu_mie; GFT = float(promote_type(eltype(radii), eltype(weights))), i.e. the caller's own choice — typically Float32 for NativeFloat32/DSEmulated). On MetalGPU(), radii/weights must be Float32 (or narrower) — Float64 throws ArgumentError immediately (Metal has no Float64 hardware support at all; the batched device arrays would be illegal), mirroring the single-ensemble _prepare_node_mie_gpu's Metal guard.

  • backend: nothing (default) derives the device-array backend from architecture (Architectures.ka_backend) whenever architecture is NON-CPU, exactly like every other GPU entry point in this module. Pass an EXPLICIT KA backend (e.g. KernelAbstractions.CPU()) to override this — the main use case is testing a GPU-style geometry (architecture non-CPU, so GFT is the caller's Float32/Float64 choice rather than the CPU path's forced Float64) on a portable backend without real GPU hardware, mirroring how compute_aerosol_optical_properties_nodes_gpu takes an explicit backend decoupled from any "architecture" concept. Whichever backend is resolved here (the override, or the architecture-derived default) is STORED on the returned geometry (MieNodeGeometry.backend) and reused verbatim by every later call, INCLUDING the top-level compute_aerosol_optical_properties_nodes_batched/ compute_aerosol_extinction_nodes_batched dispatchers — they launch kernels on this stored backend rather than re-deriving one from architecture, so an overridden geometry stays internally consistent no matter which entry point (top-level or the explicit-backend _gpu sibling) is used on it afterward. Ignored (must be left nothing) when architecture isa Architectures.CPU.

Returns

MieNodeGeometry, reusable across arbitrarily many compute_aerosol_optical_properties_nodes_batched / compute_aerosol_extinction_nodes_batched calls at different wavelengths (grouping is recomputed fresh each call — see those functions' docstrings — but the underlying node data is never re-normalized or re-uploaded).

source
vSmartMOM.Scattering.MieNodeGeometry Type
julia
MieNodeGeometry{GFT, ARCH, RV, IV}

Immutable, λ-INDEPENDENT, device-resident (for GPU architectures) caller-node size-distribution geometry backing compute_aerosol_optical_properties_nodes_batched / compute_aerosol_extinction_nodes_batched. Built once via prepare_mie_node_geometry and reused across many per-λ batched calls (and, in the intended GCHPIO usage, across many spectral bands) without re-uploading or re-normalizing anything.

Fields (all internal; use the accessors / pass the geometry object around,

do not construct or mutate directly)

  • architecture: fixed at construction; determines GFT and whether radii/ weights/node_to_ensemble are device-resident copies (GPU) or identical to the host arrays (CPU — no copy is made).

  • backend: the ACTUAL KA backend radii/weights/node_to_ensemble were built on (nothing for Architectures.CPU()) — may be an explicit override passed to prepare_mie_node_geometry's backend keyword, NOT necessarily Architectures.ka_backend(architecture). The top-level batched dispatchers (compute_aerosol_optical_properties_nodes_batched/ compute_aerosol_extinction_nodes_batched) launch kernels on THIS stored backend rather than re-deriving one from architecture, so an overridden geometry stays internally consistent no matter which entry point is used.

  • output_FT: the CALLER's ORIGINAL radii/weights promoted element type (float(promote_type(eltype(radii), eltype(weights))), captured BEFORE any internal-compute forcing) — this is what the batched CPU output narrows to (together with n_real/n_imag/wavelength's own eltypes), exactly mirroring compute_aerosol_optical_properties_nodes's own output-type formula. Distinct from GFT for Architectures.CPU() geometries, where GFT is forced to Float64 for internal compute regardless of the caller's actual radii/weights eltype (e.g. an all-Float32 CPU call must still return Float32, not Float64) — for GPU architectures output_FT === GFT trivially, since GFT is never forced there.

  • n_ensembles, n_nodes: E and total node count.

  • offsets::Vector{Int}: length E+1, ALWAYS host-resident (needed for per-λ host-side grouping and the CPU threaded loop regardless of architecture); offsets[1] == 0, offsets[end] == n_nodes.

  • radii_host::Vector{GFT}: host mirror of the concatenated, narrowed radii — needed every λ-call for the host-side per-ensemble n_max_global computation that drives GPU grouping, even when radii itself is device-resident.

  • radii::RV, weights::RV: concatenated, narrowed-to-GFT radii and PRE-NORMALIZED (per ensemble, in Float64, then narrowed — see prepare_mie_node_geometry) weights. RV === Vector{GFT} (same object as radii_host/the host weight mirror) for Architectures.CPU(); a device array type (CuArray/MtlArray) for GPU architectures.

  • node_to_ensemble::IV: length n_nodes, 1-based ensemble id (in ORIGINAL ensemble order) per node.

GFT (the geometry's own float type) is Float64 for CPU() architectures (matching the single-ensemble CPU path's IC = Float64 internal-computation convention exactly, so the pre-normalized weights need no further narrowing or widening round-trip before being fed to the shared _nai2_bulk_optics core) and float(promote_type(eltype(radii), eltype(weights))) for GPU architectures (matching the single-ensemble GPU path's kernel type FT).

source
vSmartMOM.Scattering.compute_aerosol_optical_properties_nodes_batched Function
julia
compute_aerosol_optical_properties_nodes_batched(geometry::MieNodeGeometry,
    wavelength, n_real::AbstractVector, n_imag::AbstractVector,
    polarization, truncation;
    l_max = nothing,
    precision_policy = nothing,
    workspace = nothing) -> Vector{AerosolOptics}

Batched sibling of compute_aerosol_optical_properties_nodes: given a MieNodeGeometry (E ensembles, prepared once via prepare_mie_node_geometry) and ONE wavelength, returns a Vector{AerosolOptics} of length E — element e is what compute_aerosol_optical_properties_nodes(radii_e, weights_e, n_real[e], n_imag[e], wavelength, polarization, truncation; l_max, precision_policy) would return for ensemble e's own (raw) nodes, in ORIGINAL ensemble order. Output type: AerosolOptics{FT_out} where FT_out is derived from geometry's ORIGINAL radii/weights eltype (MieNodeGeometry.output_FT, captured before any internal-compute-precision forcing) promoted with n_real/n_imag/wavelength's eltypes — an all-Float32 CPU call returns Float32, matching compute_aerosol_optical_properties_nodes exactly (this is a load-bearing contract: GCHPIO's consumer code narrows to its own working precision and must be able to rely on this without a manual post-hoc narrowing pass).

Bit-compatibility contract (precise, by backend)

  • Architectures.CPU() and KA-CPU (portable KernelAbstractions backends that aren't real CUDA/Metal hardware): bit-for-bit identical (==) to the single-ensemble path, by construction (same arithmetic order, same per-ensemble normalize-once semantics — see this file's top-of-file module docstring for the full argument).

  • Real CUDA hardware: near-machine-precision, NOT bit-for-bit (different ndrange launch shapes between the batched and single-ensemble kernel invocations can produce different, but equally valid, PTX/SASS for logically identical code — see docs/src/pages/gpu.md's "Batched caller-node Mie seam" section for the measured magnitudes).

  • Exact even on real CUDA: compute_aerosol_extinction_nodes_batched's result for ensemble e equals THIS function's .k for that same ensemble — the invariant GCHPIO's hybrid exact-τ mechanism depends on (see compute_aerosol_extinction_nodes_batched).

See test/test_mie_nodes_batched.jl for the acceptance-gate tests (bitwise batched-vs-single on CPU/KA-CPU, batching-split invariance, workspace-reuse invariance, original-order reassembly, and the dedicated strict real-CUDA .k-equality regression test).

GPU: exact-n_max_global grouping

For a GPU geometry.architecture, ensembles are grouped by EXACT n_max_global at THIS wavelength (λ-dependent — group membership may differ between wavelengths; recomputed fresh every call), and one kernel-1 launch covers each group's concatenated nodes (per-ensemble refractive index via a node_to_ensemble index — see gpu_mie_kernels.jl's _batched sibling kernels).

CPU: no grouping

For Architectures.CPU(), this is a Threads.@threads loop over ensembles (grouping is a GPU-only concern — see the module docstring for why a single sequential "CPU batch" would regress the measured 16-thread performance this proposal is built around).

Arguments

  • geometry: prepared once via prepare_mie_node_geometry; reused across as many calls (different wavelengths, n_real/n_imag, etc.) as needed without re-normalizing or re-uploading anything.

  • wavelength: ONE wavelength for this call (μm, matching geometry's radii units).

  • n_real::AbstractVector, n_imag::AbstractVector: length-E per-ensemble refractive index (convention n = n_real - i·n_imag, n_imag ≥ 0 for every ensemble), one complex value per ensemble (shared by all of that ensemble's own nodes — the same "one RI per ensemble" convention as the single-ensemble API).

  • polarization, truncation: same meaning as compute_aerosol_optical_properties_nodes; truncation (δBGE then l_max cap ordering) is applied PER ENSEMBLE on the host after the raw batched compute.

Keyword arguments

  • l_max: same meaning as the single-ensemble API, applied per ensemble.

  • precision_policy: GPU-only, see MiePrecisionPolicy; nothing auto-selects via Architectures.default_mie_precision_policy(architecture, eltype(geometry.radii_host)). Ignored on the CPU path.

  • workspace: nothing (default, allocate-and-drop) or a MieBatchWorkspace to reuse device scratch across calls.

source
vSmartMOM.Scattering.compute_aerosol_extinction_nodes_batched Function
julia
compute_aerosol_extinction_nodes_batched(geometry::MieNodeGeometry,
    wavelength, n_real::AbstractVector, n_imag::AbstractVector;
    precision_policy = nothing, workspace = nothing) -> Vector

Batched, extinction-only sibling of compute_aerosol_extinction_nodes — same geometry/grouping/workspace story as compute_aerosol_optical_properties_nodes_batched, but returns just Vector{k} (length E, one number-mean bulk extinction cross-section per ensemble, ORIGINAL order, output type from geometry.output_FT promoted with n_real/n_imag/wavelength – same load-bearing Float32-output contract as the full API, see its docstring), computing none of the amplitude matrices, phase-matrix elements, or Greek coefficients — see compute_aerosol_extinction_nodes's own docstring for the motivation (GCHPIO's hybrid exact-τ leg is the intended consumer of this batched variant).

Exact-even-on-CUDA invariant

compute_aerosol_extinction_nodes_batched's result for ensemble e equals compute_aerosol_optical_properties_nodes_batched's .k for that SAME ensemble EXACTLY (==), on EVERY backend including real CUDA, for all three precision policies — see that function's "Bit-compatibility contract" section for why (shared Kernel-1 + cross-section + segmented-weighted-sum code path). GCHPIO's hybrid exact-τ mechanism rests on this; see the dedicated strict real-CUDA regression test in test/test_mie_nodes_batched.jl.

source
vSmartMOM.Scattering.MieBatchWorkspace Type
julia
MieBatchWorkspace(backend, FT::Type, RA::Type)

Reusable, GROW-ONLY GPU scratch space for the batched caller-node Mie API, keyed on the backend VALUE (objectid(backend)) plus FT (kernel precision) and RA (reduction precision, see _mie_reduction_type) — NOT polarization: polarization does not change the raw kernel-1/reduction scratch layout, so there is no separate workspace per AbstractPolarizationType.

objectid(backend) is backend-VALUE identity, not device/stream identity: KA backend structs (KernelAbstractions.CPU(), CUDA.CUDABackend(), Metal's backend) are isbits values carrying only backend-TYPE/config fields (e.g. CUDABackend's prefer_blocks/always_inline flags) — NOT which physical GPU or CUDA stream is currently active. Two separately constructed CUDA.CUDABackend() values are objectid-EQUAL regardless of CUDA.device()/CUDA.stream() (verified: same objectid across 2 distinct GPUs on a multi-GPU host). This check therefore catches genuinely different backend TYPES/configs (a KA-CPU workspace handed a CUDA launch, or vice versa) but does NOT catch a workspace built while one CUDA device/stream was active being reused while a DIFFERENT device/stream is active — that is a real gap, not yet closed: true device/stream-aware keying (e.g. incorporating CUDA.device()/CUDA.stream() into the key) is future work. Single-device, single-stream usage (the common case) is unaffected.

Capacities are tracked independently PER SCRATCH ARRAY (grouped-node count, coefficient count, node-angle scratch, Greek-table size): a request larger than an array's current capacity reallocates JUST that array (grow-only — never shrinks, so reuse across calls with varying group sizes never thrashes). nothing (the default workspace keyword everywhere in this module) means allocate-and-drop, i.e. the pre-existing zero-workspace behavior — passing an explicit MieBatchWorkspace is a pure optimization, never a correctness requirement (see the workspace-reuse-invariance test in test/test_mie_nodes_batched.jl: warm == cold, and reuse across two different wavelengths — where group membership genuinely changes — still reproduces the fresh/cold result bit-for-bit).

NOT thread-safe. Single-owner only: one workspace per concurrent GPU stream / CPU task. The FIRST use from a given task stamps owner_task; a LATER use from a DIFFERENT task throws an ArgumentError immediately rather than silently corrupting scratch shared across two concurrent callers. This check is UNCONDITIONAL (every _ws_array! call, not gated behind bounds-checking) — it is one field compare (plus, on first use, one field write), cheap enough that making it a debug-only check would only trade a negligible cost saving for the ability to silently disable safety via @inbounds. Every use also verifies objectid(backend) === backend_id (the backend VALUE the workspace was constructed for): this catches a workspace built for one backend TYPE/config being handed to a launch against a DIFFERENT one (e.g. KA-CPU vs CUDA vs Metal, or a differently-configured CUDABackend). It does NOT distinguish CUDA devices or streams (see the caveat above) — single-device, single-stream usage only; true device/stream keying is future work.

source

Phase Matrix Reconstruction ​

vSmartMOM.Scattering.reconstruct_phase Function
julia
reconstruct_phase(greek_coefs, μ; returnLeg = false)

Reconstruct angle-space scattering-matrix elements from Greek coefficients.

Reference: Sanghavi (2014), Fourier/Greek framework around Eq. (17).

f₁₁ represents the scalar phase function and is normalized as:

Using Legendre basis matrices computed at μ, this function evaluates:

Arguments

  • greek_coefs: GreekCoefs coefficients.

  • μ: cosine of scattering angles where the matrix should be reconstructed.

  • returnLeg: if true, also return (P, P²).

Returns

  • ScatteringMatrix when returnLeg=false.

  • (ScatteringMatrix, P, P²) when returnLeg=true.

source
julia
$(FUNCTIONNAME)(greek_coefs, μ; returnLeg = false)

Returns the reconstructed elements of the 4x4 scattering matrix at positions f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄ from the greek coefficients

f₁₁ represents the phase function p for the Intensity (first Stokes Vector element) and is normalized as follows:

  • greek_coefs greek coefficients (Domke Type)

  • returnLeg if false (default), just return f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄, if true,

  • return f₁₁, f₁₂, f₂₂, f₃₃, f₃₄, f₄₄, P, P² (i.e. also the two legendre polynomials as matrices)

source

Fourier Z-Moment Assembly ​

Per-Fourier-moment phase operators on the quadrature grid: compute_Z_moments builds the diffuse Z⁺⁺/Z⁻⁺ blocks from Greek coefficients (optionally via the per-solve ZMomentTables precompute), and compute_Z_source_moments builds the rectangular direct-solar Z₀⁺/Z₀⁻ columns used by the external-solar SFI mode.

vSmartMOM.Scattering.compute_Z_moments Function
julia
compute_Z_moments(mod::AbstractPolarizationType, μ, α, β, γ, δ, ϵ, ζ, m::Int)

Compute moments of the phase matrix

source
julia
compute_Z_moments(mod::AbstractPolarizationType, μ, α, β, γ, δ, ϵ, ζ, m::Int)

Compute moments of the phase matrix.

Keywords tables/Π_pair (see ZMomentTables / make_Π_lists) supply the scatterer-independent precompute; results are identical (the same sums in the same order), only the redundant recomputation is skipped.

source
vSmartMOM.Scattering.compute_Z_source_moments Function

Forward and tangent solar-direction phase blocks without a square μ₀ node.

source
julia
compute_Z_source_moments(mod, μ_out, μ₀, greek_coefs, m; arr_type=Array)

Evaluate the two Fourier phase blocks coupling one collimated incident direction μ₀ to the diffuse output directions μ_out, without constructing an augmented square angular matrix. The angular shape is length(μ_out) × 1; after Stokes expansion the returned arrays have shape (length(μ_out)*nStokes, nStokes[, nSpec]).

Keeping the complete incident Stokes block preserves support for polarized collimated sources. For the usual unpolarized solar beam only column one is used, so its contraction is an NquadN × 1 source vector.

source
vSmartMOM.Scattering.ZMomentTables Type
julia
ZMomentTables(μ, l_max)

Scatterer-independent precompute for compute_Z_moments: the generalised spherical functions (Sanghavi 2014, eq. 15) evaluated on the stream cosines up to a global l_max (the longest Greek-coefficient expansion among all scatterers in the run). Build once per solve and pass via the tables keyword.

source
vSmartMOM.Scattering.make_Π_lists Function
julia
make_Π_lists(mod, tables, m)

Assemble the Π matrices (Sanghavi 2014, eq. 15) for Fourier moment m (0-indexed) from precomputed tables, for every l = m+1 … l_max (1-indexed). The result is shared by ALL scatterers at this m — Π does not involve the Greek coefficients. Returns (Π, Π⁻) indexed as Π[l], or nothing when m ≥ l_max (all Z moments vanish there).

source

Convenience Computations ​

Useful for cross-sections or scalar phase-function outputs without running the full RT-facing workflow.

vSmartMOM.Scattering.phase_function Function
julia
phase_function(aerosol::Aerosol, λ, r_max, nquad_radius; r_min=DEFAULT_R_MIN)
    -> (μ, w_μ, p11, C_ext, C_sca, g)

Compute scalar phase-function output for a size-distributed aerosol, with the size integral on nquad_radius log-spaced Gauss-Legendre nodes in the fixed interval [r_min, r_max] (the grid of radius_quadrature).

The returned p11(μ) is normalized by the bulk scattering cross-section and satisfies:

Asymmetry factor:

Reference: Sanghavi (2014), Mie-amplitude/cross-section setup from Eq. (1), with NAI-style angular integration discussed in Secs. 3-4.

source
julia
phase_function(r::FT, λ::FT, nᵣ::FT, nᵢ::FT) where {FT<:AbstractFloat}
    -> (μ, w_μ, p11, C_ext, C_sca, g)

Monodisperse version of phase_function, for a single radius r. Returns the same tuple as the size-distribution overload.

source
vSmartMOM.Scattering.compute_aerosol_XS Function
julia
compute_aerosol_XS(aerosol::Aerosol, λ, r_max, nquad_radius; r_min=DEFAULT_R_MIN)
    -> (XS_ext, XS_sca, Cext_eff, Csca_eff)

Compute bulk extinction/scattering cross-sections for a size-distributed aerosol without reconstructing phase matrices. The size integral uses the fixed radius interval [r_min, r_max] (the grid of radius_quadrature).

Returns

  • XS_ext: bulk extinction cross-section

  • XS_sca: bulk scattering cross-section

  • Cext_eff: area-normalized extinction

  • Csca_eff: area-normalized scattering

with

Reference: Sanghavi (2014), Eq. (1) for particle-level C_ext/C_sca.

source
vSmartMOM.Scattering.compute_ref_aerosol_extinction Function
julia
compute_ref_aerosol_extinction(model::MieModel{<:NAI2,FT}, FT_out::Type) -> Real

Compute only the bulk extinction cross-section for the aerosol model:

where is evaluated from the Mie series and is represented with quadrature over [r_min, r_max] (see radius_quadrature).

The second argument FT_out sets the return type. It is supplied by the caller and must be passed explicitly — there is no default. The internal cross-section computation always runs in Float64 for plain floats (same numerical stability rationale as compute_aerosol_optical_properties).

This helper avoids phase-function and Greek-coefficient reconstruction.

Reference: Sanghavi (2014), Eq. (1) for C_ext(r).

source

Analytic Phase Functions ​

Lightweight scattering sources for idealized tests, sensitivity experiments, and StandaloneSS validation. Not a separate RT implementation: each analytic source is projected into GreekCoefs and then consumed by the same CoreRT optical-property machinery as Mie-derived aerosols. phase_matrix_first_column evaluates the single-scatter phase-matrix column needed for an unpolarized direct solar beam at one exact sun-view geometry — the shared hook used by StandaloneSS vector paths for the sun-to-atmosphere-to-sensor term.

vSmartMOM.Scattering.AbstractAnalyticPhaseFunction Type
julia
AbstractAnalyticPhaseFunction

Analytic phase/scattering matrix source that can be converted to Greek coefficients and then used by the standard MOM optical-property path.

source
vSmartMOM.Scattering.HenyeyGreensteinPhaseFunction Type
julia
HenyeyGreensteinPhaseFunction(g)

Scalar Henyey-Greenstein phase function, (1 - g^2) / (1 + g^2 - 2g cosΘ)^(3/2), normalized so its sphere average is (1 - g^2) / (1 + g^2 - 2g cosΘ)^(3/2), normalized so its sphere average is one.

source
vSmartMOM.Scattering.SyntheticPolarizedHenyeyGreensteinPhaseFunction Type
julia
SyntheticPolarizedHenyeyGreensteinPhaseFunction(; g, polarization_fraction)

Diagnostic polarizing Henyey-Greenstein-like scattering matrix. The f11 element is standard Henyey-Greenstein; f12/f11 follows the bounded toy law polarization_fraction * (1 - cosΘ^2) / (1 + cosΘ^2). This is intended for tests and sensitivity experiments, not as a Mie substitute.

source
vSmartMOM.Scattering.greek_coefficients Function
julia
greek_coefficients(phase; l_max=64, nquad=max(2l_max+1, 64))

Project an analytic phase/scattering matrix into Greek coefficients so it can enter the same CoreRT MOM path as Mie-derived aerosol optics.

source
vSmartMOM.Scattering.analytic_aerosol_optics Function
julia
analytic_aerosol_optics(phase; single_scattering_albedo=1, extinction_cross_section=1, l_max=64)

Create an AerosolOptics object from an analytic phase/scattering matrix. The returned object is directly consumable by the existing CoreRT aerosol mixing and Fourier-moment code.

source
vSmartMOM.Scattering.phase_matrix_first_column Function
julia
phase_matrix_first_column(greek, μ₀, μv, Δϕ, Val(N))

Evaluate the first column of the scattering phase matrix for an unpolarized incident beam at one plane-parallel sun-view geometry. μ₀ and μv are positive cosines and Δϕ is the vSmartMOM relative azimuth in radians.

source

Pure-Gas Molecular Rayleigh ​

Method-explicit pure-gas molecular Rayleigh properties (no frequency redistribution). See Rayleigh and the Cabannes choice for the physics and provenance.

vSmartMOM.Scattering.AbstractMolecularRayleighMethod Type
julia
AbstractMolecularRayleighMethod

Dispatch category for independently implemented molecular-Rayleigh cross-section methods with no explicit frequency redistribution. The method is an explicit argument so a refractive-index calculation can never be confused with a polarizability- tensor calculation in provenance records or calling code.

Concrete methods:

This abstraction is part of Scattering, not Inelastic. Only UsingRefInd is independent of the Raman-module (InelasticScattering) constants. UsingPolTen deliberately uses them: it builds fresh Raman-module molecular constants on each call and never mutates shared Raman state.

source
vSmartMOM.Scattering.UsingRefInd Type
julia
UsingRefInd()
RefractiveIndexRayleigh()

Select the conventional refractive-index plus effective-King-factor Rayleigh formulation without constructing frequency-shifted Raman branches. The Habitable Worlds Observatory (HWO) exo-Earth studies use this explicit method for the Ar and CO₂ cross sections. RefractiveIndexRayleigh is a compatibility alias (RefractiveIndexRayleigh === UsingRefInd) kept so existing calling code, e.g. ExoOptics, continues to work unchanged. The polarizability-tensor implementation uses the distinct dispatch type UsingPolTen; it never silently replaces this method.

source
vSmartMOM.Scattering.UsingPolTen Type
julia
UsingPolTen(anisotropy::AbstractPolarizabilityDispersion = DispersiveAnisotropy())

Select the polarizability-tensor Rayleigh formulation, built from the mean polarizability and the anisotropy of N₂, O₂ or H₂ (Unicode aliases such as :N₂ are accepted). anisotropy selects their wavelength dependence:

  • DispersiveAnisotropy (default): and from the UsingRefInd refractivity and King factor of N₂ and O₂, with the Buldakov temperature factor on . This is the formulation of the RT model's automatic (depol: -1) Rayleigh depolarization.

  • ConstantAnisotropy: the Raman-module tensor (Buldakov , wavelength-independent ), which the Raman and Cabannes cross sections use. UsingPolTen(ConstantAnisotropy()) reproduces the UsingPolTen() results of the v2.2 pre-releases.

Both go through the Raman-module routines getMolecularConstants, compute_effective_coefficents! and compute_σ_Rayl_coeff! and never mutate shared state. See molecular_rayleigh_properties for the formulation.

source
vSmartMOM.Scattering.AbstractPolarizabilityDispersion Type
julia
AbstractPolarizabilityDispersion

Wavelength dependence of the N₂ and O₂ polarizability-tensor invariants, the mean polarizability and the anisotropy . Two formulations are kept and selected by dispatch:

  • DispersiveAnisotropy (default for UsingPolTen and for the automatic depol: -1 Rayleigh depolarization of the RT model): from the UsingRefInd refractivity through the Lorentz–Lorenz relation and from its King factor;

  • ConstantAnisotropy (default for the Raman/Cabannes cross sections): the published Raman-module tensor, Buldakov-type and a wavelength-independent .

See docs/dev_notes/molecular_rayleigh_dispersion.md for the sources and the size of the difference.

source
vSmartMOM.Scattering.DispersiveAnisotropy Type
julia
DispersiveAnisotropy()

Polarizability-tensor invariants of N₂ and O₂ that follow the refractive-index data of UsingRefInd:

where and are the refractive index and number density of the dispersion fit used at (reference temperature ), is the UsingRefInd King factor (Bates 1984 for N₂; He et al. 2021, Table 1, for O₂), and     is the Buldakov et al. temperature factor of the Raman-module mean polarizability. The anisotropy is that implied by the measured depolarization ratios behind the King-factor fits (Bates 1984, as in Bodhaine et al. 1999 and Tomasi et al. 2005), so the depolarization    equals the UsingRefInd value    at  . Validity follows the UsingRefInd fits: N₂ 254–2059 nm, O₂ 288–1542 nm measured (longer wavelengths are dispersion extrapolations, which are smooth because  ); the King factors are fitted to visible depolarization data and extrapolated into the UV and NIR.

H₂ has no refractive-index fit in UsingRefInd; for H₂ this formulation uses the Raman-module tensor of ConstantAnisotropy unchanged.

source
vSmartMOM.Scattering.ConstantAnisotropy Type
julia
ConstantAnisotropy()

The polarizability tensor of the Raman module (getMolecularConstants in src/Inelastic/src/molecular_constructors.jl): mean polarizability with the Buldakov et al. single-pole dispersion and temperature dependence,     , and the anisotropy  , independent of wavelength and temperature (0.71×10⁻²⁴ cm³ for N₂, 1.080×10⁻²⁴ cm³ for O₂). These constants produced the published Raman tables (Sanghavi 2022, JQSRT 291, 108328; the Sanghavi 2026 Table 1 checked in test/test_ray_cab_raman.jl), and they remain the default of the Raman and Cabannes cross sections. For N₂ and O₂ the constant anisotropy equals the King-factor anisotropy of DispersiveAnisotropy near 488 nm (within 0.03 % for N₂ and 0.32 % for O₂); it is 2 % (N₂) and 4 % (O₂) low at 400 nm and 4 % (N₂) and 8 % (O₂) high at 2060 nm.

source
vSmartMOM.Scattering.molecular_rayleigh_properties Function
julia
molecular_rayleigh_properties(::UsingRefInd, species, wavelength_nm)
molecular_rayleigh_properties(species, wavelength_nm)

Return no-frequency-redistribution Rayleigh properties for a pure gas at vacuum wavelength wavelength_nm (nm). Supported species are :He, :Ar, :N2, :O2, and :CO2 (Unicode aliases such as :N₂ are accepted).

The two-argument form is equivalent to passing UsingRefInd(). The machine-readable provenance record (molecular_rayleigh_provenance) currently covers Ar and CO₂ only; the He, N₂ and O₂ coefficient sets have no such record yet (their sources are listed below).

The returned named tuple contains

  • refractive_index: refractive index at the species reference number density: 2.546899×10¹⁹ cm⁻³ (288.15 K, 1013.25 hPa) for He, Ar, N₂ and CO₂ and 2.68678×10¹⁹ cm⁻³ (273.15 K, 1013.25 hPa) for O₂. A fit given at another density is converted through the Lorentz–Lorenz relation ((n²-1)/(n²+2) proportional to N). It is stored as itself, so a Float32 result resolves   only to about eps(1f0)/(n-1) relative (eps(1f0) ≈ 1.2e-7): a few times 1e-4 for N₂, O₂, Ar and CO₂, and a few times 1e-3 for He. The cross section is computed from the refractivity directly and is not affected;

  • king_factor: molecular-anisotropy correction ;

  • depolarization: depolarization ratio in the convention consumed by get_greek_rayleigh, related by    ;

  • cross_section_cm2: Rayleigh scattering cross section in cm² molecule⁻¹.

The individual dispersion and King-factor sources are returned by molecular_rayleigh_provenance for Ar and CO₂. Wilmouth and Sayres (2019) provide a direct 333–363 nm check of both gases. The cross-section expression is

where is in cm⁻¹ and N in molecule cm⁻³, the number density at which the dispersion fit gives . Every fit is used with the ideal-gas number density of its stated reference state, as in Bates (1984), Bodhaine et al. (1999), Sneep and Ubachs (2005) and He et al. (2021); gas compressibility (5–10×10⁻⁴ in   per molecule at 1 atm, 0–20 °C) is neglected. The implementation can evaluate any positive wavelength, but callers must use the distinct dispersion, King-factor, and validation evidence returned by molecular_rayleigh_provenance and mark extrapolation. The returned conventional total molecular-Rayleigh cross section uses an effective King factor; it is not a separately computed Cabannes-only cross section. No species-specific rotational or vibrational Raman redistribution is built.

N₂ and O₂ dispersion and King factors

Each gas joins two dispersion fits at the wavenumber where they agree, so   and are continuous (there is a kink in their slope):

  • N₂,   22 919.551 cm⁻¹ (  436.309 nm): Bates (1984),         at 288.15 K (He et al. 2021, Table 1; fitted over 254–468 nm). Longer wavelengths: Peck and Khanna (1966),         at 0 °C and 760 Torr (measured 467.9–2058.7 nm; checked against the refractiveindex.info N2/nk/Peck-0C.yml record). The two fits differ by 1.4×10⁻⁴ in   at 468 nm; joining them there would leave a 0.03 % jump in .

  • O₂,   21 428.071 cm⁻¹ (  466.678 nm): Bates (1984),         at 273.15 K (He et al. 2021, Table 1; fitted over 288–546 nm). Longer wavelengths: Egan (2026), Eq. (8),       ,   in µm⁻¹, at 20 °C and 100 kPa, constrained by refractivity at 0.633 and 1.542 µm (uncertainty 7×10⁻⁶ in  ), visible data and a 0.254/0.436 µm ratio. At 546 nm the two fits differ by 0.10 % in   (0.21 % in ); at 1.6 µm the extrapolated Bates fit is 0.43 % high in  .

  • King factors: N₂     (Bates 1984); O₂       (He et al. 2021, Table 1). The O₂ constant 1.09 is He et al.'s; Bates (1984) and Bodhaine et al. (1999) use 1.096, which gives a 0.54 % larger O₂ cross section (0.54 % at 400 nm, 0.55 % at 1.6 µm). Both King factors are fitted to visible depolarization data and extrapolated elsewhere.

References

source
julia
molecular_rayleigh_properties(method::UsingPolTen, species, wavelength_nm;
                              temperature_K = 288.15,
                              number_density_cm3 = 2.546899e19)

Polarizability-tensor variant. Species: :N2, :O2, :H2. For any other species an ArgumentError is thrown that names the species and suggests UsingRefInd().

Tensor invariants, selected by method.anisotropy:

  • DispersiveAnisotropy (UsingPolTen()), N₂ and O₂:    and   from the UsingRefInd Lorentz–Lorenz factor , density and King factor at   cm⁻¹, with the Buldakov temperature factor    . At   (288.15 K for N₂ below 436.3 nm, 273.15 K above; 273.15 K for O₂ below 466.7 nm, 293.15 K above) the cross section and depolarization equal the UsingRefInd values; at other temperatures only changes (by 1.8×10⁻⁶ K⁻¹ for N₂ and about 2.6×10⁻⁶ K⁻¹ for O₂ near 288 K). H₂ uses the constant-anisotropy tensor below.

  • ConstantAnisotropy (Raman-module constants of getMolecularConstants(::N₂/::O₂/::H₂, vmr)): mean polarizability with the Buldakov et al. single-pole dispersion and temperature dependence,     , and the anisotropy  , independent of wavelength and temperature.

Derived quantities (no further inputs):

  • depolarization    (the module's rho_depol_Rayl, the same convention as the RI method and get_greek_rayleigh);

  • king_factor      ;

  • cross_section_cm2  , i.e. the module's σ_Rayl_coeff · ν⁴ (cm² molecule⁻¹, in cm³, in cm⁻¹). This is the total (Cabannes plus rotational Raman) molecular Rayleigh cross section, as for the RI method;

  • refractive_index: the Lorentz–Lorenz equivalent     at number_density_cm3 (default 2.546899e19 molecule cm⁻³, 288.15 K and 1013.25 hPa). Pass number_density_cm3 = nothing to return refractive_index = nothing. The cross section, depolarization and King factor are per molecule and do not depend on number_density_cm3.

The default temperature_K = 288.15 K (kept by decision 13d, 2026-10-19) is the temperature of the default number_density_cm3, the standard state of the He, Ar, N₂ (λ < 436 nm) and CO₂ RI fits and of the Bates (1984) and Bodhaine et al. (1999) cross sections. It therefore compares the tensor with those data at their own state; the O₂ RI refractive index refers to 273.15 K and 2.68678e19 molecule cm⁻³, so compare O₂ refractive indices at that density. The temperature enters only through the small Buldakov factor on : between 200 and 300 K the cross sections change by less than 4×10⁻⁴ (N₂, O₂) and 2×10⁻³ (H₂). The RT model's own molecular-depolarization path does not use this default: it evaluates the tensor at the column-weighted mean temperature of the profile and the band-centre wavenumber. Calling this function does not change any model default.

The internal evaluation is always in Float64 (  underflows Float32); results are converted to the floating-point type of float(wavelength_nm). The method is therefore not ForwardDiff-aware (dual numbers are not propagated) and does not gain precision from BigFloat inputs. As for the RI method, a Float32 refractive_index resolves   only to about eps(1f0)/(n-1) relative (up to ~1e-3, for H₂). The Raman-module polarizability constants and their temperature coefficients come from Buldakov et al.; the O₂ and H₂ anisotropies from Asawaroengchai & Rosenblatt (1980), per the source comments. The H₂ tensor is checked against the Peck & Huang (1977) refractivity and the Raj, Hamaguchi & Witek (2018) anisotropy in test/test_molecular_rayleigh.jl.

source
vSmartMOM.Scattering.molecular_rayleigh_cross_section_ratio Function
julia
molecular_rayleigh_cross_section_ratio(method::AbstractMolecularRayleighMethod,
                                        species, wavelength_nm,
                                        reference_nm)
molecular_rayleigh_cross_section_ratio(::UsingPolTen, species,
                                        wavelength_nm, reference_nm;
                                        temperature_K = 288.15)
molecular_rayleigh_cross_section_ratio(species, wavelength_nm, reference_nm)

Return for a pure gas, with = wavelength_nm and = reference_nm. This ratio is the appropriate scaling when a comparison fixes the Rayleigh optical depth at one reference wavelength and varies molecular identity. For UsingPolTen both cross sections are evaluated at the same temperature_K; the number density cancels. The three-argument form uses UsingRefInd.

source
vSmartMOM.Scattering.molecular_rayleigh_provenance Function
julia
molecular_rayleigh_provenance(::UsingRefInd, species)

Return the coefficient provenance and separate evidence intervals for the refractive-index/King-factor coefficients of Ar and CO₂. The primary RI measurement range, adopted dispersion-formula envelope, King-factor support, and direct scattering checks are reported as distinct fields, so callers can tell measured or directly validated coverage from formula extrapolation.

Ar uses the Peck–Fisher dispersion relation and unity King factor. CO₂ uses the Bideau-Mehu dispersion relation (scaled to the stated reference density) and the Sneep–Ubachs King-factor fit to Alms–Burnham–Flygare measurements. Wilmouth and Sayres (2019) directly measured both gases at 333–363 nm; Thalman et al. (2014) checked Ar in selected regions between 350 and 660 nm; and He et al. (2021) provide a continuous 307–725 nm direct check for CO₂. This provenance API is currently limited to Ar and CO₂.

source
vSmartMOM.Scattering.get_greek_rayleigh Function
julia
get_greek_rayleigh(depol)

Returns the greek coefficients (as GreekCoefs) of the Rayleigh phase function given depolarization value.

  • depol Depolarization (best use 0 as default )
source

Wigner Utilities (PCW) ​

vSmartMOM.Scattering.compute_wigner_values Function
julia
compute_wigner_values(m_max::Integer, n_max::Integer, l_max::Integer)

Compute the Wigner 3j values for all (m, n, l) combinations up to m/n/l max, for (m1, m2, m3) = (-1, 1, 0) (wigner_A) and (-1, -1, 2) (wigner_B)

source
julia
compute_wigner_values(N_max)

Shorthand for compute_wigner_values(2N_max + 1, N_max + 1, 2N_max + 1)

source
vSmartMOM.Scattering.save_wigner_values Function
julia
save_wigner_values(filepath, wigner_A, wigner_B)

Save the Wigner A and Wigner B matrices at the given filepath

source
vSmartMOM.Scattering.load_wigner_values Function
julia
load_wigner_values(filename)

Load the Wigner A and Wigner B matrices from the given filepath

source

Core Abstract Interfaces ​

vSmartMOM.Scattering.AbstractAerosolType Type
julia
type AbstractAerosolType

Abstract aerosol type

source
vSmartMOM.Scattering.AbstractFourierDecompositionType Type
julia
type AbstractFourierDecompositionType

Abstract aerosol Fourier Decomposition computation type (NAI2 and PCW)

source
vSmartMOM.Scattering.AbstractPolarizationType Type
julia
type AbstractPolarizationType

Abstract Polarization type

source
vSmartMOM.Scattering.AbstractTruncationType Type
julia
AbstractTruncationType

Abstract supertype for phase-function truncation methods. All concrete methods are dispatched through truncate_phase and supply l_max(t) (the per-band Legendre cutoff that the RT pipeline allocates for). Subtypes:

  • NoTruncation — identity. Use when the phase function has no sharp forward peak (canopy, isotropic scattering, smooth Rayleigh).

  • δBGE — δ-BGE-fit (Sanghavi & Stephens 2015, JQSRT 159 §3); recommended for hyperspectral retrievals.

The atmospheric Δ_angle (forward exclusion half-angle) lives inside the truncation type that needs it, not as a free parameter on CoreRT.vSmartMOM_Parameters — different methods have different hyper-parameters and NoTruncation has none.

source

Output Types ​

vSmartMOM.Scattering.GreekCoefs Type
julia
GreekCoefs{FT<:Real}

Expansion coefficients of the 4×4 scattering (phase) matrix in generalised spherical functions (the "Greek" coefficients). Six independent coefficient vectors (α, β, γ, δ, ϵ, ζ) fully describe the azimuthal Fourier decomposition of the scattering matrix B for a given particle or mixture. See Eq. 16 in Sanghavi (2014) for the mapping to B elements.

For scalar (intensity-only) RT, only β (the phase-function expansion) is used.

Fields

  • α: Greek matrix coefficient α, is in B[2,2]

  • β: Greek matrix coefficient β, is in B[1,1] (only important one for scalar!)

  • γ: Greek matrix coefficient γ, is in B[2,1],B[1,2]

  • δ: Greek matrix coefficient δ, is in B[4,4]

  • ϵ: Greek matrix coefficient ϵ, is in B[3,4] and - in B[4,3]

  • ζ: Greek matrix coefficient ζ, is in B[3,3]

source
vSmartMOM.Scattering.Aerosol Type
julia
Aerosol

Aerosol microphysical properties: particle size distribution and complex refractive index. Used as input to MieModel for computing single-scattering optical properties via Lorenz-Mie theory.

The refractive index convention is n = nᵣ - i·nᵢ, where positive nᵢ indicates absorption.

The struct is parameterized by the refractive-index element type FT. FT is typically a plain AbstractFloat (Float64/Float32) but may also be a ForwardDiff.Dual so the autodiff path can track derivatives with respect to nᵣ/nᵢ; for that reason FT is intentionally not constrained to <:AbstractFloat. The outer convenience constructor Aerosol(dist, nᵣ, nᵢ) promotes nᵣ and nᵢ to a common type, so existing call sites continue to work unchanged.

Fields

  • size_distribution::ContinuousUnivariateDistribution: Particle radius distribution (e.g., LogNormal). Units: μm.

  • nᵣ::FT: Real part of the refractive index (relative to air).

  • nᵢ::FT: Imaginary part of the refractive index (absorption).

Example

julia
using Distributions
aer = Aerosol(LogNormal(log(0.3), 0.4), 1.3, 0.01)
source
vSmartMOM.Scattering.AerosolOptics Type
julia
AerosolOptics{FT<:Real}

Computed aerosol single-scattering optical properties for one aerosol type at one (or more) wavelengths. Produced by integrating the Mie solution over the particle size distribution.

Fields

  • greek_coefs: Greek matrix

  • ω̃: Single Scattering Albedo

  • k: Extinction cross-section

  • fᵗ: Truncation factor

  • derivs: Derivatives

  • phase_ν: Wavenumber nodes (cm⁻¹) at which the post-truncation phase coefficients were evaluated

  • phase_greek: Post-truncation Greek coefficients at phase_ν; retained so Z is evaluated only at these nodes

  • phase_greek_raw: UNTRUNCATED Greek coefficients at phase_ν (or at the single Mie wavelength when phase_ν === nothing). Retained when δ-BGE truncation actually removed a forward peak, so the TMS single-scattering correction can evaluate the exact phase function at the view geometries without a rebuild. nothing when no truncation occurred (exact == truncated).

source
vSmartMOM.Scattering.MieModel Type
julia
MieModel{FDT<:AbstractFourierDecompositionType, FT, ARCH}

Configuration for a Lorenz–Mie scattering computation. Specifies the aerosol (size distribution + refractive index), wavelength, polarization type, truncation strategy, integration parameters, and the compute architecture. The size integral uses nquad_radius Gauss–Legendre nodes in ln r on the fixed interval [r_min, r_max] (radius_quadrature). The computation_type selects between NAI-2 (Siewert) and PCW (Domke) Fourier decomposition algorithms.

Pre-computed Wigner symbol tables (wigner_A, wigner_B) can be supplied for PCW; they default to trivial 1×1×1 placeholders (of element type FT) when unused (NAI2 path).

FT semantics (three orthogonal precision axes)

FT is the output float type of the returned Greek coefficients and optical scalars. It is carried by λ, r_max, r_min, and aerosol (which must all share the same FT) so that the entire MieModel is FT-consistent. This is independent of two other precision choices:

  1. FT (output type) — what the user consumes downstream in the RT pipeline.

  2. The internal Dₙ continued-fraction recursion: on the CPU path the recursion is always promoted to Float64 for numerical stability regardless of FT (see _mie_dn_recursion! in mie_helper_functions.jl); Dual inputs keep native arithmetic for AD compatibility.

  3. precision_policy (GPU only) — NativeFloat64 (default on CUDA) runs the GPU Dₙ recursion in hardware FP64; DSEmulated uses Float32 double-single pairs (for Metal/L40S). This axis is inert on the CPU path.

Architecture dispatch

architecture::ARCH is either Architectures.CPU() (default) or Architectures.GPU(). compute_aerosol_optical_properties(mie_model) dispatches on it: NAI2+CPU → analytic CPU path; NAI2+GPU → the KernelAbstractions GPU pipeline; PCW+GPU → CPU fallback (no GPU PCW kernel).

Fields

  • computation_type

  • aerosol: Aerosol microphysics; element type must match the model's FT

  • λ: Wavelength [μm]; must have the same float type as r_max

  • polarization_type

  • truncation_type

  • r_max: Maximum radius [μm]

  • r_min: Minimum radius [μm]; fixed, independent of the size distribution (default DEFAULT_R_MIN)

  • nquad_radius: Number of quadrature points for integration over size distribution

  • wigner_A: Precomputed Wigner ν=0 table (PCW only; trivial 1×1×1 placeholder for NAI2)

  • wigner_B: Precomputed Wigner ν=2 table (PCW only; trivial 1×1×1 placeholder for NAI2)

  • architecture: Compute architecture (CPU() or GPU()); selects CPU vs GPU Mie path

  • precision_policy: GPU precision policy (NativeFloat64/DSEmulated); nothing = auto-select on GPU, ignored on CPU

source

Linearized Output Types ​

vSmartMOM.Scattering.linGreekCoefs Type
julia
struct GreekCoefs{FT}

A struct which holds all Greek coefficient lists (over l) in one object. See eq 16 in Sanghavi 2014 for details.

Fields

  • α̇: Greek matrix coefficient α, is in B[2,2]

  • β̇: Greek matrix coefficient β, is in B[1,1] (only important one for scalar!)

  • γ̇: Greek matrix coefficient γ, is in B[2,1],B[1,2]

  • δ̇: Greek matrix coefficient δ, is in B[4,4]

  • ϵ̇: Greek matrix coefficient ϵ, is in B[3,4] and - in B[4,3]

  • ζ̇: Greek matrix coefficient ζ, is in B[3,3]

source
vSmartMOM.Scattering.linAerosolOptics Type
julia
struct AerosolOptics

A struct which holds all computed aerosol optics

Fields

  • lin_greek_coefs: Greek matrix

  • ω̃̇: Single Scattering Albedo

  • k̇: Extinction cross-section

  • ḟᵗ: Truncation factor

  • phase_ν: Wavenumber nodes corresponding to phase_lin_greek

  • phase_lin_greek: Post-truncation Greek-coefficient tangents retained at phase nodes

source

Fourier Decomposition Modes ​

vSmartMOM.Scattering.NAI2 Type
julia
type NAI2

Perform Siewart's numerical integration method, NAI-2, to compute aerosol phase function decomposition. See: http://adsabs.harvard.edu/full/1982A%26A...109..195S

source
vSmartMOM.Scattering.PCW Type
julia
type PCW

Perform Domke's Precomputed Wigner Symbols method, PCW, to compute aerosol phase function decomposition. See: http://adsabs.harvard.edu/full/1984A%26A...131..237D

source

Polarization Types ​

vSmartMOM.Scattering.Stokes_IQUV Type
julia
struct Stokes_IQUV{FT<:AbstractFloat}

A struct which defines full Stokes Vector ([I,Q,U,V]) RT code

Fields

  • n: Number of Stokes components (int)

  • D: Vector of length n for ... (see eq in Sanghavi )

  • I₀: Incoming Stokes vector for scalar only

source
vSmartMOM.Scattering.Stokes_IQU Type
julia
struct Stokes_IQU{FT<:AbstractFloat}

A struct which defines Stokes Vector ([I,Q,U]) RT code

Fields

  • n: Number of Stokes components (int)

  • D: Vector of length n for ... (see eq in Sanghavi )

  • I₀: Incoming Stokes vector for scalar only

source
vSmartMOM.Scattering.Stokes_IQ Type
julia
struct Stokes_IQ{FT<:AbstractFloat}

A struct which defines Stokes Vector ([I,Q]) RT code.

The CoreRT solvers reject Stokes_IQ: neglecting U in multiple scattering can cause intensity errors of several percent (in the tested aerosol scene, largest in the principal plane). The type remains for the Scattering-level functions and StandaloneSS ExactSSConfig single-scattering runs. Use Stokes_IQU for RT.

Fields

  • n: Number of Stokes components (int)

  • D: Vector of length n for ... (see eq in Sanghavi )

  • I₀: Incoming Stokes vector for scalar only

source
vSmartMOM.Scattering.Stokes_I Type
julia
struct Stokes_I{FT<:AbstractFloat}

A struct which define scalar I only RT code

Fields

  • n: Number of Stokes components (int)

  • D: Vector of length n for ... (see eq in Sanghavi )

  • I₀: Incoming Stokes vector for scalar only

source

Phase-function Truncation ​

vSmartMOM.Scattering.AutoTruncation Type
julia
AutoTruncation()

Phase D — deferred-decision marker for truncation: auto in YAML (mirrors VLIDORT's DO_DELTAM_SCALING philosophy). At model_from_parameters time it resolves deterministically:

  • No aerosol scattering, or all aerosols' length(greek.β) - 1 fits within stream_l_cap ⇒ NoTruncation() (typical for Rayleigh-only or smooth-aerosol scenes).

  • Aerosols with phase_lmax > stream_l_cap ⇒ δBGE(stream_l_cap) (the retired Δ_angle is always 0).

The resolver emits an @info line stating which branch was taken so the user can always see what was applied. AutoTruncation is deliberately not threaded through the Mie/RT kernels — it is a build-time placeholder that gets replaced before any kernel sees it.

User-facing knobs in YAML:

ValueMeaning
truncation: autoThis deferred-decision mode (Phase D recommended)
truncation: NoTruncation() / nullExactly no transform; errors if coefs exceed cap
truncation: δBGE(N, Δ)Explicit; used for benchmarks / cross-validation (Δ is retired and forced to 0)
source
vSmartMOM.Scattering.NoTruncation Type
julia
NoTruncation(; l_max=typemax(Int))

Identity truncation — phase functions are passed through unchanged.

This is the correct choice for radiative transfer through media whose phase function has no sharp forward peak: canopy bi-Lambertian scattering (the f_tr → 0 limit of Sanghavi & Stephens 2015 Eq. 8 is exactly the identity), isotropic scattering, and smooth Rayleigh. For Mie aerosol or ice-cloud forward peaks use δBGE instead.

source
vSmartMOM.Scattering.δBGE Type
julia
δBGE{FT}(l_max)

δ-BGE-fit truncation, vector form (Sanghavi & Stephens 2015, JQSRT 159, §3 — extension of Hu et al. 2000 to vector RT). Fits truncated Legendre coefficients over the full angular domain and renormalises by the retained scattering fraction   . (The former forward exclusion cone Δ_angle is retired and always 0; see the field docs.)

Recommended over plain δ-m for hyperspectral retrievals because δ-m has known DSE and PTE errors near exact backscatter (Sanghavi & Stephens 2015 §2.4).

Fields

  • l_max: Truncation length for Legendre terms

  • Δ_angle: RETIRED — always 0 (design decision 2026-08-23). The δBGE least-squares fit always uses the full angular domain; the former forward-peak exclusion cone was honoured only by the forward truncate_phase, never by the linearized overload, so any non-zero value made LinMode Jacobians differentiate a different truncated phase function than the forward model evaluated. The field remains (fixed at zero) only so stored objects and printing stay compatible.

source