Surface scattering dev by m-brath · Pull Request #1194 · atmtools/arts

m-brath · 2026-10-16T09:30:04Z

Surface scattering framework

This PR introduces a full surface scattering framework to ARTS3 — a new physics module, workspace integration, Python bindings, tests, examples, and documentation. Surfaces are no longer limited to plain emissivity: they can now scatter incoming radiance diffusely and specularly, with spectrally and spatially varying properties.

What's new

Core module (src/core/surface_scattering/)

SurfaceScatteringModel base class with a property contract of mandatory BRDF and emissivity tensors (MuelmatTensor5 / StokvecTensor3)
Three concrete models:
- Lambertian — including a spatially-varying LambertianSurfaceScattererField with adaptive longitude cyclers, interpolation, and frequency extrapolation limits
- FlatScalar — angle-independent scalar scattering
- Fresnel — dielectric surface via Fresnel equations
Model blending: MapOfSurfaceScatteringModel with configurable weighting options (normalized weights, enum-controlled behavior)

Workspace integration

New workspace variable surface_models and supporting types/groups
New spectral_rad_incoming_agenda (replacing spectral_rad_observer_agenda) with predefined options, including clear-sky Rayleigh scattering
New spectral radiance methods:
- spectral_radSurfaceScatteringSpecular / spectral_radSurfaceScatteringDiffuse
- spectral_radSurfaceScatteringDirect / spectral_radSurfaceScatteringDiffuseDirect / spectral_radSurfaceScatteringSpecularDirect — direct beams derived from workspace suns with refractive LOS
- spectral_radSurfaceScatteringInit for chaining multiple scattering contributions

Sun support (src/core/sun/, src/sun_methods)

Helper to derive scattering beams from suns; suns that cannot be placed in the sky of the current planet are now rejected

Python bindings — nanobind bindings for all models, the model map with weighting options, and star/sun helpers (py_surface_scattering.cpp, py_star.cpp)

Tests — 13 new regression tests under tests/core/surf/ covering each scattering method (incl. Jacobian orientation, sun double-counting, agenda chaining, predefined agendas, Rayleigh options, interpolation/extrapolation, and map weighting)

Examples — new examples/7-surface-scattering-properties/ suite: a checkerboard scene rendered clear-sky, with atmospheric LUTs (mm-wave and IR), and a camera photograph of the scene

Documentation — three-level docs (user / concept / dev Sphinx pages) for the framework, plus doc-infrastructure fixes: broken cross-refs, case-clashing pages, and a Sphinx extension fixing inspect.getsource in plot directives

Notable fixes along the way

Emission counting, BRDF formula, and extrapolation correctness fixes
Double-counting of suns in surface scattering eliminated
Jacobian tensor orientation fixed
size_t indexing and missing-field-initializer warnings cleaned up

Notes for reviewers

spectral_radSurfaceScatteringFlatDirect / FlatDiffuse were introduced and then renamed/removed during development; the final API is the Direct/Diffuse/Specular family above.

Copilot

🟡 Changes recommended

Critical Fresnel geometry and agenda-state defects remain, alongside unresolved interpolation, serialization, blending, and solar-position issues.

6 open findings

What changed in this PR

Introduces surface scattering to ARTS3, integrating physical models with workspace radiance calculations, Python interfaces, tests, examples, and documentation.

Changes:

Adds Lambertian, Fresnel, and flat-scalar models, spatial fields, and weighted blending.
Adds scattering agendas, solar-direction helpers, and Rayleigh-only source handling.
Adds regression coverage and documentation infrastructure updates.

