indicatrix 0.5.0

//! Homogeneous Henyey-Greenstein volumetric scattering.
//!
//! The per-bounce extinction/scatter estimator ([`maybe_scatter_or_extinguish`]), its
//! bounce-loop wrapper ([`try_scatter_step`]), the phase function and its importance
//! sampler, and the scattering-point next-event-estimation contribution.

use super::{
    super::{
        NUM_CHANNELS,
        absorption::channel_absorption_alphas_assigned,
        camera::{FacetFinish, Ray},
        environment::{EnvironmentSource, sample_environment_for_nee},
        intersect::intersect_polyhedron_soa,
        refraction::{RayMaterialContext, RayWavelengthCache},
        sampling::{
            DISTANCE_SAMPLE_STREAM, NEE_ENV_DIR_U_STREAM, NEE_ENV_DIR_V_STREAM, PHASE_DIR_U_STREAM,
            PHASE_DIR_V_STREAM, hash_u32,
        },
        transport::apply_russian_roulette,
    },
    NeeContext, balance_heuristic, frosted_orthonormal_basis,
};
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_assigned`] result
/// (computed once in `try_scatter_step` below) 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_assigned`] 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`](super::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(in super::super)` (raytracer + descendants), not `pub(super)`: `transport`'s
/// bounce loop (a sibling of `scattering`, not a descendant of it) consumes this
/// directly -- the same scope this type had before the split, when it lived straight
/// in `scattering.rs` (parent: `raytracer`).
pub(in super::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.
    ///
    /// Carries the Henyey-Greenstein phase function's own value at the sampled
    /// continuation direction, paired with that continuation direction itself --
    /// `Some` only when [`NeeContext::enabled`], `None` otherwise (including every
    /// trace with NEE disabled). The caller stashes this as the
    /// pending complementary MIS weight, consumed by whichever LATER bounce-loop
    /// iteration first either escapes directly or dispatches a transmit-out event at
    /// the (necessarily convex) polyhedron's exit facet -- see `transport::dispatch_bounce`'s
    /// doc comment for why the carry must survive that intervening polished-exit
    /// refraction rather than being read only the very next iteration: a scattering
    /// point is strictly interior, so its very next iteration always hits the exit
    /// facet first, never escapes directly there. Once an eventual direct escape
    /// consumes it, that escape is weighted by the balance heuristic against
    /// [`nee_contribution_hg_scatter`]'s own light-sampling technique (evaluated at the
    /// carried INTERIOR direction, the same measure NEE sampled in), rather than
    /// double-counting the same direct-lighting contribution both ways. If a bounce
    /// happens instead while the carry is live -- a reflect, a TIR, another scatter
    /// event, or any hit while still inside the gem -- the pending weight is dropped:
    /// that path has no competing NEE sample left to weigh itself against, so it
    /// counts at full weight.
    ScatteredAndSurvived(Option<(f32, Vec3)>),
    /// 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.
///
/// `is_extraordinary` names which eigenmode this path was assigned to at its most recent
/// air->crystal entry -- see [`channel_absorption_alphas_assigned`]'s doc comment. Does
/// not need `prev_plane_normal`: the assigned-mode absorption this calls reads no
/// Stokes-frame azimuth at all.
#[expect(
    clippy::too_many_arguments,
    reason = "bundles the fixed-for-the-trace contexts (mat_ctx, cache, nee), 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; `lambdas`/`radiance` are the HDR-map NEE \
              contribution's inputs/output, threaded through rather than bundled \
              into `nee` since they vary in a way `NeeContext` (fixed for the whole \
              trace) deliberately does not; `facet_finishes` (the per-facet surface \
              finish) is threaded the same way rather than added to `NeeContext` so \
              the GPU Tier 2 harness's `NeeContext` literals stay independent of it"
)]
pub(in super::super) fn try_scatter_step(
    mat_ctx: &RayMaterialContext,
    cache: &RayWavelengthCache,
    material: &GemMaterial,
    is_extraordinary: bool,
    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],
    nee: NeeContext<'_>,
    lambdas: &[f32; NUM_CHANNELS],
    radiance: &mut [f32; NUM_CHANNELS],
    facet_finishes: &[FacetFinish],
) -> ScatterStepOutcome {
    if material.scattering_sigma_s <= 0.0 {
        return ScatterStepOutcome::NotApplicable;
    }
    let alphas = channel_absorption_alphas_assigned(mat_ctx, cache, *current_k, is_extraordinary);
    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;
    };
    let old_dir = current_ray.dir;
    let scatter_point = current_ray.origin + t_free * old_dir;

    // Direct light sample from the scattering point, MIS-weighted against
    // the phase-sampled continuation below. A no-op (no RNG draw, no accumulator touch)
    // whenever `!nee.enabled`.
    nee_contribution_hg_scatter(
        nee,
        lambdas,
        cache.n_o_ch[mat_ctx.hero_idx],
        scatter_point,
        old_dir,
        material.scattering_g,
        rng_seed,
        bounce,
        stokes,
        radiance,
        &alphas,
        material.scattering_sigma_s,
        material.absorption_path_scale,
        facet_finishes,
    );

    current_ray.origin = scatter_point;
    current_ray.dir = new_dir;
    *current_k = new_dir;
    // The phase-sampled continuation's own technique density at the direction it just
    // drew -- see `ScatterStepOutcome::ScatteredAndSurvived`'s doc comment for how the
    // caller uses this. `None` (not merely `0.0`) whenever NEE is off, so the caller can
    // tell "no competing technique" apart from "competing technique has zero density
    // here" (the latter would legitimately give the escape branch full weight too, but
    // via `balance_heuristic`'s own degenerate-input handling, not by skipping it).
    let phase_pdf_for_mis = nee.enabled.then(|| {
        (
            henyey_greenstein_phase(new_dir.dot(old_dir), material.scattering_g),
            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(phase_pdf_for_mis)
}

/// Next-event estimation at a Henyey-Greenstein scattering event.
///
/// Draws one direction from the environment's own importance distribution
/// ([`sample_environment_for_nee`]), traces a shadow ray from the scattering point to the
/// (necessarily convex) polyhedron's exit facet, applies a scalar Fresnel transmittance
/// there, refracts the sampled direction through that exit facet to look up
/// the environment radiance where the light-sampled ray actually leaves along, applies
/// the medium's own transmittance over the shadow ray's interior path length, and adds
/// the balance-heuristic-weighted contribution directly into `radiance`.
/// A no-op whenever `!nee.enabled`, whenever the environment has no importance
/// distribution to sample (`Studio`), whenever the sampled direction is totally
/// internally reflected at the exit facet (cannot reach the environment at all along that
/// direction -- not a bias: the phase-sampled continuation still has its own chance to
/// escape through a different exit, at full weight, since no NEE sample competes with it
/// there), or whenever the exit facet is [`FacetFinish::Frosted`] (a frosted
/// exit has no well-defined specular Fresnel/refraction for this shadow ray to use --
/// [`nee_contribution_frosted_exterior`](super::nee_contribution_frosted_exterior)
/// already handles NEE for a frosted exit's own diffusely-sampled surface point, so this
/// is a genuine division of labour, not a dropped case).
///
/// Consumes exactly one new 2D RNG draw
/// ([`NEE_ENV_DIR_U_STREAM`]/[`NEE_ENV_DIR_V_STREAM`]) -- a stream pair no existing draw
/// uses, so a trace with NEE disabled never even hashes against it and stays
/// bit-identical to one traced with NEE support absent altogether.
///
/// # Simplification: one achromatic exit transmittance
///
/// Like [`maybe_scatter_or_extinguish`] itself (hazard 4: one achromatic direction/
/// distance shared by every channel), this uses ONE scalar Fresnel transmittance --
/// computed from the HERO channel's own dispersion index (`n_inside_hero`) via the plain
/// isotropic Fresnel formula -- shared by every channel's own spectral radiance at the
/// sampled direction, rather than each channel's own (dispersion-shifted) exit index. A
/// real gem's dispersion at a diffuse scattering event is a second-order effect next to
/// the achromatic phase-function/free-path sampling already in place; a future revision
/// could add per-channel exit indices the way [`compute_channel_transmission`]
/// (`refraction.rs`) already does for the polished exit path.
///
/// # Medium transmittance
///
/// `alphas`/`sigma_s`/`absorption_path_scale` are the SAME quantities
/// [`maybe_scatter_or_extinguish`]'s own survive branch uses (`alphas` is the caller's
/// precomputed [`channel_absorption_alphas_assigned`] result, `sigma_s` is
/// [`GemMaterial::scattering_sigma_s`]), and `hit.t` (the shadow ray's own distance to
/// the exit facet) plays the same role that function's `hit_t` does -- so this applies
/// the identical per-channel `exp(-(alphas[k]+sigma_s)*hit.t*path_scale)` transmittance,
/// via the same [`crate::simd::exp_f32x8`] idiom, that a phase-sampled continuation
/// reaching the same boundary would have paid.
#[expect(
    clippy::too_many_arguments,
    reason = "mirrors try_scatter_step's own argument shape -- context, this bounce's own \
              geometric/RNG state, and the per-channel accumulators an NEE contribution \
              reads/writes; see that function's identical justification"
)]
pub(crate) fn nee_contribution_hg_scatter(
    nee: NeeContext<'_>,
    lambdas: &[f32; NUM_CHANNELS],
    n_inside_hero: f32,
    scatter_point: Vec3,
    scatter_dir_in: Vec3,
    g: f32,
    rng_seed: u32,
    bounce: u32,
    stokes: &[StokesVector; NUM_CHANNELS],
    radiance: &mut [f32; NUM_CHANNELS],
    alphas: &[f32; NUM_CHANNELS],
    sigma_s: f32,
    absorption_path_scale: f32,
    facet_finishes: &[FacetFinish],
) {
    if !nee.enabled {
        return;
    }
    let u0 =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ NEE_ENV_DIR_U_STREAM)) as f32) / 4_294_967_295.0;
    let u1 =
        (hash_u32(rng_seed ^ hash_u32(bounce ^ NEE_ENV_DIR_V_STREAM)) as f32) / 4_294_967_295.0;
    let Some(sample) = sample_environment_for_nee(nee.environment, u0, u1) else {
        return;
    };

    // Shadow ray from the scattering point toward the sampled direction, offset by a
    // small epsilon along it -- mirrors `try_split_exit_channel`'s identical `1e-4`
    // offset (`refraction.rs`) to clear the origin point itself before probing.
    let probe = Ray {
        origin: scatter_point + sample.dir * 1e-4,
        dir: sample.dir,
    };
    let Some(hit) = intersect_polyhedron_soa(probe, nee.plane_soa) else {
        // No exit facet found -- shouldn't happen for a point genuinely inside a closed
        // convex solid, but a degenerate/open mesh must not panic or fabricate energy.
        return;
    };
    // A frosted exit facet has no well-defined specular Fresnel/refraction
    // for this shadow ray to use -- see this function's own doc comment for why that is
    // a division of labour with `nee_contribution_frosted_exterior`, not a dropped case.
    if facet_finishes
        .get(hit.facet_idx)
        .copied()
        .unwrap_or_default()
        == FacetFinish::Frosted
    {
        return;
    }

    // Scalar (hero-index, isotropic-approximation) exit Fresnel transmittance -- see
    // this function's own "Simplification" doc section. `hit.normal` is the exit
    // facet's own OUTWARD normal (`intersect_polyhedron`'s convention: an interior
    // origin's far/exit intersection returns the plane's own stored normal, unflipped),
    // and the shadow ray travels outward through it, so `cos_i = dot(sample.dir,
    // hit.normal)` (no negation) is this interface's cosine of incidence -- the same
    // sign convention `compute_channel_transmission`'s own `n1 -> n2` exit formula uses.
    let cos_i = sample.dir.dot(hit.normal).clamp(0.0, 1.0);
    let sin2_t = (n_inside_hero * n_inside_hero * cos_i.mul_add(-cos_i, 1.0)).min(1.0);
    if sin2_t >= 1.0 {
        // Total internal reflection along this direction -- see this function's own doc
        // comment for why this is a legitimate zero, not a bias.
        return;
    }
    let cos_t = (1.0 - sin2_t).max(0.0).sqrt();
    let r_s = n_inside_hero.mul_add(cos_i, -cos_t) / n_inside_hero.mul_add(cos_i, cos_t);
    let r_p = n_inside_hero.mul_add(-cos_t, cos_i) / n_inside_hero.mul_add(cos_t, cos_i);
    let r_unpol = (0.5 * r_p.mul_add(r_p, r_s * r_s)).clamp(0.0, 1.0);
    let t_unpol = 1.0 - r_unpol;

    // The competing (phase-function) technique's own density at this SAME
    // light-sampled direction -- the balance heuristic's other half.
    let phase_cos = sample.dir.dot(scatter_dir_in);
    let phase_val = henyey_greenstein_phase(phase_cos, g);
    let mis_weight = balance_heuristic(sample.pdf, phase_val);
    if mis_weight <= 0.0 {
        return;
    }

    // The exterior direction this light sample actually leaves along --
    // Snell's law at the exit facet, `sample.dir` playing the incident-direction role
    // and (the inward-flipped) `-hit.normal` the surface normal, mirroring
    // `refraction.rs`'s own `eta*k_hat + (eta*cos_i - cos_t)*normal` vector form (see
    // this function's doc comment for the sign-convention correspondence). `sample.pdf`
    // itself stays in the INTERIOR (pre-refraction) measure `sample_environment_for_nee`
    // sampled in -- only the radiance LOOKUP moves to the refracted direction.
    let refracted_dir = (n_inside_hero * sample.dir
        - n_inside_hero.mul_add(cos_i, -cos_t) * hit.normal)
        .normalize_or_zero();
    let EnvironmentSource::HdrMap(env_map) = nee.environment else {
        // `sample` is `Some` only for `HdrMap` (see `sample_environment_for_nee`), so
        // this is unreachable in practice -- a defensive `return`, not a real branch.
        return;
    };
    let env_rgb = env_map.radiance_rgb(refracted_dir);

    // The medium transmittance a phase-sampled continuation reaching this
    // same boundary would have paid -- see this function's own "Medium transmittance"
    // doc section. Vectorized Beer-Lambert via the same `exp_f32x8` idiom
    // `maybe_scatter_or_extinguish`'s survive branch uses, so the CPU stays
    // self-consistent between the two estimators.
    let hit_t_scaled = hit.t * absorption_path_scale;
    let mut trans_args = [0f32; NUM_CHANNELS];
    for k in 0..NUM_CHANNELS {
        trans_args[k] = -(alphas[k] + sigma_s) * hit_t_scaled;
    }
    let transmittance = crate::simd::exp_f32x8(trans_args);

    let common = t_unpol * phase_val * mis_weight / sample.pdf;
    for k in 0..NUM_CHANNELS {
        let env_k = crate::renderer::env_map::rgb_to_spectral_radiance(env_rgb, lambdas[k]);
        radiance[k] =
            (stokes[k].intensity() * transmittance[k] * common * env_k).mul_add(1.0, radiance[k]);
    }
}

/// 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).
///
/// IS called by [`nee_contribution_hg_scatter`]: NEE evaluates the phase
/// function at the LIGHT-sampled direction (not the phase-sampled one), where no such
/// cancellation applies, plus at the phase-sampled continuation's own direction for the
/// complementary MIS weight -- see that function's doc comment. Also exposed so a Tier 2
/// GPU self-test can pin `phase / pdf == 1.0` at a phase-sampled direction directly.
#[must_use]
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`](super::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()
}