indicatrix 0.1.0

//! Henyey-Greenstein volumetric scattering and the frosted-facet BSDF.
//!
//! [`maybe_scatter_or_extinguish`]'s per-bounce extinction/scatter estimator, the HG
//! phase function and its importance sampler, and [`apply_frosted_bounce`]'s diffuse
//! reflect/transmit dispatch.

use super::{
    NUM_CHANNELS,
    absorption::channel_absorption_alphas,
    camera::Ray,
    refraction::{BounceRefractionGeometry, RayMaterialContext, RayWavelengthCache},
    sampling::{
        BIREFRINGENT_SPLIT_STREAM, DISTANCE_SAMPLE_STREAM, FRESNEL_BRANCH_STREAM,
        FROSTED_DIR_U_STREAM, FROSTED_DIR_V_STREAM, PHASE_DIR_U_STREAM, PHASE_DIR_V_STREAM,
        hash_u32,
    },
    transport::apply_russian_roulette,
};
use crate::optics::{materials::GemMaterial, polarization::StokesVector};
use glam::Vec3;

/// Homogeneous Henyey-Greenstein volumetric scattering, decided/sampled once for the
/// shared hero-driven geometric path and reweighted per channel -- the same structural
/// pattern [`apply_partial_fresnel_bounce`]/[`apply_refract_channel`] use for every
/// other stochastic decision in this module tree.
///
/// `alphas` is the caller's precomputed [`channel_absorption_alphas`] result (computed
/// once in `trace_spectral_ray_inner`) rather than derived internally, so
/// `shaders/transport_physics.wgsl`'s WGSL translation can share this exact signature.
///
/// # The estimator, hazard by hazard
///
/// 1. **Extinction, not absorption.** `sigma_t = sigma_a + sigma_s` (`sigma_a` is this
///    channel's chromatic [`channel_absorption_alphas`] value, unchanged from
///    [`apply_absorption`]; `sigma_s` is the material's achromatic
///    [`GemMaterial::scattering_sigma_s`]). A scattering event REDIRECTS the path
///    rather than destroying its energy -- loss happens only through `sigma_a`.
/// 2. **Every stochastic decision divides by the SAME value that drove it.** One
///    `[0, 1)` random (`DISTANCE_SAMPLE_STREAM`) inverts the HERO channel's own
///    exponential CDF (`sigma_t_hero`) into a free-path distance `t_free`. Both
///    branches divide by the density that draw has under the hero's own technique
///    (`pdf_hero` scattering; `survive_hero` surviving) -- never a per-channel value.
///    Each companion channel k's own physical weight (`tr_k`/`survive_k`, from k's own
///    `sigma_t_k`) sits in the numerator, and `path_pdf[k]` is updated with k's own
///    technique density for the final `spectral_mis_weight` combination. For
///    `k == hero_idx` both branches collapse to `sigma_s / sigma_t_hero` (the
///    single-scattering albedo) and exactly `1.0` -- the standard unbounded-free-path-
///    sampling identities. Genuinely different estimator SHAPE than `apply_absorption`'s
///    deterministic `exp(-alpha*path_len)` (unbiased only in expectation), which is why
///    this is gated behind `scattering_sigma_s > 0.0` rather than relied on to reduce
///    to the old path algebraically at `sigma_s -> 0`.
/// 3. **HG phase/pdf cancellation.** [`sample_henyey_greenstein_direction`]
///    importance-samples exactly the HG phase function's own normalized distribution,
///    so `phase / pdf == 1.0` identically -- mirroring `cosine_weighted_hemisphere`'s
///    cancellation for the frosted BSDF.
/// 4. **Achromatic direction and distance.** Exactly one `t_free` and one scattered
///    direction (from the hero's `sigma_t_hero` and the material's achromatic
///    `scattering_g`) are sampled and reused for every channel, so -- unlike the
///    polished path's `direction_matches` guard -- there is no chromatic-termination
///    check: there is only ever one direction to match.
/// 5. **New, decorrelated RNG streams** ([`DISTANCE_SAMPLE_STREAM`],
///    [`PHASE_DIR_U_STREAM`], [`PHASE_DIR_V_STREAM`]), mirrored bit-for-bit in
///    `shaders/transport_physics.wgsl`.
///
/// # Depolarization
///
/// Like [`apply_frosted_bounce`], every channel's Stokes vector collapses to
/// `StokesVector::unpolarized` at a scatter event: incoherent scattering off an
/// inclusion's internal structure scrambles the phase relationship polarization
/// depends on.
///
/// # Return value
///
/// `Some((t_free, new_dir))` if a scatter event fired strictly before `hit_t` (caller
/// must redirect `current_ray` from `t_free` along its OLD direction, not advance all
/// the way to the facet); `None` if the path survived to the facet boundary
/// (`stokes`/`path_pdf` already carry that survival's extinction weight -- the caller
/// must NOT also call [`apply_absorption`] for this same segment).
///
/// `pub(crate)`: `renderer::gpu::transport_check`'s Tier 2 ULP check calls this real
/// function directly, comparing against the WGSL translation. Visibility only.
#[expect(
    clippy::too_many_arguments,
    reason = "argument list mirrors transport_physics.wgsl's own maybe_scatter_or_extinguish \
              one-for-one; bundling into a struct would break that direct correspondence"
)]
pub(crate) fn maybe_scatter_or_extinguish(
    alphas: &[f32; NUM_CHANNELS],
    sigma_s: f32,
    g: f32,
    hero_idx: usize,
    ray_dir: Vec3,
    hit_t: f32,
    path_scale: f32,
    rng_seed: u32,
    bounce: u32,
    stokes: &mut [StokesVector; NUM_CHANNELS],
    path_pdf: &mut [f32; NUM_CHANNELS],
) -> Option<(f32, Vec3)> {
    let sigma_t_hero = alphas[hero_idx] + sigma_s;

    let dist_rand =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ DISTANCE_SAMPLE_STREAM)) as f32) / 4_294_967_295.0;
    // Unbounded exponential free-path sample (see hazard 2 above). `one_minus_u` is
    // reused as `pdf_hero` below (`exp(-sigma_t_hero * t_free) == one_minus_u`
    // algebraically) rather than recomputing the exponential. `t_free` is in
    // ABSORPTION-LENGTH units (same as `alphas`/`sigma_s`), converted back below.
    let one_minus_u = (1.0 - dist_rand).max(1e-7);
    let t_free = -(one_minus_u.ln()) / sigma_t_hero;

    // `hit_t` arrives in MODEL units; scale to absorption-length units (see
    // `GemMaterial::absorption_path_scale`). `path_scale == 1.0` (every built-in) makes
    // this an exact no-op, so `hit_t_scaled` is bit-identical to plain `hit_t`.
    let hit_t_scaled = hit_t * path_scale;

    if t_free < hit_t_scaled {
        let pdf_hero = sigma_t_hero * one_minus_u;
        // Vectorized Beer-Lambert; exp_f32x8 is a few-ULP polynomial exponential, NOT
        // bit-identical to f32::exp -- see apply_absorption's comment / src/simd.rs.
        let mut tr_args = [0f32; NUM_CHANNELS];
        for k in 0..NUM_CHANNELS {
            tr_args[k] = -(alphas[k] + sigma_s) * t_free;
        }
        let tr = crate::simd::exp_f32x8(tr_args);
        for k in 0..NUM_CHANNELS {
            let sigma_t_k = alphas[k] + sigma_s;
            let tr_k = tr[k];
            let weight = tr_k * sigma_s / pdf_hero;
            stokes[k] = StokesVector::unpolarized(stokes[k].intensity() * weight);
            path_pdf[k] *= sigma_t_k * tr_k;
        }
        let u1 =
            (hash_u32(rng_seed ^ hash_u32(bounce ^ PHASE_DIR_U_STREAM)) as f32) / 4_294_967_295.0;
        let u2 =
            (hash_u32(rng_seed ^ hash_u32(bounce ^ PHASE_DIR_V_STREAM)) as f32) / 4_294_967_295.0;
        let new_dir = sample_henyey_greenstein_direction(u1, u2, g, ray_dir);
        // Convert back to MODEL units (`ray_dir` is a unit vector in model space) before
        // the caller advances `current_ray.origin` along it. `path_scale == 1.0` makes
        // this division an exact no-op; the estimator's weights are unaffected since
        // they're computed entirely in absorption-length units above.
        let t_free_model = t_free / path_scale;
        Some((t_free_model, new_dir))
    } else {
        // Vectorized Beer-Lambert (see the scatter branch above). `survive_hero` is
        // read from this same batched result at `hero_idx` so the hero lane's ratio
        // always comes from the same exponential its companions used.
        let mut survive_args = [0f32; NUM_CHANNELS];
        for k in 0..NUM_CHANNELS {
            survive_args[k] = -(alphas[k] + sigma_s) * hit_t_scaled;
        }
        let survive = crate::simd::exp_f32x8(survive_args);
        let survive_hero = survive[hero_idx];
        for k in 0..NUM_CHANNELS {
            let survive_k = survive[k];
            stokes[k] = stokes[k].scale(survive_k / survive_hero.max(1e-30));
            path_pdf[k] *= survive_k;
        }
        None
    }
}