File	Reviewed changes and unresolved findings
`tests/core/surf/spectral_rad_surface_scattering_sun_double_count.py`	Tests sun exclusion and double-counting prevention.
`tests/core/surf/spectral_rad_surface_scattering_specular.py`	Tests specular calculations and validation.
`tests/core/surf/spectral_rad_surface_scattering_predefined_agendas.py`	Tests predefined scattering agendas.
`tests/core/surf/spectral_rad_surface_scattering_chain.py`	Tests additive chaining and initialization.
`tests/core/surf/spectral_rad_incoming_agenda_rayleigh.py`	Compares Rayleigh agenda options.
`tests/core/surf/map_weighting_option.py`	Tests model weighting.
`tests/core/surf/lambertian_scatterer_field.py`	Tests spatial Lambertian properties.
`tests/core/surf/interp_extrapolation.py`	Tests frequency extrapolation modes.
`tests/core/surf/fresnel_surface_scatterer.py`	Tests Fresnel properties.
`tests/core/surf/flat_scalar_surface_scatterer.py`	Tests scalar scattering properties.
`src/workspace_variables.cpp`	Registers models and quadrature variables.
`src/workspace_methods.cpp`	Registers scattering methods. nit, 1 vote: Flux-conservation guidance omits the required `1/pi` normalization and actual-normal incidence convention.
`src/workspace_meta_methods.cpp`	Adds a scattering-only Rayleigh chain.
`src/workspace_groups.cpp`	Registers the model map and sun invariants.
`src/workspace_group_friends.cpp`	Registers the BRDF tensor type.
`src/workspace_agendas.cpp`	Defines incoming and surface agenda contracts.
`src/workspace_agenda_creator.cpp`	Assembles predefined agendas. critical, 3 votes: Solar paths are not computed for the newly traced incoming ray, causing missing-state failures or reuse of another ray’s paths.
`src/sun_methods.h`	Declares solar-direction helpers.
`src/sun_methods.cc`	Implements solar-direction helpers. moderate, 1 vote: Planet-centred sun distance is incorrectly interpreted as ellipsoid-relative altitude.
`src/python_interface/py_star.cpp`	Exposes solar helpers to Python.
`src/python_interface/py_module.cpp`	Initializes scattering bindings.
`src/python_interface/CMakeLists.txt`	Includes scattering bindings in the build.
`src/m_sun.cc`	Validates sun placement.
`src/m_srcvec.cc`	Adds zero-source initialization.
`src/core/surface_scattering/surface_scattering_properties.h`	Defines scattering property tensors.
`src/core/surface_scattering/surface_scattering_model.h`	Defines model variants and weighted maps.
`src/core/surface_scattering/surface_scattering_model.cc`	Implements blending and serialization. moderate, 3 votes: Zero-weight models are still evaluated and can throw. moderate, 3 votes: XML omits `weighting_option`, changing blending behavior after loading.
`src/core/surface_scattering/scattering_internal.h`	Adds shared interpolation helpers.
`src/core/surface_scattering/lambertian.h`	Declares Lambertian models. nit, 1 vote: Weight documentation must specify `cos(theta_inc)*dOmega/pi`, measured from the actual surface normal.
`src/core/surface_scattering/lambertian.cc`	Implements Lambertian properties and interpolation. moderate, 3 votes: Serialization omits `interp_extrapolation`. moderate, 1 vote: Singleton spatial axes cause out-of-bounds latitude access or degenerate longitude interpolation.
`src/core/surface_scattering/fresnel.h`	Declares Fresnel models.
`src/core/surface_scattering/fresnel.cc`	Implements Fresnel tensors and interpolation. critical, 3 votes: Geodetic angles ignore the actual surface normal and polarization frame, breaking equilibrium on tilted surfaces. moderate, 1 vote: Singleton spatial axes cause out-of-bounds latitude access or degenerate longitude interpolation.
`src/core/surface_scattering/flat_scalar.h`	Declares flat-scalar models.
`src/core/surface_scattering/flat_scalar.cc`	Implements scalar properties and interpolation. moderate, 1 vote: Singleton spatial axes cause out-of-bounds latitude access or degenerate longitude interpolation.
`src/core/surface_scattering/CMakeLists.txt`	Builds the scattering module.
`src/core/sun/sun.h`	Declares sun-placement validation.
`src/core/sun/sun.cc`	Implements sun-placement validation.
`src/core/CMakeLists.txt`	Integrates the core scattering module.
`src/CMakeLists.txt`	Links workspace scattering support.
`python/doc/source/conf.py`	Enables the plot-source extension.
`python/doc/source/_ext/plot_source_code.py`	Supports source inspection in plot directives.
`python/doc/gen_pyarts_arts_autoclass_rst.py`	Updates type references and case-clashing page handling.
`examples/7-surface-scattering-properties/surface-scattering-framework.py`	Demonstrates model construction.
`examples/7-surface-scattering-properties/README.rst`	Describes the example suite.
`doc/arts/users.rst`	Links the user guide.
`doc/arts/user.surface_scattering.rst`	Documents setup and usage. nit, 1 vote: Weight arrays use sorted-grid constructors instead of Vectors, so the recipe fails. nit, 2 votes: Quadrature omits the incidence cosine and `1/pi`, violating reflectivity/emissivity balance.
`doc/arts/devs.rst`	Links the developer guide.
`doc/arts/dev.surface_scattering.rst`	Documents implementation contracts.
`doc/arts/concepts.rst`	Links conceptual documentation.
`doc/arts/concept.surface_scattering.rst`	Explains scattering physics. nit, 1 vote: The flux-conservation condition omits `1/pi`, conflicting with complementary emissivity.
`.gitignore`	Excludes local planning files.

🧠 Review effort: Balanced

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

erikssonpatrick · 2026-10-16T09:58:42Z

+      const Numeric za_eff = za_fold(za_inc_grid[zi]);
+      const auto [Rv, Rh] = fresnel(Complex{n1}, Complex{n2}, za_eff);


I'm impressed that copilote spotted this! This passed my mind at the last meeting, but I forgot to ask because we ended up discussing other things. I think this is a nice feature, especially if we want to claim that ARTS3 is truly 3D. (in a 1D view there is no surface topology, but there is in 3D).
(We consider this in v2.6. At least regarding angles, but maybe not when it comes to polarisation.)

+      break;
+    case ClearskyRayleighScattering:
+      agenda.add("ray_path_observer_agendaExecute");
+      agenda.add("spectral_radClearskyRayleighScattering");


