Skip to content

CoreRT API ​

CoreRT owns the adding-doubling solver, the RT model containers, and the linearized/Jacobian-facing model layout.

Forward and Linearized Modes ​

vSmartMOM.FwdMode Type
julia
FwdMode()

Marker selecting the forward radiative-transfer model construction path.

source
vSmartMOM.LinMode Type
julia
LinMode()

Marker selecting the linearized model construction path used for Jacobians.

source
vSmartMOM.CoreRT.model_from_parameters_lin Function
julia
model_from_parameters_lin(params; kwargs...)

Convenience alias for model_from_parameters(LinMode(), params; kwargs...). Returns (model, lin_model).

source
vSmartMOM.CoreRT.rt_run_lin Function
julia
rt_run_lin(model, lin_model, NAer, NGas, NSurf; i_band=1)

Convenience alias for the linearized rt_run overload. Equivalent to rt_run(model, lin_model, NAer, NGas, NSurf; i_band).

source

Model Types ​

vSmartMOM.CoreRT.AbstractRTModel Type

Abstract base type for all RT models

source
vSmartMOM.CoreRT.RTModel Type
julia
RTModel{ARCH, FT} <: AbstractRTModel{ARCH, FT}

Central model state for vSmartMOM radiative transfer.

Only two type parameters:

  • ARCH: compute architecture (CPU or GPU)

  • FT: floating-point precision (Float32 or Float64)

All physics sub-components are organized hierarchically:

  • solver: RT solver configuration (polarization, quadrature, truncation)

  • geometry: observation geometry (SZA, VZA, VAZ, observer altitude)

  • quad_points: precomputed quadrature points and weights

  • atmosphere: atmospheric state (profile + spectral bands)

  • optics: all precomputed optical properties (Rayleigh, aerosol, absorption)

  • surfaces: per-band surface BRDF models

Fields

  • architecture: Compute architecture (CPU/GPU)

  • solver: RT solver configuration

  • numerics: Numerical-knob parameters (doubling threshold, etc.) — see RTNumericalParameters

  • geometry: Observation geometry

  • quad_points: Quadrature points and weights

  • atmosphere: Atmospheric state (profile + spectral bands)

  • optics: Precomputed optical properties

  • surfaces: Surface models, one per spectral band

  • sources: Source-term scene (v0.6 source-term refactor). Defaults to a single SolarBeam — i.e. the historical unit-Stokes-I direct solar illumination — and can be overridden at rt_run call time via the sources= kwarg. Mirrors the abstract-typed surfaces slot.

source
vSmartMOM.CoreRT.SolverConfig Type
julia
SolverConfig{FT, PT, QT}

Immutable RT solver configuration: polarization mode, quadrature scheme, Fourier truncation, and Legendre truncation. Fixed after model construction; never differentiated.

Fields

  • polarization_type: Type of polarization (Stokes_I / IQU / IQUV)

  • quadrature_type: Quadrature type (RadauQuad / GaussLegQuad)

  • m_max_bands: Per-band Fourier loop bound, order semantics: loop runs m = 0:m_max_bands[iBand]. Equals n_fourier_moments_bands .- 1.

  • n_fourier_moments_bands: Per-band Fourier moment count: n_fourier_moments_bands[iBand] = m_max_bands[iBand] + 1. Provided so consumers don't have to remember the count↔order convention.

  • l_max: Per-band max truncated Legendre index

  • l_trunc: Legendre truncation order (user-specified)

  • Δ_angle: Exclusion angle for forward peak [deg]

  • depol: Depolarization factor

  • use_component_traits: Phase C flag: when true (current default), per-band m_max_bands are derived from component_m_max(c, ctx) traits across active components — Cox-Munk / RPV / RossLi / canopy run to their full user_l_cap instead of the historical half-truncated aggregator. Flip to false to fall back to the legacy min(ceil((l_max+1)/2), params.max_m) aggregator (bit-equal to Phase B).

source
vSmartMOM.CoreRT.Atmosphere Type
julia
Atmosphere{FT, VMR}

The atmospheric column: profile data plus spectral grid definitions.

Fields

  • profile: Atmospheric profile (T, p, q, vmr, vcd, etc.)

  • spec_bands: Spectral bands — Vector of wavenumber grids, one per band

source
vSmartMOM.CoreRT.RayleighScattering Type
julia
RayleighScattering{FT, GC}

Precomputed Rayleigh and Cabannes (elastic) scattering properties. Derived from the depolarization ratio; fixed for a given spectral band.

Fields

  • greek_rayleigh: Greek coefficients for total Rayleigh scattering (single set or per-band)

  • greek_cabannes: Greek coefficients for Cabannes (pure elastic) Rayleigh per band

  • ϖ_Cabannes: Pure elastic (Cabannes) fraction of Rayleigh scattering per band

source
vSmartMOM.CoreRT.AerosolState Type
julia
AerosolState{FT, AO}

Per-band aerosol scattering optics and optical depth profiles. Primary differentiable state for aerosol retrievals.

Fields

  • aerosol_optics: Truncated aerosol optics: aerosol_optics[iBand][iAer]

  • τ_aer: Aerosol optical depth profiles: τ_aer[iBand][iAer, nSpec, iLayer] — 3-D for both the forward and linearized paths, enabling per-wavelength aerosol optical depth within wide spectral bands.

source
vSmartMOM.CoreRT.Optics Type
julia
Optics{FT, RS, AS}

Container for all precomputed optical properties that feed into the RT solver. Groups Rayleigh scattering, aerosol scattering, absorption optical depths, and Rayleigh optical depths.

AD boundary: aerosols and τ_abs hold differentiable state; rayleigh and τ_rayl are typically fixed.

Fields

  • rayleigh: Rayleigh/Cabannes scattering properties

  • aerosols: Aerosol scattering optics and optical depths

  • τ_abs: Absorption optical depth: τ_abs[iBand][nSpec × nLayers]

  • τ_rayl: Rayleigh optical depth: τ_rayl[iBand][nSpec × nLayers]

source
vSmartMOM.CoreRT.OpticsLin Type
julia
OpticsLin{FT}

Linearized (Jacobian) counterpart of Optics. Each field stores derivatives of the corresponding forward-model field with respect to the physical state vector.

Fields

  • τ̇_abs: ∂τ_abs/∂VMR(gas,z): [NGas*Nz × nSpec × Nz], species-major

  • τ̇_aer: ∂τ_aer/∂x per band: Vector of arrays [NAer × 7 × nSpec × nLayers]

  • lin_aerosol_optics: Linearized aerosol optics per band per aerosol

source
vSmartMOM.CoreRT.vSmartMOM_Parameters Type
julia
vSmartMOM_Parameters{FT<:Real}

Top-level container for all user-specified model parameters before any derived quantities are computed. Groups radiative-transfer settings (spectral bands, BRDF, quadrature, polarization), observation geometry (SZA, VZA, VAZ), the atmospheric profile (T, p, q), and optional absorption/scattering sub-parameter structs.

Fields

  • spec_bands: Spectral bands — Vector of wavenumber grids, one per band

  • brdf: Surface (Bidirectional Reflectance Distribution Function)

  • quadrature_type: Quadrature type for RT streams (RadauQuad/GaussLegQuad)

  • polarization_type: Type of polarization (I/IQ/IQU/IQUV)

  • max_m: Hard cutoff for maximum number of Fourier moments to loop over

  • Δ_angle: Exclusion angle for forward peak in degrees (legacy — see truncation)

  • l_trunc: Truncation length for legendre terms (scalar for now, can do nBand later)

  • truncation: Phase-function truncation method. Defaults to δBGE(l_trunc, Δ_angle) for backward compatibility; set explicitly to NoTruncation() for canopy-only or smooth-phase-function runs (Sanghavi & Stephens 2015 §2 — the f_tr → 0 limit). The Δ_angle field above is then ignored.

  • depol: Depolarization factor

  • numerics: Numerical-knob parameters (doubling threshold, etc.) — defaults to RTNumericalParameters{FT}() if not specified in YAML.

  • float_type: Float type to use in the RT (Float64/Float32)

  • architecture: Architecture to use for calculations (CPU/GPU)

  • sza: Solar zenith angle [deg]

  • vza: Viewing zenith angles [deg]

  • vaz: Viewing azimuthal angles [deg]

  • obs_alt: Observer altitude(s) [km above BOA]; scalar and vector forms have distinct output semantics

  • T: Temperature Profile [K]

  • p: Pressure Profile [hPa]

  • q: Specific humidity profile

  • profile_reduction_n: Length of profile reduction

  • absorption_params: Optional struct that holds all absorption-related parameters

  • scattering_params: Optional struct that holds all aerosol scattering-related parameters

  • nstreams: Phase D — primary user-facing resolution knob: weighted streams per hemisphere. Public contract: stream_l_cap = 2·nstreams - 1. When nothing, the parser populated max_m/l_trunc from a legacy YAML; new-schema configs set this and may omit max_m/l_trunc.

  • m_max_override: Phase D — explicit per-band Fourier loop bound (order, optional). When set, clamps the trait-aggregator output. Lower hard cap on top of stream_l_cap.

  • stream_l_cap: Phase D — derived projection cap: 2·nstreams - 1 for new schema, or l_trunc for legacy. Internally consumed by the trait aggregator and the truncation-resolver. Always populated; defaults to l_trunc until nstreams is explicitly set.

  • legacy_l_cap_override: Phase D — legacy l_trunc value retained verbatim from YAML when the user explicitly set it. nothing for new-schema configs that derive stream_l_cap from nstreams instead.

  • greek_beta_cutoff: Optional post-Mie Fourier-support cutoff based specifically on the magnitude of the integrated aerosol Greek coefficient β_l. The largest degree satisfying abs(β_l) >= greek_beta_cutoff is retained for each aerosol; the per-band maximum enters the component trait aggregator. This changes Fourier-loop support after Mie integration but does not resize the previously selected quadrature. nothing disables the cutoff.

  • rayleigh_anisotropy: Wavelength dependence of the N₂/O₂ polarizability tensor used for the automatic (depol < 0) Rayleigh depolarization of greek_rayleigh and τ_rayl: DispersiveAnisotropy() (default) or ConstantAnisotropy() (the Raman-module constants with a wavelength-independent anisotropy, the v2.2 pre-release behaviour). The Cabannes fractions and Raman cross sections always use the Raman-module constants. YAML key radiative_transfer.rayleigh_anisotropy.