/// What [`try_scatter_step`] found, and what `trace_spectral_ray_inner`'s bounce loop
/// should do about it.
pub(super) enum ScatterStepOutcome {
    /// `material.scattering_sigma_s <= 0.0` -- skipped entirely (no RNG draw, no
    /// extinction applied). Caller falls through to the plain facet-dispatch code,
    /// including its own `apply_absorption` call.
    NotApplicable,
    /// No scatter event fired: the path survived to the facet boundary, and
    /// `maybe_scatter_or_extinguish` already applied this segment's extinction weight.
    /// Caller falls through to facet-dispatch but must NOT also call `apply_absorption`.
    ReachedBoundary,
    /// A scatter event fired and the path survived Russian roulette (or wasn't yet
    /// eligible). Caller should `continue` the bounce loop.
    ScatteredAndSurvived,
    /// A scatter event fired and Russian roulette terminated the path. Caller should
    /// `break` the bounce loop.
    ScatteredAndTerminated,
}

/// Attempts a Henyey-Greenstein scattering event on the segment from `current_ray`'s
/// current origin to the facet just hit (`hit_t` away), mutating
/// `current_ray`/`stokes`/`path_pdf` in place and returning what happened -- extracted
/// out of `trace_spectral_ray_inner`'s bounce loop to keep that function under
/// clippy's line-count lint.
///
/// Gated on `material.scattering_sigma_s > 0.0` -- `<= 0.0` returns
/// [`ScatterStepOutcome::NotApplicable`] immediately, before drawing from any new RNG
/// stream. See [`maybe_scatter_or_extinguish`] for the estimator itself.
///
/// Applies Russian roulette (via the same [`apply_russian_roulette`] the polished/
/// frosted facet-dispatch path calls at the end of every bounce, same `bounce > 4`
/// gate) when a scatter event fires, since it otherwise never reaches the bounce
/// loop's own trailing call.
///
/// `current_k` is the wave normal `k`, while [`maybe_scatter_or_extinguish`]'s
/// `ray_dir` argument stays `current_ray.dir` (`S`) -- different roles. A fired
/// scatter event collapses `k == S` going forward (`*current_k = new_dir`): scattering
/// already depolarizes the Stokes vector, so there's no wave-normal-vs-Poynting
/// distinction left to carry past it.
#[expect(
    clippy::too_many_arguments,
    reason = "bundles the fixed-for-the-trace contexts (mat_ctx, cache), this bounce's own \
              state, the RNG stream identity, and the per-ray state a scatter event can \
              mutate -- the same shape dispatch_bounce's reason explains in transport.rs"
)]
pub(super) fn try_scatter_step(
    mat_ctx: &RayMaterialContext,
    cache: &RayWavelengthCache,
    material: &GemMaterial,
    prev_plane_normal: Option<Vec3>,
    current_ray: &mut Ray,
    current_k: &mut Vec3,
    hit_t: f32,
    rng_seed: u32,
    bounce: u32,
    stokes: &mut [StokesVector; NUM_CHANNELS],
    path_pdf: &mut [f32; NUM_CHANNELS],
    split_radiance: &mut [f32; NUM_CHANNELS],
) -> ScatterStepOutcome {
    if material.scattering_sigma_s <= 0.0 {
        return ScatterStepOutcome::NotApplicable;
    }
    let s_axis = prev_plane_normal.unwrap_or(Vec3::ZERO);
    let alphas = channel_absorption_alphas(mat_ctx, cache, s_axis, *current_k, stokes);
    let Some((t_free, new_dir)) = maybe_scatter_or_extinguish(
        &alphas,
        material.scattering_sigma_s,
        material.scattering_g,
        mat_ctx.hero_idx,
        current_ray.dir,
        hit_t,
        material.absorption_path_scale,
        rng_seed,
        bounce,
        stokes,
        path_pdf,
    ) else {
        return ScatterStepOutcome::ReachedBoundary;
    };
    current_ray.origin += t_free * current_ray.dir;
    current_ray.dir = new_dir;
    *current_k = new_dir;
    // `split_radiance` rides along on this same survival rescale, exactly like the
    // bounce loop's own trailing Russian-roulette call -- see `apply_russian_roulette`.
    if bounce > 4 && !apply_russian_roulette(bounce, rng_seed, stokes, split_radiance) {
        return ScatterStepOutcome::ScatteredAndTerminated;
    }
    ScatterStepOutcome::ScatteredAndSurvived
}