+  xml_write_to_stream(os, x.reflectivity_spectrum, pbofs);
+
+  tag.write_to_end_stream(os);
+}
+
+void xml_io_stream<surface_scattering::LambertianSurfaceScatterer>::read(
+    std::istream& is,
+    surface_scattering::LambertianSurfaceScatterer& x,
+    bifstream* pbifs) {
+  XMLTag tag;
+  tag.read_from_stream(is);
+  tag.check_name(type_name);
+
+  xml_read_from_stream(is, x.reflectivity_spectrum, pbifs);


+  for (const auto& [key, model] : models) {
+    surface_scattering::SurfaceScatteringModelProperties model_props =
+        std::visit(visitor, model);
+
+    model_props *= weights[i];
+    bsp += model_props;
+    i++;
+  }


+  XMLTag tag(type_name, "name", name, "nelem", x.models.size());
+  tag.write_to_stream(os);


+  # Plain solid-angle hemisphere: weights sin(za) dza daa, sums to ~2*pi sr
+  za = np.arange(1, 180, 2.0)
+  aa = np.arange(0, 360, 2.0)
+  ws.zen_grid = arts.ZenGrid(za.tolist())
+  ws.az_grid  = arts.AziGrid(aa.tolist())
+  ws.zen_grid_weights = arts.ZenGrid((np.sin(np.radians(za)) * np.radians(2.0)).tolist())
+  ws.az_grid_weights  = arts.AziGrid(np.full(len(aa), np.radians(2.0)).tolist())


The total internal reflection case compared the specular BRDF against 1.0 and 0.0 with exact equality. That value is only mathematically 1.0: fresnel() reaches TIR through a complex sqrt and a complex division, so std::norm(Rv) is 1 within a few ulp and the residual depends on the platform libm. Apple's complex sqrt leaves a spurious ~8.7e-17 real part in the transmitted cosine, which perturbs the division just enough that the norm lands exactly on 1.0. glibc's correctly rounded csqrt returns exactly 0 there and yields 0.99999999999999989; MSVC differs again. The test therefore failed on every Linux and Windows CI job while passing on both macOS jobs. - Compare BRDF[0,0], BRDF[1,0] and the TIR emissivity against a 1e-12 tolerance, as the rest of the file already does - Correct the C++ comment claiming fresnel() carries a total-internal-reflection guard; it never did Resolves the only failing test in PR atmtools#1194 CI run 38880886530.

riclarsson · 2026-10-19T10:37:34Z

@m-brath I pushed some code to deal with common spectral data. Check that PR if you want it or drop these two commits.

It allows you to set a SpectralGeodeticField3 and the interpolation is done along any grids of interest. The type also uses the same approach as the Numeric and FunctionalData that Atm/Subsurf/Surf fields allow, but now one dimension is frequency. That makes it possible to pass any method def n(f, lat, lon) -> float or numeric to avoid all interpolation.

There are a few other things about this PR I would have a look at, but first I want you do judge this change is in your design target.

… properties and update related methods

…ty field

… bindings * Add surface scattering framework properties example

…hon bindings.

…tory number.

…polation

…atterer

… surface scatterer

…rfaceScattererField

…mbertian surface scatterer implementation

…ty tensors in surface scattering models

…attering models, and update Python bindings/tests accordingly

…ings with weighting options

…date Python bindings with class annotations

…enda` to support surface scattering calculations

…ty tensors and implement spectral radiance scattering calculations

…on and correct test method call

…alculate spectral radiance components for surface scattering

…a` in workspace methods and refactor variable declarations for surface scattering calculations

…or tensor initializations in surface scattering models.

StokvecTensor4 was redundantly declared in workspace_groups.cpp, shadowing its existing entry in workspace_group_friends.cpp and failing diagnostics-wsg-vs-wsv.

… remove unused import

`src/m_surface.cc` — spectral_radSurfaceScatteringDiffuse and ...Specular take a `suns` input and an `exclude_suns` gin; diffuse skips quadrature directions that hit_sun(), specular skips the incoming-agenda trace on the mirror LOS when it hits the sun (pre-zeroed radiance/Jacobians keep emission intact). `src/workspace_methods.cpp` — both registry entries add `suns` to `.in` and an `exclude_suns` gin (Index, default 0) with gin_desc, plus desc paragraphs on de-duplication when a ...Direct method is chained. `src/workspace_agendas.cpp` — spectral_rad_surface_agenda desc rewritten: non-Direct methods run with exclude_suns = 1, each sun counted once by the Direct methods, hand-written chains must set the gin themselves. `src/workspace_agenda_creator.cpp` — predefined SurfaceScatteringModel option passes SetWsv("exclude_suns", Index{1}) to the Diffuse and Specular adds. `tests/core/surf/spectral_rad_surface_scattering_sun_double_count.py` — new test (7 cases) pinning both gates against closed forms, subsolar full-chain de-duplication to factor 1, no-op when no traced direction holds the sun, and zero scattering Jacobians for pure reflectors. `tests/core/surf/spectral_rad_surface_scattering_chain.py` — adds `ws.suns` setups and a hand-written gated chain test with closed forms. `tests/core/surf/spectral_rad_surface_scattering_predefined_agendas.py` — on-sun geometry for the sum identity, method filter widened to ("_", "@") prefixes to match @exclude_suns preset inputs. `tests/core/surf/spectral_rad_surface_scattering_diffuse.py` — base workspace sets `ws.suns = []` for the new input. `tests/core/surf/spectral_rad_surface_scattering_specular.py` — base workspace sets `ws.suns = []` for the new input.

- `src/core/surface_scattering/lambertian.h` — drop erroneous /π from the BRDF formula, fix interp_extrapolation comment to refer to the frequency grid - `src/python_interface/py_surface_scattering.cpp` — emissivity docstrings now show [nf, nza_scat, naa_scat, 4, 4], interp_extrapolation docstring refers to the frequency dimension - `src/core/surface_scattering/surface_scattering_properties.h` — emissivity tensor shape comments corrected to the full Mueller-matrix dims - `src/core/surface_scattering/surface_scattering_model.h` — stale method name replaced with get_surface_scattering_model_properties - `src/workspace_groups.cpp` — same stale method name fixed in the group description - `src/workspace_variables.cpp` — remove broken :doc: link from the surface_models description - `tests/core/surf/spectral_rad_surface_scattering_specular.py` — header comment reworded, specular-capable models do exist - `src/core/surface_scattering/lambertian.cc` — delete dead commented-out emissivity declaration No functional change; comments and documentation only.

…ation - `src/m_surface.cc` — gate Direct-method sub-surface emission behind new `include_emission` gin, weight diffuse Direct scattering by cos(theta_inc) at the actual surface normal, validate quadrature weight sizes, and fix three stream-style `ARTS_USER_ERROR_IF` calls that dropped their messages. - `src/workspace_methods.cpp` — declare `include_emission` (default 1) on both Direct methods and document the single-emission contract in all four method .desc strings. - `src/core/surface_scattering/lambertian.cc` — zero-initialize isotropic BRDF/emissivity Muelmat (identity default leaked reflectance 1 into Q/U/V), adopt new extrapolation handling. - `src/core/surface_scattering/fresnel.cc` — fold za > 90 onto 180-za before Fresnel evaluation (raw angles inverted/blew up reflectance). - `src/core/surface_scattering/scattering_internal.h` — replace `frequency_extrap_limit` with `extrap_query_grid`/`extrap_postprocess`/ `grid_range` so Nearest/None/Zero actually behave as documented. - `src/workspace_agenda_creator.cpp` — set `include_emission=0` on the Direct methods in the predefined `SurfaceScatteringModel` chain. - `src/core/surface_scattering/flat_scalar.cc` — apply the new extrapolation handling to reflectivity interpolation. - `src/core/surface_scattering/surface_scattering_model.cc` — user error on negative surface masks, static_assert instead of dead `if constexpr`, host `operator+=`. - `src/core/surface_scattering/lambertian.h` — document the dimensionless-kernel BRDF convention and required cos(za)dOmega weights, explicit ctor, per-mode extrap docs. - `src/core/surface_scattering/fresnel.h` / `flat_scalar.h` — validate longitude grid in field setters, per-mode extrap docs. - `src/core/surface_scattering/surface_scattering_model.h` — accurate Weighting and model-accumulation docs. - `src/core/surface_scattering/surface_scattering_properties.h` — remove stale docs and unused include. - `src/python_interface/py_surface_scattering.cpp` — bind missing `interp_extrapolation` on plain `LambertianSurfaceScatterer`, correct Weighting enum docs. - `src/core/surface_scattering/CMakeLists.txt` — link `rtepack`. - `tests/core/surf/interp_extrapolation.py` — new test covering all four extrapolation modes. - `tests/core/surf/spectral_rad_surface_scattering_specular.py` — closed-form Fresnel mirror-direction (I+Q) and mirror-azimuth gate tests. - `tests/core/surf/spectral_rad_surface_scattering_diffuse_direct.py` — closed-form oblique-sun cos(theta_inc) tests. - `tests/core/surf/spectral_rad_surface_scattering_chain.py` — invert assertions that enshrined the doubled emission. - `tests/core/surf/spectral_rad_surface_scattering_predefined_agendas.py` / `…_sun_double_count.py` — additivity identity corrected to full == DiffuseOnly + DirectOnly - emission. - `tests/core/surf/spectral_rad_surface_scattering_diffuse.py` — weight-mismatch and unusable-sun user-error tests. - `tests/core/surf/lambertian_scatterer_field.py` — full 4x4 Mueller structure regression test for the identity-default bug. - `tests/core/surf/fresnel_surface_scatterer.py` — za > 90 fold regression test. - `.gitignore` — ignore `/.plan/`. Behavior changes: the predefined SurfaceScatteringModel agenda and chained Direct methods now add sub-surface emission exactly once (previously twice), and Nearest interp_extrapolation now edge-clamps as documented.

`src/core/surface_scattering/surface_scattering_model.h` — switch `MapOfSurfaceScatteringModel::models` from `std::map` to `std::unordered_map` for O(1) name lookup, swap the include, and note in the class doc that iteration order is now unspecified.

… of zero src/core/surface_scattering/surface_scattering_properties.h — wrap the four BRDF/emissivity members in std::optional; missing component means zero contribution for consumers src/m_surface.cc — guard each of the four RTE consumer paths and treat a missing component as a zero contribution src/core/surface_scattering/surface_scattering_model.cc — null-aware operator+= / operator*= (null + value = value); empty model map returns all-null properties src/core/surface_scattering/lambertian.cc — stop allocating zero specular tensors; leave specular components as std::nullopt src/core/surface_scattering/fresnel.cc — stop allocating zero diffuse tensors; leave diffuse components as std::nullopt src/core/surface_scattering/flat_scalar.cc — stop allocating zero diffuse tensors; leave diffuse components as std::nullopt src/python_interface/py_surface_scattering.cpp — document that unused members surface as None in Python tests/core/surf/lambertian_scatterer_field.py — assert unused specular components are None instead of zero tensors tests/core/surf/fresnel_surface_scatterer.py — assert unused diffuse components are None instead of zero tensors tests/core/surf/flat_scalar_surface_scatterer.py — assert unused diffuse components are None instead of zero tensors Breaking-Change: SurfaceScatteringModelProperties components that no model provides are now None in Python; code indexing previously-zero unused tensors must handle None.

Replace `Index` with `size_t` for loop variables to ensure compatibility with unsigned sizes in `f_grid` and `data`.

Designated initializers of SurfaceScatteringModelProperties in fresnel.cc, flat_scalar.cc, and lambertian.cc omitted the std::optional fields they do not provide, triggering clang's -Wmissing-designated-field-initializers. Initialize those fields explicitly with std::nullopt; no semantic change.

`python/doc/gen_pyarts_arts_autoclass_rst.py` — validate every attribute type marker against pyarts3, builtins and intersphinx roots so unresolvable names become generator errors instead of broken Sphinx links; keep nested enums so `MapOfSurfaceScatteringModel.Weighting` gets its own autoclass block; write module pages that clash with a same-cased sibling class as `<name>-module.rst` (and delete the stale file) so `pyarts3.arts.sun` survives on case-insensitive filesystems; reword the missing-type-info error so the `class-information` placeholder can't be pasted again. `src/python_interface/py_surface_scattering.cpp` — replace 13 literal `.. :class:`class-information`` placeholders with the real types (MuelmatTensor5/3 optionals, InterpolationExtrapolation, float, MapOfSurfaceScatteringModel.Weighting). `src/python_interface/py_star.cpp` — fix `sun.refractive_los` ws parameter type from nonexistent `~pyarts3.arts.Workspace` to `~pyarts3.workspace.Workspace`. Sphinx build warnings drop from 37 to 4 (remaining are pre-existing).

examples/7-surface-scattering-properties/1-surface-scattering-checkerboard.py -- new example driving the surface scattering framework as the lower boundary condition of a full simulation: a MapOfSurfaceScatteringModel keyed by SurfacePropertyTag masks on surf_field, two Lambertian scatterers with opposite reflectivity spectra painted into a 4 x 8 checkerboard of 45 deg cells, the sun reflected as a collimated beam through the SurfaceScatteringModelDirectOnly agenda option, and a nadir Dirac sensor at 36 000 km. It plots the masks, the reflectivity spectra with their interpolation onto f_grid, and the per-band nadir radiance, and pins the spectra plus the full 4 x 8 x 3 scan to reference arrays, so it doubles as a regression test.

- `examples/7-surface-scattering-properties/2-surface-scattering-checkerboard-atm-lut.py` — new solar-reflective example: red/blue Lambertian checkerboard seen from 36 000 km through a tropical clearsky atmosphere with absorption from an abs_lookup_data table and first-order solar Rayleigh scattering lighting the sky, at 400/550/650 nm bands, with reference-array regression asserts. - `examples/7-surface-scattering-properties/3-surface-scattering-checkerboard-atm-lut-ir.py` — new thermal-IR counterpart: no sun or Rayleigh, window wavenumbers 895/1145 /cm straddling the O3 9.6 µm band, Kirchhoff emissivity (1−r) spectra, quadratic latitudinal surface temperature (300→250 K, +10 K on red cells), Dirac sensor at 800 km with SurfaceScatteringModelDiffuseOnly agenda, linear radiance and Planck brightness-temperature plots, with regression asserts. - `examples/7-surface-scattering-properties/README.rst` — list and briefly describe both new examples.

…ctives - python/doc/source/_ext/plot_source_code.py - new extension patching matplotlib's plot_directive._run_code to compile example code with its real filename, so inspect-based ARTS decorators (arts_agenda, callback_operator) work inside plot blocks - python/doc/source/conf.py - register the plot_source_code extension

examples/7-surface-scattering-properties/3-surface-scattering-checkerboard-atm-lut-ir.py — remove O3-666 from the absorbing species, rewrite the ozone-band narrative/annotations and assertions around H2O damping, and regenerate the reference arrays examples/7-surface-scattering-properties/README.rst — drop the ozone-band wording from the IR example description

… scene `examples/7-surface-scattering-properties/4-surface-scattering-checkerboard-camera.py` — new single-band 550 nm thin-lens camera example photographing the checkerboard scene through the atmosphere (LUT absorption, Rayleigh scattering), regression-pinned `examples/7-surface-scattering-properties/README.rst` — list the new camera example

…cattering framework - doc/arts/concept.surface_scattering.rst — new page developing the BRDF lower-boundary physics: polarimetric boundary condition, Lambertian/Fresnel/FlatScalar model families, Kirchhoff consistency, and hemisphere sampling. - doc/arts/user.surface_scattering.rst — new user page on the interface: surface_models map, surf_field SurfacePropertyTag mask weighting, spectral_rad_surface_agenda wiring, with a graphviz overview diagram. - doc/arts/dev.surface_scattering.rst — new developer page mapping the implementation: source layout, model contract, and how to add a new model. - doc/arts/concepts.rst — add concept.surface_scattering to the toctree. - doc/arts/users.rst — add user.surface_scattering to the toctree. - doc/arts/devs.rst — add dev.surface_scattering to the toctree.

The total internal reflection case compared the specular BRDF against 1.0 and 0.0 with exact equality. That value is only mathematically 1.0: fresnel() reaches TIR through a complex sqrt and a complex division, so std::norm(Rv) is 1 within a few ulp and the residual depends on the platform libm. Apple's complex sqrt leaves a spurious ~8.7e-17 real part in the transmitted cosine, which perturbs the division just enough that the norm lands exactly on 1.0. glibc's correctly rounded csqrt returns exactly 0 there and yields 0.99999999999999989; MSVC differs again. The test therefore failed on every Linux and Windows CI job while passing on both macOS jobs. - Compare BRDF[0,0], BRDF[1,0] and the TIR emissivity against a 1e-12 tolerance, as the rest of the file already does - Correct the C++ comment claiming fresnel() carries a total-internal-reflection guard; it never did Resolves the only failing test in PR atmtools#1194 CI run 38880886530.

riclarsson · 2026-10-21T02:02:50Z

@m-brath If you consider something a bug in the future, please fix it rather than leave your opinion as a comment. It was a useful approximation for the common $P=\Pi_iT_i$ case but it causes too many issues elsewhere. There is no reason to put warnings and just calling it a bug all over your comments.

riclarsson · 2026-10-21T02:41:11Z

+def setup(app):
+    plot_directive._run_code = _run_code
+    return {"parallel_read_safe": True, "parallel_write_safe": True}


OK, this is extremely scary. Why are you overriding matplotlib? Is there no workaround than to claim they are wrong? It seems like a really bad idea to monkeypatch external libraries.

I reverted this, it's not needed. @m-brath Install pyqt in your mamba environment and the plots should work.

Correction: This is a real issue. It affects all examples that use the agenda decorator. I'll see if there is a better fix than the monkeypatching. Example Sphinx error:

/Users/u237002/Hacking/uhh/arts/build/python/doc/source/examples.1-basics.rst:58: WARNING: Exception occurred in plotting 5-create-an-agenda from /Users/u237002/Hacking/uhh/arts/build/python/doc/source/../../../examples/1-basics/5-create-an-agenda.py: Traceback (most recent call last): File "/Users/u237002/miniforge3/envs/py314/lib/python3.14/site-packages/matplotlib/sphinxext/plot_directive.py", line 594, in _run_code exec(code, ns) ~~~~^^^^^^^^^^ File "<string>", line 34, in <module> File "/Users/u237002/Hacking/uhh/arts/build/python/src/pyarts3/workspace/workspace.py", line 526, in arts_agenda return _arts_agenda(func, ws, fix) File "/Users/u237002/Hacking/uhh/arts/build/python/src/pyarts3/workspace/workspace.py", line 388, in _arts_agenda srccod = getsource(f) File "/Users/u237002/miniforge3/envs/py314/lib/python3.14/inspect.py", line 1166, in getsource lines, lnum = getsourcelines(object) ~~~~~~~~~~~~~~^^^^^^^^ File "/Users/u237002/miniforge3/envs/py314/lib/python3.14/inspect.py", line 1148, in getsourcelines lines, lnum = findsource(object) ~~~~~~~~~~^^^^^^^^ File "/Users/u237002/miniforge3/envs/py314/lib/python3.14/inspect.py", line 983, in findsource raise OSError('could not get source code') OSError: could not get source code [docutils]

riclarsson · 2026-10-21T02:43:41Z

+def fresnel_amplitudes(n1, n2, za_deg):
+    """Closed-form complex Fresnel amplitude coefficients (theta in degrees).
+
+    Mirrors the physics_funcs.cc pair overload.  The TIR branch is the exact
+    limit |Rv| = |Rh| = 1; the C++ reaches it only within a few ulp, via its
+    complex sqrt and division, so never bit-exactly.
+    """
+    th = math.radians(za_deg)
+    sin2 = n1 * math.sin(th) / n2
+    if abs(sin2) > 1.0:
+        return 1.0, 1.0  # total internal reflection
+    cos1 = math.cos(th)
+    cos2 = math.cos(math.asin(sin2))
+    Rv = (n2 * cos1 - n1 * cos2) / (n2 * cos1 + n1 * cos2)
+    Rh = (n1 * cos1 - n2 * cos2) / (n1 * cos1 + n2 * cos2)
+    return Rv, Rh
+
+
+def mun_fresnel_reflectance(Rv, Rh):
+    """Mirror of rtepack::fresnel_reflectance (rtepack_surface.cc)."""
+    rv = abs(Rv) ** 2
+    rh = abs(Rh) ** 2
+    rmean = 0.5 * (rv + rh)
+    rdiff = 0.5 * (rv - rh)
+    c = (Rh * Rv.conjugate()).real
+    d = 0.5 * (Rh * Rv.conjugate() - Rv * Rh.conjugate()).imag
+    M = np.zeros((4, 4))
+    M[0, 0] = rmean
+    M[1, 1] = rmean
+    M[0, 1] = rdiff
+    M[1, 0] = rdiff
+    M[2, 2] = c
+    M[3, 3] = c
+    M[2, 3] = d
+    M[3, 2] = -d
+    return M


I do not like this as a test. Copy code from C++ to python to test if the C++ code is doing the same as the copied python code seems weird. More like a prototype step than an actual implementation of the test.

riclarsson · 2026-10-21T02:49:20Z

+void spectral_radSurfaceScatteringDiffuseDirect(const Workspace&                   ws,
+                                                StokvecVector&                     spectral_rad,
+                                                StokvecMatrix&                     spectral_rad_jac,
+                                                const AscendingGrid&               freq_grid,
+                                                const AtmField&                    atm_field,
+                                                const SurfaceField&                surf_field,
+                                                const SubsurfaceField&             subsurf_field,
+                                                const MapOfSurfaceScatteringModel& surface_models,
+                                                const JacobianTargets&             jac_targets,
+                                                const PropagationPathPoint&        ray_point,
+                                                const ArrayOfSun&                  suns,
+                                                const Agenda&                      ray_path_observer_agenda,
+                                                const Agenda&                      spectral_rad_incoming_agenda,
+                                                const Agenda&                      spectral_rad_closed_surface_agenda,
+                                                const Numeric&                     angle_cut,
+                                                const Index&                       refinement,
+                                                const Index&                       include_emission) try {


My AI says this "returns r·cosθ·L_sun with no solar solid angle. With ARTS's own convention (weights cos·dΩ/π, and π·sin²α as in m_sun.cc), it should be r·cosθ·sin²α·L_sun."

riclarsson · 2026-10-21T02:53:41Z

+      for (Index as = 0; as < nas; ++as) emissivity_specular[f, zs, as] = emissivity;
+  }
+}
+