source
vSmartMOM.CoreRT.AtmosphericProfile Type
julia
AtmosphericProfile{FT, VMR}

Vertical atmospheric state on a pressure grid, ordered from TOA to BOA.

Stores temperature, pressure (full levels and half levels), humidity, dry and wet vertical column densities, and volume mixing ratios for trace gases. Constructed internally by model_from_parameters from the raw arrays in vSmartMOM_Parameters.

p_full has N levels; p_half has N+1 layer-boundary pressures. T, q, vmr_h2o, vcd_dry, vcd_h2o, and Δz are layer quantities of length N.

Fields

  • T: Temperature Profile

  • p_full: Pressure Profile (Full)

  • q: Specific humidity profile

  • p_half: Pressure Levels

  • vmr_h2o: H2O/dry-air molar-ratio profile

  • vcd_dry: Vertical Column Density (Dry)

  • vcd_h2o: Vertical Column Density (H2O)

  • vmr: Volume Mixing Ratio of Constituent Gases

  • Δz: Layer height (meters)

source
vSmartMOM.CoreRT.ObsGeometry Type

Observation geometry, including the resolved vertical-output interfaces.

source
vSmartMOM.CoreRT.QuadPoints Type
julia
QuadPoints{FT}

Gauss or Gauss-Radau quadrature points and weights used for the angular discretisation of the radiative transfer equation. Also stores the cosine of the solar zenith angle (μ₀). In the legacy representation iμ₀ indexes that direction in the diffuse grid; external-solar models instead use a zero sentinel there and iμ₀_phase for the exact phase-source column.

Fields

  • μ₀: μ₀, cos(SZA)

  • iμ₀: Index of μ₀ in the diffuse grid, or zero when the solar direction is external

  • iμ₀Nstart: First Stokes index of μ₀ in qp_μN, or zero for external-solar models

  • qp_μ: Diffuse-operator quadrature points

  • wt_μ: Weights of diffuse-operator quadrature points

  • qp_μN: Quadrature points (repeated for polarizations)

  • wt_μN: Weights of quadrature points (repeated for polarizations)

  • Nquad: Total diffuse-grid nodes (weighted streams + zero-weight VZAs, and legacy SZA)

  • Nstreams: Number of weighted streams per hemisphere (count of nonzero weights). Public contract: stream_l_cap = 2·Nstreams - 1.

  • phase_qp_μ: Angular grid used to construct phase matrices; may additionally contain the external solar direction.

  • iμ₀_phase: Index of μ₀ in phase_qp_μ.

  • iμ₀Nstart_phase: First Stokes index of μ₀ in the polarization-expanded phase grid.

  • external_solar: Whether the direct solar direction is external to the diffuse-operator quadrature grid.

source
vSmartMOM.CoreRT.CompositeLayer Type
julia
CompositeLayer{FT} <: AbstractLayer

Accumulated (composite) layer matrices produced by the interaction step of the adding/doubling RT method. Stores the reflectance (R), transmission (T), and source (J) matrices for the combined atmospheric slab from the top of atmosphere down to the current layer.

Sign convention: − = outgoing (upward, decreasing τ), + = incoming (downward, increasing τ).

Fields

  • R⁻⁺::AbstractArray{FT,3}: reflectance from incoming (+) to outgoing (−)

  • R⁺⁻::AbstractArray{FT,3}: reflectance from outgoing (−) to incoming (+)

  • T⁺⁺::AbstractArray{FT,3}: transmission in the incoming (+) direction

  • T⁻⁻::AbstractArray{FT,3}: transmission in the outgoing (−) direction

  • J₀⁺::AbstractArray{FT,3}: source in the incoming (+) direction

  • J₀⁻::AbstractArray{FT,3}: source in the outgoing (−) direction

source
vSmartMOM.CoreRT.AddedLayer Type
julia
AddedLayer{FT} <: AbstractLayer

Single homogeneous layer matrices produced by the elemental/doubling steps of the RT solver. After doubling, these represent a full atmospheric layer that is subsequently combined with the composite layer via the interaction step.

Lower-case field names (r, t, j) distinguish added-layer quantities from composite-layer quantities (R, T, J).

Fields

  • r⁻⁺::AbstractArray{FT,3}: reflectance from incoming (+) to outgoing (−)

  • t⁺⁺::AbstractArray{FT,3}: transmission in the incoming (+) direction

  • r⁺⁻::AbstractArray{FT,3}: reflectance from outgoing (−) to incoming (+)

  • t⁻⁻::AbstractArray{FT,3}: transmission in the outgoing (−) direction

  • j₀⁺::AbstractArray{FT,3}: source in the incoming (+) direction

  • j₀⁻::AbstractArray{FT,3}: source in the outgoing (−) direction

  • temp1, temp2: pre-allocated workspace arrays to avoid allocations

  • temp1_ptr, temp2_ptr: CUDA pointer arrays (ignored on CPU)

source
vSmartMOM.CoreRT.RTNumericalParameters Type
julia
RTNumericalParameters{FT}

Centralised home for tunable numerical knobs that were previously hardcoded in the RT kernels. Lives once on vSmartMOM_Parameters / RTModel, threaded through to the kernels where needed. Designed to grow as more constants are lifted (Glint Legendre expansion cap, etc.).

Fields

  • dτ_max_threshold: Doubling resolution: cap on the elemental-layer optical depth dτ_max = threshold · μ_min, where μ_min is the smallest TRUE (positive-weight) quadrature stream cosine. Smaller values produce finer initial layers (more doublings → more accurate single-scatter at the cost of more roundoff through doubling). Default 0.001 is conservative legacy; raising to ~0.01 reduces doubling iterations.

  • dτ_min_floor: Absolute floor on dτ_max (and hence on dτ_initial), expressed as a multiple of eps(FT). Prevents the doubling discretisation from collapsing below FT precision regardless of geometry. Default 1024·eps(FT) ≈ 1.2e-4 for Float32 / 2.3e-13 for Float64 — the F32 value caps ndoubl around 13–14 for τ_total≤1 (well within Float32 representability); the F64 value is far below any sensible threshold·μ_min so it never activates in practice.

  • blas_threads: BLAS thread cap applied at every rt_run invocation. nothing means leave the current BLAS thread setting alone. Why cap: the batched-GEMM call sites in cpu_batched.jl and the elemental/doubling/interaction kernels operate on small matrices (NSTREAMS·n_stokes ≈ 12-48 per side) tiled across a wide spectral batch. Multi-threaded BLAS coordination dominates the work at that matrix size and serializes against the spectral-batch parallelism. The VLIDORT baseline harness empirically picked 8 as the sweet spot on a 128-thread machine; 1 is the most defensive choice when a single process owns all spectral points and threading happens at the batched outer level. This is a process-wide, persistent setting, despite being selected from a model: it is not restored after the run and models with different values must not execute concurrently in one process. Set it once for an application via the numerics.blas_threadsYAML key, or programmatically withRTNumericalParameters{FT}(blas_threads = 8).

  • verbose: Verbose output flag. When false (default) rt_run/rt_run_lin suppress the per-call print_timer() dump at the end of the run. Set true (YAML key numerics.verbose: true or programmatically RTNumericalParameters{FT}(verbose = true)) to bring back the timing tree — useful for profiling, noisy for production loops. Per-band Computing profile for X... and AOD at band X... messages were demoted to @debug and are silent unless you set ENV["JULIA_DEBUG"] = "vSmartMOM".

  • fourier_convergence: Azimuthal Fourier-loop termination strategy. The package-wide default runs all moments. Retrieval/truth workflows may explicitly select IntensityConvergence or StokesConvergence after validating the omitted tail for their viewing geometry.

  • ss_correction: Single-scattering (first-order) truncation correction — see AbstractSingleScatteringCorrection. Default NoSSCorrection (bit-identical historical behavior). Set TMSCorrection for the Nakajima–Tanaka correction; YAML key numerics.ss_correction: tms.

source
vSmartMOM.CoreRT.AbstractFourierConvergence Type
julia
AbstractFourierConvergence

Termination strategy for the azimuthal Fourier loop. The static component traits first determine the highest order the scene can contribute; a runtime strategy may then stop below that ceiling after the accumulated directional radiances have converged. AllFourierMoments retains the complete trait-selected series, IntensityConvergence monitors Stokes I, and StokesConvergence monitors every propagated Stokes component.

source
vSmartMOM.CoreRT.AllFourierMoments Type
julia
AllFourierMoments()

Evaluate every Fourier moment through the scene's resolved m_max. This is the package default and is bit-identical to the historical full loop.

source
vSmartMOM.CoreRT.IntensityConvergence Type
julia
IntensityConvergence(tolerance; min_m=3, n_consecutive=2)

Terminate after n_consecutive eligible Fourier moments satisfy

at every monitored view and wavelength. A failing moment resets the counter. The test is disabled through min(min_m-1, m_max): the default min_m=3 prevents a structural zero at m=1 from hiding a nonzero m=2 contribution, without manufacturing moments beyond a physically shorter precomputed series. The combined analytic-linearization path follows the same forward stopping moment and propagates all Jacobian contributions through it.

source
vSmartMOM.CoreRT.StokesConvergence Type
julia
StokesConvergence(tolerance; min_m=3, n_consecutive=2)

Terminate after n_consecutive eligible Fourier moments satisfy

independently for every propagated Stokes component at every monitored view and wavelength. Thus an I/Q/U solve must converge in I, Q, and U. The same low-order guard and forward/linearized stopping contract as IntensityConvergence apply.

source

Single-Scattering Correction ​

vSmartMOM.CoreRT.AbstractSingleScatteringCorrection Type
julia
AbstractSingleScatteringCorrection