/// A stable orthonormal basis `(t, b)` perpendicular to unit vector `n` -- the same
/// branch-minimal construction `birefringence::stable_orthonormal_basis` uses, kept as
/// its own local copy rather than exposing that private helper across modules.
fn frosted_orthonormal_basis(n: Vec3) -> (Vec3, Vec3) {
    let a = if n.x.abs() > 0.9 { Vec3::Y } else { Vec3::X };
    let t = (a - n * n.dot(a)).normalize_or_zero();
    let b = n.cross(t);
    (t, b)
}

/// Cosine-weighted hemisphere direction about `n`, from two independent uniform
/// `[0, 1)` randoms (Malley's method -- polar mapping, not the concentric-disk variant,
/// for a direct WGSL trig-call translation). See `apply_frosted_bounce` for why its
/// call sites need no explicit pdf division.
///
/// `pub(crate)`: `renderer::gpu::transport_check`'s Tier 2 ULP check compares this real
/// function against `shaders/transport_physics.wgsl`'s translation. Visibility only.
pub(crate) fn cosine_weighted_hemisphere(u1: f32, u2: f32, n: Vec3) -> Vec3 {
    let r = u1.sqrt();
    let theta = 2.0 * std::f32::consts::PI * u2;
    let (sin_t, cos_t) = theta.sin_cos();
    let (t, b) = frosted_orthonormal_basis(n);
    (t * (r * cos_t) + b * (r * sin_t) + n * (1.0 - u1).max(0.0).sqrt()).normalize_or_zero()
}

/// The Henyey-Greenstein phase function, normalized to integrate to `1.0` over the
/// full sphere (4*pi steradians).
///
/// `cos_theta` is the cosine of the angle between the SCATTERED direction and the ray's
/// ORIGINAL propagation direction (not "look back toward the source") -- so `g > 0`
/// means forward scattering, matching
/// [`crate::optics::materials::GemMaterial::scattering_g`]. `g == 0.0` reduces to the
/// isotropic constant `1 / (4*pi)`.
///
/// Not called by [`maybe_scatter_or_extinguish`]'s own estimator (hazard 3:
/// [`sample_henyey_greenstein_direction`] importance-samples this exact distribution,
/// so `phase / pdf` cancels to `1.0` and is never evaluated at a real scattering
/// event) -- exposed so a Tier 2 GPU self-test can pin `phase / pdf == 1.0` directly.
#[must_use]
// Evaluated only by tests and the GPU Tier 2 equivalence check; kept as the reference
// the sampler is checked against, not something to "wire up" into the transport path.
#[cfg(any(test, feature = "gpu"))]
pub(crate) fn henyey_greenstein_phase(cos_theta: f32, g: f32) -> f32 {
    let g2 = g * g;
    let denom = (2.0 * g).mul_add(-cos_theta, 1.0 + g2).max(1e-6).powf(1.5);
    (1.0 - g2) / (4.0 * std::f32::consts::PI * denom)
}

/// Importance-samples a direction from the Henyey-Greenstein phase function about
/// `forward` (the ray's current propagation direction), from two independent uniform
/// `[0, 1)` randoms.
///
/// # Derivation
///
/// The marginal density in `mu = cos_theta` is `p(mu) = (1 - g^2) / (2 * (1 + g^2 -
/// 2*g*mu)^1.5)`. Inverting its CDF (under [`henyey_greenstein_phase`]'s sign
/// convention) gives, for `g != 0`, `mu = (1 + g^2 - ((1 - g^2) / (1 - g +
/// 2*g*u1))^2) / (2*g)`, with the isotropic special case `mu = 1 - 2*u1` near `g == 0`
/// (the general formula divides by `g` twice and is unstable there). `phi` is uniform
/// (no azimuthal dependence); the direction reuses [`frosted_orthonormal_basis`], the
/// same basis [`cosine_weighted_hemisphere`] uses, for a line-for-line WGSL translation.
///
/// Because this samples exactly [`henyey_greenstein_phase`]'s own distribution,
/// `henyey_greenstein_phase(dot(result, forward), g) / pdf(result) == 1.0` for every
/// `(u1, u2, g)` -- see [`maybe_scatter_or_extinguish`] hazard 3.
#[must_use]
pub(crate) fn sample_henyey_greenstein_direction(u1: f32, u2: f32, g: f32, forward: Vec3) -> Vec3 {
    let cos_theta = if g.abs() < 1e-3 {
        2.0f32.mul_add(-u1, 1.0)
    } else {
        let one_minus_g2 = g.mul_add(-g, 1.0);
        let denom = (2.0 * g).mul_add(u1, 1.0 - g);
        let sq = one_minus_g2 / denom;
        sq.mul_add(-sq, g.mul_add(g, 1.0)) / (2.0 * g)
    };
    let cos_theta = cos_theta.clamp(-1.0, 1.0);
    let sin_theta = f32::mul_add(cos_theta, -cos_theta, 1.0).max(0.0).sqrt();
    let phi = 2.0 * std::f32::consts::PI * u2;
    let (sin_p, cos_p) = phi.sin_cos();
    let (t, b) = frosted_orthonormal_basis(forward);
    (t * (sin_theta * cos_p) + b * (sin_theta * sin_p) + forward * cos_theta).normalize_or_zero()
}