I am not sure we are following this convention. I think that brdf should not be Muelmat{r} but mirror(Muelmat{r}). (That means that U and V should have negative signs.)

Can you confirm this is correct?

riclarsson · 2026-10-21T03:03:59Z

+    // Write the variant model
+    xml_write_to_stream(os, model, pbofs);
+  }
+


Why is this necessary? If both keys and values are already supported, unordered maps should just work? I think it would be good to fix the bug if that is not the case. Besides that, if you make the enum an arts_option, you can just claim the entire struct is an xml_io_stream_aggregate to write none of this code at all. I feel it is overly complicated to write this all manually.

riclarsson · 2026-10-21T03:09:39Z

+        try:
+            lambertian_brdf_i(sc, f_query)
+        except RuntimeError as e:
+            assert NONE_ERROR in str(e), f"Unexpected error message: {e}"


I really do not like this way of testing. It is something AI agents do all the time. Please be aware that this tests says "You are never allowed to change the error message to anything more user-friendly, this is the most user-friendly message possible" and "please use our error messages logic in your own scripts". I think neither claim should ever be made.

riclarsson · 2026-10-21T03:09:59Z

+    try:
+        lambertian_brdf_i(sc, f_query)
+    except RuntimeError as e:
+        assert NONE_ERROR in str(e), f"Unexpected error message: {e}"


