indicatrix 0.1.0

//! The spectral transport loop itself.
//!
//! [`trace_spectral_ray`] and [`trace_spectral_ray_with_finish`] (the public entry
//! points), [`trace_spectral_ray_with_finish_soa`] (S4: the same, but with the
//! `PlanesSoA32` intersection arena built once by the caller rather than once per
//! sample -- see its own doc comment), and [`trace_spectral_ray_inner`] (their shared
//! bounce loop), plus the per-bounce facet dispatch and
//! Russian-roulette/mode-coupling machinery they drive.
//!
//! P6 (2026-09-07): [`trace_spectral_ray_inner`] also builds the per-trace
//! [`ExitSplitCtx`] and threads it (alongside its own `radiance`/`path_pdf`
//! accumulators -- exit-event splits fold directly into those, no separate accumulator)
//! down into [`dispatch_bounce`] -- see `refraction`'s own top-of-file "Exit-event
//! spectral splitting" doc comment for the estimator this supports.
//! `enable_exit_splitting` mirrors `enable_internal_mode_coupling`'s established
//! pattern exactly: hardcoded `true` at every public entry point below, exercised as an
//! explicit A/B parameter only by this file's own `#[cfg(test)]` module (see
//! `exit_splitting_tests`).

use super::{
    NUM_CHANNELS,
    absorption::{apply_absorption, rotate_stokes_to_plane_of_incidence},
    camera::{FacetFinish, HitRecord, Ray},
    color::{
        apply_von_kries_white_balance, integrate_channels_to_xyz,
        integrate_channels_to_xyz_families,
    },
    environment::{EnvironmentSource, environment_white_balance, sample_environment_channel},
    intersect::{build_plane_soa, intersect_polyhedron_soa, shading_normal_near_edge},
    refraction::{
        BounceRefractionGeometry, ExitSplitCtx, RayMaterialContext, RayWavelengthCache,
        apply_partial_fresnel_bounce, apply_tir_bounce, build_ray_wavelength_cache,
        compute_bounce_refraction_geometry, entry_eigenmode_selection, poynting_dir_for_mode,
    },
    sampling::{MODE_COUPLING_STREAM, RUSSIAN_ROULETTE_STREAM, hash_u32},
    scattering::{ScatterStepOutcome, apply_frosted_bounce, try_scatter_step},
};
use crate::{
    geometry::plane::GpuFacetPlane,
    optics::{
        materials::{CrystalSystem, GemMaterial},
        polarization::StokesVector,
    },
};
use glam::Vec3;

/// Fix G (Part 1): builds the `N`-member wrapped hero-wavelength comb from a single
/// `hero_rand` in `[0, 1)`.
///
/// Per `GEMSTONE_RENDERING_BLUEPRINT.md`'s Algorithm 2, line 5
/// (`wavelengths[i] = 380 + fmod((lambda_hero - 380) + i*(400/N), 400)`, with
/// `lambda_hero = 380 + hero_rand*400` drawn over the FULL visible range rather than
/// just the first sub-band). Pulled out to a standalone, directly-testable function
/// (rather than left inlined in `trace_spectral_ray`) so the wraparound arithmetic
/// itself -- that every generated wavelength stays within `[380, 780]`, and that the
/// hero always lands at index 0 by construction -- can be unit tested in isolation,
/// without needing a full ray/material/planes setup. Uses a const-generic array
/// (`std::array::from_fn`) rather than returning a `Vec`, so calling it from
/// `trace_spectral_ray`'s hot per-sample path costs no heap allocation.
#[must_use]
pub fn wrapped_hero_wavelengths<const N: usize>(hero_rand: f32) -> [f32; N] {
    const SPECTRUM_SPAN: f32 = 780.0 - 380.0;
    let channel_width = SPECTRUM_SPAN / N as f32;
    let lambda_hero = hero_rand.mul_add(SPECTRUM_SPAN, 380.0);
    std::array::from_fn(|k| {
        let offset = (k as f32).mul_add(channel_width, lambda_hero - 380.0);
        380.0 + offset.rem_euclid(SPECTRUM_SPAN)
    })
}

/// Accumulates the environment radiance for a ray that exited or missed the gemstone
/// entirely, terminating `trace_spectral_ray`'s bounce loop. A direct extraction of the
/// pre-extraction inline block: the same per-channel `f32::mul_add` chain in the same
/// order, called exactly once (this is the loop's break path), so extracting it changes
/// nothing about the computed bits.
///
/// `StudioRig::new` (~40 sin/cos plus 16 vector normalizes) is constant across
/// an entire ray -- `light_yaw`/`light_pitch` don't vary per spectral channel -- so it
/// is built ONCE here, before the `NUM_CHANNELS` loop, rather than once per channel
/// inside `sample_studio_environment` as it used to be (measured: 204ns per
/// `StudioRig::new`, x8 redundant rebuilds per ray ~= the 1373ns/ray gap this fix
/// closes). Bit-identical: same rig values, computed once instead of eight times.
fn accumulate_miss_radiance(
    environment: EnvironmentSource<'_>,
    ray_dir: Vec3,
    lambdas: &[f32; NUM_CHANNELS],
    stokes: &[StokesVector; NUM_CHANNELS],
    radiance: &mut [f32; NUM_CHANNELS],
) {
    let studio_rig = match environment {
        EnvironmentSource::Studio {
            light_yaw,
            light_pitch,
            ..
        } => Some(crate::optics::studio_rig::StudioRig::new(
            light_yaw,
            light_pitch,
        )),
        EnvironmentSource::HdrMap(_) => None,
    };
    for k in 0..NUM_CHANNELS {
        let env_spectral =
            sample_environment_channel(environment, ray_dir, lambdas[k], studio_rig.as_ref());
        // P7 (GPU/CPU parity at the environment lookup): `StokesVector::intensity`
        // clamps `I` to >= 0 (`self.i.max(0.0)`) before it reaches the environment
        // sample here, matching `spectral_transport.wgsl`'s own `max(stokes[k].x,
        // 0.0)` at its equivalent miss/environment-lookup site -- negative `I` is
        // unphysical on either side. Kept as an explicit cross-reference rather than a
        // behavioural change: this call already went through `.intensity()`, not the
        // raw `.i` field, before this comment was added.
        radiance[k] = f32::mul_add(stokes[k].intensity(), env_spectral, radiance[k]);
    }
}

/// Dispatches ONE bounce's TIR / partial-reflect / refract event -- or, for a
/// `FacetFinish::Frosted` facet, `apply_frosted_bounce` replacement -- and then
/// applies internal mode re-coupling if that event turned out to be an internal
/// reflection. This is `trace_spectral_ray_inner`'s entire per-bounce dispatch, extracted
/// into one function (mutating `current_ray`/`inside_gem`/`is_extraordinary` through
/// their `&mut` references) purely to keep that loop -- and clippy's function-length
/// lint -- manageable; every floating-point operation below is exactly what a bare
/// inline three-way dispatch would perform, just written once instead of at two call
/// sites. `was_internal_reflection`'s formula is the SAME one already established
/// for the polished reflect arm (`inside_gem && is_extraordinary_update.is_none() &&
/// new_inside_gem == inside_gem`, evaluated on the PRE-bounce `inside_gem`): it is
/// unconditionally `true` whenever a TIR bounce (polished or frosted) fires, since TIR is
/// only ever geometrically reachable with `inside_gem == true` to begin with (`n1 > n2`
/// for the hero channel requires it -- see `compute_bounce_refraction_geometry`), and it
/// is `false` for any transmit/refract event, since those always flip `inside_gem`.
///
/// P2: `current_k` is the wave normal `k` (see `refraction`'s own "wave normal vs
/// Poynting direction" design note) -- read here as the INCOMING `k` for this bounce's
/// TIR/Fresnel dispatch, and updated in place to the OUTGOING `k'` alongside
/// `current_ray.dir` (`S'`). A `FacetFinish::Frosted` bounce collapses `k' == S'`: a
/// diffuse (cosine-weighted-hemisphere) bounce already depolarizes exactly like a
/// Henyey-Greenstein scatter event (see `try_scatter_step`'s own identical treatment),
/// so there is no coherent wave-normal-vs-Poynting distinction left to track past it --
/// this is a documented simplification, not something rule 3 of the P2 design note
/// covers (that rule is specifically about the POLISHED TIR/reflect/refract dispatch).
#[expect(
    clippy::too_many_arguments,
    reason = "already bundles every value fixed for the trace (ctx) or the bounce \
              (cache, geo) into the established context structs; the rest are this \
              bounce's own scalar inputs (hit_point, normal, current_plane_normal, \
              finish), the RNG stream identity (rng_seed, bounce), and the FOUR separate \
              pieces of per-ray state a bounce can mutate (stokes, path_pdf, \
              current_ray, current_k, inside_gem, is_extraordinary) -- collapsing those \
              into one struct would hide which of six independently-updated fields each \
              bounce kind actually touches behind a single opaque &mut, which is a worse \
              trade than the argument count; P6 adds `exit` (the per-trace \
              `ExitSplitCtx`, holding its own staging `split_radiance` accumulator --\
              see that struct's doc comment), used only by the \
              `apply_partial_fresnel_bounce` branch below -- the frosted/TIR branches \
              never reach an exit event"
)]
fn dispatch_bounce(
    ctx: &RayMaterialContext,
    cache: &RayWavelengthCache,
    geo: &BounceRefractionGeometry,
    hit_point: Vec3,
    normal: Vec3,
    current_plane_normal: Vec3,
    finish: FacetFinish,
    rng_seed: u32,
    bounce: u32,
    stokes: &mut [StokesVector; NUM_CHANNELS],
    path_pdf: &mut [f32; NUM_CHANNELS],
    exit: &mut ExitSplitCtx<'_>,
    current_ray: &mut Ray,
    current_k: &mut Vec3,
    inside_gem: &mut bool,
    is_extraordinary: &mut bool,
) {
    let (new_k, new_s, new_inside_gem, is_extraordinary_update, exact_p_o) =
        if finish == FacetFinish::Frosted {
            let (new_dir, new_inside_gem, is_extraordinary_update) = apply_frosted_bounce(
                ctx,
                geo,
                normal,
                *inside_gem,
                *is_extraordinary,
                rng_seed,
                bounce,
                stokes,
                path_pdf,
            );
            (
                new_dir,
                new_dir,
                new_inside_gem,
                is_extraordinary_update,
                None,
            )
        } else if geo.sin2_t > 1.0 {
            let (k_prime, exact_p_o) = apply_tir_bounce(
                ctx,
                geo,
                *is_extraordinary,
                *current_k,
                normal,
                stokes,
                path_pdf,
            );
            let s_prime = poynting_dir_for_mode(ctx, geo, k_prime, *inside_gem, *is_extraordinary);
            (k_prime, s_prime, *inside_gem, None, exact_p_o)
        } else {
            apply_partial_fresnel_bounce(
                ctx,
                cache,
                geo,
                *current_k,
                normal,
                current_plane_normal,
                *inside_gem,
                *is_extraordinary,
                rng_seed,
                bounce,
                stokes,
                path_pdf,
                exit,
                hit_point,
            )
        };
    let was_internal_reflection =
        *inside_gem && is_extraordinary_update.is_none() && new_inside_gem == *inside_gem;
    current_ray.dir = new_s;
    current_ray.origin = hit_point + new_s * 1e-4;
    *current_k = new_k;
    *inside_gem = new_inside_gem;
    if let Some(updated) = is_extraordinary_update {
        *is_extraordinary = updated;
    }
    // P1: `current_k` (just updated above to the OUTGOING `k'`) and `stokes` (already
    // carrying this bounce's own Mueller matrix) are the CURRENT bounce's own
    // post-event state -- exactly what the internal re-coupling draw below needs to
    // weight itself by the light's actual polarization relative to the NEW wave
    // normal's eigenbasis, still expressed in `current_plane_normal` (unrotated until
    // the next facet hit) -- see `maybe_apply_internal_mode_coupling`'s doc comment.
    // P2 exit wiring (2026-09-07): `exact_p_o` is the SAME closed-form `R_o/(R_o+R_e)`
    // split `apply_tir_bounce`/`apply_uniaxial_internal_bounce` already computed for
    // this bounce's own reflected energy (`Some` only when this bounce was a uniaxial
    // internal reflection; `None` for a frosted bounce, a biaxial material, or a
    // transmit/exit event, in which case `maybe_apply_internal_mode_coupling` falls
    // back to its pre-existing polarization-projection heuristic).
    maybe_apply_internal_mode_coupling(
        ctx,
        geo.is_biaxial,
        current_plane_normal,
        *current_k,
        stokes,
        was_internal_reflection,
        exact_p_o,
        rng_seed,
        bounce,
        is_extraordinary,
    );
}

