VANvonZHANG · Aug 10, 2026 · Aug 10, 2026 · Aug 11, 2026 · Aug 11, 2026 · Aug 11, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 Highlights: new `pyradtran.data` data layer with bundled libRadtran data, a LEGO-style "blocks" aerosol API, OPAC folding with real Mie phase functions, and a publication-ready visualization/workflow suite. Contains **breaking changes** to the aerosol API (see Removed).
 
+#### Example self-validation + analysis expansion
+- `compute_energy_budget` + `assert_energy_conservation` (`pyradtran.core.postprocess`): column energy identity `F_inc = eup_TOA + (1-a)(edir+edn)_surf + F_abs_atm` with hard physical-bound assertions.
+- `parse_heating_ascii` (`pyradtran.core.output_parser`): parses libRadtran's wide heating-rate format (K/day per zout per wavelength). `Runner.execute` auto-detects heating mode (scene has `heating_rate` but no `output_user` quantities) and dispatches to it.
+- Viz: `plot_size_distributions`, `plot_phase_functions`, `plot_legendre_decay`, `plot_block_spectral_optics`, `plot_drf_spectral` (direct radiative forcing), `plot_spectral_attribution`.
+- `examples/multicomponent_viz`: expanded from 10 to 23 plots (per-block spectral optics, scattering phase functions, Legendre decay, heating-rate profile, DRF via no-aerosol baseline, spectral attribution) with a hard energy-conservation assertion and a column-τ-stacking check. Heating rates require a second uvspec invocation (libRadtran's `heating_rate` mode replaces flux output rather than appending a column); the demo runs both and merges.
+- `examples/multicomponent_viz/canonical.py`: shared scene config incl. `build_scene_heating` (used by both the demo and the regression test so they cannot drift).
+- `tests/fixtures/multicomponent_baseline.json` + `tests/test_multicomponent_regression.py` (slow/uvspec-gated): end-to-end regression baseline for the canonical 3-block scene.
+- `scripts/regen_baseline.py` to regenerate the baseline.
+
 ### Added
 
 #### Data layer (`pyradtran.data`)
@@ -49,7 +58,20 @@ Highlights: new `pyradtran.data` data layer with bundled libRadtran data, a LEGO
 - Comprehensive viz + workflow demo: `MieSpecies` → DISORT → all plots + attribution.
 - Sphinx user guide, expanded API reference, and rewritten README.
 
+#### YAML configuration front-end (`pyradtran.config`)
+- `config_version: 1` YAML schema (strict: scene / aerosol / analysis; block kinds `mie`/`bulk`/`opac_preset`/`explicit_layer`; placements `od_inversion`/`mass`/`exponential`/`tabulated`), `load_config` / `export_config`, and the `run_config` orchestrator (energy-conservation assertion, heating second run, DRF baseline, leave-one-out attribution, plot registry, NetCDF export).
+- CLI `pyradtran run|validate|export-config` (also `python -m pyradtran`): `validate` checks a config without invoking uvspec; `export-config` writes canonical YAML.
+- Hard round-trip guarantee, enforced by `tests/test_config_roundtrip.py`: a YAML config and its API-built twin emit byte-identical uvspec input and byte-identical `.master` layer files; `export_config` output re-loads identically.
+- `examples/multicomponent_viz/canonical.yaml` — YAML twin of the canonical scene, generated from `canonical.py` by `make_yaml.py`.
+#### Benchmarks (`pyradtran.benchmarks`)
+- Randles et al. (2013) AeroCom shortwave RT intercomparison replication: `run_randles2013` (3 cases × 2 AFGL atmospheres × 2 SZAs × 2 bands = 24 uvspec runs; Å=1 power law, 0–2 km linear taper, HG phase function as PMOM `beta_l = g**l`, LBL-median normalization constants).
+- Bundled reference `reference/randles2013_lbl.json`: 64 tabulated LBL values plus a `_meta` block (thresholds: fluxes ±8% PASS / ±12% WARN; RF ≤15% or ≤1.5 W/m²).
+- `compare_benchmark` / `format_report` / `write_report` / `plot_benchmark_overlay`: per-row status classification, Markdown + CSV report, LBL-overlay PNG.
+- Ready-to-run example `examples/randles2013_benchmark/run_benchmark.py`: matrix -> report + overlay PNG in one command; `--from-json` regenerates reports offline; non-zero exit on any FAIL (CI gate).
+- Regression test `tests/test_benchmark_randles.py` (slow, uvspec-gated): real run vs libRadtran 2.0.6 — 68 comparison rows, 64 PASS / 0 WARN / 0 FAIL (+4 report-only n/a).
+- Ratified deviation: benchmark results persist as JSON, not NetCDF (nested ragged dict).
 ### Changed
+- `evaluate_blocks_on_grid` (`pyradtran.core.postprocess`) now also returns per-block `ssa` and `g` (was `tau`/`rho_kg_m3` only).
 - `OpacPreset` / `OpacCustom` are now `PlacedBlock` factories; OPAC convenience functions expose `output_dir`.
 - Runner resolves the libRadtran data path through `DataResolver`.
 - **Default `Output.format` is now `"ascii"`** (was `"netcdf"`): uvspec's NetCDF output is broken in many libRadtran builds (a 0-byte `.nc` from a libnetcdf ABI mismatch), so the previous default crashed every simulation that omitted `format=`. ASCII works everywhere and yields an equivalent `xarray.Dataset`; NetCDF remains available via `format="netcdf"`. The convenience functions no longer hardcode `"netcdf"` — they inherit the ASCII default.
@@ -62,6 +84,9 @@ Highlights: new `pyradtran.data` data layer with bundled libRadtran data, a LEGO
 - **Breaking:** removed `ExternalFile` / `ExternalAerosol` and the `run_with_aerosol` convenience function — the only non-LEGO aerosol path. Pre-computed explicit files now go through `DirectLayerOpticsBlock` mixed into a `CompositeAerosol`. Also removed the unused `AerosolModel.set_tau_at_wvl` / `king_byrne` fields (never set by any code path).
 
 ### Fixed
+- ASCII output parser names the 8th uvspec column `heating_rate` (was `col_7`) so `HEATING_RATE_COLUMN` resolves when libRadtran emits heating rates alongside the standard 7 flux columns.
+- `_trapz` guard added to `tests/test_composite_aerosol_unit.py` normalization tests (numpy <2.0 compat; production `mie.py` already had the guard).
+- `run_satellite` test now skips cleanly when libRadtran's `MPS` satellite-geometry file is absent (not bundled).
 - libRadtran integration issues for explicit aerosol coupling.
 - ASCII column-name mapping and aerosol-layer file cache collision.
 - Mixing-rule HG fallback now uses `g_l` (PMOM) Legendre moments consistently across species.
@@ -70,7 +95,10 @@ Highlights: new `pyradtran.data` data layer with bundled libRadtran data, a LEGO
 - CI branch triggers, lint, and test assertions realigned for green CI.
 - Convenience functions no longer crash on the NetCDF parse path: they now use ASCII output (see **Changed**). `_parse_netcdf` raises a clear error pointing at `format='ascii'` when uvspec produces an empty/missing NetCDF file instead of xarray's cryptic backend-mismatch message.
 - `run_with_opac_preset` / `run_with_opac_custom` now set `disort_intcor="moments"` on their DISORT solver. The OPAC folding produces Legendre-moment phase functions, which DISORT rejects without this flag (`you need to specify 'disort_intcor moments'`); these functions previously crashed at runtime.
+- **i550 wavelength-index bug** in the multicomponent demo and its regression test: `argmin(|wavelength - 0.55|)` compared an nm grid against 0.55 µm, silently selecting 401 nm everywhere a 550 nm scalar was labeled (energy log, heating log, baseline fixture). Also removed a spurious `× 1000` on the DRF wavelength axis (same µm/nm confusion) that made `np.interp(550, ...)` clamp to the 401 nm endpoint. The committed baseline now holds true 550 nm values (edir_surf 798.75, eup_toa 194.79 W/m², F_abs_atm 336.94 W/m²).
 
+- `Runner.execute_many` no longer swallows worker failures: failing scenes were silently returned to callers as `(idx, exception)` result entries; the first failure now raises `RuntimeError(f"scene {i} failed")` chained to the original exception, after cancelling outstanding futures.
+- `Runner.execute_many` uses a thread pool instead of a process pool: scenes carrying `CompositeAerosol` can hold unpicklable payloads (aerosol3D bulk size distributions store local closures), which crashed `ProcessPoolExecutor` pickling; uvspec runs as a subprocess and releases the GIL while waiting, so threads parallelize the invocations just as well and nothing needs pickling.
 ## [0.1.0] - 2026-05-12
 
 ### Added

diff --git a/docs/api.rst b/docs/api.rst
@@ -253,3 +253,67 @@ Workflow
    :members:
    :undoc-members:
    :show-inheritance:
+
+YAML Configuration Front-end
+----------------------------
+
+.. automodule:: pyradtran.config
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Configuration Schema
+^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: pyradtran.config.schema
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Configuration Loader
+^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: pyradtran.config.loader
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Config Orchestrator
+^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: pyradtran.config.orchestrator
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Command Line
+^^^^^^^^^^^^
+
+.. automodule:: pyradtran.config.cli
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Benchmarks
+----------
+
+.. automodule:: pyradtran.benchmarks
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Randles et al. (2013)
+^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: pyradtran.benchmarks.randles2013
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Benchmark Comparison and Reporting
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: pyradtran.benchmarks.compare
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/guide/benchmarks.rst b/docs/guide/benchmarks.rst
@@ -0,0 +1,209 @@
+Benchmarks
+==========
+
+``pyradtran.benchmarks`` replicates published radiative-transfer
+intercomparisons end-to-end with real uvspec runs and compares the results
+against the published line-by-line (LBL) reference values. It currently ships
+one benchmark: the shortwave aerosol experiment of Randles et al. (2013).
+
+Randles et al. (2013) protocol
+------------------------------
+
+Randles, C. A. et al., "Intercomparison of shortwave radiative transfer
+schemes to improve the aerosol and surface albedo treatment in climate
+models", Atmos. Chem. Phys. 13, 2347-2362 (2013) — the AeroCom RT experiment,
+Tables 3, 5, A3 and A4. Three cases:
+
+.. list-table::
+   :header-rows: 1
+
+   * - Case
+     - Aerosol
+     - AOD550
+     - Ångström
+     - SSA
+     - g
+   * - ``case1``
+     - none (Rayleigh only)
+     - —
+     - —
+     - —
+     - —
+   * - ``case2a``
+     - lowest 2 km
+     - 0.2
+     - 1
+     - 1.0
+     - 0.7
+   * - ``case2b``
+     - as 2a
+     - 0.2
+     - 1
+     - 0.8
+     - 0.7
+
+Fixed conditions: surface albedo 0.2 (Lambert), SZA 30° and 75°, AFGL subarctic
+winter (``saw``) and tropics (``trop``) atmospheres, no clouds, and two bands —
+broadband 0.2–4.0 µm (``bb``) and UV-VIS 0.2–0.7 µm (``uvvis``). The full
+matrix, 3 cases × 2 atmospheres × 2 SZAs × 2 bands, is 24 uvspec invocations
+(~20–60 s). The solver setup is DISORT with 16 streams, ``disort_intcor
+moments``, pseudospherical geometry and reptran correlated-k with
+band-integrated output; the phase function is supplied as 32 Legendre moments.
+
+Replication assumptions
+-----------------------
+
+The paper leaves details to each participating model; this replication fixes
+them as:
+
+1. Built-in AFGL profiles are used (the participants used their
+   individual-model profiles; the comparison tolerance absorbs this).
+2. "Aerosol linear in the lowest 2 km" is a linear taper from the surface to
+   zero at 2 km, ``w(z) = max(0, 1 - z / 2 km)``.
+3. ``g = 0.7`` Henyey–Greenstein phase function represented through PMOM
+   moments ``beta_l = g**l`` (the DISORT/libRadtran coefficient form).
+4. Ångström exponent 1 means ``tau(lambda) = 0.2 * (lambda / 0.55 um)**-1``.
+
+Normalization
+-------------
+
+Each result is divided by the same run's band TOA incident flux and multiplied
+by the LBL median constants, which removes the solar-constant choice.
+Formally ``F_norm = F_raw / edir_toa_raw(same run) * C[band][sza]`` with:
+
+.. list-table::
+   :header-rows: 1
+
+   * - Constant (W/m²)
+     - SZA 30°
+     - SZA 75°
+   * - Broadband ``bb``
+     - 1189.28
+     - 355.43
+   * - UV-VIS ``uvvis``
+     - 563.38
+     - 168.37
+
+Reference and thresholds
+------------------------
+
+``pyradtran/benchmarks/reference/randles2013_lbl.json`` bundles the 64
+tabulated LBL values (per case and quantity, for the four configs
+``saw30`` / ``saw75`` / ``trop30`` / ``trop75``) plus a ``_meta`` block holding
+the citation and the comparison thresholds:
+
+.. list-table::
+   :header-rows: 1
+
+   * - Quantity
+     - PASS
+     - WARN
+     - FAIL
+   * - Fluxes (everything except ``*_rf``)
+     - ``|rel| <= 8%``
+     - ``|rel| <= 12%``
+     - beyond
+   * - Radiative forcing (``toa_rf`` / ``sfc_rf``)
+     - ``|rel| <= 15%`` **or** ``|abs| <= 1.5 W/m²``
+     - ``|rel| <= 25%`` or ``|abs| <= 3 W/m²``
+     - beyond both
+
+Case-2a atmospheric forcing (``toa_rf - sfc_rf``) is below 1 W/m² and
+RSD-invalid per the paper; those four rows are classified ``n/a`` and are
+report-only.
+
+Running the benchmark
+---------------------
+
+.. code-block:: python
+
+   from pyradtran.benchmarks import (
+       CASES,
+       compare_benchmark,
+       format_report,
+       load_reference,
+       plot_benchmark_overlay,
+       run_randles2013,
+       write_report,
+   )
+
+   reference = load_reference()  # bundled LBL values + _meta
+   print(CASES, sorted(k for k in reference if not k.startswith("_")))
+
+   def full_benchmark(out_dir, uvspec_exe=None, data_path=None):
+       """Run the matrix (needs uvspec), compare against LBL, emit report + plot."""
+       results = run_randles2013(out_dir, uvspec_exe=uvspec_exe, data_path=data_path)
+       rows = compare_benchmark(results)  # one status row per case/quantity/config
+       md, csv = write_report(rows, f"{out_dir}/randles2013.md")
+       png = plot_benchmark_overlay(rows, f"{out_dir}/randles2013_overlay.png")
+       return format_report(rows), md, csv, png
+
+:func:`~pyradtran.benchmarks.randles2013.run_randles2013` runs the whole
+matrix (pass ``cases=`` to restrict it) and returns a nested results dict.
+
+The regression suite runs the real matrix. It is marked ``slow`` and skips
+automatically when uvspec or the libRadtran data directory are not found:
+
+.. code-block:: console
+
+   $ python -m pytest tests/test_benchmark_randles.py -v -s   # ~25 s, real uvspec runs
+   $ python -m pytest -q -m "not slow"                        # the rest of the suite
+
+A ready-to-run driver for the whole workflow lives at
+``examples/randles2013_benchmark/run_benchmark.py``: it runs the matrix,
+compares against the reference and writes the Markdown/CSV report and the
+overlay PNG to its ``output/`` directory. ``--from-json`` regenerates the
+report from an earlier results JSON without re-running uvspec; the exit
+code is non-zero if any row FAILs.
+
+Outputs
+-------
+
+- ``compare_benchmark(results)`` — one row per (case, quantity, config) with
+  ``lbl``, ``ours``, ``abs_diff``, ``rel_diff`` and ``status`` (``PASS`` /
+  ``WARN`` / ``FAIL`` / ``n/a``).
+- ``format_report(rows)`` — a Markdown report string, grouped by case with a
+  status-count summary.
+- ``write_report(rows, path)`` — writes the Markdown report plus a CSV twin
+  next to it; returns both paths.
+- ``plot_benchmark_overlay(rows, path)`` — a PNG overlaying this run's values
+  on the LBL reference.
+
+.. note::
+
+   ``run_randles2013`` persists its results as **JSON**, not NetCDF: the
+   matrix is a nested, ragged dict (per-case quantities differ), which NetCDF
+   cannot represent. The JSON additionally carries a ``results_path`` key
+   pointing at the artifact.
+
+Validation result
+-----------------
+
+Real run against libRadtran 2.0.6 (uvspec, 2026-08-16): the regression suite
+passes — 28 tests in 22.84 s — and the full-matrix comparison yields
+**68 rows: 64 PASS / 0 WARN / 0 FAIL** (plus the 4 report-only ``n/a`` rows).
+The three largest deviations, all comfortably inside the 8% PASS band:
+
+.. list-table::
+   :header-rows: 1
+
+   * - Row (case / quantity / config)
+     - LBL
+     - pyRadtran
+     - rel.
+   * - case1 / absorptance / trop75
+     - 0.307
+     - 0.2858
+     - -6.91%
+   * - case1 / nir_sfc_down / trop75
+     - 101.1
+     - 104.95
+     - +3.81%
+   * - case1 / direct_bb_sfc_down / trop75
+     - 179.6
+     - 184.26
+     - +2.59%
+
+The paper's central case-2b result reproduces: positive TOA forcing at SZA 30
+(ours +11.4 / +10.3 vs LBL +11.6 / +10.3 W/m² for ``saw`` / ``trop``) and
+negative surface forcing in every configuration.