riclarsson · 2026-10-21T03:10:13Z

+    try:
+        specular_i(sc_none, f_query)
+    except RuntimeError as e:
+        assert NONE_ERROR in str(e), f"Unexpected error message: {e}"


OK, last on in this file, but see above

riclarsson · 2026-10-21T03:14:17Z

+    h = constants.h
+    c = constants.c
+    k_B = constants.k
+    return (2 * h * f**3 / c**2) / (np.exp(h * f / (k_B * T)) - 1)


Again, I do not like using non-C++ copies of things. Besides that, this spreads a bad version of Planck function. You should use np.expm1 instead since it is accurate for low falues.

riclarsson · 2026-10-21T03:21:49Z

+    for (Index zi = 0; zi < nzi; ++zi)
+      for (Index ai = 0; ai < nai; ++ai)
+        for (Index zs = 0; zs < nzs; ++zs)
+          for (Index as = 0; as < nas; ++as) brdf_specular[f, zi, ai, zs, as] = brdf;


I will leave this here as well because that makes the comments dense and easy to find.

Why do you prefer all of these loops over something like brdf_specular[f] = brdf? The latter is so much easier to read in my opinion.

riclarsson

I have a few places I think should be changed. Most of my comments are about style. (I did some of the comments via vscode and some via browser because of convenience, so do not trust github here.)