/// The replacement for the specular TIR/partial-reflect/refract dispatch at a
/// `FacetFinish::Frosted` facet.
///
/// # The model, honestly
///
/// A real bruted/ground surface's reflectance doesn't carry the polished interface's
/// sharp per-wavelength Fresnel structure -- roughness averages over local micro-facet
/// angles, washing out the fine structure -- so this uses ONE broadband reflect/
/// transmit fraction (`r_unpol`, hero channel only) for every spectral channel, unlike
/// the polished path's per-channel `r_unpol_k`. Every channel shares both the sampled
/// direction (roughness scattering isn't meaningfully dispersive) and the reflect/
/// transmit split probability (a modeling simplification).
///
/// Direction is drawn from the cosine-weighted hemisphere about the correct macroscopic
/// normal (`normal` reflect, `-normal` transmit) -- standard importance sampling for a
/// Lambertian BRDF/BTDF, so `f * cos(theta) / pdf` is the constant albedo (`1.0`,
/// already spent on `r_unpol`/`1 - r_unpol`) -- folded into the `1.0 / r_unpol` (or
/// `t_unpol`) throughput scale below, with no separate pdf division needed.
///
/// Depolarizes every channel (`StokesVector::unpolarized`): multiple internal
/// micro-scattering events scramble the coherent phase relationship polarization
/// depends on.
///
/// # Estimator composition -- why this needs no chromatic termination
///
/// The polished path drops a companion channel to exactly zero when its specular
/// direction misses the hero-driven direction, because a delta BSDF has zero density
/// elsewhere. Neither premise holds here: direction is shared and the BSDF is a smooth
/// hemisphere, so `path_pdf[k] *= r_unpol` (or `t_unpol`) uses the same factor for
/// every channel.
///
/// Reuses `FRESNEL_BRANCH_STREAM` for the reflect-vs-transmit branch (divides by its
/// own selection probability, disjoint energy). `BIREFRINGENT_SPLIT_STREAM` on an
/// air->crystal entry picks which ~0.5 energy SHARE this path becomes, so it is NOT
/// divided by its own 0.5 (same reasoning as `apply_refract_bounce`). Mode
/// re-coupling is applied by the caller.
// `pub(crate)`: `renderer::gpu::transport_check`'s Tier 2 ULP check compares this real
// function against `shaders/transport_physics.wgsl`'s translation. Visibility only.
#[expect(
    clippy::too_many_arguments,
    reason = "bundles ctx/geo, the same contexts apply_partial_fresnel_bounce and \
              apply_refract_bounce use in refraction.rs; the rest matches \
              transport_physics.wgsl's own apply_frosted_bounce parameter-for-parameter"
)]
pub(crate) fn apply_frosted_bounce(
    ctx: &RayMaterialContext,
    geo: &BounceRefractionGeometry,
    normal: Vec3,
    inside_gem: bool,
    is_extraordinary: bool,
    rng_seed: u32,
    bounce: u32,
    stokes: &mut [StokesVector; NUM_CHANNELS],
    path_pdf: &mut [f32; NUM_CHANNELS],
) -> (Vec3, bool, Option<bool>) {
    let u1 =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ FROSTED_DIR_U_STREAM)) as f32) / 4_294_967_295.0;
    let u2 =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ FROSTED_DIR_V_STREAM)) as f32) / 4_294_967_295.0;

    if geo.sin2_t > 1.0 {
        // Forced reflect (TIR), probability 1 -- no draw, no pdf division, mirroring
        // `apply_tir_bounce` for the polished path.
        let new_dir = cosine_weighted_hemisphere(u1, u2, normal);
        for s in &mut *stokes {
            *s = StokesVector::unpolarized(s.intensity());
        }
        return (new_dir, inside_gem, None);
    }

    let cos_t = (1.0 - geo.sin2_t).max(0.0).sqrt();
    let r_s = f32::mul_add(geo.n2, -cos_t, geo.n1 * geo.cos_i)
        / f32::mul_add(geo.n2, cos_t, geo.n1 * geo.cos_i);
    let r_p = f32::mul_add(geo.n1, -cos_t, geo.n2 * geo.cos_i)
        / f32::mul_add(geo.n1, cos_t, geo.n2 * geo.cos_i);
    let r_unpol = (0.5 * r_p.mul_add(r_p, r_s * r_s)).clamp(1e-4, 1.0 - 1e-4);
    let rng_bounce =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ FRESNEL_BRANCH_STREAM)) as f32) / 4_294_967_295.0;

    if rng_bounce < r_unpol {
        let new_dir = cosine_weighted_hemisphere(u1, u2, normal);
        for k in 0..NUM_CHANNELS {
            stokes[k] = StokesVector::unpolarized(stokes[k].intensity() / r_unpol);
            path_pdf[k] *= r_unpol;
        }
        (new_dir, inside_gem, None)
    } else {
        let new_dir = cosine_weighted_hemisphere(u1, u2, -normal);
        let entering_anisotropic = !inside_gem && ctx.is_anisotropic;
        // Mode SELECTION is still a stochastic 50/50 draw -- no throughput weighting
        // accompanies it (a `split_pdf` divisor would estimate energy no interface can
        // deliver; see `apply_refract_bounce` in `refraction.rs`).
        let use_extraordinary = if entering_anisotropic {
            let split_rand = (hash_u32(rng_seed ^ hash_u32(bounce ^ BIREFRINGENT_SPLIT_STREAM))
                as f32)
                / 4_294_967_295.0;
            split_rand < 0.5
        } else {
            is_extraordinary
        };
        let t_unpol = 1.0 - r_unpol;
        for k in 0..NUM_CHANNELS {
            // No `/ split_pdf` here (see `apply_refract_channel` in `refraction.rs`):
            // the selected mode already carries only its ~0.5 energy share, so dividing
            // by the 0.5 selection probability on top would double-count it.
            stokes[k] = StokesVector::unpolarized(stokes[k].intensity() / t_unpol);
            // No `* split_pdf` either -- `spectral_mis_weight` is scale-invariant under
            // multiplying every channel's `path_pdf` by the same uniform factor, so this
            // was a no-op on the actual MIS weight, just one risking underflow.
            path_pdf[k] *= t_unpol;
        }
        (
            new_dir,
            !inside_gem,
            entering_anisotropic.then_some(use_extraordinary),
        )
    }
}