/// Russian Roulette termination with weighted survival (Fix D). The previous hard
/// cutoff (`if max_intensity < 0.02 { break }`) simply discarded the remaining energy
/// outright once a path got dim, which biases the estimator dark -- worst on the long
/// internal bounce trains inside high-index stones. Instead, survive with probability
/// `q` (derived from the path's maximum spectral throughput) and, on survival, divide
/// every Stokes vector by `q` so the estimator stays unbiased in expectation
/// (`E[survive] * (1/q) = 1`). Returns `false` if the path should terminate (the
/// pre-extraction code's `break`), `true` if it survives (Stokes vectors already
/// rescaled in place, exactly as before).
///
/// P6 (Bias A fix, 2026-09-07): `split_radiance` is a THIRD radiometric quantity this
/// same survival draw governs, alongside `stokes` -- a channel that split off at an
/// earlier bounce (this bounce index or any prior one, `bounce > 4`) already folded its
/// fully-resolved contribution into `exit.split_radiance[k]` from the UNSCALED `stokes`
/// at that exit bounce, and `trace_spectral_ray_inner` commits it into `radiance` only
/// if the hero path goes on to survive every remaining Russian-roulette draw and
/// eventually escape. That commit happens with probability `q` per surviving draw --
/// exactly like every other quantity riding along on this same path -- so leaving
/// `split_radiance` un-rescaled while `stokes` gets `1/q` means a split contribution is
/// committed with probability `q` but never compensated by the matching `1/q`,
/// understating it by a factor of `q` at every exit that occurred at bounce index >= 5
/// (`bounce > 4` is when this function first becomes reachable). Previously measured as
/// a real, sample-count-growing bias (this crate's own `exit_splitting_tests` -- see
/// this task's report for the before/after numbers): scaling `split_radiance` by the
/// identical `1/q` on survival, in the identical branch, restores
/// `E[commit] * (1/q) = 1` for it too.
pub(super) fn apply_russian_roulette(
    bounce: u32,
    rng_seed: u32,
    stokes: &mut [StokesVector; NUM_CHANNELS],
    split_radiance: &mut [f32; NUM_CHANNELS],
) -> bool {
    let max_intensity = stokes.iter().fold(0.0f32, |a, s| a.max(s.intensity()));
    let q = max_intensity.clamp(0.05, 1.0);
    let rr_rand =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ RUSSIAN_ROULETTE_STREAM)) as f32) / 4_294_967_295.0;
    if rr_rand > q {
        return false;
    }
    for s in &mut *stokes {
        *s = s.scale(1.0 / q);
    }
    for r in &mut *split_radiance {
        *r /= q;
    }
    true
}

/// Stochastic o<->e (uniaxial) / mode-A<->mode-B (biaxial)
/// re-coupling at an INTERNAL reflection inside an anisotropic crystal.
///
/// # What this models, and what it does not
///
/// Today's model assigns a path's eigenmode ONCE, at the air->crystal entry
/// (`BIREFRINGENT_SPLIT_STREAM`), and every subsequent internal bounce reuses that same
/// mode's scalar index (`n_medium_ch` in `compute_bounce_refraction_geometry`) via an
/// ISOTROPIC Fresnel calculation. That is wrong in a specific way: a uniaxial crystal's
/// o/e eigenbasis is defined relative to the LOCAL wave-normal direction at each facet
/// (`ordinary_eigen_polarization`/`extraordinary_eigen_polarization`, both functions of
/// `wave_normal`), so a facet with a different orientation relative to `c_axis` has a
/// DIFFERENT eigenbasis than the facet the path last bounced off. A polarization state
/// that was purely in the old facet's o (or e) eigenmode is therefore, in general, a
/// SUPERPOSITION of the new facet's o' and e' eigenmodes -- real light partially converts
/// between the two labels at every internal bounce, which is a genuine contributor to
/// the "doubling" complexity high-birefringence stones (zircon, moissanite) show.
///
/// This function models that conversion as a fresh, UNCONDITIONAL 50/50 coin flip of
/// which mode's index governs the path from here on -- deliberately the same blanket
/// 50/50 the entry split already uses, not a per-bounce probability derived from the
/// actual angle between the old and new eigenbases (that would require projecting the
/// carried Stokes/eigenmode state onto the new facet's eigenbasis, i.e. genuine
/// anisotropic Fresnel coefficients at an anisotropic interface -- explicitly out of
/// scope, see this module's callers' doc comments). What this buys, honestly: internal
/// bounces stop being frozen to whichever mode the path entered as, so a
/// many-internal-bounce path (exactly the pavilion TIR trains responsible for a cut
/// stone's brilliance) now samples a genuine MIX of o-like and e-like propagation across
/// its life, instead of committing to one for its entire interior traversal. What it does
/// NOT buy: a physically-derived conversion FRACTION per bounce, or coherent
/// superposition -- this is still scalar Fresnel with an effective index at an
/// anisotropic interface, exactly as the entry split already was, just re-rolled more
/// often.
///
/// # This is a RELABELING, not a SPLIT -- and neither draw needs a `1/p` division
///
/// It is tempting to reach for a `1 / selection_probability` reweighting for this
/// draw, since it is keyed off an unconditional 50/50 coin flip just like the entry
/// split's own mode-selection draw in `apply_refract_bounce`. Neither draw actually
/// needs one, though -- for two structurally different reasons.
///
/// At the entry split, one incident ray genuinely becomes TWO physically distinct rays
/// (the o-ray and the e-ray), each carrying only its own ~0.5 share of the incident
/// energy. The selected mode's own full Fresnel transmittance is evaluated against the
/// FULL incident Stokes vector -- i.e. it always computes "what if ALL the incident
/// energy took this mode's index," which overstates that mode's true contribution by
/// exactly the 2x its ~0.5 energy share corrects for. Since the mode is drawn with
/// probability 0.5 -- the SAME 0.5 as its true energy share -- the plain, unscaled
/// sample is already an unbiased estimator of the average `0.5*T_o + 0.5*T_e`:
/// `E[g(X)] = 0.5*g(o) + 0.5*g(e)` IS that average directly, no `1/p` correction
/// needed, precisely because the target quantity is itself a probability-weighted
/// average and the sampling distribution already matches those weights. (A `1/p`
/// correction is the right tool when the target is instead a SUM of disjoint
/// contributions sampled at non-matching probabilities -- e.g. the reflect/refract
/// branch's `r_unpol`/`1 - r_unpol` split, where the two branches are mutually
/// exclusive fates for the FULL incident beam, not simultaneous energy-sharing
/// alternatives, so recovering the full per-branch magnitude genuinely does require
/// dividing back out by the branch's own selection probability.)
///
/// Here, nothing new is created. The path already has exactly one Stokes vector and one
/// `path_pdf` -- this function does not touch either of them beyond the label they are
/// filed under. All it does is re-roll which of the two eigenmode LABELS ("o" or "e")
/// governs the refractive index used for the path's NEXT bounce (see
/// `compute_bounce_refraction_geometry`'s `n_medium_ch` lookup). That is a
/// redistribution of an existing quantity's bookkeeping, not an estimate of a sum of two
/// quantities. The energy carried by the path before this call and after it is, and must
/// be, identical -- otherwise a lossless closed system (a furnace) would gain or lose
/// energy on every single internal reflection, purely from re-labeling which index
/// governs the next bounce, which is nonsensical: the physical light did not get
/// brighter or dimmer because this code changed its mind about which eigenmode name to
/// call it.
///
/// The reason the entry split's unscaled draw and this function's unscaled draw are
/// BOTH correct despite looking so similar: the entry split's `0.5` selection
/// probability is deliberately set to match the SAME `0.5` blanket energy-fraction
/// assumption used for each mode's true contribution, so the direct sample already
/// IS the true per-mode contribution to the average sum -- see the section above, no
/// division needed. Here there is no sum and no per-mode energy fraction to recover --
/// the `0.5` is purely which label applies to the one quantity already in hand, so the
/// matching, unbiased scale factor is `1.0`, i.e. no scaling at all -- the same `1.0`
/// conclusion as the entry split, reached for a different reason.
///
/// # Unbiasedness -- and why P1 can weight this draw too, unlike the entry split
///
/// The draw is independent of every channel (one shared coin flip drives the one shared
/// geometric path, exactly like `BIREFRINGENT_SPLIT_STREAM`/`FRESNEL_BRANCH_STREAM`), so
/// every channel's own hypothetical technique would draw the identical outcome from this
/// same `(rng_seed, bounce)`-keyed stream -- but since this is a relabeling and not a
/// split (see above), neither `stokes` (the actual radiometric quantity) nor `path_pdf`
/// (each channel's own technique-density bookkeeping for the final spectral MIS weight)
/// needs any adjustment: neither `stokes` (the actual radiometric quantity) nor
/// `path_pdf` (each channel's own technique-density bookkeeping for the final spectral
/// MIS weight) is touched here -- unlike the entry split, this function takes no
/// mutable `stokes`/`path_pdf` parameters at all, since it has nothing to do to them.
/// This has a second consequence P1 relies on: since NEITHER accumulator is touched
/// regardless of which label this draw picks, the draw's own SELECTION PROBABILITY
/// cannot introduce bias no matter what it is -- unlike the entry split (where `p_sel`
/// must equal each mode's true physical energy fraction for the unscaled sample to stay
/// unbiased, see [`crate::optics::raytracer::refraction::apply_partial_fresnel_bounce`]'s
/// own doc comment), there is no matching constraint here to violate. So P1 weighting
/// this draw by [`crate::optics::raytracer::refraction`]'s
/// `entry_eigenmode_selection` (same formula as the entry split, evaluated at the
/// bounce's OWN post-reflection wave normal and Stokes state) can only improve which
/// index the NEXT bounce actually uses -- it cannot make the estimator less unbiased,
/// since nothing here was ever unbiased-BECAUSE-OF a specific probability to begin
/// with. Biaxial materials keep the pre-P1 blanket 50/50 (`ordinary_eigen_polarization`
/// is uniaxial-only -- see `entry_eigenmode_selection`'s own doc comment).
///
/// Leaving `path_pdf` untouched is also a robustness win over the old `*= 0.5`: across
/// the long internal-bounce trains a TIR-heavy pavilion produces (tens to ~100 internal
/// reflections is not unusual), `0.5^n` underflows to exactly `0.0` well before `n`
/// reaches 100, which would have driven `sum(path_pdf) == 0` and tripped
/// `spectral_mis_weight`'s "should not happen" guard. Returns the freshly-selected
/// `is_extraordinary` value; the caller applies it only when the internal reflection it
/// just dispatched actually happened while `inside_gem && is_anisotropic`.
// P2 exit wiring (2026-09-07): this function's relabeling probability now uses the
// EXACT closed-form `internal_solve` Poynting-weighted `R_o/(R_o+R_e)` split whenever
// the caller has one to give (`exact_p_o: Some(..)`, threaded up from the SAME
// `internal_solve` call `apply_tir_bounce`/`apply_uniaxial_internal_bounce` already ran
// for this bounce's own reflected energy -- no extra boundary-value solve needed here).
// The pre-existing `entry_eigenmode_selection` polarization-projection heuristic below
// is now reached ONLY as the `exact_p_o == None` fallback: a biaxial material (out of
// this task's scope -- `BiaxialIndicatrix` has no uniaxial closed-form solve to draw
// an exact split from) or any other caller that has no per-bounce solve result on hand.
// Since this function's own doc comment above already establishes the draw's
// SELECTION PROBABILITY cannot introduce bias no matter what it is (neither `stokes`
// nor `path_pdf` is touched here), swapping in the exact value can only improve which
// label the next bounce uses -- see `mode_coupling_tests` for the re-derived
// expectations this exact split produces.
#[expect(
    clippy::too_many_arguments,
    reason = "thin re-labeling draw: ctx/is_biaxial/current_plane_normal/new_k/stokes \
              are P1's own inputs (see the doc comment above), exact_p_o is P2's \
              closed-form override of the same probability, plus the RNG stream \
              identity -- see the doc comment above for why each is needed"
)]
fn apply_internal_mode_coupling(
    ctx: &RayMaterialContext,
    is_biaxial: bool,
    current_plane_normal: Vec3,
    new_k: Vec3,
    stokes: &[StokesVector; NUM_CHANNELS],
    exact_p_o: Option<f32>,
    rng_seed: u32,
    bounce: u32,
) -> bool {
    let p_o = exact_p_o.unwrap_or_else(|| {
        if is_biaxial {
            None
        } else {
            entry_eigenmode_selection(
                ctx.c_axis,
                current_plane_normal,
                new_k,
                stokes[ctx.hero_idx],
            )
        }
        .map_or(0.5, |(p, ..)| p)
    });
    let split_rand =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ MODE_COUPLING_STREAM)) as f32) / 4_294_967_295.0;
    split_rand >= p_o
}