There are some big things.

Several places there are C++ code that seems to have been copied to become python code that is run to check that the copy-paste operation worked and the two produce the same results. I think these tests are nonsense tests.
Several python tests makes the current error messages canon. I think this is quite silly. It does not matter what you write in the error message. You might want to change it in the future if you find it bad for some reason. Please just check the error type if you want these negative tests.
While it is good python practice to write scripts as you do in the examples, I think it is bad from a tutorial POV. I cannot easily inspect what happens inside build_workspace(). And it is called only once in an otherwise easy to introspect script.
The examples also contain quite a lot of "handwritten" constants instead of referring to ARTS. Same with basic functions like nm2hz. If you need to repeat these that often, make them pyarts3 or pyarts3.arts callables, or even real methods (like measure_nadir - a very simply method that serves a clear purpose). Several other examples contain callables but they are either there to simplify special logic, to name things, or to plot them. Or because the examples/ are growing and we couldn't do the new way when the example was created.

riclarsson · 2026-10-21T03:23:26Z

+def planck(f, T):
+    """Planck function [W m-2 sr-1 Hz-1]."""
+    h = constants.h
+    c = constants.c
+    k_B = constants.k
+    return (2 * h * f**3 / c**2) / (np.exp(h * f / (k_B * T)) - 1)


