Backward Monte Carlo Sky Radiance Spectral radiative transfer for a planetary atmosphere, rendered from a spacecraft — current target: Mars from ~100,000 km.

This code renders the full sky (or a planetary disk) as seen by an observer, by tracing light backward from the camera through a spherical, scattering atmosphere and off the ground. It is an independent implementation of the MYSTIC-style method — the backward Monte Carlo technique with next-event ("local") estimation of Mayer (2009), the algorithm behind MYSTIC / libRadtran — with several variance-reduction and modeling choices of its own (see Relationship to MYSTIC). Per color channel it produces an absolute radiance image that is then gamma-mapped to RGB.

Camera ray→ Forced first scatter→ Random walk (natural free paths)→ Local estimate to Sun at each vertex→ + deterministic direct terms→ Pixel radiance

1Geometry & observer

The marcher is fully spherical and planet-aware: a planet radius is selected per body (Mars is the default; Earth, Titan and Pluto blocks exist). The terrain datum sits at htmin = 0 and the top of atmosphere at htmax = 100 km; the observer height is a command-line override.

Space observer. A vacuum fast-skip advances each ray to the atmosphere entry (htmsl = htmax) before any sampling. The disk then spans nadir (alt = −90°) out to the limb (≈ −88.06° at 100,000 km).
Ground observer. Rays start in-atmosphere at the observer height — the original surface-based configuration this core was first validated against.

2The backward random walk

nray photons are launched per pixel from the observer along the line of sight. The walk is sun-independent: it depends only on the atmosphere along the path, never on where the Sun is. The Sun enters only through the local estimate (§3).

Forced first collision. A uniform random number selects the optical depth of the first scatter, drawn from an exponential free-path law but truncated to the total line-of-sight optical depth od_slant — so the scatter is forced to fall somewhere inside the atmospheric column rather than passing straight through (τ₁ = −ln(1 − ξ · prob_1st_ext), with ξ ∈ [0,1) and prob_1st_ext = 1 − trans(od_slant) the fraction of the column that scatters). The photon weight w₀ = prob_1st_ext unbiases that forcing. The geometric location is then found by marching a scalar optical-depth accumulator tausum out to τ₁, trapezoidally in the extinction coefficient.
Subsequent scatters use natural free-path sampling — the photon escapes thin air on its own, so multiple scattering terminates by itself. The walk continues while photon_wt > prob_1st_ext · 1e−3.
Scattering direction at each vertex is set by three uniform random draws, with the phase function supplying all the angular structure (importance sampling):
- One draw picks the scatterer — gas, coarse aerosol, or fine aerosol — in proportion to each species' local scattering contribution.
- One draw sets the deflection angle θ by inverse-CDF sampling: the phase function is pre-integrated over the sphere into a cumulative distribution and inverted into a lookup table, so a flat random number ξ ∈ [0,1) comes out distributed exactly as the phase function — mostly forward for dust, no rejection looping. (Rayleigh uses the same idea in closed form.)
- One draw sets the azimuth φ uniformly in [0°, 360°), since scattering is rotationally symmetric about the incoming direction.
The (θ, φ) pair is then rotated from the photon frame into the world frame to give the new direction. The walk samples the delta-scaled phase function while the local estimate (§3) deposits with the full one — the TMS correction (§5).
Weight decay. Each scatter multiplies the weight by the single-scattering albedo ω₀; each ground bounce multiplies by the surface albedo. No sun-path or shadow factor enters the walk.

3The local estimate (next-event)

At every scattering vertex, the deterministic contribution of light arriving directly from the Sun is added — this is the only place solar geometry enters, and it keeps variance low (no extra rays are spawned toward the Sun):

rad_local_est = 3e9 · weighted_pf · sun_trans      (atmospheric scatter)
rad_local_est = 6e9 · weighted_pf                  (ground bounce)

summed as   Σ  rad_local_est · photon_wt

sun_trans = trans(od_2_sun) is the Sun-path transmission, read from a precomputed grazing-path table (§6). A blocked Sun (shadow sentinel) contributes 0.
weighted_pf = frac_a·aero_pf + frac_g·rayleigh_pf mixes the aerosol and gas phase functions by their local extinction fractions. The aerosol term uses the full two-term Henyey–Greenstein phase function even on a delta-scaled walk (the TMS correction — see §5).

4Deterministic direct terms (noise-free)

Two contributions are computed analytically rather than sampled, so they carry no Monte Carlo noise:

Direct surface: prob_1st_trans · rad_direct_sfc — the unscattered camera ray reaching the ground, reflecting toward the Sun (12e9 · sin(alt_2sun) · trans(sun-path) · albedo · f_BRDF). Gives subsolar cosine-brightening and surface terrain detail; f_BRDF is the simplified-Hapke surface shape factor (see §7).
Direct sky: prob_1st_trans · rad_sky_dir — the unscattered straight-through term.

The final per-pixel radiance combines the direct terms with the averaged walk contribution:

rad_ray = prob_1st_trans·rad_sky_dir  +  rad_dir_g  +  rad_sky_sum / nray

5Phase function & delta-scaling

Aerosol scattering uses a per-color two-term Henyey–Greenstein function (forward/back blend with asymmetries g₁, g₂, the sharp forward lobe reaching g₂ ≈ 0.91–0.945). A sampling CDF is built once and reused.

Delta-scaling tames the runaway multiple-forward-scattering that the sharp coarse lobe would otherwise produce. A per-color truncation angle (the coarse/fine-lobe crossover, ~15°, parameter-free) removes the forward cone from the walk's sampling CDF and reassigns the removed fraction f to the unscattered beam (ω₀′ = ω₀(1−f)/(1−ω₀f), ext′ = ext·(1−ω₀f)). Crucially, the local-estimate deposit keeps the full phase function (TMS correction), preserving the bright forward limb while stabilizing the walk.

