Overview and notation

The problem

We solve the azimuthally averaged radiative transfer equation for a plane-parallel atmosphere composed of a stack of homogeneous layers. Each layer \(l\) is characterised by

• an optical depth \(\tau_l\) (delta_tau[l]),

• a single-scattering albedo \(\omega_l\) (single_scat_albedo[l]),

• a set of Legendre phase-function moments \(\chi_\ell\) (phase_function_moments[l]),

together with, for thermal problems, the Planck emission at the bounding levels.

The atmosphere is discretised in the polar-angle cosine \(\mu\) using Gauss–Legendre quadrature on the half-range \([0, 1]\) with \(M\) nodes \(\mu_i\) and weights \(w_i\) (num_quadrature). A radiation field is represented by two length-\(M\) vectors: the upward intensities \(I^{\uparrow}_i \equiv I(+\mu_i)\) and the downward intensities \(I^{\downarrow}_i \equiv I(-\mu_i)\).

The reflection / transmission operators

Each layer, and every composite of layers, is described by four operators and a pair of source vectors:

Symbol	Meaning
\(\mat{R}_{ab}\)	Reflection for illumination from above (top interface \(a\)).
\(\mat{R}_{ba}\)	Reflection for illumination from below (bottom interface \(b\)).
\(\mat{T}_{ab}\)	Transmission for illumination from above.
\(\mat{T}_{ba}\)	Transmission for illumination from below.
\(\vect{s}^{+}\)	Internal source emerging upward from the top of the (composite) layer.
\(\vect{s}^{-}\)	Internal source emerging downward from the bottom.

For a single homogeneous layer the operators are symmetric under a flip of the illumination direction: \(\mat{R}_{ab} = \mat{R}_{ba}\) and \(\mat{T}_{ab} = \mat{T}_{ba}\). This symmetry is broken as soon as two different layers are combined by adding, because a composite of inhomogeneous sub-layers is itself inhomogeneous; the code therefore tracks all four matrices from the outset.

The two-step algorithm

The solve proceeds in two stages:

Doubling (Doubling: a single homogeneous layer)

builds the \(\mat{R}\), \(\mat{T}\) and source vectors of one homogeneous layer of arbitrary optical depth by starting from an optically thin sublayer and repeatedly doubling it.

Adding (Adding: combining layers)

combines two (possibly inhomogeneous) layers into one composite by summing the geometric series of inter-reflections between them.

The full atmosphere is assembled by adding all layers together. To recover the intensity at internal interfaces, the code builds composites both from the top down (rtop) and from the bottom up (rbase) and matches them at each interface (Boundary conditions and intensity reconstruction).

Index conventions

Output arrays are indexed by interface (level) number:

index 0            = top of atmosphere (above layer 0)
index num_layers   = bottom of atmosphere (below the last layer / surface)

Layer \(l\) (for \(l = 0, \dots, N_{\text{lyr}}-1\)) is bounded above by level \(l\) and below by level \(l+1\). By default the user supplies arrays ordered top-of-atmosphere first; setting index_from_bottom lets you supply and receive them bottom-of-atmosphere first (the solver reverses internally).