Again, see the other file

riclarsson · 2026-10-21T03:26:52Z

+def fresnel_R(theta_deg, n1, n2):
+    """Unpolarized Fresnel reflectance mean and difference at incidence angle.
+
+    Mirrors physics_funcs::fresnel (Rv the p-form (n2 c1 - n1 c2)/(n2 c1 +
+    n1 c2), Rh the s-form) and rtepack::fresnel_reflectance, returning
+    (rmean, rdiff) = (0.5(|Rv|^2+|Rh|^2), 0.5(|Rv|^2-|Rh|^2)).
+    """
+    ti = np.radians(theta_deg)
+    ci = np.cos(ti)
+    ct = np.sqrt(1.0 - (n1 * np.sin(ti) / n2) ** 2)
+    Rv = (n2 * ci - n1 * ct) / (n2 * ci + n1 * ct)
+    Rh = (n1 * ci - n2 * ct) / (n1 * ci + n2 * ct)
+    rv, rh = abs(Rv) ** 2, abs(Rh) ** 2
+    return 0.5 * (rv + rh), 0.5 * (rv - rh)


Look, the comment even states that this is a copy. So this was never even a prototype

riclarsson · 2026-10-21T03:32:23Z

+def nm2hz(wavelength_nm):
+    """Wavelength [nm] to frequency [Hz]"""
+    return SPEED_OF_LIGHT / (np.asarray(wavelength_nm, dtype=float) * 1e-9)


I examples I am very strongly against using custom calls for existing ARTS infra. Please change this. If you strongly believe an wavlennm2hz or nm2hz is needed, add it in py_convert

riclarsson · 2026-10-21T03:33:08Z

+def nm2hz(wavelength_nm):
+    """Wavelength [nm] to frequency [Hz]"""
+    return SPEED_OF_LIGHT / (np.asarray(wavelength_nm, dtype=float) * 1e-9)


It seems you do strongly believe it is needed

riclarsson · 2026-10-21T03:39:17Z