Strategy for correcting the single-scattering (first-order) part of the radiance for phase-function truncation. δ-M/δ-BGE truncation removes the aerosol forward peak from the Fourier-expanded optics — harmless for the multiple-scattered field but visibly wrong for the first-order term at the view geometries (aureole, rainbow/glory structure). Concrete strategies:

  • NoSSCorrection — historical behavior, bit-identical to the uncorrected solver (default).

  • TMSCorrection — Nakajima & Tanaka (1988) TMS: the truncated direct-beam source is zeroed at the zero-weight view nodes in every Fourier mode, and the exact (untruncated) single scattering is added in real space after the Fourier sum, using the solver's own δ-scaled transmittances.

Selected via RTNumericalParameters.ss_correction (YAML: numerics.ss_correction: none | tms). The hierarchy is open for future variants (e.g. an IMS-augmented correction, canopy hot-spot or sun-glint corrections) — add a subtype, never a boolean flag.

Forward/linearized parity policy: a strategy must either be implemented in BOTH the forward and linearized drivers or be explicitly rejected by the one that lacks it — silent divergence between the two is a bug by definition (see docs/dev_notes/forward_lin_parity_policy.md).

source
vSmartMOM.CoreRT.NoSSCorrection Type
julia
NoSSCorrection()

No single-scattering correction — the truncated Fourier solution is returned as-is. Bit-identical to the historical solver. Default.

source
vSmartMOM.CoreRT.TMSCorrection Type
julia
TMSCorrection()