/// Applies [`apply_internal_mode_coupling`] in place on `is_extraordinary` exactly when
/// the bounce dispatch that just ran was an internal reflection inside an anisotropic
/// crystal; a no-op otherwise. `was_internal_reflection` is `true` unconditionally for a
/// caller reporting a TIR bounce (`apply_tir_bounce` is reachable only with `n1 > n2` for
/// the hero channel, i.e. only from inside the gem -- see
/// `compute_bounce_refraction_geometry`: `n1 == 1.0` outside it), and
/// `is_extraordinary_update.is_none() && new_inside_gem == inside_gem_before_bounce` for
/// a caller reporting the outcome of `apply_partial_fresnel_bounce` (its reflect arm
/// returns `is_extraordinary_update == None` and never changes `inside_gem`, unlike its
/// refract arm, which always flips `inside_gem` and always returns `Some`). Extracted --
/// and made void/in-place rather than returning `Option<bool>` -- purely to keep both
/// this function's own call sites and `trace_spectral_ray_inner` itself compact enough
/// for clippy's function-length lint.
#[expect(
    clippy::too_many_arguments,
    reason = "thin dispatcher for apply_internal_mode_coupling's own P1 inputs \
              (is_biaxial, current_plane_normal, new_k, stokes) plus the pre-existing \
              was_internal_reflection gate and RNG stream identity -- see that \
              function's own doc comment for why each is needed"
)]
fn maybe_apply_internal_mode_coupling(
    ctx: &RayMaterialContext,
    is_biaxial: bool,
    current_plane_normal: Vec3,
    new_k: Vec3,
    stokes: &[StokesVector; NUM_CHANNELS],
    was_internal_reflection: bool,
    exact_p_o: Option<f32>,
    rng_seed: u32,
    bounce: u32,
    is_extraordinary: &mut bool,
) {
    if ctx.enable_internal_mode_coupling && ctx.is_anisotropic && was_internal_reflection {
        *is_extraordinary = apply_internal_mode_coupling(
            ctx,
            is_biaxial,
            current_plane_normal,
            new_k,
            stokes,
            exact_p_o,
            rng_seed,
            bounce,
        );
    }
}

/// Why a traced path's bounce loop stopped.
///
/// Bookkeeping the `bounce_cost` benchmark harness (`examples/bounce_cost.rs`) needs to
/// answer "what fraction of paths hit the bounce cap", not something any existing caller
/// of [`trace_spectral_ray_inner`] cares about. Maps 1:1 onto the loop's four
/// `break`/fall-through sites in [`trace_spectral_ray_inner`]:
/// - [`Self::Escaped`]: the ray missed every facet (or exited the polyhedron) and
///   [`accumulate_miss_radiance`] sampled the environment.
/// - [`Self::ScatterAbsorbed`]: [`try_scatter_step`] returned
///   [`ScatterStepOutcome::ScatteredAndTerminated`] -- only reachable for a material with
///   nonzero `scattering_sigma_s`; plain Beer-Lambert absorption in a non-scattering
///   material never terminates a path on its own; it only dims `stokes` until
///   [`apply_russian_roulette`] eventually kills it (recorded as [`Self::RussianRoulette`]
///   below, not this variant).
/// - [`Self::RussianRoulette`]: [`apply_russian_roulette`] drew a kill.
/// - [`Self::HitCap`]: none of the above fired before the `for bounce in 0..max_bounces`
///   loop ran out of iterations.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PathTermination {
    /// The ray missed the gemstone (or exited it) and the environment was sampled.
    Escaped,
    /// A Henyey-Greenstein scattering event extinguished the path (absorbed).
    ScatterAbsorbed,
    /// Russian roulette's stochastic kill fired.
    RussianRoulette,
    /// The loop ran out of bounces (`max_bounces`) before any of the above happened.
    HitCap,
}

/// Traces a single spectral ray through the 3D gemstone polyhedron with:
/// 1. 8-Channel Stratified Hero Wavelength Spectral Sampling (HWSS)
/// 2. Full 4D Stokes-Mueller Polarized Wave Tracking (with TIR phase shift and Brewster extinction)
/// 3. Birefringent Extraordinary Ray Walk-Off & Optical Doubling
/// 4. Directional Pleochroic Beer-Lambert Absorption Tensors
///
/// # `hero_rand`
///
/// The hero-wavelength draw is now an explicit `[0, 1)` parameter rather than derived
/// internally from `rng_seed` via `hash_u32` -- callers that want it STRATIFIED across a
/// pixel's own sample sequence (every real production call site: the worker, the
/// viewer, and the exporter) compute it via `low_discrepancy_base2` +
/// `cranley_patterson_rotate` from the absolute sample index and pixel index
/// separately, information a single opaque `rng_seed` cannot carry (a hash's whole job
/// is destroying input structure, so a caller cannot recover a stratifiable value from
/// `rng_seed` alone once pixel and sample have already been mixed together into it).
/// Callers that don't care about stratification (most tests, and the Tier 2/3 GPU
/// self-tests that intentionally reproduce the OLD unstratified formula to stay
/// comparable to their own historical baselines) can simply pass
/// `(hash_u32(rng_seed) as f32) / 4_294_967_295.0`, reproducing the pre-Fix-4 behaviour
/// exactly. `rng_seed` itself is unchanged: it still seeds every per-bounce draw
/// (Fresnel branch, Russian roulette, birefringent split) exactly as before
/// only touches the hero-wavelength and (in its callers) pixel-jitter draws.
#[expect(
    clippy::too_many_arguments,
    reason = "the crate's public raytracing entry point: each parameter is one \
              independent, named input a caller supplies once per traced ray (which \
              ray, which geometry, which material, the bounce cap, the lighting \
              environment, the RNG seed, the hero-wavelength draw, an optional debug \
              output) -- unlike the per-bounce internal functions in this module tree, \
              this is called once per sample, not threaded through every bounce, and \
              wrapping it in a context struct would only make every existing caller \
              (apps, examples, GPU/CPU parity tests) build a struct instead of calling a \
              function for no reduction in what they actually supply"
)]
#[must_use]
pub fn trace_spectral_ray(
    initial_ray: Ray,
    planes: &[GpuFacetPlane],
    material: &GemMaterial,
    max_bounces: u32,
    environment: EnvironmentSource<'_>,
    rng_seed: u32,
    hero_rand: f32,
    primary_hit_out: Option<&mut Option<HitRecord>>,
) -> Vec3 {
    // Internal o<->e mode re-coupling is always on for the
    // public entry point -- see `trace_spectral_ray_inner`'s doc comment for the
    // mechanism. The `bool` parameter exists so this file's own `#[cfg(test)]` module
    // can render the SAME material/ray/seed with it forced off, as a same-crate,
    // same-code-path A/B comparison proving the effect is real (see
    // `mode_coupling_tests::internal_mode_coupling_changes_zircon_render`) -- it is
    // deliberately not part of the public signature above it. Girdle finish:
    // `&[]` for `facet_finishes` -- every facet index looks up
    // `FacetFinish::default() == Polished` (see that lookup at the top of the bounce
    // loop below), so this reproduces the exact pre-Task-2 all-polished behaviour for
    // every one of this function's existing callers, unchanged.
    //
    // S4: builds the `PlanesSoA32` arena fresh, exactly as this function always has --
    // a single traced ray is not a hot per-sample batch loop, so there is nothing to
    // amortize here. See [`trace_spectral_ray_with_finish_soa`] for the entry point a
    // caller tracing MANY samples against the same fixed `planes` should use instead.
    let plane_soa = build_plane_soa(planes);
    trace_spectral_ray_inner(
        initial_ray,
        planes,
        &plane_soa,
        &[],
        material,
        max_bounces,
        environment,
        rng_seed,
        hero_rand,
        primary_hit_out,
        true,
        true,
        None,
    )
}

/// [`trace_spectral_ray`] with an explicit
/// per-facet `finish`.
///
/// See [`FacetFinish`]'s doc comment for why this is a separate function (not a new
/// parameter on `trace_spectral_ray` itself) and for its CPU-only status.
/// `facet_finishes[i]` is looked up for `planes[i]`; a facet index with no corresponding
/// entry (a shorter slice, or `&[]`) defaults to `FacetFinish::Polished`, so passing
/// `&[]` here is exactly equivalent to calling `trace_spectral_ray`.
///
/// S4: a thin wrapper that builds the `PlanesSoA32` arena fresh and delegates to
/// [`trace_spectral_ray_with_finish_soa`] -- see that function's doc comment for the
/// entry point a caller tracing many samples against the same fixed `planes` (a
/// batch/frame) should use instead, to build the arena once rather than once per
/// sample.
#[expect(
    clippy::too_many_arguments,
    reason = "see trace_spectral_ray's own reason -- same public entry point, plus the \
              one additional per-facet finish slice this variant exists for"
)]
#[must_use]
pub fn trace_spectral_ray_with_finish(
    initial_ray: Ray,
    planes: &[GpuFacetPlane],
    facet_finishes: &[FacetFinish],
    material: &GemMaterial,
    max_bounces: u32,
    environment: EnvironmentSource<'_>,
    rng_seed: u32,
    hero_rand: f32,
    primary_hit_out: Option<&mut Option<HitRecord>>,
) -> Vec3 {
    let plane_soa = build_plane_soa(planes);
    trace_spectral_ray_with_finish_soa(
        initial_ray,
        planes,
        &plane_soa,
        facet_finishes,
        material,
        max_bounces,
        environment,
        rng_seed,
        hero_rand,
        primary_hit_out,
    )
}

/// S4: [`trace_spectral_ray_with_finish`] with the [`crate::simd::PlanesSoA32`] arena
/// built by the CALLER and passed in, rather than rebuilt from `planes` on every call.
///
/// For a batch/frame tracing many samples against the same fixed `planes` (measured:
/// ~10-20% of per-sample cost on a 57-plane cut), build `plane_soa` via
/// [`build_plane_soa`] ONCE outside the sample loop and pass the SAME reference to every
/// call. `planes` must be the exact slice `plane_soa` was built from -- this function
/// does not, and cannot cheaply, verify that itself. Results are bit-identical to
/// [`trace_spectral_ray_with_finish`]'s own (rebuild-every-call) behaviour, since
/// [`build_plane_soa`] is a pure function of `planes`' contents.
#[expect(
    clippy::too_many_arguments,
    reason = "see trace_spectral_ray's own reason -- same public entry point, plus the \
              per-facet finish slice and (S4) the caller-supplied `plane_soa` arena \
              alongside the `planes` slice it was built from"
)]
#[must_use]
pub fn trace_spectral_ray_with_finish_soa(
    initial_ray: Ray,
    planes: &[GpuFacetPlane],
    plane_soa: &crate::simd::PlanesSoA32,
    facet_finishes: &[FacetFinish],
    material: &GemMaterial,
    max_bounces: u32,
    environment: EnvironmentSource<'_>,
    rng_seed: u32,
    hero_rand: f32,
    primary_hit_out: Option<&mut Option<HitRecord>>,
) -> Vec3 {
    trace_spectral_ray_inner(
        initial_ray,
        planes,
        plane_soa,
        facet_finishes,
        material,
        max_bounces,
        environment,
        rng_seed,
        hero_rand,
        primary_hit_out,
        true,
        true,
        None,
    )
}

