2 · Vector RTE & Discretization

For: readers building the equation-level mental model — atmospheric scientists, theorists, retrieval developers who want to understand what the solver actually solves.
Prev: 1 · The Problem & MOM Thesis · Next: 3 · Layer Optical Properties

This page writes the equation that vSmartMOM solves, then discretizes it. By the end, the diffuse per-layer problem reduces to square operators of shape (NquadN, NquadN, nSpec). An opt-in source-function path may evaluate the phase matrix on one additional, external solar direction without adding that direction to the diffuse operators.

The vector RTE

For a plane-parallel atmosphere illuminated by direct solar flux along , the Stokes vector for diffuse light obeys (Sanghavi 2014, Eq. 2):

Reading the right-hand side term by term:

— extinction along the path.
— thermal emission (zero in solar bands).
The single-scatter direct-beam source, exponentially attenuated from TOA.
The diffuse multiple-scatter integral over all incoming directions .

Sign convention: is downward, is upward. is the single-scattering albedo. is the 4×4 phase matrix (Sanghavi 2014, Eq. 8). The diffuse solution adds to the direct delta-source to give the total radiation field; we solve only for the diffuse part to avoid the singularity.

Boundary conditions

At the top of atmosphere (Sanghavi 2014, Eq. 6):

At the surface, an emissivity and a polarized BRDF close the system (Sanghavi 2014, Eqs. 33–37). Concepts/05 covers the surface in detail.

Fourier decomposition in azimuth

Because the problem is azimuthally symmetric about the solar direction, can be expanded as a finite Fourier series in (Sanghavi 2014, Eq. 8):

is block-diagonal (mixing intensity with linear polarization ); is block-off-diagonal (mixing with the rest). Each Fourier moment solves an independent RTE — the -loop in rt_run.jl is therefore embarrassingly parallel.

Quadrature in μ

After Fourier decomposition, each per-moment integral over is replaced by a finite sum. vSmartMOM offers two quadrature schemes (src/CoreRT/tools/rt_set_streams.jl):

GaussLegQuad — half-space Gauss-Legendre on . By default (external_solar=false), the solar direction is appended to the diffuse operator as a zero-weight node alongside the viewing directions. With the opt-in external_solar=true (TOA-only), only viewing directions are appended; exact remains a scalar source direction and a phase-evaluation column.
RadauQuad — block-Radau composite on with as a true quadrature node carrying non-zero weight (Sanghavi 2014 App. B, Eqs. B.1–B.2). This is the historical DNI-oriented solution to off-quadrature solar geometry. It retains embedded and does not support external-solar mode.

Historical DNI design and the current external-solar SFI path Sanghavi 2014 §2.2 deliberately excluded direct-solar single scattering from . That historical formulation uses the shorter matrix-operator recurrences (Eqs. 23–28) and Block-Radau (App. B) to place on the quadrature. This explains the DNI/Radau lineage, but it is not a statement about today's SFI kernel.

The current solver computes direct-beam source-function integration (SFI). Its backward-compatible default keeps in the operator grid as a zero-weight node because several solver paths still depend on that layout. The opt-in external-solar form (external_solar=true, a TOA-only path consumed through rt_run_toa) removes this representational relic: is not a diffuse stream. The elemental solver receives it as a scalar propagation direction. Exact rectangular columns generate finite-layer operators, which are then contracted into . Diffuse R/T operators remain square on the diffuse grid.

This path is restricted to Gauss quadrature, elastic forward/linearized noRS, forward rotational Raman RRS, Lambertian surfaces, and endpoint TOA output. Call rt_run_toa(model) or rt_run_toa(rs,model); BOA, HDR, and BHR arrays are not allocated. Raman Jacobians, VRS, single-scatter-only, non-Lambertian, and interior-sensor calculations retain embedded . :::

The supermatrix form (per Fourier moment)

For each Fourier moment , the RTE becomes (Sanghavi 2014, Eq. 12):

After quadrature in , the upwelling and downwelling streams are stacked into supervectors of length (Sanghavi 2014, Eqs. 14–15):