6Atmosphere & aerosol profile

Three components per color: gas (Rayleigh) extinction, aerosol (dust) extinction, and ozone optical depth, each with its own single-scattering albedo.
Dust vertical structure = a surface exponential (scale height ~11 km on Mars, a well-mixed troposphere) tapered linearly to zero up to ~64 km, plus a detached trapezoidal slab. A horizontal modulation adds large-scale azimuthal structure.
Shared profile module. Both the random walk and the Sun-path lookup table read the density profile from one place, so they cannot drift apart.
Sun-path table. A precomputed, color-independent table tabulates the grazing Sun-path optical-depth integral over (altitude, solar elevation); blocked geometry returns a sentinel. This is what makes the next-event estimate cheap.

7Surface albedo

Scalar default (Mars): albedo = 0.29 / 0.17 / 0.09 for R/G/B.
Spatial map. An optional Mars albedo image overrides the scalar per ground point. It is gamma-decoded, then per-channel area-weighted mean-normalized so each channel's map mean equals the scalar target — preserving spatial pattern while anchoring overall level and the Mars color. A common ceiling renders saturated polar frost neutral rather than red. The map is absent ⇒ off (safe).
Sub-observer registration. A chosen Mars latitude/longitude is placed at the disk center via command-line flags; the map center column is 0°E with east positive.
Directional reflectance (simplified Hapke). Rather than a flat Lambertian surface, the reflection carries a dimensionless shape factor f_BRDF(μ₀, μ, g) that reproduces the two regolith effects Mars shows: an opposition surge (brightening toward zero phase angle, g→0) and sunward darkening (shadow-hiding dimming toward the Sun, g→180°). It combines a Lommel–Seeliger limb term, a backscattering Henyey–Greenstein phase function, and Hapke's opposition term B(g) = B₀ / (1 + tan(g/2)/h). The incidence μ₀, emergence μ, and phase angle g are all read off at the reflection point. The factor is self-normalized to ≈1 at a reference geometry, so the albedo map stays the energy scale and the absolute calibration is closely preserved. It is applied to both the deterministic direct-surface term (§4) and the ground-bounce local estimate.

8Absolute calibration

The model carries one geometry-independent constant: 3e9 = F_☉/4π for atmospheric scattering, 6e9 for ground reflection, and 12e9 = F_☉/π for the direct surface term (Lambertian normalization, reshaped by the surface BRDF factor). All solar-geometry and path dependence is carried by the transport itself; a single re-tune of this constant rescales the absolute brightness if a physics change shifts the overall level.

9Summary of the data flow

For each pixel and color channel:

Launch nray photons backward from the observer.
Force the first scatter inside the line-of-sight optical depth; set the photon weight.
Random-walk by natural free paths, decaying the weight by ω₀ (and surface albedo at bounces) until it falls below threshold.
At every vertex, add the next-event estimate of direct sunlight, weighted by the full phase function and the tabulated Sun-path transmission.
Add the two deterministic direct terms (sky + surface).
Average over rays, scale by the calibration constant, and gamma-map the three channels to RGB.

10Relationship to MYSTIC — and what's original here

This solver follows the same algorithm family as MYSTIC (the libRadtran Monte Carlo solver), as formulated by Mayer (2009), but it is an independent implementation with its own variance-reduction structure, an optimization for the space-observer geometry, and application-specific Mars physics. In some respects — full 3-D geometry, measured Mie scattering — it is also deliberately simpler than MYSTIC.

Shared with MYSTIC (the standard backward-MC core)

Backward photon tracing from the observer, with the local (next-event) estimate as the only point where the Sun enters.
Natural free-path sampling on subsequent scatters; the walk terminates by photon-weight threshold, weight decaying as ∏ω₀ · ∏albedo (Mayer's w₀).
Phase-function importance sampling by inverted CDF, and delta-scaling / TMS (sample the truncated phase function, deposit the full one).

Original choices here

Forced first collision with analytic truncation (τ₁ = −ln(1 − ξ·prob_1st_ext), weight prob_1st_ext), applied only to the first scatter — a hybrid forcing/natural scheme that guarantees every ray contributes without wasting clean-transit rays.
Deterministic direct terms split out as noise-free analytic contributions: the direct ray→ground→Sun term (12e9 · sinα · trans · albedo · f_BRDF) and the direct-sky term.
Simplified-Hapke surface BRDF — a self-normalized directional shape factor giving the Mars opposition surge and sunward shadow-hiding, applied to both the direct-surface term and the ground-bounce local estimate (see §7), rather than a flat Lambertian surface.
Precomputed grazing Sun-path table (mc_lookup) over (altitude, solar elevation) — an optimization for rendering a planet from space, where MYSTIC would trace those Sun paths directly.
Application physics not part of any general solver: a detached dust slab plus horizontal/latitude-band aerosol modulation, a high-altitude aerosol taper, a two-term Henyey–Greenstein dust phase function, a spatial Mars albedo map with sub-observer registration, and a single F/4π solar-relative calibration constant.

Where it is simpler than MYSTIC

Surface BRDF in the local estimate only. The simplified-Hapke factor reshapes the radiance seen along the view direction, but the continuation of a reflected photon is still sampled with a cosine (Lambertian) bounce — so multiply-scattered surface light is not yet fully BRDF-consistent. A BRDF-sampled reflected walk is the natural next step.
Some plane-parallel shortcuts remain in the slant-optical-depth bookkeeping even though the marcher itself is spherical; MYSTIC is fully 3-D spherical.
Henyey–Greenstein-family aerosol phase function rather than full Mie / measured scattering.