/// Instrumented variant of [`trace_spectral_ray_with_finish`] for the `bounce_cost`
/// benchmark harness.
///
/// Identical computation (same code path, `enable_internal_mode_coupling = true`, same
/// call to [`trace_spectral_ray_inner`]), plus how many bounces the path actually took
/// and why it stopped -- see [`PathTermination`]'s doc comment for the four possible
/// reasons and how "bounces taken" is counted (the `for bounce in 0..max_bounces` loop's
/// iteration index at the point of termination; `max_bounces` itself for
/// [`PathTermination::HitCap`]).
#[expect(
    clippy::too_many_arguments,
    reason = "see trace_spectral_ray's own reason -- same public entry point, plus the \
              per-facet finish slice and the `bounce_cost` harness's own instrumentation \
              this variant exists for"
)]
#[must_use]
pub fn trace_spectral_ray_with_finish_instrumented(
    initial_ray: Ray,
    planes: &[GpuFacetPlane],
    facet_finishes: &[FacetFinish],
    material: &GemMaterial,
    max_bounces: u32,
    environment: EnvironmentSource<'_>,
    rng_seed: u32,
    hero_rand: f32,
) -> (Vec3, u32, PathTermination) {
    let plane_soa = build_plane_soa(planes);
    let mut termination = (max_bounces, PathTermination::HitCap);
    let radiance = trace_spectral_ray_inner(
        initial_ray,
        planes,
        &plane_soa,
        facet_finishes,
        material,
        max_bounces,
        environment,
        rng_seed,
        hero_rand,
        None,
        true,
        true,
        Some(&mut termination),
    );
    (radiance, termination.0, termination.1)
}

/// Records the primary ray's first hit for the denoiser's guide buffers -- a
/// direct extraction from `trace_spectral_ray_inner` (same operations, same
/// order; out of line purely for clippy's line-count lint). Only bounce 0
/// writes: the ray at that point is still exactly the camera ray, and the
/// single traced geometric path is always driven by the hero channel, so this
/// unambiguously IS the hero channel's own first hit.
fn capture_primary_hit(
    primary_hit_out: &mut Option<&mut Option<HitRecord>>,
    bounce: u32,
    hit: Option<HitRecord>,
) {
    if bounce == 0
        && let Some(slot) = primary_hit_out.as_deref_mut()
    {
        *slot = hit;
    }
}

/// The facet's finish, with `Polished` for any index `facet_finishes` doesn't
/// cover -- a direct extraction from `trace_spectral_ray_inner` (same
/// operations, same order; out of line purely for clippy's line-count lint).
fn facet_finish_for(facet_finishes: &[FacetFinish], facet_idx: usize) -> FacetFinish {
    facet_finishes.get(facet_idx).copied().unwrap_or_default()
}

/// Interior-side handling of one facet hit -- a direct extraction from
/// `trace_spectral_ray_inner` (same operations in the same order, out of line
/// purely for clippy's line-count lint): flips the geometric normal to face
/// the interior ray, then applies the segment's plain Beer-Lambert absorption.
///
/// A scattering-active material's extinction for this segment was
/// already applied by `try_scatter_step`'s no-scatter branch -- calling the
/// plain `apply_absorption` here too would charge this segment's absorption
/// TWICE. `scattering_sigma_s <= 0.0` is exactly the pre-Task-1 case
/// (including every existing scene), where that branch never ran and this is
/// the ONLY absorption application, byte for byte the same call this module tree
/// always made. Does nothing while the ray is outside the gem.
///
/// `ray_ctx` bundles `(&RayMaterialContext, &RayWavelengthCache)` into a single
/// parameter (R3: rather than two separate ones) purely to keep this function's
/// argument count within clippy's `too_many_arguments` limit without a new `#[allow]`.
///
/// P2: `k_hat` is the WAVE NORMAL `k` (`current_k` in the caller), not the Poynting
/// direction `S` -- `apply_absorption`'s `channel_absorption_alphas` call derives its
/// eigen-polarizations from this direction, which is a property of `k`, not `S`; `hit_t`
/// (the path LENGTH) stays geometric, unaffected. See `refraction`'s own design note,
/// rule 6.
fn apply_interior_segment(
    ray_ctx: (&RayMaterialContext, &RayWavelengthCache),
    current_plane_normal: Vec3,
    k_hat: Vec3,
    hit_t: f32,
    inside_gem: bool,
    normal: &mut Vec3,
    stokes: &mut [StokesVector; NUM_CHANNELS],
) {
    let (ctx, cache) = ray_ctx;
    if !inside_gem {
        return;
    }
    *normal = -*normal;
    if ctx.material.scattering_sigma_s <= 0.0 {
        // See `apply_absorption`'s doc comment: a direct extraction of this
        // same block, same floating-point operations in the same order.
        apply_absorption(ctx, cache, current_plane_normal, k_hat, hit_t, stokes);
    }
}

/// Builds the per-sample [`RayMaterialContext`] -- a direct extraction from
/// `trace_spectral_ray_inner` (same operations in the same order, kept out of
/// line purely for clippy's line-count lint). The anisotropy gate and c-axis
/// lookup are byte-for-byte the code that used to sit inline.
fn build_ray_material_context(
    material: &GemMaterial,
    lambdas: [f32; NUM_CHANNELS],
    hero_idx: usize,
    enable_internal_mode_coupling: bool,
) -> RayMaterialContext<'_> {
    // Optical c-axis for anisotropic birefringence, per-material (
    // previously hard-coded to Vec3::Y for every material).
    let c_axis = material.c_axis;
    let is_anisotropic = material.crystal_system != CrystalSystem::Cubic
        && material.birefringence_delta.abs() > 1e-4;
    RayMaterialContext {
        material,
        lambdas,
        hero_idx,
        c_axis,
        is_anisotropic,
        enable_internal_mode_coupling,
    }
}

/// R3: bundles [`build_ray_material_context`] and [`build_ray_wavelength_cache`] into one
/// call (out of line purely for clippy's line-count lint on
/// `trace_spectral_ray_inner`) -- `RayWavelengthCache`'s own doc comment covers what it
/// caches and why that is bit-identical to the pre-hoist per-bounce recomputation.
fn build_ray_context(
    material: &GemMaterial,
    lambdas: [f32; NUM_CHANNELS],
    hero_idx: usize,
    enable_internal_mode_coupling: bool,
) -> (RayMaterialContext<'_>, RayWavelengthCache) {
    let mat_ctx =
        build_ray_material_context(material, lambdas, hero_idx, enable_internal_mode_coupling);
    let wavelength_cache = build_ray_wavelength_cache(&mat_ctx);
    (mat_ctx, wavelength_cache)
}

/// Records why/where the bounce loop stopped into the caller's `termination_out` slot, a
/// no-op if it is `None` -- extracted purely to keep each of
/// [`trace_spectral_ray_inner`]'s three break sites down to one line, for clippy's
/// function-length lint.
fn record_termination(
    termination_out: &mut Option<&mut (u32, PathTermination)>,
    bounce: u32,
    reason: PathTermination,
) {
    if let Some(out) = termination_out.as_deref_mut() {
        *out = (bounce, reason);
    }
}