/// Tests for [`henyey_greenstein_phase`], [`sample_henyey_greenstein_direction`],
/// [`maybe_scatter_or_extinguish`], and the `scattering_sigma_s` gate in
/// `trace_spectral_ray_inner`'s bounce loop.
#[cfg(test)]
mod scattering_tests {
    use super::{
        super::{
            camera::Camera,
            color::cie_1931_cmf,
            environment::{EnvironmentSource, LightingPreset},
            transport::trace_spectral_ray,
        },
        *,
    };
    use crate::{
        geometry::cuts::StandardGemCuts,
        renderer::env_map::{EnvironmentMap, rgb_to_spectral_radiance},
    };

    /// [`henyey_greenstein_phase`] must integrate to `1.0` over the full sphere for any
    /// `g`, via dense quadrature over `cos_theta` (no azimuthal dependence, so the
    /// solid-angle integral reduces to `2*pi * integral_{-1}^{1} phase(mu) dmu`).
    #[test]
    fn hg_phase_integrates_to_one_over_the_sphere() {
        for g in [-0.8f32, -0.3, 0.0, 0.3, 0.7, 0.95] {
            let steps = 200_000u32;
            let dmu = 2.0 / f64::from(steps);
            let mut integral = 0.0f64;
            for i in 0..steps {
                let mu = f64::mul_add(f64::from(i) + 0.5, dmu, -1.0);
                integral = f64::mul_add(
                    f64::from(henyey_greenstein_phase(mu as f32, g)),
                    dmu,
                    integral,
                );
            }
            integral *= 2.0 * std::f64::consts::PI;
            assert!(
                (integral - 1.0).abs() < 1e-3,
                "HG phase function should integrate to 1.0 over the sphere for g={g}, got {integral}"
            );
        }
    }

    /// `henyey_greenstein_phase(1.0, g)` (exact forward direction) must strictly
    /// increase with `g` for `g > 0` -- pins the sign convention ("positive g
    /// forward-scatters") against the phase function itself, independent of the sampler.
    #[test]
    fn hg_phase_is_more_forward_peaked_for_larger_positive_g() {
        let p0 = henyey_greenstein_phase(1.0, 0.0);
        let p5 = henyey_greenstein_phase(1.0, 0.5);
        let p9 = henyey_greenstein_phase(1.0, 0.9);
        assert!(
            p0 < p5 && p5 < p9,
            "forward-direction phase value should strictly increase with g: p(g=0)={p0}, \
             p(g=0.5)={p5}, p(g=0.9)={p9}"
        );
    }

    /// The decisive check for hazard 3: [`sample_henyey_greenstein_direction`]'s mean
    /// `cos_theta` must converge to `g` exactly -- the closed-form mean of the HG
    /// distribution. A sign error, swapped `u1`/`u2`, or wrong CDF inversion would
    /// converge to the wrong mean and be caught here.
    #[test]
    fn hg_sampling_mean_cos_theta_converges_to_g() {
        let forward = Vec3::new(0.2, 0.6, 0.77).normalize();
        for g in [-0.7f32, -0.2, 0.0, 0.4, 0.85] {
            let trials = 400_000u32;
            let mut sum = 0.0f64;
            for i in 0..trials {
                let u1 = (hash_u32(i ^ 0x1111_1111) as f32) / 4_294_967_295.0;
                let u2 = (hash_u32(i ^ 0x2222_2222) as f32) / 4_294_967_295.0;
                let dir = sample_henyey_greenstein_direction(u1, u2, g, forward);
                sum += f64::from(dir.dot(forward));
            }
            let mean = (sum / f64::from(trials)) as f32;
            assert!(
                (mean - g).abs() < 0.01,
                "mean cos_theta for g={g} should converge to g itself, got {mean} over {trials} trials"
            );
        }
    }

    /// Sampled directions must stay unit-length and finite across the full `g` range,
    /// including at the isotropic branch threshold.
    #[test]
    fn hg_sampling_always_produces_finite_unit_directions() {
        let forward = Vec3::new(-0.3, 0.1, 0.94).normalize();
        for g in [-0.999f32, -0.5, -1e-3, 0.0, 1e-3, 0.5, 0.999] {
            for i in 0..2000u32 {
                let u1 = (hash_u32(i ^ 0x3333_3333) as f32) / 4_294_967_295.0;
                let u2 = (hash_u32(i ^ 0x4444_4444) as f32) / 4_294_967_295.0;
                let dir = sample_henyey_greenstein_direction(u1, u2, g, forward);
                assert!(
                    dir.is_finite(),
                    "non-finite direction for g={g}, i={i}: {dir:?}"
                );
                assert!(
                    (dir.length() - 1.0).abs() < 1e-4,
                    "non-unit direction for g={g}, i={i}: {dir:?} (len={})",
                    dir.length()
                );
            }
        }
    }

