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
FwdMode()Marker selecting the forward radiative-transfer model construction path.
sourcevSmartMOM.LinMode Type
LinMode()Marker selecting the linearized model construction path used for Jacobians.
sourcevSmartMOM.CoreRT.model_from_parameters_lin Function
model_from_parameters_lin(params; kwargs...)Convenience alias for model_from_parameters(LinMode(), params; kwargs...). Returns (model, lin_model).
vSmartMOM.CoreRT.rt_run_lin Function
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).
Model Types
vSmartMOM.CoreRT.RTModel Type
RTModel{ARCH, FT} <: AbstractRTModel{ARCH, FT}Central model state for vSmartMOM radiative transfer.
Only two type parameters:
ARCH: compute architecture (CPUorGPU)FT: floating-point precision (Float32orFloat64)
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 weightsatmosphere: 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 configurationnumerics: Numerical-knob parameters (doubling threshold, etc.) — seeRTNumericalParametersgeometry: Observation geometryquad_points: Quadrature points and weightsatmosphere: Atmospheric state (profile + spectral bands)optics: Precomputed optical propertiessurfaces: Surface models, one per spectral bandsources: Source-term scene (v0.6 source-term refactor). Defaults to a singleSolarBeam— i.e. the historical unit-Stokes-I direct solar illumination — and can be overridden atrt_runcall time via thesources=kwarg. Mirrors the abstract-typedsurfacesslot.
vSmartMOM.CoreRT.SolverConfig Type
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 runsm = 0:m_max_bands[iBand]. Equalsn_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 indexl_trunc: Legendre truncation order (user-specified)Δ_angle: Exclusion angle for forward peak [deg]depol: Depolarization factoruse_component_traits: Phase C flag: whentrue(current default), per-bandm_max_bandsare derived fromcomponent_m_max(c, ctx)traits across active components — Cox-Munk / RPV / RossLi / canopy run to their fulluser_l_capinstead of the historical half-truncated aggregator. Flip tofalseto fall back to the legacymin(ceil((l_max+1)/2), params.max_m)aggregator (bit-equal to Phase B).
vSmartMOM.CoreRT.Atmosphere Type
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
vSmartMOM.CoreRT.RayleighScattering Type
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
vSmartMOM.CoreRT.AerosolState Type
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.
vSmartMOM.CoreRT.Optics Type
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 propertiesaerosols: Aerosol scattering optics and optical depthsτ_abs: Absorption optical depth: τ_abs[iBand][nSpec × nLayers]τ_rayl: Rayleigh optical depth: τ_rayl[iBand][nSpec × nLayers]
vSmartMOM.CoreRT.OpticsLin Type
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
vSmartMOM.CoreRT.vSmartMOM_Parameters Type
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 bandbrdf: 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 — seetruncation)l_trunc: Truncation length for legendre terms (scalar for now, can donBandlater)truncation: Phase-function truncation method. Defaults toδBGE(l_trunc, Δ_angle)for backward compatibility; set explicitly toNoTruncation()for canopy-only or smooth-phase-function runs (Sanghavi & Stephens 2015 §2 — thef_tr → 0limit). TheΔ_anglefield above is then ignored.depol: Depolarization factornumerics: Numerical-knob parameters (doubling threshold, etc.) — defaults toRTNumericalParameters{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 semanticsT: Temperature Profile [K]p: Pressure Profile [hPa]q: Specific humidity profileprofile_reduction_n: Length of profile reductionabsorption_params: Optional struct that holds all absorption-related parametersscattering_params: Optional struct that holds all aerosol scattering-related parametersnstreams: Phase D — primary user-facing resolution knob: weighted streams per hemisphere. Public contract:stream_l_cap = 2·nstreams - 1. Whennothing, the parser populatedmax_m/l_truncfrom a legacy YAML; new-schema configs set this and may omitmax_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 ofstream_l_cap.stream_l_cap: Phase D — derived projection cap:2·nstreams - 1for new schema, orl_truncfor legacy. Internally consumed by the trait aggregator and the truncation-resolver. Always populated; defaults tol_truncuntilnstreamsis explicitly set.legacy_l_cap_override: Phase D — legacyl_truncvalue retained verbatim from YAML when the user explicitly set it.nothingfor new-schema configs that derivestream_l_capfromnstreamsinstead.greek_beta_cutoff: Optional post-Mie Fourier-support cutoff based specifically on the magnitude of the integrated aerosol Greek coefficientβ_l. The largest degree satisfyingabs(β_l) >= greek_beta_cutoffis 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.nothingdisables the cutoff.rayleigh_anisotropy: Wavelength dependence of the N₂/O₂ polarizability tensor used for the automatic (depol < 0) Rayleigh depolarization ofgreek_rayleighandτ_rayl:DispersiveAnisotropy()(default) orConstantAnisotropy()(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 keyradiative_transfer.rayleigh_anisotropy.
vSmartMOM.CoreRT.AtmosphericProfile Type
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 Profilep_full: Pressure Profile (Full)q: Specific humidity profilep_half: Pressure Levelsvmr_h2o: H2O/dry-air molar-ratio profilevcd_dry: Vertical Column Density (Dry)vcd_h2o: Vertical Column Density (H2O)vmr: Volume Mixing Ratio of Constituent GasesΔz: Layer height (meters)
vSmartMOM.CoreRT.ObsGeometry Type
Observation geometry, including the resolved vertical-output interfaces.
sourcevSmartMOM.CoreRT.QuadPoints Type
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 externaliμ₀Nstart: First Stokes index of μ₀ in qp_μN, or zero for external-solar modelsqp_μ: Diffuse-operator quadrature pointswt_μ: Weights of diffuse-operator quadrature pointsqp_μ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 μ₀ inphase_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.
vSmartMOM.CoreRT.CompositeLayer Type
CompositeLayer{FT} <: AbstractLayerAccumulated (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 (+) directionT⁻⁻::AbstractArray{FT,3}: transmission in the outgoing (−) directionJ₀⁺::AbstractArray{FT,3}: source in the incoming (+) directionJ₀⁻::AbstractArray{FT,3}: source in the outgoing (−) direction
vSmartMOM.CoreRT.AddedLayer Type
AddedLayer{FT} <: AbstractLayerSingle 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 (+) directionr⁺⁻::AbstractArray{FT,3}: reflectance from outgoing (−) to incoming (+)t⁻⁻::AbstractArray{FT,3}: transmission in the outgoing (−) directionj₀⁺::AbstractArray{FT,3}: source in the incoming (+) directionj₀⁻::AbstractArray{FT,3}: source in the outgoing (−) directiontemp1,temp2: pre-allocated workspace arrays to avoid allocationstemp1_ptr,temp2_ptr: CUDA pointer arrays (ignored on CPU)
vSmartMOM.CoreRT.RTNumericalParameters Type
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 depthdτ_max = threshold · μ_min, whereμ_minis 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). Default0.001is conservative legacy; raising to~0.01reduces doubling iterations.dτ_min_floor: Absolute floor ondτ_max(and hence ondτ_initial), expressed as a multiple ofeps(FT). Prevents the doubling discretisation from collapsing below FT precision regardless of geometry. Default1024·eps(FT)≈ 1.2e-4 for Float32 / 2.3e-13 for Float64 — the F32 value capsndoublaround 13–14 for τ_total≤1 (well within Float32 representability); the F64 value is far below any sensiblethreshold·μ_minso it never activates in practice.blas_threads: BLAS thread cap applied at everyrt_runinvocation.nothingmeans leave the current BLAS thread setting alone. Why cap: the batched-GEMM call sites incpu_batched.jland 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;1is 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 thenumerics.blas_threadsYAML key, or programmatically withRTNumericalParameters{FT}(blas_threads = 8).verbose: Verbose output flag. Whenfalse(default)rt_run/rt_run_linsuppress the per-callprint_timer()dump at the end of the run. Settrue(YAML keynumerics.verbose: trueor programmaticallyRTNumericalParameters{FT}(verbose = true)) to bring back the timing tree — useful for profiling, noisy for production loops. Per-bandComputing profile for X...andAOD at band X...messages were demoted to@debugand are silent unless you setENV["JULIA_DEBUG"] = "vSmartMOM".fourier_convergence: Azimuthal Fourier-loop termination strategy. The package-wide default runs all moments. Retrieval/truth workflows may explicitly selectIntensityConvergenceorStokesConvergenceafter validating the omitted tail for their viewing geometry.ss_correction: Single-scattering (first-order) truncation correction — seeAbstractSingleScatteringCorrection. DefaultNoSSCorrection(bit-identical historical behavior). SetTMSCorrectionfor the Nakajima–Tanaka correction; YAML keynumerics.ss_correction: tms.
vSmartMOM.CoreRT.AbstractFourierConvergence Type
AbstractFourierConvergenceTermination 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.
vSmartMOM.CoreRT.AllFourierMoments Type
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.
vSmartMOM.CoreRT.IntensityConvergence Type
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.
vSmartMOM.CoreRT.StokesConvergence Type
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.
Single-Scattering Correction
vSmartMOM.CoreRT.AbstractSingleScatteringCorrection Type
AbstractSingleScatteringCorrectionStrategy 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).
sourcevSmartMOM.CoreRT.NoSSCorrection Type
NoSSCorrection()No single-scattering correction — the truncated Fourier solution is returned as-is. Bit-identical to the historical solver. Default.
sourcevSmartMOM.CoreRT.TMSCorrection Type
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.
vSmartMOM.CoreRT.rt_run_ss_exact Function
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.
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
component_m_max(component, ctx) -> IntMaximum 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.
vSmartMOM.CoreRT._aggregate_m_max Function
_aggregate_m_max(active_components, ctx) -> IntAggregate 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.
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
_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 firstSolarBeaminsources, a user presetRS_type.F₀(nonzero, shaped(nStokes, nSpec)), thenmodel.sources; it is routed into theRS_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 explicitsources=applies to its call only.streams_callback: fires with the per-moment stream state (used byrt_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), andm_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.
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
ObserverRTResultHeight-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.
vSmartMOM.CoreRT.LevelRadiance Type
LevelRadianceRadiances 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.
vSmartMOM.CoreRT.ObserverRTResultLin Type
ObserverRTResultLinHeight-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).
vSmartMOM.CoreRT.LevelRadianceLin Type
LevelRadianceLinForward 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.
vSmartMOM.CoreRT.total_downwelling Function
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.
sourceReturn diffuse plus unscattered downwelling radiance at an interior level.
sourcevSmartMOM.CoreRT.total_downwelling_jacobian Function
Return the analytic Jacobian of total_downwelling.
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
BatchContextPre-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
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 ...
endWhat may change per scene
Temperature profile
T[K]Pressure half-levels
p_half[hPa]Specific humidity
q[kg/kg]Trace-gas VMRs via the
vmrkeyword dict (same keys as the originalparams.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: theRTModelthat is updated in place byupdate_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: cachedAtmosphericAbsorption.LineByLineModelobjects, one per(band, species)pair — eliminates HITRAN re-parsing per sceneh2o_models: cached H₂OLineByLineModelper band (nothingonly when the configuration has noabsorption_paramsat all; otherwise the model is cached unconditionally so a scene whoseqchanges 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 thatupdate_model!was last called with (initialised fromparamsat construction).nothingkeyword arguments toupdate_model!fall back to these stored values — never toparams— 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 frommodel.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λ_refand vertical distribution), initialised fromparamsand written by bothupdate_aerosol_loading!andupdate_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 bookmarksprofile_reduction_n: reduction target passed to the anchoredprepare_observer_profilepath;-1means no reduction
vSmartMOM.CoreRT.update_model! Function
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:
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], lengthNz(full levels, TOA to BOA). Ifprofile_reduction_n != -1, pass the unreduced profile length (same as the originalparams.T).p_half::AbstractVector: Half-level pressures [hPa], lengthNz + 1(orN_orig + 1before 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 eitherNzlayer-center values orNz + 1interface 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 toparams.absorption_params.vmrand rebuild the context instead). Values may be scalars, length-Nzlayer vectors, or length-Nz+1interface 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
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)
endvSmartMOM.CoreRT.update_aerosol_loading! Function
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. Whennothing, the value stored inparams.scattering_params.rt_aerosols[i_aer].τ_refis kept (only the distribution is updated).profile_dist::Distributions.Distribution: New vertical pressure distribution (aDistributions.jlDistributionobject — the same type used inRT_Aerosol.profile). Whennothing, 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,SolverConfigFourier bounds, gas absorption, or Rayleigh properties.
Example
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)vSmartMOM.CoreRT.update_aerosol_microphysics! Function
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 bandctx.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: TheBatchContextto update.i_aer: 1-based index of the aerosol species to replace.aerosol: NewAerosol(fromScattering.Aerosol), carrying the newsize_distribution,nᵣ, andnᵢ.
Keyword arguments
τ_ref::Real: New column optical depth atλ_ref. Whennothing, the current scene valuectx.current_τ_ref[i_aer]is kept (i.e. the most recently applied loading, not the originalparamsvalue).
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
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)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
rt_run_atmosphere(model::RTModel;
target_brdfs = nothing,
i_band::Integer = 1,
sources = nothing,
cache_mode::Symbol = :auto) -> AtmosphereRTCacheRun 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 viart_run_surface. The cache's Fourier loop bound is sized formax(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 tomodel.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 inrt_run.cache_mode::Symbol = :auto—:full(all six composite blocks at every Fourier moment),:slim(proposal §4: six blocks only atm=0;J₀⁺/J₀⁻/per-source slots only atm>0— valid only when everytarget_brdfsentry is Lambertian-family), or:auto(picks:slimiff every entry inbrdfs—target_brdfs, ormodel.surfaceswhen left at its default — hascomponent_m_max(b, ctx) == 0, else:full). See the memory-footprint note onAtmosphereRTCache.
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.
vSmartMOM.CoreRT.rt_run_surface Function
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.
vSmartMOM.CoreRT.rt_run_multi_surface Function
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.
vSmartMOM.CoreRT.AtmosphereRTCache Type
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_devarray the monolithic surface step consumes, so the replay preserves rt_run's host/device contract).AT2— concrete 2-D array type forsurface_F₀((pol_n, nSpec)).AT— the array-type constructor (Array/CuArray); applied asarr_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; alwaysnoRS{FT}for this cache (the only supported type — see the module-level scope note inrt_run_split.jl).PS—prepared_sourcesconcrete type.IS—I_staticconcrete type (Diagonal{FT, AT1}).SCI <: AbstractScatteringInterface— concrete scattering-interface type used at the surface step (uniform acrossmfor a given scene).
Every field is concrete so the hot replay loop in rt_run_surface stays type-stable end-to-end.
Analytic Lambertian Albedo Closure
vSmartMOM.CoreRT.LambertianClosure Type
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 ata = 0((n_vza, pol_n, nSpec)), i.e.rt_run_surface(cache, LambertianSurfaceScalar(0))[1]. Includes every Fourier moment (onlym=0depends ona;m>0is albedo-independent and folded in here). When the cache was built withTMSCorrection, the surface-independent exact single-scattering term is already folded in here too (viart_run_surface), so every albedo of the sweep carries the TMS correction at no extra cost.t_up :: Array{FT,3}—T⁻⁻urestricted to the view-direction rows and pre-weighted by the exactm=0azimuthal/Stokes factor thatpostprocessing_vza!applies (w₀ · Diagonal([1,1,0,0][1:pol_n]), i.e.cos(0)=1on I/Q,sin(0)=0on 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— them=0Fourier weight (0.5/π); stored for reference — already folded intot_up.
See lambertian_closure (constructor), the functor (c::LambertianClosure)(a), albedo_jacobian, and invert_albedo.
vSmartMOM.CoreRT.lambertian_closure Function
lambertian_closure(cache::AtmosphereRTCache) -> LambertianClosureBuild 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.
vSmartMOM.CoreRT.albedo_jacobian Function
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.
vSmartMOM.CoreRT.invert_albedo Function
invert_albedo(c::LambertianClosure, R_meas; i_vza::Int = 1) -> VectorClosed-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).
Scenario Sweeps
vSmartMOM.CoreRT.ScenarioSweep Type
ScenarioSweep{FT}Declarative cross product (SZA × view-pair × BRDF) over a single scene.
sza :: Vector{FT}— crossedview_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.
vSmartMOM.CoreRT.SweepResult Type
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.
vSmartMOM.CoreRT.run_sweep Function
run_sweep(model::RTModel, params::vSmartMOM_Parameters, sweep::ScenarioSweep;
i_band::Integer = 1, sources = nothing) -> SweepResultExecute 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).
vSmartMOM.CoreRT.remake_geometry Function
remake_geometry(model::RTModel, params::vSmartMOM_Parameters;
sza = model.geometry.sza,
vza = model.geometry.vza,
vaz = model.geometry.vaz) -> RTModelRebuild 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.
Stream-Level RT Output
vSmartMOM.CoreRT.StreamRTResult Type
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 (lengthNquad).iμ₀ :: Int— index intoqp_μ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, whereNquadN = 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
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
endA 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_streamsrejects it with anArgumentError.the Cox–Munk exact direct-solar glint residual (
apply_ss_correction!): the streams of aCoxMunkSurfacemodel hold the Fourier-truncated glint only, so a reconstruction differs fromrt_runnear the specular direction. Add the residual for each reconstructed geometry if it is needed.
vSmartMOM.CoreRT.rt_run_streams Method
rt_run_streams(model; i_band=1, sources=nothing) -> StreamRTResultRun 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.
Jacobian Parameter Layout
vSmartMOM.CoreRT.AbstractParameterLayout Type
Common interface for full and retrieval-selected Jacobian layouts.
sourcevSmartMOM.CoreRT.ParameterLayout Type
ParameterLayoutDescribes 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
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) # 15vSmartMOM.CoreRT.ParameterKey Type
ParameterKeyStable, 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.
vSmartMOM.CoreRT.ActiveParameterLayout Type
ActiveParameterLayoutCompact 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.
vSmartMOM.CoreRT.gas_profile_range Function
Index range for all vertical layers of gas species igas (species-major ordering).
vSmartMOM.CoreRT.gas_layer_index Function
Jacobian index for gas species igas in atmospheric layer iz.
vSmartMOM.CoreRT.surface_index Function
Index of a specific surface parameter (1-based within the surface block).
sourcevSmartMOM.CoreRT.surface_parameter_count Function
Number of retrieval parameters contributed by a surface model.
sourcevSmartMOM.CoreRT.supports_linearized_surface Function
supports_linearized_surface(surface) -> BoolWhether 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.
vSmartMOM.CoreRT.sif_range Function
Index range for SIF parameters: per-wavenumber radiance at 760 nm, then dSIF/dν.
sourcevSmartMOM.CoreRT.atmosphere_range Function
Atmospheric-state range; column 1 is surface pressure.
sourcevSmartMOM.CoreRT.n_layer_params Function
Number of layer-level parameters (aerosol + gas, excluding surface and canopy).
sourcevSmartMOM.CoreRT.canopy_range Function
Index range for canopy parameters (LAI, leaf_R, leaf_T, ...).
sourcevSmartMOM.sif_reference_state Function
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.
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.
sourcevSmartMOM.CoreRT.OCO_RRS_synth Type
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.
vSmartMOM.CoreRT.JacobianPlan Type
JacobianPlanA 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.
vSmartMOM.CoreRT.PlannedRTModelLin Type
PlannedRTModelLinBackward-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.
vSmartMOM.CoreRT.jacobian_plan Function
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.
vSmartMOM.CoreRT.requires_aerosol_microphysics_jacobians Function
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.
vSmartMOM.CoreRT.requires_h2o_jacobians Function
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.
vSmartMOM.CoreRT.requires_pressure_jacobians Function
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.
vSmartMOM.CoreRT.globalize_jacobian Function
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.
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
AbstractVerticalGridConventionHow 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 boundaryp_half[end] = p_surfmoves.FixedAltitudeGrid: every interface is fixed in geometric altitude above the surface; a change inp_surfrescales 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:
surface_pressure_interfacesrebuilds the interface pressures for a trialp_surf;model_from_parameters(LinMode(), params; vertical_grid=...)selects the matching surface-pressure tangents (psurf_profile_tangents).
vSmartMOM.CoreRT.FixedPressureGrid Type
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.
vSmartMOM.CoreRT.FixedAltitudeGrid Type
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.
vSmartMOM.CoreRT.surface_pressure_interfaces Function
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.
Quadrature Types
vSmartMOM.CoreRT.RadauQuad Type
struct RadauQuadUse Gauss Radau quadrature scheme, which includes the SZA as point (see Sanghavi vSmartMOM)
sourcevSmartMOM.CoreRT.GaussLegQuad Type
struct GaussLegQuadHalf-space Gauss-Legendre quadrature on [0, 1] for the upper hemisphere. The default for plane-parallel RT.
Quadrature Construction
vSmartMOM.CoreRT.rt_set_streams Function
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.
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
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.
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.