/// The actual body of [`trace_spectral_ray`]/[`trace_spectral_ray_with_finish`] -- see
/// those functions' doc comments for the full parameter list.
/// `enable_internal_mode_coupling` is addition: `true` (what both public
/// wrappers always pass) applies `apply_internal_mode_coupling` at every internal
/// reflection inside an anisotropic crystal; `false` reproduces the exact pre-Task-1
/// behaviour (mode fixed at entry for the whole interior traversal), used only by this
/// file's own tests to measure the difference the mechanism makes.
///
/// S4: `plane_soa` is the caller-supplied [`crate::simd::PlanesSoA32`] arena
/// `intersect_polyhedron_soa` scans -- built ONCE by the caller (a batch/frame of many
/// samples all tracing against the same fixed `planes`) rather than once per sample as
/// this function used to do internally. `planes` (the raw slice) is still taken
/// separately -- `shading_normal_near_edge` needs the per-facet plane records
/// themselves, not just the `SoA` arena -- so this is an ADDITIONAL parameter, not a
/// replacement; every caller passes `planes` plus a `&PlanesSoA32` built from that
/// exact same slice. [`build_plane_soa`] is a pure, deterministic function of
/// `planes`' contents, so building it once outside the per-sample hot loop instead of
/// inside it changes nothing about the bits any individual sample computes.
#[expect(
    clippy::too_many_arguments,
    reason = "the shared bounce loop every public entry point above delegates to -- see \
              trace_spectral_ray's own reason for the public-facing parameters, plus the \
              two debug/instrumentation output hooks (primary_hit_out, termination_out) \
              only `bounce_cost` and this file's own tests populate, plus (S4) the \
              caller-supplied `plane_soa` arena alongside the `planes` slice it was \
              built from"
)]
#[expect(
    clippy::too_many_lines,
    reason = "already at the pedantic line-length threshold pre-instrumentation; the \
              `bounce_cost` harness's termination bookkeeping (one `record_termination` \
              call per break site) pushes it a few lines over -- see `record_termination`'s \
              own doc comment for why that extraction, not further splitting this already \
              heavily-decomposed loop, is the right amount of surgery"
)]
fn trace_spectral_ray_inner(
    initial_ray: Ray,
    planes: &[GpuFacetPlane],
    plane_soa: &crate::simd::PlanesSoA32,
    facet_finishes: &[FacetFinish],
    material: &GemMaterial,
    max_bounces: u32,
    environment: EnvironmentSource<'_>,
    rng_seed: u32,
    hero_rand: f32,
    mut primary_hit_out: Option<&mut Option<HitRecord>>,
    enable_internal_mode_coupling: bool,
    enable_exit_splitting: bool,
    mut termination_out: Option<&mut (u32, PathTermination)>,
) -> Vec3 {
    // Fix G (Part 1): hero is now drawn over the FULL visible range [380, 780), not
    // just the first channel_width sub-band, and the remaining channels wrap around
    // instead of running off the top past 780nm -- see `wrapped_hero_wavelengths`'s
    // doc comment for the formula (Algorithm 2, line 5 of
    // GEMSTONE_RENDERING_BLUEPRINT.md) and why it matters for MIS.
    //
    // The property this buys, that the old construction (hero confined to
    // [380, 430), no wraparound) did NOT have: a hero draw `h` and `h + channel_width`
    // now generate the exact SAME 8-member wavelength comb, just cyclically rotated --
    // so for a FIXED comb, each of its N members is equally likely (probability
    // channel_width / 400 = 1/N) to be the one drawn as hero on any given invocation.
    // That is the premise spectral MIS (`spectral_mis_weight` below) requires: across
    // the ensemble of independent ray samples a render accumulates, every wavelength
    // must be reachable as the driving ("hero") channel at positive, uniform
    // probability. Under the OLD construction it was not (channel 0 was *always* the
    // hero, i.e. there was exactly one technique, which is why `mis_weighted_radiance`
    // used to be a hard-coded identity function).
    let lambdas: [f32; NUM_CHANNELS] = wrapped_hero_wavelengths(hero_rand);

    let mut stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
    let mut radiance = [0.0f32; NUM_CHANNELS];

    // `hero_idx` names which slot of `lambdas` (and every other per-channel array
    // below) holds the wavelength actually driving the shared geometric path, so call
    // sites read `lambdas[hero_idx]` etc. instead of a bare magic-number `[0]`. Under
    // this exact construction (`lambdas[k] = wrap(lambda_hero + k*channel_width)`, see
    // `wrapped_hero_wavelengths`), substituting k=0 shows `lambdas[0] == lambda_hero`
    // identically on every call -- so `hero_idx` is provably 0 for every invocation,
    // not something that varies at runtime; what varies from call to call (Part 1's
    // whole point) is which physical WAVELENGTH VALUE lambda_hero itself is, now
    // spanning the full spectrum instead of being confined to the first sub-band.
    // `hero_idx` is still threaded through explicitly (rather than hardcoding `[0]` at
    // each of the several call sites below) so the code documents *which* channel is
    // playing the hero role, and so a future change to this construction (e.g.
    // rotating which array slot the hero lands in) has one place to update instead of
    // several silently-wrong `[0]`s.
    let hero_idx: usize = 0;

    // Fix G (Part 2): per-channel running density of "technique k (channel k as hero)
    // would have generated this exact realized path" -- see the TIR/reflect/refract
    // branches below for how each factor is derived, and `spectral_mis_weight`'s doc
    // comment for the final combination. Starts at 1.0 (the multiplicative identity)
    // for every channel: before any interaction, every technique is equally capable
    // of having produced the (as yet nonexistent) path.
    let mut path_pdf = [1.0f32; NUM_CHANNELS];

    let mut current_ray = initial_ray;
    // P2: the wave normal `k`, tracked ALONGSIDE `current_ray.dir` (the Poynting/energy
    // direction `S`) -- see `refraction`'s own "wave normal vs Poynting direction"
    // design note. Starts equal to the initial ray's direction: outside the gem (air)
    // `k == S` always (isotropic medium, no walk-off), so this is exactly `current_ray.dir`
    // until the first air->crystal entry.
    let mut current_k = initial_ray.dir;
    let mut inside_gem = false;
    // `None` until the first well-defined plane of incidence is recorded, so no
    // spurious frame rotation is applied at the very first surface (there was
    // previously no real "previous" plane of incidence to rotate from).
    let mut prev_plane_normal: Option<Vec3> = None;
    // Which eigenmode the ray currently inside the crystal was stochastically
    // assigned to at its most recent air->crystal entry. Meaningless while
    // `!inside_gem`; carried across internal bounces so a path that entered as (say)
    // the extraordinary ray keeps using the extraordinary index for every subsequent
    // internal reflection until it exits.
    //
    // For a UNIAXIAL material this is exactly the ordinary/extraordinary
    // split it always was. For a genuinely BIAXIAL material (Alexandrite, Topaz,
    // Tanzanite) there is no "ordinary" ray -- both eigenmodes are direction-dependent
    // and both walk off -- so this flag is generalized to a plain two-valued mode
    // selector: `false` selects "mode A" (the faster, lower-index root of
    // `BiaxialIndicatrix::wave_indices`), `true` selects "mode B" (the slower,
    // higher-index root). The name is kept (rather than renamed throughout) since the
    // uniaxial case, which is still the overwhelming majority of anisotropic built-ins,
    // keeps its exact original meaning and every read site below already documents the
    // biaxial generalization locally.
    let mut is_extraordinary = false;

    // Fixed across every bounce below -- see `RayMaterialContext`/`RayWavelengthCache`'s
    // own doc comments for why each is bundled rather than passed field by field.
    let (mat_ctx, wavelength_cache) =
        build_ray_context(material, lambdas, hero_idx, enable_internal_mode_coupling);

    // P6: precomputed once per trace, exactly like `accumulate_miss_radiance`'s own
    // `studio_rig` (same rationale -- see `ExitSplitCtx::studio_rig`'s doc comment).
    let exit_split_studio_rig = match environment {
        EnvironmentSource::Studio {
            light_yaw,
            light_pitch,
            ..
        } => Some(crate::optics::studio_rig::StudioRig::new(
            light_yaw,
            light_pitch,
        )),
        EnvironmentSource::HdrMap(_) => None,
    };
    // P6: exit-event splitting's own STAGING accumulator -- see
    // `ExitSplitCtx::split_radiance`'s doc comment for why every split channel's
    // contribution lands here first, folded into `radiance` below ONLY if the
    // shared/hero path itself terminates via `PathTermination::Escaped`.
    let mut split_radiance = [0.0f32; NUM_CHANNELS];
    // P6: per-trace context for exit-event spectral splitting -- see `refraction`'s own
    // top-of-file "Exit-event spectral splitting" doc comment.
    let mut exit_split_ctx = ExitSplitCtx {
        plane_soa,
        environment,
        studio_rig: exit_split_studio_rig,
        split_radiance: &mut split_radiance,
        enabled: enable_exit_splitting,
        compat: [u8::MAX; NUM_CHANNELS],
    };
    // P6: set alongside `record_termination(.., PathTermination::Escaped)` below, the
    // ONLY point `accumulate_miss_radiance` (and hence any REAL `radiance` for this
    // trace) is ever populated -- see `ExitSplitCtx::split_radiance`'s doc comment for
    // why `split_radiance` must only be committed when this ends up `true`.
    let mut path_escaped = false;

    for bounce in 0..max_bounces {
        let hit = intersect_polyhedron_soa(current_ray, plane_soa);

        // Denoiser wiring: the primary ray's first-hit depth/normal/facet
        // index feed the À-Trous denoiser's guide buffers (see `renderer::denoise`).
        // Captured only at bounce 0 -- `current_ray` at this point is still exactly
        // `initial_ray`, and since the single traced geometric path is always driven by
        // the HERO channel (see `hero_idx`'s doc comment above), this unambiguously *is*
        // the hero channel's own first hit, not an average or a guess across the 8
        // dispersed channels a companion wavelength might have refracted toward.
        capture_primary_hit(&mut primary_hit_out, bounce, hit);

        let Some(hit_rec) = hit else {
            // Ray exited or missed the gemstone -> sample the environment source. See
            // `accumulate_miss_radiance`'s doc comment.
            accumulate_miss_radiance(
                environment,
                current_ray.dir,
                &lambdas,
                &stokes,
                &mut radiance,
            );
            path_escaped = true;
            record_termination(&mut termination_out, bounce, PathTermination::Escaped);
            break;
        };
        // Attempt a
        // Henyey-Greenstein scattering event somewhere along this segment (from the
        // current ray origin to the facet just hit, `hit_rec.t` away) BEFORE processing
        // that facet at all -- see `try_scatter_step`'s doc comment for the full
        // estimator/control-flow rationale (extracted into its own function purely to
        // keep this already-long function under clippy's line-count lint; every
        // floating-point operation is exactly what the inline version performed).
        if inside_gem {
            match try_scatter_step(
                &mat_ctx,
                &wavelength_cache,
                material,
                prev_plane_normal,
                &mut current_ray,
                &mut current_k,
                hit_rec.t,
                rng_seed,
                bounce,
                &mut stokes,
                &mut path_pdf,
                exit_split_ctx.split_radiance,
            ) {
                ScatterStepOutcome::NotApplicable | ScatterStepOutcome::ReachedBoundary => {}
                ScatterStepOutcome::ScatteredAndSurvived => {
                    // The scattered Stokes vectors are already depolarized (see
                    // `maybe_scatter_or_extinguish`'s doc comment), so the previous
                    // plane of incidence is no longer physically meaningful -- reset it
                    // exactly like the pre-first-bounce `None`, so the NEXT real facet
                    // hit applies no spurious frame rotation (mirroring
                    // `prev_plane_normal`'s own doc comment above).
                    prev_plane_normal = None;
                    continue;
                }
                ScatterStepOutcome::ScatteredAndTerminated => {
                    record_termination(
                        &mut termination_out,
                        bounce,
                        PathTermination::ScatterAbsorbed,
                    );
                    break;
                }
            }
        }

        let hit_point = current_ray.origin + hit_rec.t * current_ray.dir;
        // Facet edge rounding: see `shading_normal_near_edge`'s doc comment.
        let mut normal = shading_normal_near_edge(
            planes,
            hit_point,
            hit_rec.facet_idx,
            hit_rec.normal,
            material.edge_rounding_radius,
        );

        // See `rotate_stokes_to_plane_of_incidence`'s doc comment: a direct extraction
        // of this same block, same floating-point operations in the same order. P2: the
        // Stokes plane-of-incidence frame is defined by the WAVE NORMAL `k`, not the
        // Poynting direction `S` -- see `refraction`'s own design note.
        let current_plane_normal =
            rotate_stokes_to_plane_of_incidence(current_k, normal, prev_plane_normal, &mut stokes);
        prev_plane_normal = Some(current_plane_normal);

        // See `apply_interior_segment`'s doc comment: a direct extraction of
        // this same block, same floating-point operations in the same order. P2: the
        // propagation direction feeding pleochroic absorption's eigen-polarizations is
        // `k`, not `S` -- `hit_rec.t` (the path LENGTH) stays geometric/`S`-based,
        // unaffected -- see `refraction`'s own design note, rule 6.
        apply_interior_segment(
            (&mat_ctx, &wavelength_cache),
            current_plane_normal,
            current_k,
            hit_rec.t,
            inside_gem,
            &mut normal,
            &mut stokes,
        );

        // See `compute_bounce_refraction_geometry`'s doc comment: a direct
        // extraction of this same block, same floating-point operations in the same
        // order, just packaged as a function. P2: index lookups (`theta_c`, biaxial
        // `wave_indices`), `cos_i`/`sin_i`, and every downstream Snell/Fresnel
        // evaluation use the WAVE NORMAL `k`, not `S` -- see `refraction`'s own design
        // note.
        let geo = compute_bounce_refraction_geometry(
            &mat_ctx,
            &wavelength_cache,
            normal,
            current_k,
            inside_gem,
            is_extraordinary,
        );

        // Girdle finish: which specular/diffuse treatment this facet gets --
        // `Polished` (the default for any index `facet_finishes` doesn't cover) takes
        // the exact pre-Task-2 dispatch; `Frosted` takes `apply_frosted_bounce` instead.
        // See `dispatch_bounce`'s doc comment for the unified call this feeds.
        let finish = facet_finish_for(facet_finishes, hit_rec.facet_idx);
        dispatch_bounce(
            &mat_ctx,
            &wavelength_cache,
            &geo,
            hit_point,
            normal,
            current_plane_normal,
            finish,
            rng_seed,
            bounce,
            &mut stokes,
            &mut path_pdf,
            &mut exit_split_ctx,
            &mut current_ray,
            &mut current_k,
            &mut inside_gem,
            &mut is_extraordinary,
        );

        // See `apply_russian_roulette`'s doc comment: a direct extraction of this same
        // block, same floating-point operations in the same order. P6 (Bias A fix):
        // `split_radiance` now rides along on the same `1/q` survival rescale as
        // `stokes` -- see that function's own doc comment. Reborrowed through
        // `exit_split_ctx.split_radiance` (rather than the outer `split_radiance`
        // local directly) since `exit_split_ctx` already holds that local borrowed
        // mutably for the whole loop -- see `ExitSplitCtx::split_radiance`'s doc
        // comment.
        if bounce > 4
            && !apply_russian_roulette(bounce, rng_seed, &mut stokes, exit_split_ctx.split_radiance)
        {
            record_termination(
                &mut termination_out,
                bounce,
                PathTermination::RussianRoulette,
            );
            break;
        }
    }

    // P6: each channel's MIS family, final after the last interior dispersive event.
    let compat = exit_split_ctx.compat;

    // P6: commit the staged exit-split contributions into the REAL `radiance` array
    // ONLY if the shared/hero path itself reached its own environment lookup -- see
    // `ExitSplitCtx::split_radiance`'s doc comment for why this all-or-nothing gate is
    // required (mirrors `accumulate_miss_radiance`'s own all-or-nothing behaviour for
    // a plain matching channel). A no-op whenever nothing ever split, or splitting was
    // disabled (`split_radiance` stays all-zero in both cases). Every channel -- hero,
    // matching companion and split companion alike -- is then integrated under the ONE
    // shared `spectral_mis_weight`; see `refraction`'s top-of-file "Exit-event spectral
    // splitting" doc comment for why that single weight is the correct combination for
    // all of them.
    if path_escaped {
        for k in 0..NUM_CHANNELS {
            radiance[k] += split_radiance[k];
        }
    }

    // See `integrate_channels_to_xyz`'s doc comment: a direct extraction of this same
    // block, same floating-point operations in the same order. P6: with splitting
    // enabled every channel is weighted over its OWN family instead -- see
    // `integrate_channels_to_xyz_families`'s doc comment for why the shared weight is
    // biased once companions contribute.
    let xyz = if enable_exit_splitting {
        integrate_channels_to_xyz_families(&radiance, &lambdas, &path_pdf, hero_idx, compat)
    } else {
        integrate_channels_to_xyz(&radiance, &lambdas, &path_pdf, hero_idx)
    };

    // Von Kries white-balance (diagonalised in Bradford LMS, not raw XYZ -- see
    // `compute_illuminant_white_balance`'s doc comment) so the chosen illuminant itself
    // renders as neutral white, rather than baking its colour temperature into every
    // pixel.
    //
    // P6 (2026-09-07): `.max(Vec3::ZERO)` clamps the white-balanced result the same way
    // `StokesVector::intensity`/`cie_1931_cmf(_x8)` already clamp every OTHER physical
    // radiance quantity in this crate to non-negative -- pre-white-balance `xyz` above is
    // provably non-negative (`integrate_channels_to_xyz` sums non-negative
    // `cmf * radiance * mis_weight * norm_factor` terms only), but the Bradford LMS
    // chromatic-adaptation matrices (`BRADFORD_XYZ_TO_LMS`/`BRADFORD_LMS_TO_XYZ`, both
    // with negative off-diagonal entries) are a linear transform that does not itself
    // preserve non-negativity for a sufficiently saturated/spectrally-narrow input --
    // a PRE-EXISTING property of that untouched transform, not something P6 introduces.
    // Exit-event splitting simply makes more previously chromatically-terminated
    // companion channels survive to contribute their own genuine spectral content,
    // which can shift a dispersive material's raw XYZ into that narrow-input regime for
    // some (material, ray, seed) combinations where it previously did not -- confirmed
    // by feeding a representative saturated XYZ directly into
    // `apply_von_kries_white_balance` with P6 entirely uninvolved and observing the same
    // small negative component. Clamping here (once, at the one place every entry point
    // returns through) is the same physical non-negativity floor already applied
    // elsewhere in this file, not a change to the estimator itself.
    apply_von_kries_white_balance(xyz, environment_white_balance(environment)).max(Vec3::ZERO)
}