    fn round_brilliant_colourless_scattering_material(sigma_s: f32, g: f32) -> GemMaterial {
        // Colourless, non-dispersive, cubic -- isolates the scattering estimator from
        // the chromatic absorption/dispersion machinery.
        GemMaterial::new_custom("scattering furnace probe", 1.5, 0.0, 0.0, [0.0, 0.0, 0.0])
            .with_scattering(sigma_s, g)
    }

    /// Non-negotiable regression guard: a material with `scattering_sigma_s == 0.0`
    /// reached via `GemMaterial`'s plain default must trace bit-identically to the
    /// same material with scattering explicitly set to `(0.0, g)` for several `g` --
    /// proving the `scattering_sigma_s > 0.0` gate disables the new code path, since
    /// the estimator is not guaranteed to reduce to the old one algebraically at
    /// `sigma_s=0` (see [`maybe_scatter_or_extinguish`] hazard 2).
    #[test]
    fn default_off_scattering_is_bit_identical_regardless_of_g() {
        let planes = StandardGemCuts::standard_round_brilliant();
        let material_default = GemMaterial::ruby();
        assert!(
            material_default.scattering_sigma_s <= 0.0,
            "test premise: Ruby's own built-in scattering_sigma_s must be 0.0"
        );
        let camera = Camera::new(0.35, 0.28, 5.0, 18.0);
        let env = || LightingPreset::Daylight.studio(1.0, 0.4, 0.35);

        for g in [-0.6f32, 0.0, 0.4, 0.95] {
            let material_explicit_zero = material_default.clone().with_scattering(0.0, g);
            for iy in 0..6usize {
                for ix in 0..6usize {
                    let ray = camera.generate_ray(ix as f32, iy as f32, 6.0, 6.0, 0.5, 0.5);
                    for s in 0..4u32 {
                        let pixel_id = (iy as u32) * 6 + (ix as u32);
                        let seed = hash_u32(pixel_id ^ hash_u32(s ^ 0x1357_9BDF));
                        let hero_rand = (hash_u32(seed) as f32) / 4_294_967_295.0;
                        let a = trace_spectral_ray(
                            ray,
                            &planes,
                            &material_default,
                            10,
                            env(),
                            seed,
                            hero_rand,
                            None,
                        );
                        let b = trace_spectral_ray(
                            ray,
                            &planes,
                            &material_explicit_zero,
                            10,
                            env(),
                            seed,
                            hero_rand,
                            None,
                        );
                        assert_eq!(
                            a.to_array(),
                            b.to_array(),
                            "sigma_s=0.0 (default) vs sigma_s=0.0 (explicit, g={g}) must be \
                             BIT identical at pixel ({ix},{iy}) sample {s}"
                        );
                    }
                }
            }
        }
    }

    /// The decisive check: a lossless scattering medium (`sigma_a == 0`, `sigma_s > 0`)
    /// immersed in a spatially uniform environment must still render at exactly that
    /// environment's own radiance -- scattering redirects energy, it must not create
    /// or destroy it. Mirrors `tests/raytracer_tests.rs`'s
    /// `frosted_girdle_white_furnace_energy_conservation_still_holds`.
    #[test]
    fn lossless_scattering_white_furnace_energy_conservation_holds() {
        const L0: f32 = 2.5;
        const SAMPLES_PER_PIXEL: u32 = 96;
        const GRID: usize = 12;
        const TOLERANCE: f32 = 0.08; // generous: CPU-only unit test sample budget

        let planes = StandardGemCuts::standard_round_brilliant();
        let material = round_brilliant_colourless_scattering_material(1.2, 0.4);
        let env_map = EnvironmentMap::uniform(1, 1, [L0, L0, L0]);

        let camera = Camera::new(0.35, 0.28, 5.0, 18.0);
        let mut sum = Vec3::ZERO;
        let mut count = 0u32;
        for iy in 0..GRID {
            for ix in 0..GRID {
                let ray =
                    camera.generate_ray(ix as f32, iy as f32, GRID as f32, GRID as f32, 0.5, 0.5);
                for s in 0..SAMPLES_PER_PIXEL {
                    let pixel_id = (iy as u32) * (GRID as u32) + (ix as u32);
                    let seed = hash_u32(pixel_id ^ hash_u32(s ^ 0x5CA1_AB1E));
                    sum += trace_spectral_ray(
                        ray,
                        &planes,
                        &material,
                        16,
                        EnvironmentSource::HdrMap(&env_map),
                        seed,
                        (hash_u32(seed) as f32) / 4_294_967_295.0,
                        None,
                    );
                    count += 1;
                }
            }
        }
        let mean = sum / count as f32;

        let mut target = Vec3::ZERO;
        for step in 0..=(780 - 380) {
            let lambda = 380.0f32 + step as f32;
            let spec = rgb_to_spectral_radiance([L0, L0, L0], lambda);
            target += cie_1931_cmf(lambda) * spec;
        }
        target /= 106.856;

        let rel_err = |v: f32, t: f32| (v - t).abs() / t.abs().max(1e-6);
        let (ex, ey, ez) = (
            rel_err(mean.x, target.x),
            rel_err(mean.y, target.y),
            rel_err(mean.z, target.z),
        );
        println!(
            "[lossless-scattering furnace] mean={mean:?} target={target:?} rel_err=({ex:.4}, {ey:.4}, {ez:.4}) over {count} samples"
        );
        assert!(
            ex <= TOLERANCE && ey <= TOLERANCE && ez <= TOLERANCE,
            "a lossless scattering medium should still converge to the uniform \
             environment's own radiance (mean={mean:?}, target={target:?}, \
             rel_err=({ex}, {ey}, {ez}), tolerance={TOLERANCE})"
        );
    }