Nakajima–Tanaka TMS single-scattering correction ("never generate, then add"): zero the truncated beam source at zero-weight view nodes in every Fourier mode; add the exact single scattering post-Fourier-sum from the untruncated phase functions, with the N–T scaling ω/(1−ωf) paired with δ-scaled optical depths (implemented as β_scat/Δτ′, which is algebraically identical and reuses the solver's own scaled τ_sum). Requires an unpolarized solar source; currently forward-elastic (noRS) only — the linearized driver and Raman paths reject it explicitly, as do interior sensor heights and per-moment stream exports (rt_run_streams). Combined-band solves evaluate the correction per band.

source
vSmartMOM.CoreRT.rt_run_ss_exact Function
julia
rt_run_ss_exact(model; i_band=1, sources=nothing) -> Array (nvza × nStokes × nSpec)

Standalone exact (untruncated) single-scattering TOA radiance — the analytic first-order reference. Uses the same engine as the TMS correction with f = 0 semantics: UNSCALED optical depths and the untruncated phase functions, no Fourier loop at all (O(nlayers × ngeom × nSpec)). This is the exact-SS reference the Fourier-space rt_run_ss (which runs on truncated per-mode optics) cannot be: the azimuth expansion of a peaked kernel converges too slowly.

Uses model.sources unless overridden by sources. Supports one unpolarized SolarBeam with its supplied irradiance, or NoSource for zero radiance. This atmosphere-only reference includes no surface reflection or emission; nonzero/retrievable SurfaceSIF, thermal sources, and polarized incident beams throw ArgumentError.

source

Fourier Loop Bound ​

The azimuthal Fourier loop bound is derived per band at model build: each active component (surface, sources, Rayleigh, aerosol modes) declares the highest Fourier order it can contribute via the component_m_max trait, and the band bound is their maximum, clamped by the user's max_m. At run time, the IntensityConvergence strategy above may terminate the loop below that ceiling once the accumulated Stokes-I contribution at the user view angles stabilises — the two mechanisms compose.

vSmartMOM.CoreRT.component_m_max Function
julia
component_m_max(component, ctx) -> Int

Maximum Fourier order this component contributes to. Used by _aggregate_m_max when SolverConfig.use_component_traits == true. The ctx NamedTuple carries user_l_cap, stream_l_cap, truncation, and optionally greek_beta_cutoff, so each component can declare a tighter bound. The aerosol cutoff applies only to abs(β_l); the other Greek families are not Fourier-support tests.

source
vSmartMOM.CoreRT._aggregate_m_max Function
julia
_aggregate_m_max(active_components, ctx) -> Int

Aggregate the per-component Fourier supports into a single per-band loop bound (order). Internal helper used by model_from_parameters when SolverConfig.use_component_traits == true.

active_components is an iterable of components participating in the band — typically (rayleigh, aerosols..., surface, sources). The aggregator takes the maximum Fourier support, then clamps to ctx.user_l_cap and (if non-nothing) ctx.m_max_override.

Note the naming collision with rt_run's m_max_override keyword, which has the OPPOSITE sense: here (model build) the override NARROWS the trait ceiling; there (runtime, used by rt_run_atmosphere to size surface-replay caches) it only WIDENS the loop bound. They are independent knobs that happen to share a name.

source

Column Core and TOA Fast Path ​

Every public forward entry point funnels into the internal single-column solver _rt_run_column; rt_run_toa is the lean TOA-only wrapper used by retrieval loops that never consume BOA or interior fields (documented under the top-level API).

vSmartMOM.CoreRT._rt_run_column Function
julia
_rt_run_column(RS_type, model::RTModel, iBand; kwargs...)

Internal single-column solver: the Fourier loop over elemental → doubling → interaction → surface for one full atmosphere/surface column. Every public entry point (rt_run, rt_run_toa, rt_run_atmosphere) funnels into this function; user code should call those instead — they own the observer conventions and result wrapping.

Returns the raw SFI tuple (R, T, ieR, ieT, hdr, bhr_uw, bhr_dw) (or (R, T) without SFI; only R — elastic and inelastic — under toa_only).

Keyword arguments

  • sources::AbstractSource: optional source-term set (v0.6 source-term refactor, Phase 2). The solar F₀ of the solve is, in order of precedence, the F₀ of the first SolarBeam in sources, a user preset RS_type.F₀ (nonzero, shaped (nStokes, nSpec)), then model.sources; it is routed into the RS_type.F₀ carrier ahead of the SFI kernel call. This core does not restore the carrier; the public entry points do (_with_call_scoped_F₀), so an explicit sources= applies to its call only.

  • streams_callback: fires with the per-moment stream state (used by rt_run_streams).

  • Atmosphere/surface-split hooks consumed by rt_run_atmosphere: atm_snapshot_callback (fired once per Fourier moment right after the layer loop, before the surface step), stop_after_atmosphere (skip the surface step + postprocessing), and m_max_override (widen, never narrow, the Fourier loop bound). All three default to a no-op / bit-exact-compatible value.

  • toa_only::Bool: skip BOA/transmittance postprocessing and return only the TOA field (multisensor fast path).

The Fourier loop bound comes from the model build (component_m_max traits aggregated per band); within the loop, model.numerics.fourier_convergence (an AbstractFourierConvergence strategy) may terminate the series early once the observed Stokes-I contribution stabilises. The two mechanisms compose: traits set the ceiling, convergence stops below it.

source

Height-Aware Forward and Linearized Output ​

rt_run returns a named ObserverRTResult. Endpoint fields are nothing when not selected by geometry.obs_alt; strict-interior outputs are stored as LevelRadiance records in result.levels. See the geometry schema for the scalar/vector selection convention.

The analytic linearized solver returns ObserverRTResultLin. It stores endpoint Jacobians alongside the endpoint radiances and uses LevelRadianceLin records for co-located forward fields and Jacobians at each strict-interior height. Both result types retain their historical tuple iteration conventions.

vSmartMOM.CoreRT.ObserverRTResult Type
julia
ObserverRTResult

Height-aware forward-RT result. toa is upwelling at TOA, boa is downwelling at BOA, and levels contains co-located diffuse up/down radiances plus a separate unscattered solar component at requested interior heights. Unrequested endpoints are nothing.

The result remains iterable/indexable over the historical seven forward slots (toa, boa, inelastic_toa, inelastic_boa, hdr, bhr_uw, bhr_dw) so the default obs_alt: [0] preserves historical destructuring and indexing.

source
vSmartMOM.CoreRT.LevelRadiance Type
julia
LevelRadiance

Radiances at one strict-interior atmospheric interface. upwelling and downwelling are the diffuse elastic fields solved by MOM. unscattered_downwelling is the separately reported attenuated collimated solar beam at the solar ordinate (zero for other viewing zenith angles). Use total_downwelling to add the diffuse and unscattered components.

Over a CoxMunkSurface, upwelling carries the direct-solar glint as resolved by the Fourier series: the exact-glint correction (apply_ss_correction!) is applied to the TOA radiance only.

source
vSmartMOM.CoreRT.ObserverRTResultLin Type
julia
ObserverRTResultLin

Height-aware result from the analytic tangent-linear RT solver. toa and boa contain the requested endpoint radiances, their matching Jacobians are stored in toa_jacobian and boa_jacobian, and levels contains co-located up/down results at every requested strict-interior height.

For backward compatibility the result is iterable and indexable as the historical four slots (toa, boa, toa_jacobian, boa_jacobian).

source
vSmartMOM.CoreRT.LevelRadianceLin Type
julia
LevelRadianceLin

Forward radiances and analytic Jacobians at one strict-interior atmospheric interface. The three radiance fields mirror LevelRadiance; each *_jacobian field appends the ParameterLayout dimension as its last axis.

The unscattered solar beam is kept separate from the diffuse MOM field. Use total_downwelling and total_downwelling_jacobian when the sum of both components is desired.

As in the forward LevelRadiance, a CoxMunkSurface glint enters upwelling and upwelling_jacobian through the Fourier series only; the exact-glint correction is applied at TOA only.

source
vSmartMOM.CoreRT.total_downwelling Function
julia
total_downwelling(level::LevelRadiance)

Return the total elastic downwelling radiance at an interior level: the MOM diffuse field plus the separately reported unscattered solar beam.

source

Return diffuse plus unscattered downwelling radiance at an interior level.

source
vSmartMOM.CoreRT.total_downwelling_jacobian Function

Return the analytic Jacobian of total_downwelling.

source

Source-term types (SolarBeam, BlackbodySource, SurfaceSIF, SourceSet, NoSource, AD-mode traits, prepare_source / surface_source_contribute! dispatchers) are documented on the Source terms page.

Batch Scene Processing ​

BatchContext enables efficient multi-scene loop workflows (e.g. ensemble retrievals, parameter sweeps) by caching expensive one-time setup work (Mie, Fourier decomposition, HITRAN parsing) and exposing cheap per-scene update functions.

vSmartMOM.CoreRT.BatchContext Type
julia
BatchContext

Pre-built context for efficient batch radiative-transfer over many atmospheric scenes that share the same spectral configuration, geometry, surface BRDF, aerosol microphysics, and quadrature setup but differ in temperature/pressure profiles, humidity, and/or trace-gas VMRs.

Construct once from a vSmartMOM_Parameters object; then call update_model! before each rt_run invocation.

Example usage

julia
params = parameters_from_yaml("my_config.yaml")
ctx    = BatchContext(params)

for scene in scenes
    update_model!(ctx;
        T      = scene.T,        # length-Nz temperature profile [K]
        p_half = scene.p_half,   # length-(Nz+1) half-level pressures [hPa]
        q      = scene.q,        # length-Nz or -(Nz+1) specific humidity [kg/kg]
        vmr    = scene.vmr)      # Dict{String,Any} — same keys as params.absorption_params.vmr
    R, T = rt_run(ctx.model)
    # ... process R, T ...
end

What may change per scene

  • Temperature profile T [K]

  • Pressure half-levels p_half [hPa]

  • Specific humidity q [kg/kg]

  • Trace-gas VMRs via the vmr keyword dict (same keys as the original params.absorption_params.vmr; adding new species requires a full rebuild)

What requires a full rebuild (new BatchContext)

  • Spectral bands (params.spec_bands)

  • Geometry (SZA, VZA, VAZ, observer altitude)

  • Surface BRDF

  • Polarization type, quadrature type, truncation

  • A different number of aerosols or gas species

Aerosol changes do not require a rebuild: use update_aerosol_loading! for τ_ref / vertical-placement changes (cheap, no Mie recomputation) and update_aerosol_microphysics! for size-distribution / refractive-index changes (reruns Mie for that aerosol and re-derives the Fourier-loop bounds).

Thread safety

One BatchContext per thread. update_model! and rt_run(ctx.model) mutate shared arrays and must not be called concurrently on the same ctx. If you run parallel scenes, create one BatchContext per worker thread (each BatchContext owns its own RTModel with independent mutable arrays).

Derived caches do not survive updates

An AtmosphereRTCache built from ctx.model (via rt_run_atmosphere) snapshots the optics at build time — every update_model! / update_aerosol_*! call invalidates it (and any LambertianClosure derived from it). Rebuild such caches after each scene update, never reuse them across updates. See the Fast Re-runs & Batch Processing guide.

Fields

  • model: the RTModel that is updated in place by update_model!

  • params: an owned parameter copy for molecule lists, configuration and float type. Loaded LUT payloads are shared read-only with the caller.

  • absorption_models: cached AtmosphericAbsorption.LineByLineModel objects, one per (band, species) pair — eliminates HITRAN re-parsing per scene

  • h2o_models: cached H₂O LineByLineModel per band (nothing only when the configuration has no absorption_params at all; otherwise the model is cached unconditionally so a scene whose q changes from all-zero to non-zero can apply H₂O line absorption without a rebuild)

  • k_ref: reference-wavelength aerosol extinction coefficients, one per aerosol species (needed for Phase 2 aerosol-loading updates)

  • current_T, current_p_half, current_q, current_vmr: the unreduced scene state that update_model! was last called with (initialised from params at construction). nothing keyword arguments to update_model! fall back to these stored values — never to params — so successive partial updates compose incrementally. Stored already FT-converted so repeated updates do not re-convert or drift. The unreduced state cannot be recovered from model.atmosphere.profile (which is reduced), which is why it lives here.

  • current_τ_ref, current_profile_dist: the current per-aerosol loading state (column optical depth at λ_ref and vertical distribution), initialised from params and written by both update_aerosol_loading! and update_aerosol_microphysics!. update_model!'s τ_aer redistribution reads these so a prior loading update is not silently wiped.

  • input_sources: the user-supplied source tree on the unreframed input grid; vertically resolved fields are re-interpolated from this stable copy on each scene update (avoids cumulative interpolation drift).

  • scratch_τ_abs, scratch_τ_rayl, scratch_τ_aer: context-owned trial optical-depth arrays. update_model! prepares a complete candidate scene in these reusable buffers before committing it to the live model.

  • n_bands, n_aerosols, Nz: scene-invariant dimension bookmarks

  • profile_reduction_n: reduction target passed to the anchored prepare_observer_profile path; -1 means no reduction

source
vSmartMOM.CoreRT.update_model! Function
julia
update_model!(ctx::BatchContext;
              T      = nothing,
              p_half = nothing,
              q      = nothing,
              vmr    = nothing)

Update ctx.model in place for a new atmospheric scene.

Only the arguments you supply are updated. nothing for a field means "keep the current scene value" — i.e. the value from the most recent successful update_model! call (or, if none yet, the value the BatchContext was built with). It does not fall back to the original params. This makes successive partial updates compose incrementally:

julia
update_model!(ctx; vmr = B)   # vmr = B, everything else still original
update_model!(ctx; T   = C)   # T = C AND vmr is still B (not reset to original)

The merge is per-field: T, p_half, q, and the vmr override dict are each remembered on the context (ctx.current_T, ctx.current_p_half, ctx.current_q, ctx.current_vmr) and overwritten only by the fields you pass. The remembered state is stored unreduced and already converted to the model's float type (params.float_type), so repeated updates neither re-convert nor drift.

For vmr, the supplied keys are overlaid onto the current override dict (which itself started as params.absorption_params.vmr); keys you do not pass keep their last value. The merged overrides are then layered over the configured defaults exactly as model_from_parameters does.

After this call rt_run(ctx.model) will produce radiances for the new scene. Candidate profiles and optical properties are prepared in detached storage. If validation, interpolation, spectroscopy, or continuum loading throws, the live model and remembered scene remain at the last successful state and are safe to use or retry. Caller-owned arrays are copied; mutating them after a successful call does not mutate the context.

Keyword arguments

  • T::AbstractVector: Temperature profile [K], length Nz (full levels, TOA to BOA). If profile_reduction_n != -1, pass the unreduced profile length (same as the original params.T).

  • p_half::AbstractVector: Half-level pressures [hPa], length Nz + 1 (or N_orig + 1 before reduction). Must be strictly positive and monotonically increasing (TOA to BOA, i.e. p_half[end] is surface pressure).

  • q::AbstractVector: Specific humidity [kg/kg], with either Nz layer-center values or Nz + 1 interface values on the unreframed input grid.

  • vmr::Dict{String,Any}: Trace-gas volume mixing ratios. Keys must be a subset of the configured molecules — supplying an unknown species key raises an error (add it to params.absorption_params.vmr and rebuild the context instead). Values may be scalars, length-Nz layer vectors, or length-Nz+1 interface vectors.

Thread safety

update_model! and rt_run(ctx.model) share mutable arrays. Do not call them concurrently on the same ctx. Use one BatchContext per thread for parallelism.

Example

julia
params = parameters_from_yaml("config.yaml")
ctx    = BatchContext(params)

for scene in scenes
    update_model!(ctx;
        T      = scene.T,
        p_half = scene.p_half,
        q      = scene.q,
        vmr    = Dict("O2" => 0.21, "CH4" => scene.ch4_vmr))
    R, T_trans = rt_run(ctx.model)
end
source
vSmartMOM.CoreRT.update_aerosol_loading! Function
julia
update_aerosol_loading!(ctx::BatchContext, i_aer::Int;
                        τ_ref        = nothing,
                        profile_dist = nothing)

Cheap path — update only the column-integrated optical depth (τ_ref) and/or the vertical pressure distribution (profile_dist) for aerosol species i_aer, using the cached aerosol optics (aerosol_optics) and reference extinction coefficient (ctx.k_ref[i_aer]). No Mie recomputation is performed.

After this call rt_run(ctx.model) produces radiances for the new aerosol loading. All other model state (gas absorption, profile, surface) is unchanged.

All bands are staged before committing. A failed calculation leaves live optics and remembered loading unchanged; scratch buffers may contain trial data.

Keyword arguments

  • τ_ref::Real: New column optical depth at λ_ref. When nothing, the value stored in params.scattering_params.rt_aerosols[i_aer].τ_ref is kept (only the distribution is updated).

  • profile_dist::Distributions.Distribution: New vertical pressure distribution (a Distributions.jl Distribution object — the same type used in RT_Aerosol.profile). When nothing, the existing distribution is kept.

Cost

O(Nz) per band — evaluates the distribution PDF on the pressure grid and scales the optical depth profile. No Mie, no HITRAN.

What it invalidates

  • model.optics.aerosols.τ_aer[i_band][i_aer, :, :] — all bands (3-D: nSpec × nLayers).

  • Does NOT change aerosol_optics, k_ref, SolverConfig Fourier bounds, gas absorption, or Rayleigh properties.

Example

julia
ctx = BatchContext(params)
# Change τ_ref for aerosol 1 to 0.3, keep vertical distribution
update_aerosol_loading!(ctx, 1; τ_ref = 0.3)
R, T = rt_run(ctx.model)
source
vSmartMOM.CoreRT.update_aerosol_microphysics! Function
julia
update_aerosol_microphysics!(ctx::BatchContext, i_aer::Int, aerosol::Aerosol;
                              τ_ref = nothing)

Expensive path — replace the microphysics of aerosol species i_aer with a new Aerosol (new size distribution and/or refractive index) and recompute all derived quantities via the same Mie code path that model_from_parameters uses.

This updates:

  • model.optics.aerosols.aerosol_optics[i_band][i_aer] for every band

  • ctx.k_ref[i_aer] (the new reference-wavelength extinction coefficient)

  • model.solver.m_max_bands, n_fourier_moments_bands, l_max (in-place via .=) — the critical Fourier-loop re-derivation that prevents silent wrong results when the new particle size changes the length of the Greek- coefficient series.

  • model.optics.aerosols.τ_aer[i_band][i_aer, :, :] for every band (3-D: nSpec × nLayers).

After this call rt_run(ctx.model) gives the same result as building a fresh model with the new aerosol.

Mie optics, loading and Fourier bounds are prepared before committing. A failed calculation leaves the live model and remembered state unchanged.

Arguments

  • ctx: The BatchContext to update.

  • i_aer: 1-based index of the aerosol species to replace.

  • aerosol: New Aerosol (from Scattering.Aerosol), carrying the new size_distribution, nᵣ, and nᵢ.

Keyword arguments

  • τ_ref::Real: New column optical depth at λ_ref. When nothing, the current scene value ctx.current_τ_ref[i_aer] is kept (i.e. the most recently applied loading, not the original params value).

SolverConfig mutability note

SolverConfig is an immutable struct, but its three Vector fields (m_max_bands, n_fourier_moments_bands, l_max) are mutable Julia arrays and are updated in-place with .=. The scalar fields (l_trunc, Δ_angle, depol, polarization_type, quadrature_type) cannot be changed in-place; they are guaranteed to be unaffected by aerosol microphysics changes since they derive from params, not from per-aerosol optics. If a caller somehow requires a new l_trunc or Δ_angle, a full model_from_parameters rebuild is necessary.

Cost

Full Mie per band (O(n_bands × nquad_radius × size_distribution_points)), i.e. the same cost as the aerosol section of model_from_parameters. The radius nodes are those of ctx.params.scattering_params (r_min, r_max, nquad_radius); they do not move with the new size distribution.

What it invalidates

  • aerosol_optics[i_band][i_aer] — all bands.

  • ctx.k_ref[i_aer].

  • solver.m_max_bands, solver.n_fourier_moments_bands, solver.l_max.

  • τ_aer[i_band][i_aer, :, :] — all bands (3-D: nSpec × nLayers).

Example

julia
ctx = BatchContext(params)
# Save initial state for comparison
R_A, _ = rt_run(ctx.model)

# Replace aerosol 1 with larger particles (reff ~2 μm)
new_aerosol = Aerosol(LogNormal(log(2.0), 0.4), 1.3, 1e-8)
update_aerosol_microphysics!(ctx, 1, new_aerosol; τ_ref = 0.1)
R_B, _ = rt_run(ctx.model)
source

Atmosphere/Surface Split ​

Run the atmosphere phase once, replay the surface phase per BRDF — see the Fast Re-runs & Batch Processing guide.

vSmartMOM.CoreRT.rt_run_atmosphere Function
julia
rt_run_atmosphere(model::RTModel;
                  target_brdfs = nothing,
                  i_band::Integer = 1,
                  sources = nothing,
                  cache_mode::Symbol = :auto) -> AtmosphereRTCache

Run the atmosphere phase of rt_run once and cache the per-Fourier-moment composite-layer snapshot. The returned AtmosphereRTCache can be consumed by rt_run_surface repeatedly to swap surfaces over the same atmosphere without paying the m-loop layer-accumulation cost again.

Keyword arguments

  • target_brdfs — surfaces the caller intends to swap in via rt_run_surface. The cache's Fourier loop bound is sized for max(model's own m_max, maximum(component_m_max(b, ctx) for b in target_brdfs)), so e.g. a Lambertian model can still serve a Cox-Munk replay. Defaults to model.surfaces (least surprising — the cache always covers the model's own surface at minimum).

  • i_band::Integer = 1 — spectral band index.

  • sources — source set, as in rt_run.

  • cache_mode::Symbol = :auto — :full (all six composite blocks at every Fourier moment), :slim (proposal §4: six blocks only at m=0; J₀⁺/ J₀⁻/per-source slots only at m>0 — valid only when every target_brdfs entry is Lambertian-family), or :auto (picks :slim iff every entry in brdfs — target_brdfs, or model.surfaces when left at its default — has component_m_max(b, ctx) == 0, else :full). See the memory-footprint note on AtmosphereRTCache.

Only noRS is supported (see the module-level scope note above); CanopySurface in target_brdfs (or the model's own surface, when target_brdfs is left at its default) raises ArgumentError.

Caches do not survive scene updates

The cache snapshots the model's optics at build time. Any subsequent update_model!, update_aerosol_loading! / update_aerosol_microphysics!, or direct mutation of model.optics silently invalidates it (and any LambertianClosure derived from it) — the replay would keep answering for the old scene. Per-scenario order must always be: update the scene → rt_run_atmosphere → sweep surfaces. See the Fast Re-runs & Batch Processing guide.

source
vSmartMOM.CoreRT.rt_run_surface Function
julia
rt_run_surface(cache::AtmosphereRTCache, brdf;
               verbose::Bool = false)
    -> (R_SFI, T_SFI, ieR_SFI, ieT_SFI, hdr, bhr_uw, bhr_dw)

Replay the surface phase of rt_run against a cached atmosphere. Returns the same tuple shape as rt_run(model) under SFI (cache.SFI is always true — rt_run hardcodes SFI).

Mirrors rt_run.jl's surface block (create_surface_layer! → surface_source_contribute! → interaction! → interaction_hdrf! → postprocessing_vza!/_hdrf!) call-for-call, so results are bit-exact with a monolithic rt_run call using the same brdf and the same Fourier loop bound.

Throws ArgumentError for CanopySurface (not yet supported — see the module-level scope note in this file), and when component_m_max(brdf, ctx) — evaluated against the model's Fourier-support cap (cache.user_l_cap), not the cache's own width (cache.m_max); see AtmosphereRTCache — exceeds cache.m_max: rebuild the cache with this BRDF (or one of the same family) in target_brdfs. Also throws ArgumentError when cache.cache_mode == :slim and brdf is not Lambertian-family (component_m_max(brdf, ctx) != 0) — the cache's m>0 blocks were discarded on the assumption that only a Lambertian-family surface would ever be replayed; rebuild with cache_mode = :full.

source
vSmartMOM.CoreRT.rt_run_multi_surface Function
julia
rt_run_multi_surface(model::RTModel, brdfs::AbstractVector;
                     i_band::Integer = 1, sources = nothing)

Sweep a set of BRDFs over a single atmosphere. The atmosphere is built once via rt_run_atmosphere (with target_brdfs = brdfs, so the cache is sized for every BRDF in the sweep); each BRDF is then run through rt_run_surface.

Returns a Vector of rt_run-style tuples, one per BRDF in brdfs, in the same order.

Acceptance

rt_run_multi_surface(model, [model.surfaces[1]])[1] is bit-exact with rt_run(model) — same kernels, same values, no separate code path for the "just replay the model's own surface" case.

source
vSmartMOM.CoreRT.AtmosphereRTCache Type
julia
AtmosphereRTCache{FT, AT3, AT1, AT2, AT, PT, QT, ARCH, RS, PS, IS, SCI}

Per-Fourier-moment snapshot of the atmosphere-only composite layer, plus the context needed to drive the surface phase of rt_run. Built by rt_run_atmosphere and consumed by rt_run_surface; the two together replace one full rt_run call when iterating BRDFs over the same atmosphere — the m-loop layer accumulation (elemental → doubling → interaction over Nz layers) runs once per scene, not once per surface.

Memory footprint

:full mode (cache_mode): device-resident (m_max+1) copies of six (NquadN, NquadN, nSpec)-ish blocks (R⁻⁺, R⁺⁻, T⁺⁺, T⁻⁻ are (NquadN, NquadN, nSpec); J₀⁺, J₀⁻ are (NquadN, 1, nSpec), budgeted the same order for a rough estimate):

bytes ≈ (m_max+1) · sizeof(FT) · nSpec · NquadN · (4·NquadN + 2)

E.g. Float64, NquadN=60, nSpec=10⁴, m_max=20: ≈ 24 GB — fine for small nSpec / scalar sweeps, prohibitive for full spectral bands on GPU.

:slim mode (proposal §4, Lambertian-family targets only): the four big (NquadN, NquadN, nSpec) blocks are retained only at m=0 (needed by lambertian_closure and the a=0 replay); at m>0 a Lambertian surface layer is exactly zero (_zero_surface_layer!), so the final interaction is the identity on J₀⁺/J₀⁻ and only those two (NquadN, 1, nSpec) blocks (+ per-source slots) are kept. Footprint drops to 4·NquadN²·nSpec + m_max·NquadN·nSpec (Float64/NquadN=60/nSpec=10⁴/m_max=20: 24 GB → ≈1.2 GB). rt_run_atmosphere(...; cache_mode = :auto)(the default) picks:slimiff every target BRDF hascomponent_m_max(b, ctx) == 0(the Lambertian family), else:full.rt_run_surfacethrowsArgumentErrorif asked to replay a non-Lambertian-family BRDF against a:slimcache — rebuild withcache_mode = :full.

Type parameters

  • FT — working float type (Float64 / Float32).

  • AT3 — concrete 3-D array type for R/T/J snapshots (Array{FT,3} on CPU, CuArray{FT,3} on GPU).

  • AT1 — concrete 1-D array type for τ_sum_surf (device-resident, same as the τ_sum_end_dev array the monolithic surface step consumes, so the replay preserves rt_run's host/device contract).

  • AT2 — concrete 2-D array type for surface_F₀ ((pol_n, nSpec)).

  • AT — the array-type constructor (Array / CuArray); applied as arr_type(x) when replaying against a different architecture than the one the cache was built on (not currently exercised — CPU→CPU / GPU→GPU only).

  • PT — polarization type (e.g. Stokes_IQU{FT}).

  • QT — QuadPoints{FT} instance type.

  • ARCH — architecture (CPU / GPU).

  • RS — Raman type; always noRS{FT} for this cache (the only supported type — see the module-level scope note in rt_run_split.jl).

  • PS — prepared_sources concrete type.

  • IS — I_static concrete type (Diagonal{FT, AT1}).

  • SCI <: AbstractScatteringInterface — concrete scattering-interface type used at the surface step (uniform across m for a given scene).

Every field is concrete so the hot replay loop in rt_run_surface stays type-stable end-to-end.

source

Analytic Lambertian Albedo Closure ​

vSmartMOM.CoreRT.LambertianClosure Type
julia
LambertianClosure{FT}

Analytic closure for the TOA Lambertian-surface reflectance as a function of albedo a, built once from an AtmosphereRTCache's m=0 blocks. Evaluating c(a) for a new albedo (scalar or per-spectral-point vector) is an O(n_vza · pol_n · nSpec) broadcast — no surface build, no matrix inversion, no Fourier loop.

Fields

  • R_black :: Array{FT,3} — TOA reflectance at a = 0 ((n_vza, pol_n, nSpec)), i.e. rt_run_surface(cache, LambertianSurfaceScalar(0))[1]. Includes every Fourier moment (only m=0 depends on a; m>0 is albedo-independent and folded in here). When the cache was built with TMSCorrection, the surface-independent exact single-scattering term is already folded in here too (via rt_run_surface), so every albedo of the sweep carries the TMS correction at no extra cost.

  • t_up :: Array{FT,3} — T⁻⁻u restricted to the view-direction rows and pre-weighted by the exact m=0 azimuthal/Stokes factor that postprocessing_vza! applies (w₀ · Diagonal([1,1,0,0][1:pol_n]), i.e. cos(0)=1 on I/Q, sin(0)=0 on U/V) — so evaluating the closure is a plain broadcast against this array, not a matrix multiply.

  • S̄ :: Vector{FT} — atmosphere's spherical albedo viewed from below, 2·vᵀR⁺⁻u, one value per spectral point.

  • E_dw :: Vector{FT} — total (diffuse + direct) downward flux at BOA, vᵀJ₀⁺ + μ₀F₀ᴵe^(−τ_surf/μ₀).

  • w₀ :: FT — the m=0 Fourier weight (0.5/π); stored for reference — already folded into t_up.

See lambertian_closure (constructor), the functor (c::LambertianClosure)(a), albedo_jacobian, and invert_albedo.

source
vSmartMOM.CoreRT.lambertian_closure Function
julia
lambertian_closure(cache::AtmosphereRTCache) -> LambertianClosure

Build the analytic Lambertian-albedo closure (proposal §3.1) from an AtmosphereRTCache. Works on both :full and :slim caches (§4) — only the m=0 blocks are read, and both cache modes always keep those in full.

Solar-only: throws ArgumentError if the cache carries per-source J₀ slots (e.g. thermal emission, which declares a source_key and therefore gets its own composite slot) — the closure term for those sources is a follow-up (§3 notes this is analogous to "no closure exists in general" for non-Lambertian BRDFs). Also throws for an active SurfaceSIF source: SIF is injected directly into the surface's legacy j₀⁻ slot (not a per-source composite slot, so the check above wouldn't catch it), but it still bounces between surface and atmosphere in an albedo-dependent way that (†) does not account for.

GPU note: the cache's device-resident arrays are read exactly once here (Array(...)-converted up front), so the closure itself is built on the host regardless of the cache's architecture; the returned LambertianClosure is always host-resident, and evaluating it is plain CPU broadcasting.

Like the cache it is built from, the closure describes the scene as of rt_run_atmosphere — any subsequent update_model! / update_aerosol_*! scene change requires a fresh cache and a fresh closure.

source
vSmartMOM.CoreRT.albedo_jacobian Function
julia
albedo_jacobian(c::LambertianClosure, a) -> Array{FT,3}

Exact analytic ∂R_SFI/∂a at albedo a (same shape/broadcast rules as the functor): w₀ · 2·E_dw/(1−a·S̄)² · t_up (proposal §3.2). Differentiating (†) directly — R_black doesn't depend on a, so it drops out.

source
vSmartMOM.CoreRT.invert_albedo Function
julia
invert_albedo(c::LambertianClosure, R_meas; i_vza::Int = 1) -> Vector

Closed-form albedo retrieval from a measured TOA Stokes-I reflectance R_meas (same (n_vza, pol_n, nSpec) shape as the functor's output — e.g. from a real instrument, or from rt_run/rt_run_surface) at view index i_vza. No iteration, no linearized RT run: invert (†) for a per spectral point,

y = (R_meas − R_black) / t_upᴵ;  a = y / (2·E_dw + S̄·y)

where t_upᴵ = c.t_up[i_vza, 1, :] already carries the w₀/Stokes-I weighting.

When |t_upᴵ[s]| or |E_dw[s]| is at the machine-epsilon floor (the Stokes-I upward transmission or the downward flux effectively vanishes — an ill-conditioned spectral point, e.g. a fully opaque path), the inversion is undefined and a[s] = NaN (rather than an arbitrarily large or garbage value from dividing by ~0).

source

Scenario Sweeps ​

vSmartMOM.CoreRT.ScenarioSweep Type
julia
ScenarioSweep{FT}

Declarative cross product (SZA × view-pair × BRDF) over a single scene.

  • sza :: Vector{FT} — crossed

  • view_pairs :: Vector{Tuple{FT,FT}} — paired (vza, vaz). postprocessing_vza! produces one result per pair, not a Cartesian product.

  • brdfs :: Vector{<:AbstractSurfaceType} — crossed

Build with the keyword constructor ScenarioSweep(; sza, vza, vaz, brdfs), which accepts either paired (view_mode = :paired) or crossed (view_mode = :cross) (vza, vaz) inputs, or an explicit view_pairs list.

source
vSmartMOM.CoreRT.SweepResult Type
julia
SweepResult{FT}

Output of run_sweep. R and T are 5-D arrays of shape (N_sza, N_view_pair, N_brdf, pol_n, nSpec). The sweep field carries the input axes for downstream introspection.

source
vSmartMOM.CoreRT.run_sweep Function
julia
run_sweep(model::RTModel, params::vSmartMOM_Parameters, sweep::ScenarioSweep;
         i_band::Integer = 1, sources = nothing) -> SweepResult

Execute the (SZA × view-pair × BRDF) sweep declared by sweep over the scene carried by model/params. For each SZA, remake_geometry builds a cheap per-geometry RTModel — sharing atmosphere/optics/ solver/surfaces with model, no HITRAN/Mie re-run — then rt_run_multi_surface amortises the atmosphere phase (elemental → doubling → interaction over Nz layers) across every BRDF in sweep.brdfs within that SZA.

model and params must satisfy remake_geometry's contract: model was built via model_from_parameters(params) (or is otherwise geometry-consistent with it — same quadrature_type/l_trunc/ polarization_type/architecture). sweep.view_pairs supplies ONE shared (vza, vaz) list reused at every SZA — rt_run_multi_surface runs every view pair in a single model call (they are just entries in the model's vza/vaz vectors), so per-SZA cost is O(N_brdf) surface replays, not O(N_view_pair · N_brdf).

Returns a SweepResult with R/T shaped (N_sza, N_view_pair, N_brdf, pol_n, nSpec).

source
vSmartMOM.CoreRT.remake_geometry Function
julia
remake_geometry(model::RTModel, params::vSmartMOM_Parameters;
               sza = model.geometry.sza,
               vza = model.geometry.vza,
               vaz = model.geometry.vaz) -> RTModel

Rebuild an RTModel at a different observation geometry WITHOUT re-running the atmosphere build (model_from_parameters's HITRAN absorption, Mie aerosol optics, and per-band Fourier-support derivation) — the expensive part of model_from_parameters for a spectral-band scene.

Only geometry::ObsGeometry and quad_points::QuadPoints depend on (sza, vza, vaz): quad_points is built by rt_set_streams, which appends SZA/VZA cosines as (mostly) zero-weight output nodes onto a fixed set of weighted quadrature nodes. Everything else — atmosphere (profile, τ_abs, τ_rayl), optics (aerosol Mie/greek coefs), solver (m_max_bands, truncation) and surfaces — is geometry-independent, so this function shares them by reference with model and only constructs fresh ObsGeometry/QuadPoints objects, THE SAME WAY model_from_parameters does (identical rt_set_streams call). The observer heights resolved at build time (interior sensor levels, TOA/BOA outputs) and the external_solar representation are kept from model.

params must be the SAME parameters object (or an unmutated deepcopy) used to build model in the first place — it supplies the quadrature scheme, l_trunc, polarization type and architecture that model.solver was derived from (params.quadrature_type, params.l_trunc, params.polarization_type, params.architecture); this function does NOT re-derive solver.m_max_bands from these, it reuses model.solver verbatim (see the bit-exactness note below for why that's safe). Mutate sza/vza/vaz via the keyword arguments here, not by mutating params — a params with anything OTHER than geometry changed (e.g. a different l_trunc) would silently desync model.solver from the quadrature this function builds.

Bit-exactness

rt_run(remake_geometry(model, params; sza=θ)) is bit-exact with rt_run(model_from_parameters(p′)) for p′ = deepcopy(params) with p′.sza = θ (and same vza/vaz) — see test/test_scenario_sweep.jl. This relies on quad_points.Nstreams (hence solver.m_max_bands, which is NOT recomputed here) being independent of sza/vza: rt_set_streams builds the weighted Gauss-Legendre/Radau nodes first (always nonzero weight) and only afterwards appends the SZA/VZA cosines as additional nodes, padded with zero weight; Nstreams = count(!iszero, wt_μ) therefore only changes if a SZA/VZA cosine happens to land exactly on an existing weighted node in floating point (it would then be dropped by unique instead of appended) — a pre-existing knife-edge property of rt_set_streams itself (shared with the monolithic model_from_parameters path), not something this function introduces. quad_points.Nquad (the total node count, including zero-weight SZA/VZA duplicates) MAY legitimately differ between geometries — that's expected and harmless, since rt_run sizes its per-call NquadN from model.quad_points.Nquad fresh on every call, not from anything cached at model-build time.

See docs/dev_notes/proposals/surface_split_albedo_sweep.md §6/§7 (PR 3) for the design context — this is the seam run_sweep uses to amortise the per-SZA HITRAN/Mie rebuild across a scenario sweep.

source

Stream-Level RT Output ​

vSmartMOM.CoreRT.StreamRTResult Type
julia
StreamRTResult{FT}

Per-Fourier-moment radiative-transfer result at all internal quadrature streams — the output of rt_run_streams. Downstream consumers (e.g. ExoOptics disk integration) reconstruct the Stokes vector at arbitrary (μ_v, μ_0, Δφ) by interpolating over (μ_v, μ_0) and Fourier-summing over m with cos(m·Δφ) / sin(m·Δφ) weights, instead of one rt_run per pixel.

Fields

  • qp_μ :: Vector{FT} — full quadrature mu nodes (length Nquad).

  • iμ₀ :: Int — index into qp_μ of the SZA stream (the closest to μ₀).

  • μ₀ :: FT — cosine of the chosen incident-beam SZA.

  • pol_n :: Int — number of Stokes components (pol_type.n).

  • weight :: Vector{FT} — Fourier weight (0.5/π for m=0; 1/π otherwise) used for each moment by the internal post-processor.

  • R⁻⁺_per_m :: Vector{Array{FT, 3}} — (NquadN, NquadN, nSpec) per Fourier moment, where NquadN = Nquad · pol_n. This is the full-stream reflection matrix (stokes-out blocks × stokes-in blocks).

  • J⁻_per_m :: Vector{Array{FT, 3}} — (NquadN, 1, nSpec) per Fourier moment, the combined per-source SFI Stokes at all output streams (legacy slot + all per-source slots summed). Use this when you want the SFI-style output without re-running RT.

  • J⁺_per_m :: Vector{Array{FT, 3}} — analog for BOA transmittance.

Recovering rt_run's output

julia
streams = rt_run_streams(model; sources = sources)
# At a single (vza, vaz):
iμ = nearest_point(streams.qp_μ, cosd(vza))
_, istart, iend = get_indices(iμ, pol_type)
R_recovered = zero(rt_run(model; sources)[1][1, :, :])  # (n_stokes, nSpec)
for (mi, m) in enumerate(0:length(streams.J⁻_per_m)-1)
    cosmφ = cosd(m * vaz); sinmφ = sind(m * vaz)
    w_stokes = pol_n == 1 ? streams.weight[mi] * cosmφ :
               streams.weight[mi] *
               Diagonal([cosmφ, cosmφ, sinmφ, sinmφ][1:pol_n])
    for s in 1:nSpec
        R_recovered[:, s] .+= w_stokes * streams.J⁻_per_m[mi][istart:iend, 1, s]
    end
end

A unit test pins the bit-exact agreement of this reconstruction against rt_run for one published-figure geometry.

Two rt_run terms are applied after the Fourier sum at the requested views and therefore cannot appear in per-moment streams:

  • TMSCorrection (exact single scattering): rt_run_streams rejects it with an ArgumentError.

  • the Cox–Munk exact direct-solar glint residual (apply_ss_correction!): the streams of a CoxMunkSurface model hold the Fourier-truncated glint only, so a reconstruction differs from rt_run near the specular direction. Add the residual for each reconstructed geometry if it is needed.

source
vSmartMOM.CoreRT.rt_run_streams Method
julia
rt_run_streams(model; i_band=1, sources=nothing) -> StreamRTResult

Run the RT solver and return per-Fourier-moment Stokes matrices at all quadrature streams instead of post-processed (vza, vaz) outputs. See StreamRTResult for the data layout and a worked recovery example.

Internally calls the same production full-column kernel as rt_run with a callback that copies composite_layer.R⁻⁺, composite_layer.J₀⁺/⁻, and the combined per-source-slot SFI contributions out of the live layer accumulators once per Fourier moment.

Stream output is intrinsically full-column. geometry.obs_alt still affects the atmospheric grid (interior heights are exact interfaces), but it does not filter the TOA/BOA stream matrices stored in StreamRTResult.

TMSCorrection models are rejected, and the Cox–Munk glint residual is not included; see StreamRTResult.

source

Jacobian Parameter Layout ​

vSmartMOM.CoreRT.AbstractParameterLayout Type

Common interface for full and retrieval-selected Jacobian layouts.

source
vSmartMOM.CoreRT.ParameterLayout Type
julia
ParameterLayout

Describes the ordering of physical parameters in the Jacobian derivative dimension.

Instead of hardcoding 1 + 7*NAer + NGas + NSurf throughout the codebase, all index arithmetic goes through this struct. Each aerosol carries aerosol_params sub-parameters (currently 7: τ_ref, nᵣ, nᵢ, μ_logr, σ_logr, profile location, profile width), followed by the layer-resolved gas VMR columns and surface parameters. The native size columns differentiate LogNormal(μ_logr, σ_logr): the YAML median radius and geometric width are exp(μ_logr) and exp(σ_logr). Convert to derivatives with respect to those YAML values by dividing the native columns by the median radius and geometric width, respectively. Gas columns are flattened species-major: all TOA-to-BOA layers for gas 1, then all layers for gas 2, and so on. Each is the native partial ∂τ/∂VMR = σ·N_dry, plus ABSCO H₂O-broadener responses in the q-driven H₂O columns, as described by AbundanceTangentContract. Surface parameters are followed by the optional two-column SIF block [SIF760, mSIF] ([SIF755, slope] for the legacy SurfaceSIF interface), then canopy parameters. The final two columns are (p₀, σ_p) for a pressure-form Normal profile and (z₀, σ₀) for an altitude-form LogNormal profile.

Example

julia
layout = ParameterLayout(n_aerosols=1, n_gases=6, n_surface=1) # 2 gases × 3 layers
psurf_index(layout)        # 1
aerosol_range(layout, 1)   # 2:8
gas_profile_range(layout, 1, 3) # 9:11
gas_layer_index(layout, 2, 3, 3) # 14
surface_range(layout)      # 15:15
n_total(layout)            # 15
source
vSmartMOM.CoreRT.ParameterKey Type
julia
ParameterKey

Stable, retrieval-independent identity for one physical Jacobian parameter. band == 0 denotes a parameter shared by all bands and layer == 0 denotes a parameter without a vertical-layer index. Retrieval flavours use these keys to compile small band-local tangent spaces into one global state vector.

kind identifies the parameter family; component identifies an aerosol, gas, surface, or source within that family; field identifies the physical quantity; and layer/band locate it vertically and spectrally. The key is metadata only: unit and coordinate transformations remain the responsibility of the retrieval boundary.

source
vSmartMOM.CoreRT.ActiveParameterLayout Type
julia
ActiveParameterLayout

Compact parameter layout used by one spectral band of a selective linearized calculation. keys are ordered exactly as the trailing tangent dimension propagated by the RT kernels. local_to_global scatters those columns into global_keys; native_layer_columns selects the corresponding columns from the historical full physical-optics layout before the expensive adding-doubling propagation begins.

Layer parameters must precede surface/source parameters. Native layer columns must be positive, unique, sorted, and remain in the historical component order (pressure, aerosols, gas); local_to_global supplies any desired retrieval reordering. Surface and SIF columns must each be contiguous because their analytic seed routines accept unit ranges.

source
vSmartMOM.CoreRT.n_total Function

Total number of retrieval parameters.

source
vSmartMOM.CoreRT.aerosol_range Function

Index range for aerosol iaer (1-based).

source
vSmartMOM.CoreRT.gas_range Function

Index range for all gas VMR parameters.

source
vSmartMOM.CoreRT.gas_profile_range Function

Index range for all vertical layers of gas species igas (species-major ordering).

source
vSmartMOM.CoreRT.gas_layer_index Function

Jacobian index for gas species igas in atmospheric layer iz.

source
vSmartMOM.CoreRT.surface_range Function

Index range for surface parameters.

source
vSmartMOM.CoreRT.surface_index Function

Index of a specific surface parameter (1-based within the surface block).

source
vSmartMOM.CoreRT.surface_parameter_count Function

Number of retrieval parameters contributed by a surface model.

source
vSmartMOM.CoreRT.supports_linearized_surface Function
julia
supports_linearized_surface(surface) -> Bool

Whether the linearized driver can build the surface layer and its tangent (a linearized create_surface_layer! method). Checked at the linearized solve boundary so unsupported surfaces fail with an ArgumentError naming the surface, not a deep MethodError. A surface that adds a linearized layer builder should also specialize this trait. The default is false.

source
vSmartMOM.CoreRT.sif_range Function

Index range for SIF parameters: per-wavenumber radiance at 760 nm, then dSIF/dν.

source
vSmartMOM.CoreRT.atmosphere_range Function

Atmospheric-state range; column 1 is surface pressure.

source
vSmartMOM.CoreRT.n_layer_params Function

Number of layer-level parameters (aerosol + gas, excluding surface and canopy).

source
vSmartMOM.CoreRT.canopy_range Function

Index range for canopy parameters (LAI, leaf_R, leaf_T, ...).

source
vSmartMOM.sif_reference_state Function
julia
sif_reference_state(; total_sif=0.5, reference_wavelength_nm=760,
                    column=:SIF_OLD, path=sif_data_path("sif-spectra.csv"))

Normalize the supplied SIF shape so that its wavelength-space integral is total_sif (mW m⁻² sr⁻¹), convert it to spectral radiance per cm⁻¹, and return the local wavenumber-linear state at reference_wavelength_nm:

(SIF760=SIF(ν₇₆₀), mSIF=dSIF/dν|₇₆₀, ν_ref=ν₇₆₀, ν, spectrum).

The SIF760 property retains that name for the retrieval convention even if a different diagnostic reference wavelength is requested. The slope is a centered secant through the nearest tabulated samples bracketing the reference.

source

Retrieval-Selected Jacobian Plans ​

The policy types below compile a retrieval definition into a small band-local tangent space before the layer optical properties and MOM operators are allocated. See Compute Jacobians for the OCO_RRS_synth state ordering, units, and chain rules.

vSmartMOM.CoreRT.AbstractJacobianFlavor Type

Base type for retrieval-specific Jacobian selection policies.

source
vSmartMOM.CoreRT.OCO_RRS_synth Type
julia
OCO_RRS_synth()

Selective Jacobian flavour for the synthetic three-band OCO RRS/XCO2 retrieval. It selects the native physical columns corresponding to surface pressure, the non-fixed CO2 layers, aerosol reference optical depth and profile location, three Legendre surface coefficients per band, and the two O2-A-band SIF parameters. Conversion from native aerosol tau_ref/z0 to the retrieval's log(AOD760)/log(z0) coordinates is intentionally applied at the retrieval boundary. Aerosol microphysics and profile widths, H2O, and the fixed upper CO2 layers remain forward-model inputs only.

source
vSmartMOM.CoreRT.JacobianPlan Type
julia
JacobianPlan

A compiled global named physical basis plus one compact ActiveParameterLayout per spectral band. global_keys determine cross-band identity and order; bands[i] determines the derivative dimension actually allocated and propagated by band i. Retrieval-specific policy is fully compiled into this object before the RT kernels are entered.

source
vSmartMOM.CoreRT.PlannedRTModelLin Type
julia
PlannedRTModelLin

Backward-compatible wrapper pairing the ordinary physical-derivative model with a compiled JacobianPlan. Unknown properties are forwarded to base, so diagnostic code written for RTModelLin continues to work.

source
vSmartMOM.CoreRT.jacobian_plan Function
julia
jacobian_plan(flavor, params, model, lin_model)

Compile a retrieval-specific Jacobian flavour into a global named state and compact band-local tangent layouts. New retrievals extend this function by multiple dispatch; the radiative-transfer kernels consume only the resulting generic JacobianPlan.

source
vSmartMOM.CoreRT.requires_aerosol_microphysics_jacobians Function
julia
requires_aerosol_microphysics_jacobians(flavor)

Upstream-work trait. Return false only when refractive index and size distribution are fixed: forward Mie optics are still computed, but their ForwardDiff tangents and Fourier phase derivatives are omitted.

source
vSmartMOM.CoreRT.requires_h2o_jacobians Function
julia
requires_h2o_jacobians(flavor)

Upstream-work trait controlling the q-driven H2O tangent block. Returning false retains the complete forward H2O absorption while leaving its native gas derivative columns zero and therefore eligible for omission by the plan.

source
vSmartMOM.CoreRT.requires_pressure_jacobians Function
julia
requires_pressure_jacobians(flavor)

Upstream-work trait. Return false when surface pressure is fixed. All forward absorption is retained, including absorbers without a pressure derivative implementation. A plan may not select pressure when disabled.

source
vSmartMOM.CoreRT.globalize_jacobian Function
julia
globalize_jacobian(jacobian, layout)

Scatter a band-local Jacobian into the plan's global named physical basis. The trailing dimension is replaced by the global dimension, inactive parameters are initialized to exact zero, and active columns are assigned through layout.local_to_global. This operation reorders/zero-fills only; it does not apply unit, AOD-reference, or logarithmic-coordinate chain rules.

source

Surface-Pressure Vertical-Grid Conventions ​

How the vertical grid responds to a change in surface pressure, both when a retrieval rebuilds the profile at a trial pressure and in the analytic surface-pressure Jacobian (model_from_parameters(LinMode(), params; vertical_grid=...)). FixedPressureGrid() is the default.

vSmartMOM.CoreRT.AbstractVerticalGridConvention Type
julia
AbstractVerticalGridConvention

How the model's vertical grid responds to a change in surface pressure p_surf while the layer temperatures, specific humidities, and trace-gas mixing ratios are held fixed.

Concrete conventions:

  • FixedPressureGrid (default): every pressure interface above the surface is fixed; only the bottom boundary p_half[end] = p_surf moves.

  • FixedAltitudeGrid: every interface is fixed in geometric altitude above the surface; a change in p_surf rescales every interface pressure.

A convention enters in two places, which must agree for the analytic Jacobian to be the derivative of the forward model actually evaluated:

  1. surface_pressure_interfaces rebuilds the interface pressures for a trial p_surf;

  2. model_from_parameters(LinMode(), params; vertical_grid=...) selects the matching surface-pressure tangents (psurf_profile_tangents).

source
vSmartMOM.CoreRT.FixedPressureGrid Type
julia
FixedPressureGrid()

Default convention. The interface pressures p_half[1:end-1] are fixed and only the bottom boundary p_half[end] = p_surf moves. Consequently only the bottom layer changes: its dry and H₂O columns, its geometric thickness, and its mid-layer pressure p_full[end], with dp_full[end]/dp_surf = 1/2. Geometric altitudes above the bottom layer shift with the bottom-layer thickness.

The tangent describes the final model grid. With profile_reduction > 0 the reduced interfaces are spaced from the top pressure to p_surf, so a rebuilt reduced grid moves more than its bottom interface; the bottom-only tangent is then the derivative of a different (fixed-final-grid) forward model.

source
vSmartMOM.CoreRT.FixedAltitudeGrid Type
julia
FixedAltitudeGrid()

Every interface is fixed in geometric altitude above the surface, with the layer temperatures T, specific humidities q, and mixing ratios fixed on those altitude layers. This is the natural convention for a fixed-height observer (flux tower, aircraft at a fixed altitude), around which the layering must not move when p_surf changes.

vSmartMOM's layer hydrostatics (compute_atmos_profile_fields) give each layer the thickness

Δzᵢ = (R Tᵢ / (g₀ Mᵢ)) · ln(p_half[i+1] / p_half[i]),

with constant g₀ and a moist-air molar mass Mᵢ fixed by qᵢ. Holding every Δzᵢ, Tᵢ, and Mᵢ fixed therefore fixes every pressure ratio p_half[i+1]/p_half[i], and the whole grid scales with the surface pressure:

p_half[i](p_surf) = p_half_ref[i] · p_surf / p_surf_ref.

This is exact for the model's isothermal-layer, constant-gravity hydrostatics; it is not a first-order approximation. It follows that, for every layer i,

dp_half[i]/dp_surf = p_half[i]/p_surf,   dp_full[i]/dp_surf = p_full[i]/p_surf,
dN_dry,i/dp_surf   = N_dry,i/p_surf,     dN_H₂O,i/dp_surf   = N_H₂O,i/p_surf,
dΔzᵢ/dp_surf       = 0,

because each hydrostatic column is proportional to its pressure thickness. Every layer's cross sections therefore respond to p_surf, not only the bottom layer's.

Aerosol profiles defined in geometric altitude above the surface (LogNormal(log(z₀), σ₀)) keep their layer allocation, because no interface moves in altitude. Aerosol profiles defined in pressure (Normal(p₀, σp)) keep p₀ and σp fixed in pressure and are redistributed as the layer pressures scale.

source
vSmartMOM.CoreRT.surface_pressure_interfaces Function
julia
surface_pressure_interfaces(convention, p_half, p_surf)

Return the interface pressures [hPa, TOA to BOA] for a trial surface pressure p_surf, starting from the reference interfaces p_half (whose last entry is the reference surface pressure). The input is not modified.

At p_surf == p_half[end] both conventions return an exact copy of p_half, so the reference atmosphere is identical bit for bit.

source

Quadrature Types ​

vSmartMOM.CoreRT.AbstractQuadratureType Type

Quadrature Types for RT streams

source
vSmartMOM.CoreRT.RadauQuad Type
julia
struct RadauQuad

Use Gauss Radau quadrature scheme, which includes the SZA as point (see Sanghavi vSmartMOM)

source
vSmartMOM.CoreRT.GaussLegQuad Type
julia
struct GaussLegQuad

Half-space Gauss-Legendre quadrature on [0, 1] for the upper hemisphere. The default for plane-parallel RT.

source

Quadrature Construction ​

vSmartMOM.CoreRT.rt_set_streams Function
julia
rt_set_streams(::GaussLegQuad, Ltrunc, obs_geom, pol_type, arr_type;
               external_solar=false)

Half-space Gauss-Legendre quadrature on [0, 1] for the upper hemisphere. Builds (Ltrunc + 2) ÷ 2 weighted nodes (Sanghavi: the +2 avoids collapsing to zero streams when Ltrunc = 0). All VZAs are appended as zero-weight output nodes. The default (external_solar=false) also appends SZA to the square operator grid — the historical representation.

With external_solar=true, the collimated solar direction is instead a scalar source geometry plus the exact Zₘ(μᵢ,μ₀) column on phase_qp_μ; it is not a diffuse stream. In this mode iμ₀ == iμ₀Nstart == 0 are deliberate sentinels, and source code must use iμ₀_phase. For example, five weighted streams plus one distinct VZA with Stokes_IQU gives Nquad=6 and an 18×18 diffuse operator, rather than the legacy 21×21 operator that also embeds SZA.

source

rt_set_streams(::RadauQuad, Ltrunc::Int, obs_geom::ObsGeometry, pol_type, arr_type; external_solar=false)

Computes hemispheric quadrature points with Gauss-Radau method in two intervals between [0,SZA] and [SZA,1] to include SZA as full quadrature point for DNI SZA included as full weighted node, VZAs included with 0 weight external_solar=true is rejected because Radau's defining construction embeds the solar direction as a weighted ordinate. Returns computed quadrature points as QuadPoints

source
julia
rt_set_streams(::GaussLegQuad, obs_geom, pol_type, arr_type;
               nstreams::Int, external_solar=false)

Build a half-space Gauss-Legendre quadrature with exactly nstreams weighted nodes per hemisphere. Equivalent to the positional form with Ltrunc = 2·nstreams - 2 (so (Ltrunc+2)÷2 == nstreams matches the inverse of Sanghavi's even-Ltrunc formula).

By default, the direct solar direction remains outside the diffuse operator while retaining its exact phase-matrix source column. Pass external_solar=false explicitly for the legacy embedded-μ₀ representation.

source
julia
rt_set_streams(::RadauQuad, obs_geom, pol_type, arr_type;
               nstreams::Int, external_solar=false)

Build a Gauss-Radau quadrature whose public-contract stream_l_cap is 2·nstreams - 1. The actual weighted-stream count qp.Nstreams may exceed nstreams in the typical SZA-not-on-node branch (doubled around μ₀); the public contract stream_l_cap = 2·nstreams - 1 still holds. Per Sanghavi (2026-05-07), RadauQuad is supported but not recommended — GaussLegQuad is the default and the cheaper choice for the same resolving power. external_solar=true is unsupported; use GaussLegQuad for that representation.

source