/// O<->e mode re-coupling at internal reflections.
#[cfg(test)]
mod mode_coupling_tests {
    use super::{
        super::{
            environment::LightingPreset,
            uniaxial_fresnel::{self, UniaxialFrame},
        },
        *,
    };
    use crate::geometry::cuts::StandardGemCuts;

    /// Builds the same `(k_hat, normal, frame)` shape `uniaxial_fresnel`'s own tests
    /// use (`normal = -Z`, incidence in the XZ plane) so a `UniaxialFrame` can be built
    /// directly from an incidence angle and optic axis, without a full bounce/geometry
    /// setup.
    fn frame_for(ang_deg: f32, c_axis: Vec3) -> (Vec3, Vec3, UniaxialFrame) {
        let theta = ang_deg.to_radians();
        let cos_i = theta.cos();
        let sin_i = theta.sin();
        let k_hat = Vec3::new(sin_i, 0.0, cos_i);
        let normal = Vec3::new(0.0, 0.0, -1.0);
        (
            k_hat,
            normal,
            UniaxialFrame::build(k_hat, normal, c_axis, cos_i, sin_i),
        )
    }

    /// P2 exit wiring (2026-09-07): `apply_internal_mode_coupling`'s relabeling
    /// probability is now the EXACT closed-form `uniaxial_fresnel::internal_solve`
    /// Poynting-weighted split `p_o = R_o / (R_o + R_e)` (see that function's own
    /// updated doc comment) whenever the caller supplies one, rather than the old
    /// polarization-projection heuristic -- these two tests replace the pre-P2 ones
    /// (which exercised the heuristic via `current_plane_normal`/`new_k` alone) and are
    /// re-derived directly from the closed form instead of from Malus's law.
    ///
    /// Mechanism-level check: for a moderate, non-degenerate uniaxial internal-
    /// reflection geometry (25 degrees incidence, c-axis tilted off both the facet
    /// normal and the incidence plane -- genuine 3D o<->e coupling, the same shape the
    /// walk-off slab test below uses), compute `p_o` directly from `internal_solve`'s
    /// own `R_o`/`R_e`, then check `apply_internal_mode_coupling`'s empirical draw
    /// frequency over many seeds reproduces THAT value (not a hard-coded 0.5) to
    /// within the binomial sampling error at this trial count -- proving the RNG
    /// plumbing correctly implements whatever probability the closed-form solve hands
    /// it, not stuck at a constant.
    #[test]
    fn internal_mode_coupling_draw_matches_exact_closed_form_split() {
        let lambdas = [589.3f32; NUM_CHANNELS];
        let zircon = GemMaterial::by_name("Zircon").expect("Zircon must be a built-in material");
        let ctx = build_ray_material_context(&zircon, lambdas, 0, true);
        let n_o = zircon.dispersion.evaluate(589.3);
        let n_e = zircon.extraordinary_index_at(589.3, n_o);

        let c_axis = Vec3::new(0.4, 0.5, 0.767_2).normalize();
        let (_, _, frame) = frame_for(25.0, c_axis);
        let sol = uniaxial_fresnel::internal_solve(n_o, n_o, n_e, c_axis, &frame, true);
        let ro_pow = sol.r_o.norm_sqr() * sol.flux_ro;
        let re_pow = sol.r_e.norm_sqr() * sol.flux_re;
        let p_o_exact = ro_pow / (ro_pow + re_pow).max(1e-12);

        let stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
        let mut extraordinary_count = 0u32;
        let trials = 20_000u32;
        for seed in 0..trials {
            if apply_internal_mode_coupling(
                &ctx,
                false,
                Vec3::X,
                Vec3::NEG_Z,
                &stokes,
                Some(p_o_exact),
                seed,
                3,
            ) {
                extraordinary_count += 1;
            }
        }
        let frac = f64::from(extraordinary_count) / f64::from(trials);
        let expected = f64::from(1.0 - p_o_exact);
        // 4-sigma binomial tolerance plus a small floor for the f32-vs-f64 conversion
        // noise in `p_o_exact` itself.
        let tol = 4.0f64.mul_add(
            (expected * (1.0 - expected) / f64::from(trials)).sqrt(),
            0.005,
        );
        assert!(
            (frac - expected).abs() < tol,
            "draw frequency should match the exact closed-form split p_o_exact={p_o_exact} \
             (expected extraordinary frac={expected}), got {frac} over {trials} trials \
             (tol={tol})"
        );
    }

    /// The exact split's real physical payoff over the old polarization heuristic: at
    /// an incidence angle between the ordinary and extraordinary modes' own critical
    /// angles, one backward mode is in genuine TIR (`|r|==1`, its own transmitted exit
    /// wave evanescent) while the other is still only partially reflecting -- the
    /// sharpest asymmetry `R_o`/`R_e` can produce, and a case the OLD heuristic (which
    /// never looked at either mode's own reflectance at all, only the incident
    /// polarization's axis alignment) could not represent. Found by a angle sweep
    /// (rather than a hand-derived critical angle) so the test stays correct against
    /// Zircon's tabulated dispersion without re-deriving the crossover by hand.
    #[test]
    fn internal_mode_coupling_draw_is_strongly_biased_near_a_modes_own_critical_angle() {
        let lambdas = [589.3f32; NUM_CHANNELS];
        let zircon = GemMaterial::by_name("Zircon").expect("Zircon must be a built-in material");
        let ctx = build_ray_material_context(&zircon, lambdas, 0, true);
        let n_o = zircon.dispersion.evaluate(589.3);
        let n_e = zircon.extraordinary_index_at(589.3, n_o);
        let c_axis = Vec3::new(0.4, 0.5, 0.767_2).normalize();

        let mut best_angle = 0.0f32;
        let mut best_p_o = 0.5f32;
        let mut best_skew = 0.0f32;
        let steps: u16 = 4000;
        for i in 0..=steps {
            let ang = 20.0f32.mul_add(f32::from(i) / f32::from(steps), 20.0);
            let (_, _, frame) = frame_for(ang, c_axis);
            let sol = uniaxial_fresnel::internal_solve(n_o, n_o, n_e, c_axis, &frame, true);
            let ro_pow = sol.r_o.norm_sqr() * sol.flux_ro;
            let re_pow = sol.r_e.norm_sqr() * sol.flux_re;
            let p_o = ro_pow / (ro_pow + re_pow).max(1e-12);
            let skew = (p_o - 0.5).abs();
            if skew > best_skew {
                best_skew = skew;
                best_angle = ang;
                best_p_o = p_o;
            }
        }
        assert!(
            best_skew > 0.45,
            "test premise: sweep should find a strongly asymmetric angle between the \
             two modes' own critical angles -- best found: angle={best_angle} \
             p_o={best_p_o} skew={best_skew}"
        );

        let stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
        let mut extraordinary_count = 0u32;
        let trials = 20_000u32;
        for seed in 0..trials {
            if apply_internal_mode_coupling(
                &ctx,
                false,
                Vec3::X,
                Vec3::NEG_Z,
                &stokes,
                Some(best_p_o),
                seed,
                7,
            ) {
                extraordinary_count += 1;
            }
        }
        let frac = f64::from(extraordinary_count) / f64::from(trials);
        let expected = f64::from(1.0 - best_p_o);
        let tol = 4.0f64.mul_add(
            (expected * (1.0 - expected) / f64::from(trials)).sqrt(),
            0.005,
        );
        assert!(
            (frac - expected).abs() < tol,
            "at angle={best_angle} (p_o_exact={best_p_o}), draw frequency should match \
             expected={expected}, got {frac} over {trials} trials (tol={tol})"
        );
    }

    /// The decisive measurement: render the SAME real Zircon material (real
    /// `birefringence_delta = +0.0590`, the largest in the built-in catalogue) at the
    /// SAME ray and seeds, averaged over many samples to beat down unrelated sampling
    /// noise, with internal mode coupling forced on vs. forced off
    /// (`trace_spectral_ray_inner`'s `enable_internal_mode_coupling` flag). Any nonzero,
    /// reproducible difference here is caused by exactly one thing: whether the eigenmode
    /// governing internal TIR/reflect bounces is re-rolled after entry -- not by any
    /// other change, since every other code path and RNG draw is identical between the
    /// two calls (same seeds, same ray, same material, same `max_bounces`).
    #[test]
    fn internal_mode_coupling_changes_zircon_render() {
        let planes = StandardGemCuts::standard_round_brilliant();
        // S4: built once, reused across every `trace_spectral_ray_inner` call below --
        // `planes` is fixed for this whole test.
        let plane_soa = build_plane_soa(&planes);
        let zircon = GemMaterial::by_name("Zircon").expect("Zircon must be a built-in material");
        assert!(
            zircon.birefringence_delta.abs() > 0.01,
            "test assumes Zircon is strongly birefringent"
        );

        // Oblique ray (not aligned with c_axis == Y), 12 max_bounces -- the same ray
        // shape `non_biaxial_materials_render_bit_identical_to_pre_chapter_04_golden_values`
        // uses, chosen there specifically because it produces real internal TIR trains.
        let ray = Ray {
            origin: Vec3::new(0.0, 2.5, 0.0),
            dir: Vec3::new(0.18, -1.0, 0.07).normalize(),
        };
        let env = || LightingPreset::RingLights.studio(1.0, 0.85, 0.95);
        let samples = 256u32;

        let mut sum_with = Vec3::ZERO;
        let mut sum_without = Vec3::ZERO;
        for i in 0..samples {
            let seed = 1000 + i;
            let hero_rand = (hash_u32(seed) as f32) / 4_294_967_295.0;
            sum_with += trace_spectral_ray_inner(
                ray,
                &planes,
                &plane_soa,
                &[],
                &zircon,
                12,
                env(),
                seed,
                hero_rand,
                None,
                true,
                true,
                None,
            );
            sum_without += trace_spectral_ray_inner(
                ray,
                &planes,
                &plane_soa,
                &[],
                &zircon,
                12,
                env(),
                seed,
                hero_rand,
                None,
                false,
                true,
                None,
            );
        }
        let mean_with = sum_with / samples as f32;
        let mean_without = sum_without / samples as f32;
        let diff = (mean_with - mean_without).length();
        assert!(
            diff > 1e-4,
            "internal mode coupling should measurably change Zircon's rendered XYZ \
             (mean_with={mean_with:?}, mean_without={mean_without:?}, diff={diff})"
        );

        // Null check: the identical comparison for a CUBIC (isotropic, `is_anisotropic
        // == false`) material must show EXACTLY zero difference -- proving this is the
        // `is_anisotropic` gate doing its job, not sampling noise leaking through a loose
        // threshold above.
        let diamond = GemMaterial::by_name("Diamond").expect("Diamond must be a built-in material");
        let mut sum_with_diamond = Vec3::ZERO;
        let mut sum_without_diamond = Vec3::ZERO;
        for i in 0..samples {
            let seed = 1000 + i;
            let hero_rand = (hash_u32(seed) as f32) / 4_294_967_295.0;
            sum_with_diamond += trace_spectral_ray_inner(
                ray,
                &planes,
                &plane_soa,
                &[],
                &diamond,
                12,
                env(),
                seed,
                hero_rand,
                None,
                true,
                true,
                None,
            );
            sum_without_diamond += trace_spectral_ray_inner(
                ray,
                &planes,
                &plane_soa,
                &[],
                &diamond,
                12,
                env(),
                seed,
                hero_rand,
                None,
                false,
                true,
                None,
            );
        }
        assert_eq!(
            sum_with_diamond, sum_without_diamond,
            "Diamond (cubic, is_anisotropic == false) must be bit-identical with the flag \
             either way -- the mechanism must not leak into isotropic materials"
        );
    }
}