    /// The scattering estimator must actually redistribute energy directionally, not
    /// just pass an energy-conservation check by coincidence: a strongly forward-biased
    /// scattering medium should measurably change a transmissive scene's face-up
    /// appearance relative to `sigma_s = 0`.
    #[test]
    fn scattering_measurably_changes_face_up_appearance() {
        const SAMPLES_PER_PIXEL: u32 = 96;
        const GRID: usize = 14;

        let planes = StandardGemCuts::standard_round_brilliant();
        let clear = GemMaterial::diamond();
        assert!(clear.scattering_sigma_s <= 0.0);
        let hazy = clear.clone().with_scattering(1.5, 0.3);
        let camera = Camera::new(0.35, 0.28, 5.0, 18.0);
        let env = || LightingPreset::RingLights.studio(1.0, 0.85, 0.95);

        let mut sum_clear = Vec3::ZERO;
        let mut sum_hazy = Vec3::ZERO;
        let mut count = 0u32;
        for iy in 0..GRID {
            for ix in 0..GRID {
                let ray =
                    camera.generate_ray(ix as f32, iy as f32, GRID as f32, GRID as f32, 0.5, 0.5);
                for s in 0..SAMPLES_PER_PIXEL {
                    let pixel_id = (iy as u32) * (GRID as u32) + (ix as u32);
                    let seed = hash_u32(pixel_id ^ hash_u32(s ^ 0x0BAD_F00D));
                    let hero_rand = (hash_u32(seed) as f32) / 4_294_967_295.0;
                    sum_clear +=
                        trace_spectral_ray(ray, &planes, &clear, 12, env(), seed, hero_rand, None);
                    sum_hazy +=
                        trace_spectral_ray(ray, &planes, &hazy, 12, env(), seed, hero_rand, None);
                    count += 1;
                }
            }
        }
        let mean_clear = sum_clear / count as f32;
        let mean_hazy = sum_hazy / count as f32;
        let delta_y = (mean_hazy.y - mean_clear.y).abs();
        let relative_change = delta_y / mean_clear.y.max(1e-6);
        println!(
            "[scattering face-up] clear Y={:.5} hazy Y={:.5} delta_y={:.5} ({:.2}%) over {count} samples",
            mean_clear.y,
            mean_hazy.y,
            delta_y,
            100.0 * relative_change
        );
        assert!(
            relative_change > 0.01,
            "a forward-scattering inclusion medium should measurably change face-up \
             brightness (>1%), not render identically to a clear stone -- got {:.4}% \
             (clear Y={:.5}, hazy Y={:.5})",
            100.0 * relative_change,
            mean_clear.y,
            mean_hazy.y
        );
    }

    /// [`maybe_scatter_or_extinguish`]'s hero channel must reproduce the textbook
    /// single-scattering-albedo/unity-survival identities (hazard 2) exactly, checked
    /// directly rather than only inferred from the furnace test above.
    #[test]
    fn hero_channel_weight_matches_albedo_and_unity_survival_identities() {
        let material = round_brilliant_colourless_scattering_material(0.8, 0.0);
        let alphas = [0.0f32; NUM_CHANNELS]; // sigma_a == 0 for this material, every channel
        let sigma_t_hero = material.scattering_sigma_s; // sigma_a == 0 for this material
        let albedo = material.scattering_sigma_s / sigma_t_hero; // == 1.0 (lossless)

        let mut scatter_trials = 0u32;
        let mut survive_trials = 0u32;
        for seed in 0..4000u32 {
            let mut stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
            let mut path_pdf = [1.0f32; NUM_CHANNELS];
            let hit_t = 1.0f32; // fixed boundary distance
            let outcome = maybe_scatter_or_extinguish(
                &alphas,
                material.scattering_sigma_s,
                material.scattering_g,
                0,
                Vec3::Z,
                hit_t,
                1.0,
                seed,
                0,
                &mut stokes,
                &mut path_pdf,
            );
            if outcome.is_some() {
                scatter_trials += 1;
                assert!(
                    (stokes[0].intensity() - albedo).abs() < 1e-4,
                    "hero channel's scatter-branch weight should equal the single-scattering \
                     albedo sigma_s/sigma_t = {albedo} exactly (lossless), got {}",
                    stokes[0].intensity()
                );
            } else {
                survive_trials += 1;
                assert!(
                    (stokes[0].intensity() - 1.0).abs() < 1e-5,
                    "hero channel's no-scatter-branch weight should be exactly 1.0 \
                     (unity-survival identity), got {}",
                    stokes[0].intensity()
                );
            }
        }
        assert!(
            scatter_trials > 100 && survive_trials > 100,
            "test premise: both branches should fire many times at sigma_t*hit_t = {} \
             (scatter={scatter_trials}, survive={survive_trials})",
            sigma_t_hero * 1.0
        );
    }