+Created on Wed May 20 10:44:35 2026
+
+@author: u242031
+"""


Please do not leave this kind of information around. Especially for an example. I think in general it is a bad idea for all of pyarts3/arts3 to have this stuff.

riclarsson · 2026-10-21T03:52:04Z

+void flat_scalar_tensors(const auto&     r_data,
+                         Index           nf,
+                         Index           nzi,
+                         Index           nai,
+                         Index           nzs,
+                         Index           nas,
+                         MuelmatTensor5& brdf_specular,
+                         MuelmatTensor3& emissivity_specular) {


This breaks a implicit thing we have been doing in ARTS for a long while and I sort of want it fixed. In most other place (all?) ARTS tend to put inout and out as the first parameters of the function call.

Beyond that, I think it also should probably be using views instead of references to the tensors as these are quite handy when size does not change.

riclarsson · 2026-10-21T04:02:36Z

+
+  /// Refractive index of medium 1 (the medium the radiation propagates in).
+  /// Must be > 0.
+  Numeric n1{1.0};


This n1. Should it not also be a specular gridded value? Because it seems to me that n1=1 is a very space oriented view for what might sometimes be the atmosphere, and even in the atmosphere it is not constant by frequency.

With the new "Data" type, you can still keep at 1.0 as a default and the "interpolation" cost should be very small.

riclarsson · 2026-10-21T04:05:55Z

+  Muelmat isotropic_emissivity{0.0};
+
+  for (Index f = 0; f < nf; ++f) {
+    const Numeric r = std::clamp(r_data[f], Numeric{0}, Numeric{1});


Future tips for these weird templates. Either std::clamp<Numeric> or just 0.0 and 1.0 would work. Just as a nicer way to keep it all short for yourself.

Edit, github removed the template parameter...

riclarsson · 2026-10-21T04:16:52Z

  ${ARTS_WORKSPACE_INTERNALS_GROUPS}
 )
-target_link_libraries(workspace_internals PUBLIC workspace_core artscore scattering)
+target_link_libraries(workspace_internals PUBLIC workspace_core artscore scattering surface_scattering)


Why are we doing this here rather than linking these to artscore again? They both live in core/. This is of course beyond this PR, but it brings it to light

Reverts 64fd88e. It's not needed.

m-brath requested review from riclarsson and a balanced review from Copilot October 16, 2026 09:30

Copilot started reviewing on behalf of m-brath October 16, 2026 09:30

View session

Copilot AI reviewed Oct 16, 2026

View reviewed changes

riclarsson mentioned this pull request Oct 20, 2026

Rtepack scattering #1197

Merged

m-brath added 23 commits October 21, 2026 10:23

Add surface scattering module with Lambertian model and Python bindings

0799137

Rename bulk surface scattering properties to surface scattering model…

b605113

… properties and update related methods

Extend Lambertian surface scattering model to use spectral reflectivi…

cdf1144

…ty field

Add spatially-varying Lambertian surface scattering model with Python…

80659c6

… bindings * Add surface scattering framework properties example

Remove SurfacePropertyTag from Lambertian surface scatterer and Pyt…

788e1df

…hon bindings.

Rename example file for surface scattering framework to correct direc…

60e8771

…tory number.

Add adaptive longitude cyclers for Lambertian surface scatterer inter…

270bcd3

…polation

Add interpolation and extrapolation support for Lambertian surface sc…

873225f

…atterer

Add longitude grid validation and simplify cycler usage in Lambertian…

bc51b3d

… surface scatterer

Add frequency extrapolation limit handling and tests for LambertianSu…

9924cf0

…rfaceScattererField

Rename identity to grid_identity and mark unused parameters in La…

6a46656

…mbertian surface scatterer implementation

Update to MuelmatTensor5 and StokvecTensor3 for BRDF and emissivi…

3a0a201

…ty tensors in surface scattering models

Add weighting options and associated methods for combining surface sc…

7378900

…attering models, and update Python bindings/tests accordingly

Normalize weights in surface scattering models and extend Python bind…

af12ed1

…ings with weighting options

Add workspace variables for spectral radiance and angle grids, and up…

3a16a52

…date Python bindings with class annotations

Add surface_models workspace variable and `spectral_rad_incoming_ag…

21ad4bf

…enda` to support surface scattering calculations

Refactor surface scattering properties to use mandatory BRDF/emissivi…

95312ca

…ty tensors and implement spectral radiance scattering calculations

Fix initialization issue in surface scattering model property summati…

1a8335d

…on and correct test method call

Add spectral_radSurfaceScatteringFlatDiffuse workspace method and c…

e779e7c

…alculate spectral radiance components for surface scattering

Update spectral_rad_observer_agenda to `spectral_rad_incoming_agend…

106ce2d

…a` in workspace methods and refactor variable declarations for surface scattering calculations

Remove unused tensor alias, streamline output constraints, and refact…

00aa138

…or tensor initializations in surface scattering models.

Remove StokvecTensor4 group declaration

e65a85a

StokvecTensor4 was redundantly declared in workspace_groups.cpp, shadowing its existing entry in workspace_group_friends.cpp and failing diagnostics-wsg-vs-wsv.

Fix emissivity tensor indexing in lambertian scatterer field test and…

7734463

… remove unused import

m-brath and others added 18 commits October 21, 2026 10:23

fix(surface): use size_t for indexing in scattering_internal.h

9f04ed9

Replace `Index` with `size_t` for loop variables to ensure compatibility with unsigned sizes in `f_grid` and `data`.

Allow common spectral data

151c08f

Apply formatting

d97f30f

Fix stale text

b39e73c

riclarsson force-pushed the surface_scattering_dev branch from 0888e69 to b39e73c Compare October 21, 2026 01:58

riclarsson reviewed Oct 21, 2026

View reviewed changes

riclarsson requested changes Oct 21, 2026

View reviewed changes

olemke and others added 5 commits October 21, 2026 09:03

docs: drop local plot directive source patch

f3b7ed5

Reverts 64fd88e. It's not needed.

Fix bug in circular interpolation

8dbdde9

Update IGRF to use custom data

9260f30

Not using internals

66a45d1

Small fixes

a6b6034

		const Numeric za_eff = za_fold(za_inc_grid[zi]);
		const auto [Rv, Rh] = fresnel(Complex{n1}, Complex{n2}, za_eff);

		XMLTag tag(type_name, "name", name, "nelem", x.models.size());
		tag.write_to_stream(os);

Conversation

m-brath commented Oct 16, 2026

Surface scattering framework

What's new

Notable fixes along the way

Notes for reviewers

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Uh oh!

Choose a reason for hiding this comment

Uh oh!

riclarsson commented Oct 19, 2026

Uh oh!

riclarsson commented Oct 21, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

olemke Oct 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

riclarsson left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

riclarsson Oct 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

olemke Oct 21, 2026 •

edited

Loading

riclarsson Oct 21, 2026 •

edited

Loading