/// P6 (2026-09-07): requirement 1's empirical unbiasedness check plus requirement 2's
/// variance-ratio measurement, both via `trace_spectral_ray_inner`'s own
/// `enable_exit_splitting` A/B switch (the SAME same-crate, same-code-path pattern
/// `mode_coupling_tests` already established for `enable_internal_mode_coupling` --
/// see `trace_spectral_ray`'s own doc comment). Two INDEPENDENT sample sets (disjoint
/// seed ranges, not a paired same-seed comparison) so the two-sample z-test below is
/// the standard, slightly-conservative unpaired form.
#[cfg(test)]
mod exit_splitting_tests {
    use super::*;
    use crate::{geometry::cuts::StandardGemCuts, optics::LightingPreset};

    /// Mean and (population) variance of `trace_spectral_ray_inner`'s XYZ output over
    /// `samples` independent draws (seeds `seed_base..seed_base+samples`), for one
    /// material/ray/`enable_exit_splitting` setting. `f64` accumulation for `sumsq`
    /// avoids catastrophic cancellation in `E[X^2] - E[X]^2` at these sample counts.
    struct RenderStats {
        mean: Vec3,
        var: Vec3,
        n: f64,
    }

    #[expect(
        clippy::too_many_arguments,
        reason = "test-only helper bundling exactly the inputs `trace_spectral_ray_inner` \
                  itself needs plus this sweep's own sample count/seed base/on-off switch"
    )]
    fn render_stats(
        material: &GemMaterial,
        ray: Ray,
        planes: &[GpuFacetPlane],
        plane_soa: &crate::simd::PlanesSoA32,
        env: EnvironmentSource<'_>,
        samples: u32,
        seed_base: u32,
        enable_exit_splitting: bool,
    ) -> RenderStats {
        let mut sum = Vec3::ZERO;
        let mut sumsq_x = 0.0f64;
        let mut sumsq_y = 0.0f64;
        let mut sumsq_z = 0.0f64;
        for i in 0..samples {
            let seed = seed_base.wrapping_add(i);
            let hero_rand = (hash_u32(seed) as f32) / 4_294_967_295.0;
            let xyz = trace_spectral_ray_inner(
                ray,
                planes,
                plane_soa,
                &[],
                material,
                12,
                env,
                seed,
                hero_rand,
                None,
                true,
                enable_exit_splitting,
                None,
            );
            sum += xyz;
            let (xd, yd, zd) = (f64::from(xyz.x), f64::from(xyz.y), f64::from(xyz.z));
            sumsq_x = xd.mul_add(xd, sumsq_x);
            sumsq_y = yd.mul_add(yd, sumsq_y);
            sumsq_z = zd.mul_add(zd, sumsq_z);
        }
        let n = f64::from(samples);
        let mean = sum / samples as f32;
        let (mx, my, mz) = (f64::from(mean.x), f64::from(mean.y), f64::from(mean.z));
        let var = Vec3::new(
            mx.mul_add(-mx, sumsq_x / n).max(0.0) as f32,
            my.mul_add(-my, sumsq_y / n).max(0.0) as f32,
            mz.mul_add(-mz, sumsq_z / n).max(0.0) as f32,
        );
        RenderStats { mean, var, n }
    }

    /// Two-sample z-score per XYZ component: `(mean_a - mean_b) / sqrt(var_a/n_a +
    /// var_b/n_b)`. Matches the pooling `renderer::gpu::estimator_check`'s own Tier 3
    /// image-comparison harness uses (see its "pooled z-score" / "image-aggregate mean
    /// z" diagnostics) -- a magnitude of a few units is ordinary sampling noise at
    /// these trial counts, not evidence of bias; the harness's own real runs (see
    /// project history) have observed a handful of |z| > 3 pixels out of thousands
    /// purely from binomial chance.
    fn z_scores(a: &RenderStats, b: &RenderStats) -> Vec3 {
        let se =
            |va: f32, na: f64, vb: f32, nb: f64| (f64::from(va) / na + f64::from(vb) / nb).sqrt();
        let sx = se(a.var.x, a.n, b.var.x, b.n).max(1e-20);
        let sy = se(a.var.y, a.n, b.var.y, b.n).max(1e-20);
        let sz = se(a.var.z, a.n, b.var.z, b.n).max(1e-20);
        Vec3::new(
            (f64::from(a.mean.x - b.mean.x) / sx) as f32,
            (f64::from(a.mean.y - b.mean.y) / sy) as f32,
            (f64::from(a.mean.z - b.mean.z) / sz) as f32,
        )
    }

    /// Requirement 1 (unbiasedness) and requirement 2 (variance reduction) together:
    /// for each of Zircon, Rutile (this crate's most strongly dispersive built-in --
    /// see `GemMaterial::all_materials`'s own comment) and Synthetic Moissanite, render
    /// the SAME oblique, refraction-heavy ray (the same ray/geometry
    /// `internal_mode_coupling_changes_zircon_render` above uses, chosen there
    /// specifically because it produces real internal TIR/refraction trains) at 4096
    /// independent samples with exit-event splitting OFF (today's pre-P6 chromatic
    /// termination) and ON (this pass), from two disjoint seed ranges.
    ///
    /// If splitting changed the estimator's EXPECTATION (a bug, not the intended pure
    /// variance reduction), `z_scores` would show it: `mean_on` and `mean_off` are two
    /// independent estimates of the exact same underlying integral, so under the null
    /// hypothesis (unbiased) their difference divided by its own standard error is
    /// asymptotically standard normal, and stays a handful of units at most across
    /// three materials x three colour channels at this sample count purely from
    /// sampling noise (matches the existing GPU equivalence harness's own tolerance --
    /// see [`z_scores`]'s doc comment). The variance ratio (`var_off / var_on`) is
    /// reported (not asserted, since the ray/material-dependent magnitude varies) --
    /// see this task's own report for the measured numbers.
    /// Regression guard for the bounce-budget truncation asymmetry this pass found
    /// and fixed (see [`super::ExitSplitCtx::split_radiance`]'s doc comment): a split
    /// channel resolves inline (one bounded intersect test) while the shared/hero path
    /// needs a WHOLE EXTRA loop iteration to discover its own miss, so at a LOW
    /// `max_bounces` cap (more likely to cut the hero off before that extra iteration)
    /// splitting-on/off means diverged sharply -- up to 2.6x for Synthetic Moissanite
    /// at `max_bounces = 4` before the `split_radiance` staging + conditional-commit
    /// fix landed. Pins the ratio close to 1 across a range of caps, catching a
    /// regression of that same asymmetry.
    #[test]
    fn exit_splitting_respects_bounce_budget_truncation() {
        const SAMPLES: u32 = 6_000;
        let planes = StandardGemCuts::standard_round_brilliant();
        let plane_soa = build_plane_soa(&planes);
        let ray = Ray {
            origin: Vec3::new(0.0, 2.5, 0.0),
            dir: Vec3::new(0.18, -1.0, 0.07).normalize(),
        };
        let env = LightingPreset::RingLights.studio(1.0, 0.85, 0.95);
        let material = GemMaterial::by_name("Synthetic Moissanite")
            .expect("Synthetic Moissanite must be a built-in material");
        for max_bounces in [2u32, 4, 6, 8, 12] {
            let mut sum_on = Vec3::ZERO;
            let mut sum_off = Vec3::ZERO;
            for i in 0..SAMPLES {
                let seed_on = 10_000 + i;
                let hero_rand_on = (hash_u32(seed_on) as f32) / 4_294_967_295.0;
                sum_on += trace_spectral_ray_inner(
                    ray,
                    &planes,
                    &plane_soa,
                    &[],
                    &material,
                    max_bounces,
                    env,
                    seed_on,
                    hero_rand_on,
                    None,
                    true,
                    true,
                    None,
                );
                let seed_off = 90_000 + i;
                let hero_rand_off = (hash_u32(seed_off) as f32) / 4_294_967_295.0;
                sum_off += trace_spectral_ray_inner(
                    ray,
                    &planes,
                    &plane_soa,
                    &[],
                    &material,
                    max_bounces,
                    env,
                    seed_off,
                    hero_rand_off,
                    None,
                    true,
                    false,
                    None,
                );
            }
            let mean_on = sum_on / SAMPLES as f32;
            let mean_off = sum_off / SAMPLES as f32;
            let ratio_x = mean_on.x / mean_off.x;
            println!(
                "max_bounces={max_bounces}: mean_on={mean_on:?} mean_off={mean_off:?} ratio_x={ratio_x:.4}"
            );
            assert!(
                (0.97..1.03).contains(&ratio_x),
                "max_bounces={max_bounces}: splitting-on/off diverged sharply \
                 (ratio_x={ratio_x:.4}) -- likely a bounce-budget truncation \
                 regression (mean_on={mean_on:?}, mean_off={mean_off:?})"
            );
        }
    }

    #[test]
    fn exit_splitting_is_unbiased_and_reduces_variance() {
        const SAMPLES: u32 = 16_384;
        let planes = StandardGemCuts::standard_round_brilliant();
        let plane_soa = build_plane_soa(&planes);
        let ray = Ray {
            origin: Vec3::new(0.0, 2.5, 0.0),
            dir: Vec3::new(0.18, -1.0, 0.07).normalize(),
        };
        let env = LightingPreset::RingLights.studio(1.0, 0.85, 0.95);

        for name in ["Zircon", "Rutile", "Synthetic Moissanite"] {
            let material = GemMaterial::by_name(name)
                .unwrap_or_else(|| panic!("{name} must be a built-in material"));
            let on = render_stats(
                &material, ray, &planes, &plane_soa, env, SAMPLES, 10_000, true,
            );
            let off = render_stats(
                &material, ray, &planes, &plane_soa, env, SAMPLES, 90_000, false,
            );
            let z = z_scores(&on, &off);
            let var_ratio = Vec3::new(
                f64::from(off.var.x / on.var.x.max(1e-12)) as f32,
                f64::from(off.var.y / on.var.y.max(1e-12)) as f32,
                f64::from(off.var.z / on.var.z.max(1e-12)) as f32,
            );
            println!(
                "P6 unbiasedness/variance -- {name}: mean_on={:?} mean_off={:?} z={z:?} \
                 var_on={:?} var_off={:?} var_ratio(off/on)={var_ratio:?}",
                on.mean, off.mean, on.var, off.var
            );
            assert!(
                z.x.abs() < 3.5 && z.y.abs() < 3.5 && z.z.abs() < 3.5,
                "{name}: exit-event splitting changed the estimator's expectation \
                 (z={z:?}, mean_on={:?}, mean_off={:?}) -- should be unbiased",
                on.mean,
                off.mean
            );
        }
    }
}

/// P2 (wave normal vs Poynting direction): a plane-parallel uniaxial slab, e-mode
/// forced. Drives `compute_bounce_refraction_geometry`/`apply_partial_fresnel_bounce`
/// directly (bypassing the full stochastic bounce loop, which only returns integrated
/// XYZ radiance with no directional output) at exactly two facets -- the slab's entry
/// and exit faces -- so the resulting `k`/`S` at each step can be inspected directly.
#[cfg(test)]
mod p2_wave_normal_tests {
    use super::{super::intersect::intersect_polyhedron, *};
    use crate::optics::materials::CrystalSystem;