    /// A companion (non-hero) channel's weight, in both branches, against an
    /// independently-derived closed form -- not merely re-checking the implementation
    /// against itself. Uses two channels with deliberately different chromatic
    /// absorption (Tourmaline's o-ray vs. e-ray at 550nm/430nm) so the hero-vs-companion
    /// divergence actually exercises the chromatic term, unlike a colourless material
    /// where `sigma_t_k == sigma_t_hero` for every k would hide a bug that swaps the
    /// shared hero-pdf denominator for a per-channel one.
    ///
    /// Closed forms (independently re-derived, not copied from the implementation):
    /// - no-scatter branch: `weight_k = exp(-(sigma_t_k - sigma_t_hero) * hit_t)`
    /// - scatter branch: `weight_k = (sigma_s / sigma_t_hero) * exp(-(sigma_t_k -
    ///   sigma_t_hero) * t_free)`, recovered by re-deriving `t_free` from the SAME
    ///   `dist_rand` the function drew (`t_free = -ln(1 - dist_rand) / sigma_t_hero`)
    ///   using a fixed `rng_seed`/`bounce` so the test can predict the exact draw.
    #[test]
    fn companion_channel_weight_matches_independently_derived_closed_form() {
        let material = GemMaterial::by_name("Tourmaline")
            .expect("\"Tourmaline\" must be a built-in material")
            .with_scattering(0.8, 0.0);
        let lambda_hero = 550.0f32;
        let lambda_companion = 430.0f32; // strongly dichroic vs. 550nm, see module doc comment
        let mat_ctx = RayMaterialContext {
            material: &material,
            lambdas: [
                lambda_hero,
                lambda_companion,
                lambda_hero,
                lambda_hero,
                lambda_hero,
                lambda_hero,
                lambda_hero,
                lambda_hero,
            ],
            hero_idx: 0,
            c_axis: material.c_axis,
            is_anisotropic: true,
            enable_internal_mode_coupling: true,
        };
        let cache = super::super::refraction::build_ray_wavelength_cache(&mat_ctx);
        let ray_dir = Vec3::new(0.3, 0.9, 0.3).normalize();
        let sigma_s = material.scattering_sigma_s;

        // Uses the real `channel_absorption_alphas` for alpha_hero/alpha_companion --
        // legitimate here since this test checks the downstream weight formula, not
        // pleochroic absorption itself.
        let stokes_probe = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
        let alphas =
            channel_absorption_alphas(&mat_ctx, &cache, Vec3::ZERO, ray_dir, &stokes_probe);
        assert!(
            (alphas[0] - alphas[1]).abs() > 1e-4,
            "test premise: hero (550nm) and companion (430nm) must have genuinely \
             different alpha (got {} vs {})",
            alphas[0],
            alphas[1]
        );
        let sigma_t_hero = alphas[0] + sigma_s;
        let sigma_t_companion = alphas[1] + sigma_s;

        let mut no_scatter_checked = false;
        let mut scatter_checked = false;
        for seed in 0..3000u32 {
            let hit_t = 0.9f32;
            let mut stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
            let mut path_pdf = [1.0f32; NUM_CHANNELS];
            let outcome = maybe_scatter_or_extinguish(
                &alphas,
                sigma_s,
                material.scattering_g,
                0,
                ray_dir,
                hit_t,
                1.0,
                seed,
                0,
                &mut stokes,
                &mut path_pdf,
            );
            let dist_rand =
                (hash_u32(seed ^ hash_u32(DISTANCE_SAMPLE_STREAM)) as f32) / 4_294_967_295.0;
            let one_minus_u = (1.0 - dist_rand).max(1e-7);
            let t_free = -(one_minus_u.ln()) / sigma_t_hero;

            if outcome.is_some() {
                if t_free >= hit_t {
                    continue; // outcome disagreed with our recomputed t_free; skip (guard only)
                }
                scatter_checked = true;
                let expected =
                    (sigma_s / sigma_t_hero) * (-(sigma_t_companion - sigma_t_hero) * t_free).exp();
                let actual = stokes[1].intensity();
                assert!(
                    (actual - expected).abs() < 1e-3 * expected.abs().max(1.0),
                    "scatter-branch companion weight mismatch at seed={seed}: expected \
                     {expected} (closed form), got {actual}"
                );
            } else {
                no_scatter_checked = true;
                let expected = (-(sigma_t_companion - sigma_t_hero) * hit_t).exp();
                let actual = stokes[1].intensity();
                assert!(
                    (actual - expected).abs() < 1e-3 * expected.abs().max(1.0),
                    "no-scatter-branch companion weight mismatch at seed={seed}: expected \
                     {expected} (closed form), got {actual}"
                );
            }
        }
        assert!(
            no_scatter_checked && scatter_checked,
            "test premise: both branches should be exercised and checked \
             (no_scatter_checked={no_scatter_checked}, scatter_checked={scatter_checked})"
        );
    }

    /// Negative control: if a companion channel's `path_pdf` update used the hero's own
    /// technique density instead of its own, `spectral_mis_weight` would silently
    /// collapse toward `N` instead of reflecting genuine per-channel agreement. Pins
    /// the correct behaviour: two channels with different absorption must end up with
    /// different `path_pdf` values after a scatter event.
    #[test]
    fn scatter_event_path_pdf_is_genuinely_per_channel_when_absorption_is_chromatic() {
        // Real built-in Tourmaline (strongly dichroic o-ray vs. e-ray absorption),
        // opted into scattering.
        let material = GemMaterial::by_name("Tourmaline")
            .expect("\"Tourmaline\" must be a built-in material")
            .with_scattering(0.8, 0.0);
        let mat_ctx = RayMaterialContext {
            material: &material,
            lambdas: [550.0, 430.0, 550.0, 550.0, 550.0, 550.0, 550.0, 550.0],
            hero_idx: 0,
            c_axis: material.c_axis,
            is_anisotropic: true,
            enable_internal_mode_coupling: true,
        };
        let cache = super::super::refraction::build_ray_wavelength_cache(&mat_ctx);
        let ray_dir = Vec3::new(0.3, 0.9, 0.3).normalize();
        let probe_stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
        let alphas =
            channel_absorption_alphas(&mat_ctx, &cache, Vec3::ZERO, ray_dir, &probe_stokes);
        let mut found_scatter = false;
        for seed in 0..2000u32 {
            let mut stokes = [StokesVector::unpolarized(1.0); NUM_CHANNELS];
            let mut path_pdf = [1.0f32; NUM_CHANNELS];
            let outcome = maybe_scatter_or_extinguish(
                &alphas,
                material.scattering_sigma_s,
                material.scattering_g,
                0,
                ray_dir,
                1.0,
                1.0,
                seed,
                0,
                &mut stokes,
                &mut path_pdf,
            );
            if outcome.is_some() {
                found_scatter = true;
                assert!(
                    (path_pdf[0] - path_pdf[1]).abs() > 1e-6,
                    "channel 0 (550nm) and channel 1 (430nm, strongly dichroic) must get \
                     DIFFERENT path_pdf values from a scatter event when the material's \
                     absorption is genuinely chromatic -- got path_pdf[0]={}, path_pdf[1]={}",
                    path_pdf[0],
                    path_pdf[1]
                );
            }
        }
        assert!(
            found_scatter,
            "test premise: at least one trial should hit the scatter branch"
        );
    }
}