Numerical verification against R
==================================

Every function in ``toolsrtm`` and ``scopeinpython`` has been ported
from, and verified directly against, the corresponding ``ToolsRTM``/
``SCOPEinR`` R function. Each has a ``pytest`` regression test comparing
its output against reference values generated by running the R package
itself (``python/scratch/scratch_export.R``) and saving the results to
CSV. Reference CSVs are bundled under each package's ``tests/refdata/``.

.. code-block:: bash

   cd python/toolsrtm && python -m pytest tests -q   # 32 tests
   cd python/scopeinpython && python -m pytest tests -q   # 18 tests

All 50 tests pass, and every ported function matches its R counterpart
to floating-point precision or better:

.. list-table::
   :header-rows: 1
   :widths: 30 40 30

   * - Function
     - Reference call
     - Max relative diff (Python vs R)
   * - ``calctav``
     - ``ToolsRTM::calctav``
     - 0 (exact)
   * - ``prospect_d`` (refl/tran)
     - ``ToolsRTM:::prospect_DB``
     - 5.8e-14
   * - ``prospect_pro`` (refl/tran)
     - ``ToolsRTM::prospect_PRO``
     - ~1e-14
   * - ``dladgen``, ``campbell``, ``volscatt``
     - ``ToolsRTM::dladgen/campbell/volscatt``
     - 0 (exact)
   * - ``foursail``
     - ``ToolsRTM::foursail`` (PROSPECT-PRO)
     - 3.4e-14
   * - ``foursail2``
     - ``ToolsRTM::foursail2`` (two-layer green/brown canopy)
     - ~1e-15
   * - ``inform`` (TypeLidf=1 and 2)
     - ``ToolsRTM::inform`` (PROSPECT-D)
     - ~1e-15
   * - ``spart_toc`` + its BSM soil step
     - ``ToolsRTM::SPART`` (TOC portion) / ``getBSM.toolsRTM``
     - ~1e-15
   * - ``get_smac``, all 9 atmospheric terms, Sentinel-2A (13 bands)
     - ``ToolsRTM::get.smac``
     - < 1e-6
   * - ``spart_toa``, TOA/TOC reflectance and radiance, Sentinel-2A (13 bands)
     - ``ToolsRTM::SPART`` (full pipeline)
     - < 1e-6
   * - ``marmit2`` (via ``get_marmit_rsoil(version='marmit2')``)
     - ``ToolsRTM::get.marmit.rsoil`` (``version='marmit2'``)
     - ~5e-7 (E1 exponential-integral: R integrates numerically, ``scipy.special.exp1`` uses a rational approximation)
   * - ``get_indices``, all 3 domains
     - ``ToolsRTM::getIndices``
     - ~5e-10 (worst: ``CR.red.nir.2``, an R-vs-scipy linear-extrapolation rounding difference; most indices agree to ~1e-14)
   * - ``fluspect_b`` (refl/tran/kChlrel/MbI/MbII/MfI/MfII)
     - ``ToolsRTM::getFluspect.B``
     - ~1e-15 (refl/tran), ~1e-19 (fluorescence matrices)
   * - ``fluspect_cx`` (refl/tran/kChlrel/kCarrel/Mb/Mf)
     - ``ToolsRTM::getFluspect.Cx``
     - ~1e-15 (refl/tran), ~1e-19 (fluorescence matrices)
   * - ``liberty`` (refl/tran/RR)
     - ``ToolsRTM::liberty``
     - 6.9e-12
   * - ``foursail``/``foursail2``/``inform`` wired with Liberty, Fluspect-B, Fluspect-B-Cx
     - same R functions with those ``LeafModel``\ s
     - ~5e-12 (Liberty), ~1e-15 (Fluspect)
   * - ``get_bsm`` / ``soilwat``
     - ``SCOPEinR::getBSM``
     - 4.2e-15
   * - ``run_rtmo``
     - ``SCOPEinR:::getRTMo``
     - 2.7e-14
   * - ``get_biochemical``, C3 and C4 (both temperature-corrected)
     - ``SCOPEinR::get.biochemical``
     - ~1e-13 (worst field ``rcw``; most ~1e-15 to 1e-18)
   * - ``get_fluspect_cx_scope``, step=5 and step=1
     - ``SCOPEinR::getFluspect.Cx.SCOPE``
     - ~1e-15 (refl/tran/kChlrel/kCarrel), ~1e-19 (Mb/Mf)
   * - ``fluspect_mscope`` (3-profile-layer, 10-canopy-layer case)
     - ``SCOPEinR::get.fluspect_mSCOPE``
     - ~1e-15 (refl/tran/kChlrel/kCarrel/phiI/phiII), ~1e-19 (Mb/Mf)
   * - ``rtmf``, native-grid (pre-spline) quantities
     - ``SCOPEinR::get.RTMf``
     - ~1e-6 (loop-accumulated over 30 layers)
   * - ``rtmf``, interpolated ``LoF_``/``EoutF_`` (211-pt display grid)
     - ``SCOPEinR::get.RTMf``
     - ~1.5e-3 absolute worst case (documented spline approximation, see :doc:`scopeinpython/rtmf`)
   * - ``rtmz``, ``rso``/``rdo``/``Eout_`` corrections (500-600nm)
     - ``SCOPEinR::get.RTMz``
     - ~1e-5
   * - ``get_scope``, exact-formula outputs (``LAIsunlit``/``Pnsun_Cab``/``Pnsha_Cab``/``Pntot_Cab``/TOC ``refl``)
     - ``SCOPEinR::get.SCOPE`` (package's own bundled example LUT)
     - ~1e-9 to ~1e-13
   * - ``get_scope``, iterative-convergence outputs (temperatures, energy-balance totals, ``Ja``/``PNPQ``/``fqe``)
     - ``SCOPEinR::get.SCOPE``
     - ~1-2% (~5% worst case, ``fqe``) -- small floating-point divergence compounding over the per-layer iterative solve, not a functional discrepancy

Two things worth knowing about, not discrepancies:

1. ``run_rtmo``'s ``rso``/``refl`` legitimately produce ``NaN``/``inf`` at
   far-thermal wavelengths where ``Esun_`` underflows to ~0 -- this
   reproduces R's own behaviour at the same wavelengths.
2. ``soilwat``'s use of ``scipy.stats.poisson.pmf(round(mu), k)`` mirrors
   the original R call ``stats::dpois(round(mu), k)`` exactly, and is
   covered by the BSM regression test.

See ``python/README.md`` at the repo root for the full narrative version
of this table.