    /// Two real slab faces (normal `+-Y`, `y` in `[-HALF_THICKNESS, HALF_THICKNESS]`)
    /// plus four "far blank" planes (`+-X`/`+-Z` at `+-1000`) purely so
    /// `intersect_polyhedron` sees a legitimate bounded polyhedron -- the walk-off
    /// lateral displacement this test measures is many orders of magnitude smaller than
    /// 1000 model units, so the blanks are never actually reached by the traced ray.
    const HALF_THICKNESS: f32 = 1.0;
    fn slab_planes() -> Vec<GpuFacetPlane> {
        vec![
            GpuFacetPlane::new(Vec3::Y, -HALF_THICKNESS),
            GpuFacetPlane::new(Vec3::NEG_Y, -HALF_THICKNESS),
            GpuFacetPlane::new(Vec3::X, -1000.0),
            GpuFacetPlane::new(Vec3::NEG_X, -1000.0),
            GpuFacetPlane::new(Vec3::Z, -1000.0),
            GpuFacetPlane::new(Vec3::NEG_Z, -1000.0),
        ]
    }

    /// Drives the entry facet (bounce 0, `!inside_gem`) and the exit facet (bounce 1,
    /// `inside_gem`) directly via [`compute_bounce_refraction_geometry`]/
    /// [`apply_partial_fresnel_bounce`], searching `rng_seed` for the first value that
    /// (a) selects the EXTRAORDINARY eigenmode at entry and (b) transmits (never
    /// reflects/TIRs) at BOTH facets -- the "plane-parallel slab, straight through"
    /// scenario this test needs. Returns `(entry_k, entry_s, exit_k, exit_s,
    /// lateral_displacement)`.
    ///
    /// # Panics
    ///
    /// Panics if no seed in `0..SEED_SEARCH_LIMIT` satisfies every condition -- a test
    /// premise, not an expected runtime failure (each condition individually has a
    /// probability of a large fraction of 1, so the combined search is expected to
    /// succeed within a handful of trials).
    #[expect(
        clippy::too_many_lines,
        reason = "P6 added the fixed (`enabled: false`) `ExitSplitCtx` both \
                  `apply_partial_fresnel_bounce` calls below now require (this test's \
                  own assertions only read the hero-direction outputs, unaffected by \
                  splitting -- see the setup's own comment); the entry/exit search loop \
                  body itself is unchanged"
    )]
    fn trace_forced_extraordinary_slab(
        material: &GemMaterial,
        incident_origin: Vec3,
        incident_dir: Vec3,
    ) -> (Vec3, Vec3, Vec3, Vec3, f32) {
        const SEED_SEARCH_LIMIT: u32 = 20_000;
        let planes = slab_planes();
        let lambdas = [589.3f32; NUM_CHANNELS]; // sodium D line, every channel (direction-only test)
        let mat_ctx = build_ray_material_context(material, lambdas, 0, false);
        let cache = build_ray_wavelength_cache(&mat_ctx);

        // The entry hit point depends only on the fixed incident ray/geometry, not on
        // `seed` -- computed once, outside the search loop.
        let entry_hit = intersect_polyhedron(
            Ray {
                origin: incident_origin,
                dir: incident_dir,
            },
            &planes,
        )
        .expect("incident ray must hit the slab's top face");
        let entry_point = incident_origin + entry_hit.t * incident_dir;

        // P6: this test's own assertions only read `k1`/`s1`/`k2`/`s2` (the hero's own
        // direction, computed the same way regardless of splitting) and
        // `lateral_displacement` derived from them -- `enabled: false` keeps every
        // stokes/path_pdf/split_radiance side effect exactly pre-P6, avoiding any
        // dependence on this synthetic slab's `EnvironmentSource` choice.
        let plane_soa = build_plane_soa(&planes);
        let mut split_radiance = [0.0f32; NUM_CHANNELS];
        let mut exit_ctx = ExitSplitCtx {
            plane_soa: &plane_soa,
            environment: EnvironmentSource::Studio {
                preset: crate::optics::raytracer::environment::LightingPreset::RingLights,
                exposure: 1.0,
                light_yaw: 0.0,
                light_pitch: 0.85,
            },
            studio_rig: None,
            split_radiance: &mut split_radiance,
            enabled: false,
            compat: [u8::MAX; NUM_CHANNELS],
        };

        for seed in 0..SEED_SEARCH_LIMIT {
            let mut stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
            let mut path_pdf = [1.0f32; NUM_CHANNELS];

            // Entry facet: normal +Y (top of the slab), current_k == current_ray.dir ==
            // incident_dir (air, isotropic -- k == S trivially before any interface).
            let entry_normal = Vec3::Y;
            let geo0 = compute_bounce_refraction_geometry(
                &mat_ctx,
                &cache,
                entry_normal,
                incident_dir,
                false,
                false,
            );
            if geo0.sin2_t > 1.0 {
                continue; // TIR is geometrically unreachable entering from air (n1==1), but stay defensive.
            }
            // P1: `stokes` here is unpolarized (`Q == U == 0`), so
            // `entry_eigenmode_selection` returns `None` regardless of this value
            // (see its own degenerate-case guard) -- computed properly anyway
            // (mirroring `rotate_stokes_to_plane_of_incidence`'s own formula) rather
            // than an arbitrary placeholder, purely for clarity.
            let entry_plane_normal = incident_dir.cross(entry_normal).normalize_or_zero();
            let (k1, s1, inside_gem_1, is_extraordinary_update, _) = apply_partial_fresnel_bounce(
                &mat_ctx,
                &cache,
                &geo0,
                incident_dir,
                entry_normal,
                entry_plane_normal,
                false,
                false,
                seed,
                0,
                &mut stokes,
                &mut path_pdf,
                &mut exit_ctx,
                entry_point,
            );
            let entered_extraordinary = is_extraordinary_update == Some(true);
            if !inside_gem_1 || !entered_extraordinary {
                continue; // reflected off the top face, or entered the ORDINARY mode -- keep searching.
            }

            // Advance geometrically along S (k1/s1 as computed) to the exit facet.
            let post_entry_origin = entry_point + s1 * 1e-4;
            let exit_ray = Ray {
                origin: post_entry_origin,
                dir: s1,
            };
            let Some(exit_hit) = intersect_polyhedron(exit_ray, &planes) else {
                continue;
            };
            let exit_point = post_entry_origin + exit_hit.t * s1;
            // `intersect_polyhedron`'s hit normal is the OUTWARD-facing plane normal
            // (here, -Y, the bottom face's own defining normal); `apply_interior_segment`
            // flips it to face the interior ray before geometry/dispatch -- see
            // `trace_spectral_ray_inner`'s own `*normal = -*normal` under `inside_gem`.
            let exit_normal = -exit_hit.normal;

            let geo1 =
                compute_bounce_refraction_geometry(&mat_ctx, &cache, exit_normal, k1, true, true);
            if geo1.sin2_t > 1.0 {
                continue; // Would TIR back into the slab -- keep searching for a seed that transmits.
            }
            // Exiting the crystal (`inside_gem == true`), so `entering_anisotropic`
            // is false at this call regardless -- this value is never consulted.
            let exit_plane_normal = k1.cross(exit_normal).normalize_or_zero();
            let (k2, s2, inside_gem_2, _, _) = apply_partial_fresnel_bounce(
                &mat_ctx,
                &cache,
                &geo1,
                k1,
                exit_normal,
                exit_plane_normal,
                true,
                true,
                seed,
                1,
                &mut stokes,
                &mut path_pdf,
                &mut exit_ctx,
                exit_point,
            );
            if inside_gem_2 {
                continue; // Reflected back into the slab instead of transmitting out -- keep searching.
            }

            // Lateral displacement: the perpendicular distance from `exit_point` to the
            // infinite line through `entry_point` along `incident_dir` -- nonzero iff the
            // walk-off genuinely displaced the exit point sideways relative to where a
            // straight (non-birefringent) path would have exited.
            let to_exit = exit_point - entry_point;
            let along = to_exit.dot(incident_dir);
            let perp = to_exit - along * incident_dir;
            let lateral_displacement = perp.length();

            return (k1, s1, k2, s2, lateral_displacement);
        }
        panic!(
            "no seed in 0..{SEED_SEARCH_LIMIT} entered the extraordinary mode and \
             transmitted cleanly through both slab faces -- test premise violated"
        );
    }

    /// The decisive P2 correctness check: for a plane-parallel uniaxial slab with a
    /// tilted c-axis, the EXTRAORDINARY ray's exit wave normal `k` (and, once back in
    /// isotropic air, `S == k`) must come out parallel to the incident ray within 1e-5
    /// -- exactly the classical "parallel slab" result Snell's law gives for ANY single
    /// constant index, applied twice (once per interface) to the SAME `k` throughout the
    /// slab's interior (see `refraction.rs`'s design note, rule 4: exit refraction uses
    /// `k`, not the walked-off `S`). The exit point must ALSO be laterally displaced by a
    /// nonzero amount from the straight-through path -- the walk-off did something real,
    /// it just didn't change the OUTGOING direction.
    #[test]
    fn plane_parallel_uniaxial_slab_extraordinary_ray_exits_parallel_and_displaced() {
        let mut material = GemMaterial::by_name("Zircon")
            .expect("\"Zircon\" is a built-in uniaxial material in GemMaterial::all_materials()");
        assert_eq!(material.crystal_system, CrystalSystem::Tetragonal);
        assert!(
            material.birefringence_delta.abs() > 0.01,
            "test premise: Zircon must be strongly birefringent"
        );
        // c-axis deliberately tilted away from BOTH the slab normal (Y) and the
        // incidence plane (XY) -- genuine 3D walk-off, not a coincidental in-plane one.
        material.c_axis = Vec3::new(0.3, 0.8, 0.5).normalize();

        let incident_origin = Vec3::new(0.0, 5.0, 0.0);
        // ~16.7 degrees off normal incidence -- comfortably sub-critical for n~1.9, and
        // oblique enough that Snell's law genuinely bends k (unlike normal incidence,
        // where k passes straight through regardless of index and the k/S distinction
        // this task fixes could never show up in the exit direction at all).
        let incident_dir = Vec3::new(0.3, -1.0, 0.0).normalize();

        let (entry_k, entry_s, exit_k, exit_s, lateral_displacement) =
            trace_forced_extraordinary_slab(&material, incident_origin, incident_dir);

        // Sanity: the extraordinary mode's walk-off must have actually fired (entry_k
        // != entry_s) -- otherwise this test would be silently checking the degenerate
        // ordinary-mode-equivalent case instead of what it claims to.
        assert!(
            (entry_k - entry_s).length() > 1e-4,
            "test premise: the extraordinary mode's walk-off should visibly separate k \
             from S at entry (entry_k={entry_k:?}, entry_s={entry_s:?})"
        );

        let cos_parallel_k = exit_k.dot(incident_dir).clamp(-1.0, 1.0);
        let cos_parallel_s = exit_s.dot(incident_dir).clamp(-1.0, 1.0);
        assert!(
            (1.0 - cos_parallel_k).abs() < 1e-5,
            "exit wave normal k must be parallel to the incident ray within 1e-5 \
             (exit_k={exit_k:?}, incident_dir={incident_dir:?}, 1-cos={})",
            1.0 - cos_parallel_k
        );
        // Once back in air, S == k exactly (isotropic medium) -- both must agree.
        assert!(
            (1.0 - cos_parallel_s).abs() < 1e-5,
            "exit Poynting direction S (== k in air) must be parallel to the incident ray \
             within 1e-5 (exit_s={exit_s:?}, incident_dir={incident_dir:?}, 1-cos={})",
            1.0 - cos_parallel_s
        );
        assert!(
            (exit_k - exit_s).length() < 1e-6,
            "S must equal k exactly once back in isotropic air (exit_k={exit_k:?}, \
             exit_s={exit_s:?})"
        );
        assert!(
            lateral_displacement > 1e-4,
            "the exit point must be laterally displaced from the straight-through path \
             by a nonzero amount (the walk-off's real, physical effect) -- got {lateral_displacement}"
        );
    }
}