and are diagonal supermatrices of quadrature angles and weights. The four blocks are the directional pairings of the phase matrix.

The / superscripts on that you see throughout the source code (r⁻⁺, t⁺⁺, R⁻⁺, etc.) are exactly this convention: is downward, is upward.

The four arrays each layer reduces to

For each Fourier moment , every atmospheric layer is captured by two spectral vectors and two phase-matrix arrays. Define for the diffuse operator and for phase evaluation:

Symbol	Meaning	Shape
	optical depth (per spectral point)	`(nSpec,)`
	single-scattering albedo	`(nSpec,)`
	forward (downward → downward) Fourier-moment phase matrix	`(NphaseN, NphaseN, nSpec)`
	backscatter (downward → upward) Fourier-moment phase matrix	`(NphaseN, NphaseN, nSpec)`

The third dimension is spectral. That is the whole batched-matmul story — see Concepts/07. Normally NphaseN == NquadN. With external-solar SFI, phase evaluation may have one additional direction; the diffuse elemental/doubling/adding operators remain (NquadN, NquadN, nSpec), and only the exact solar phase column is consumed outside that block. Five weighted streams plus one distinct VZA in Stokes_IQU therefore produces an diffuse operator, not .

Polarization is a type, not a runtime branch

Stokes_I (, scalar intensity), Stokes_IQ (), Stokes_IQU (), Stokes_IQUV () each subtype AbstractPolarizationType in src/Scattering/types.jl:142–215. Each carries the polarization symmetry vector (used in Concepts/04) and the unit solar Stokes vector .

The kernels specialize on the type at compile time. There is no if pol_type.n == 4 branch inside the inner loops — Julia's multiple dispatch emits separate code paths for each polarization at JIT time. Adding Stokes_IQ took adding one struct definition; the kernels picked it up automatically. (The RT solvers now reject Stokes_IQ: dropping U from multiple scattering can cause intensity errors of several percent; in the tested aerosol scene they were largest in the principal plane. The type remains for the Scattering-level functions and StandaloneSS ExactSSConfig runs.)

Code anchors

Concept	Source
Top-level RT loop	`src/CoreRT/rt_run.jl` (`rt_run` → `_rt_run_column`)
Fourier-moment dispatch	`src/CoreRT/rt_run.jl::_rt_run_column` (`for m = 0:m_max`)
Quadrature construction	`src/CoreRT/tools/rt_set_streams.jl:23–256`
`QuadPoints` struct	`src/CoreRT/types.jl::QuadPoints`
External direct-solar phase coupling	`src/CoreRT/CoreKernel/elemental.jl::get_elem_rt_SFI!`
Lean TOA-only entry point	`src/CoreRT/rt_run.jl::rt_run_toa`
Polarization types	`src/Scattering/types.jl:142–215`
Phase matrix Z from Greek	`src/Scattering/compute_Z_matrices.jl::compute_Z_moments`
Supermatrix layer types	`src/CoreRT/types.jl::CompositeLayer/AddedLayer`

Hands-on tutorials

Runnable examples with Plotly figures:

References

Sanghavi et al. (2014), vSmartMOM, JQSRT 133:412–433, doi:10.1016/j.jqsrt.2013.09.004. Eqs. (2), (4), (6), (8)–(10), (12), (14)–(15); App. B (Block-Radau). Primary reference for everything on this page.
Chandrasekhar (1950), Radiative Transfer. (Original vector RTE in this form.)
Hovenier (1971), Symmetry relationships for scattering of polarized light…, J. Atmos. Sci. 28:488. (Block diagonal/anti-diagonal structure of .)
Full crib sheet: docs/dev_notes/theory_references.md §A, §E.

2 · Vector RTE & Discretization ​

The vector RTE ​

Boundary conditions ​

Fourier decomposition in azimuth ​

Quadrature in μ ​

The supermatrix form (per Fourier moment) ​

The four arrays each layer reduces to ​

Polarization is a type, not a runtime branch ​

Code anchors ​

Hands-on tutorials ​

References ​