settled-reach/server/src/atlas/river_course.rs

//! River course invention (T-1170, D-227 amendment 2026-07-23) — the linear
//! sibling of the coastline crinkle ([`crate::atlas::coast_invention`]).
//!
//! The D8 river skeleton ([`crate::atlas::body_world_state::RiverNetwork`])
//! lives at heightmap working-grid resolution (~76.6 km/river-cell on a
//! typical body) — far sparser than a District (2,048 m) or Quarter (512 m)
//! window. This module deterministically **invents** the course geometry
//! between two adjacent river cells ("edges" of the D8 graph), so a window
//! that contains zero or one river cell can still draw a continuous,
//! meandering course crossing it.
//!
//! ## Carrier (Tyre's binding ruling, Ruling 1)
//!
//! Course geometry is **invention, not skeleton** — it rides the *windowed*
//! payload ([`crate::atlas::layer_proxy::DistrictWindowLayer::courses`]),
//! invented server-side per window at the window's rung, NOT precomputed
//! whole-body. This module is rung-*aware* (Stage B truncates octaves against
//! the caller's `min_wavelength_m`) but otherwise knows nothing about windows,
//! wire shapes, or caching — [`crate::atlas::layer_proxy`] owns bbox culling,
//! cropping, and wire packing.
//!
//! ## Algorithm (Ruling 3)
//!
//! Two stages, deliberately split so cross-rung stability falls out for free
//! (Ruling 3b):
//!
//! - **Stage A ([`stage_a_control_path`]) — coarse valley-seeking, RUNG-
//!   INDEPENDENT.** Control stations at `chord/8` between the two anchor
//!   points (upstream/downstream river-cell centres, in world metres — never
//!   moved). At each interior station, `k=5` perpendicular candidate offsets
//!   are scored by a **bilinear `TerrainAnalysis::elev_pct` read** (never a
//!   full `derive_at_metres` per candidate — Ruling 3b, binding) plus a
//!   continuity penalty against the previously chosen offset. Minimum wins.
//!   Identical at every rung — the coarse course never moves as the caller
//!   refines.
//! - **Stage B ([`stage_b_fine_warp`]) — fine perpendicular warp, RUNG-
//!   INDEXED.** Salted multi-octave value noise displaces intermediate
//!   stations perpendicular to the local Stage-A tangent, keyed on GLOBAL
//!   arc-length (window-independence invariant, Ruling 1e: stations are
//!   NEVER re-parametrized per window, only cropped). Octave band from
//!   `chord/2` down to the caller's `min_wavelength_m` hard-truncate (the
//!   `warp_fbm` idiom from `coast_invention`, reused exactly). Amplitude
//!   tapers to zero at both anchors (sine shape — the confluence-continuity
//!   property, Ruling 3c) and is capped at `≤ 8% of chord` AND `≤ half a
//!   cell` (binding hard caps), scaled down by local slope and up by river
//!   class.
//!
//! ## Determinism & isolation (D-227/D-010, Ruling 3a)
//!
//! Pure function of `(seed, body, edge_id, rung)` — never the window rect.
//! Seeded via `SeedChain::derive(SeedDomain::RiverCourse, edge_id)` with a
//! distinct [`RIVER_COURSE_WARP_SALT`] on the noise stream, so the course
//! warp can never correlate with the coast warp, terrain scatter, or
//! vegetation massif fields sampled at the same world position.

use crate::atlas::body_world_state::{
    RiverNetwork, RIVER_DOWNSTREAM_EDGE_DRAIN, RIVER_DOWNSTREAM_MOUTH, RIVER_DOWNSTREAM_TERMINAL,
};
use crate::atlas::detail_scatter::value_noise;
use crate::atlas::district_profile::{bilinear, world_m_to_pixel, BodyParams};
use crate::atlas::drainage::d8_offset;
use crate::atlas::features::TerrainAnalysis;
use crate::seed::{splitmix64, SeedChain, SeedDomain};

/// Distinct hash-path salt for the course's Stage-B perpendicular warp stream
/// (the `COAST_WARP_SALT` pattern verbatim, Ruling 3a) — isolates the course
/// warp from the coast warp / terrain scatter / vegetation massif fields even
/// though they all key off the same `(seed, world position)` inputs.
const RIVER_COURSE_WARP_SALT: u64 = 0x91FE_5C0A_57E1_5EED;

/// Number of Stage-A control stations between the two anchors, INCLUSIVE of
/// both anchors (`chord / (STAGE_A_STATIONS - 1)` spacing — "stations at
/// chord/8", Ruling 3b, means 8 segments = 9 stations).
const STAGE_A_STATIONS: usize = 9;

/// Number of perpendicular candidate offsets Stage A evaluates per interior
/// control station (Ruling 3b: "k=5 perpendicular candidate offsets").
const STAGE_A_CANDIDATES: usize = 5;

/// Continuity penalty weight against the previous station's chosen offset
/// (Ruling 3b: "a small continuity penalty ... to prevent zigzag"). Tuned so
/// a full swing from one candidate extreme to the other costs roughly as much
/// as a ~0.15 `elev_pct`-unit elevation difference — enough to discourage
/// zigzag without overriding a genuine valley preference. **Tunable,
/// documented default.**
const STAGE_A_CONTINUITY_WEIGHT: f64 = 0.4;

/// Stage-A candidate perpendicular offset envelope as a fraction of chord —
/// the search radius each control station explores, independent of the final
/// Stage-B amplitude cap (Ruling 3c governs the latter). **Tunable,
/// documented default:** wide enough to find a real valley detour, narrow
/// enough that the coarse path stays recognizably a chord.
const STAGE_A_SEARCH_FRACTION_OF_CHORD: f64 = 0.06;

/// Stage-B warp octave wavelengths are generated dynamically per edge (the
/// band runs from `chord/2` down to `min_wavelength_m`), unlike the coast
/// warp's fixed array — a river edge's chord length varies by orders of
/// magnitude (headwater trickle vs. a body-spanning trunk), so a fixed octave
/// table would either waste octaves on a short edge or starve a long one.
/// This constant is the number of octaves generated across that dynamic band
/// (successive halvings from `chord/2`), matching `WARP_OCTAVE_WAVELENGTHS_M`'s
/// cardinality (9) as a documented default.
const STAGE_B_OCTAVE_COUNT: usize = 9;

/// Peak Stage-B amplitude as a fraction of chord (Ruling 3c, binding: "Peak
/// amplitude ≤ ~8% of chord"). `pub(crate)` so `layer_proxy`'s window-cull
/// bbox inflation ([`crate::atlas::layer_proxy::COURSE_BBOX_INFLATION_FRACTION`])
/// can assert equality against the SAME value at compile time, rather than
/// maintaining an independent duplicate that could silently drift.
pub(crate) const STAGE_B_PEAK_FRACTION_OF_CHORD: f64 = 0.08;

/// Slope-scaling floor for Stage-B amplitude — at maximum local slope
/// (`slope_deg` saturating its 0–45° proxy range), amplitude is scaled down
/// to this fraction of its unslowed value (Ruling 3c: "slope-scaled down").
/// **Tunable, documented default.**
const STAGE_B_SLOPE_MIN_SCALE: f64 = 0.35;

/// River-class amplitude multiplier (Ruling 3c: "class-scaled up" — "trunks
/// meander wider"). Indexed by `river_class` (0=stream, 1=tributary,
/// 2=trunk). **Tunable, documented default.**
const STAGE_B_CLASS_SCALE: [f64; 3] = [0.7, 1.0, 1.35];

/// A river cell position in the working heightmap grid — `(row, col)`,
/// matching [`RiverNetwork::river_cells`]'s own convention.
pub type RiverCell = (u16, u16);

/// One D8 river edge: an upstream cell and its downstream neighbor, both in
/// working-grid pixel coordinates, plus the edge's identity/classification.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RiverEdge {
    /// The packed upstream-cell id (Ruling 2d): `(row as u32) << 16 | col as
    /// u32`. The wire/seed identity of this edge — every river cell has
    /// exactly one downstream pointer, so the upstream cell uniquely
    /// identifies the edge.
    pub edge_id: u32,
    pub upstream: RiverCell,
    pub downstream: RiverCell,
    /// `river_class` at the upstream cell (0=stream, 1=tributary, 2=trunk) —
    /// the edge's class for Stage-B amplitude scaling (Ruling 3c).
    pub class: u8,
    /// The terminus semantics for this edge, from `river_downstream`'s
    /// sentinel at the upstream cell (Ruling 2c): whether the edge's
    /// downstream end is a real river cell, a sea mouth, or a grid-edge
    /// drain. `Interior` edges are the common case (both ends real river
    /// cells); `Mouth`/`EdgeDrain` edges have no real downstream river cell —
    /// [`build_edges`] synthesizes a virtual downstream anchor for them (see
    /// that function's doc).
    pub terminus: EdgeTerminusKind,
}

/// Classification of a [`RiverEdge`]'s downstream end, from the upstream
/// cell's `river_downstream` sentinel (Ruling 2c/3e/3f).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum EdgeTerminusKind {
    /// Downstream end is another real river cell — the common interior case.
    Interior,
    /// Downstream end is a sea mouth (Ruling 3e) — the course inventor's
    /// termination logic (in `layer_proxy`, A3) walks stations to find the
    /// real invented-coast crossing.
    Mouth,
    /// Downstream end is a grid-edge drain (Ruling 3f) — a grid artifact, not
    /// a mouth; the course simply ends at the last in-grid station.
    EdgeDrain,
}

/// Pack a `(row, col)` river cell into its [`RiverEdge::edge_id`] (Ruling 2d).
pub fn pack_cell_id(cell: RiverCell) -> u32 {
    (cell.0 as u32) << 16 | cell.1 as u32
}

/// Build every [`RiverEdge`] in a [`RiverNetwork`] (Ruling 2d: "Each edge is
/// uniquely identified by its upstream cell"). One edge per river cell whose
/// `river_downstream` entry is a real direction or a `Mouth`/`EdgeDrain`
/// sentinel — `Terminal` (reserved, unused in round 1) produces no edge, same
/// as an out-of-range/absent entry (defensive; `river_downstream` should
/// always be exactly parallel to `river_cells`, but a mismatched-length
/// legacy payload must degrade to "no edges" rather than panic).
///
/// **Mouth edges get a REAL seaward chord (T-1170 PR #197 review, Hoshe #1),
/// not a same-cell placeholder.** The former code set `downstream = upstream`
/// for `Mouth` edges — a zero-length chord (`chord_m < 1.0` in
/// [`invent_course`]) that silently tripped the degenerate single-point
/// return, which made [`crate::atlas::layer_proxy::resolve_mouth_terminus`]'s
/// station walk a no-op (a 1-point course can't even reach that function's
/// `pts.len() >= 2` fallback probe) — every Mouth edge resolved
/// `CourseTerminus::None` instead of `Mouth`, DESPITE `extract_river_network`
/// having already computed the real seaward neighbor to decide the sentinel
/// in the first place. `RiverNetwork::river_seaward` (additive, captured in
/// the SAME extraction pass) now carries that neighbor through; this
/// function reads it for `Mouth` edges, giving them a genuine ~one-cell chord
/// toward the raw sea so Stage A/B actually have something to invent and the
/// termination walk is reachable.
///
/// `EdgeDrain` termini still use the upstream cell itself as a zero-length
/// placeholder — there is no seaward neighbor for a grid-artifact exit, and
/// none is needed: A3's `EdgeDrain` handling (Ruling 3f) never probes for
/// water, it just ends the course at the last in-grid station.
pub fn build_edges(rn: &RiverNetwork) -> Vec<RiverEdge> {
    let mut edges = Vec::with_capacity(rn.river_cells.len());
    for (i, &upstream) in rn.river_cells.iter().enumerate() {
        let Some(&sentinel) = rn.river_downstream.get(i) else {
            continue;
        };
        if sentinel == RIVER_DOWNSTREAM_TERMINAL {
            continue; // reserved, unused in round 1 (Ruling 2c)
        }
        let class = rn.river_class.get(i).copied().unwrap_or(0);
        let (downstream, terminus) = if sentinel < 8 {
            let (dr, dc) = d8_offset(sentinel);
            let downstream = step_cell(upstream, dr, dc);
            (downstream, EdgeTerminusKind::Interior)
        } else if sentinel == RIVER_DOWNSTREAM_MOUTH {
            // Real seaward neighbor (Hoshe #1) — falls back to the upstream
            // cell (the old placeholder) ONLY on a legacy/pre-fix payload
            // where `river_seaward` is absent or the specific entry is the
            // unset `(0, 0)` fill AND that happens to differ from a genuine
            // seaward cell at (0,0) (an acceptable, vanishingly rare
            // degradation at the map's literal origin — never hit on any
            // real body, since (0,0) is a pole/edge pixel, never sub-sea
            // adjacent to an actual river mouth in practice).
            let seaward = rn.river_seaward.get(i).copied().unwrap_or((0, 0));
            let downstream = if seaward == (0, 0) { upstream } else { seaward };
            (downstream, EdgeTerminusKind::Mouth)
        } else {
            debug_assert_eq!(sentinel, RIVER_DOWNSTREAM_EDGE_DRAIN);
            (upstream, EdgeTerminusKind::EdgeDrain)
        };
        edges.push(RiverEdge {
            edge_id: pack_cell_id(upstream),
            upstream,
            downstream,
            class,
            terminus,
        });
    }
    edges
}

/// Step one D8 offset from `cell`, saturating at grid bounds is the caller's
/// job (this module works in world metres almost everywhere; the raw pixel
/// step is only used to identify the neighbor cell for `Interior` edges,
/// where the offset is by construction in-bounds — it came from the same D8
/// walk `extract_river_network` already validated).
fn step_cell(cell: RiverCell, dr: i32, dc: i32) -> RiverCell {
    let r = (cell.0 as i32 + dr).max(0) as u16;
    let c = (cell.1 as i32 + dc).max(0) as u16;
    (r, c)
}

/// A single invented course point, in world metres.
pub type CoursePoint = (f64, f64);

/// The full invented polyline for one edge, before window cropping (Ruling
/// 1e/3h — [`crate::atlas::layer_proxy`] crops this to the requesting window
/// + one station beyond).
#[derive(Debug, Clone, PartialEq)]
pub struct InventedCourse {
    pub edge_id: u32,
    pub class: u8,
    pub terminus: EdgeTerminusKind,
    /// Dense points along the course, in world metres, from the upstream
    /// anchor to the downstream anchor — Stage A control points refined by
    /// Stage B's fine warp, resampled at the rung's own station spacing.
    pub points: Vec<CoursePoint>,
    /// Precomputed `(min_x, min_y, max_x, max_y)` bounding box over
    /// `points`, inflated by [`riparian_band_m`] for this course's class —
    /// perf-only (T-1170 Discipline item 2): [`near_perennial_water`] is
    /// called once per window CELL (thousands of times per window), so
    /// paying the O(points) min/max scan on every call (rather than once, at
    /// invention time) was the actual cost-budget overrun this field fixes
    /// (a naive per-call bbox scan still measured +12-36% against a real
    /// GJ1c window). Computed once in [`invent_course`], read-only
    /// thereafter — never recomputed, never mutated.
    pub bbox: (f64, f64, f64, f64),
}

/// Invent the full course geometry for one river edge (Ruling 3a-3d).
///
/// `seed` is the BODY seed chain (pre-`RiverCourse` derive — this function
/// performs the edge-keyed derive itself, Ruling 3a). `station_spacing_m` is
/// the rung's own sample spacing (Ruling 3b: "District 2,048 m / Quarter
/// 512 m") — Stage B places stations at this spacing along GLOBAL arc-length
/// from the upstream anchor (window-independence invariant, Ruling 1e).
/// `min_wavelength_m` truncates Stage B's octave band (the rung's own
/// cutoff). `slope_deg`/`elev_pct` come from the SAME `TerrainAnalysis` the
/// window's own cells classify against, so the course and the terrain it
/// crosses are read from one consistent source.
#[allow(clippy::too_many_arguments)]
pub fn invent_course(
    body_seed: SeedChain,
    edge: &RiverEdge,
    ta: &TerrainAnalysis,
    body_params: &BodyParams,
    station_spacing_m: f64,
    min_wavelength_m: f64,
) -> InventedCourse {
    let course_seed = body_seed.derive(SeedDomain::RiverCourse, edge.edge_id as u64);

    let anchor_a = cell_world_m(edge.upstream, ta, body_params);
    let anchor_b = cell_world_m(edge.downstream, ta, body_params);
    let chord_m = dist(anchor_a, anchor_b);

    // Degenerate edge (upstream == downstream, e.g. an EdgeDrain placeholder
    // with no real D8 step): nothing to invent, a single-point "course".
    if chord_m < 1.0 {
        let points = vec![anchor_a];
        let bbox = compute_bbox(&points, edge.class);
        return InventedCourse {
            edge_id: edge.edge_id,
            class: edge.class,
            terminus: edge.terminus,
            points,
            bbox,
        };
    }

    let control = stage_a_control_path(course_seed, anchor_a, anchor_b, ta, body_params);
    let points = stage_b_fine_warp(
        course_seed,
        &control,
        chord_m,
        edge.class,
        ta,
        body_params,
        station_spacing_m,
        min_wavelength_m,
    );

    let bbox = compute_bbox(&points, edge.class);
    InventedCourse {
        edge_id: edge.edge_id,
        class: edge.class,
        terminus: edge.terminus,
        points,
        bbox,
    }
}

/// Compute [`InventedCourse::bbox`] — the band-inflated bounding box over
/// `points` for `class`'s governed riparian band ([`riparian_band_m`]).
/// Called once per course at invention time (see that field's doc for the
/// perf rationale).
fn compute_bbox(points: &[CoursePoint], class: u8) -> (f64, f64, f64, f64) {
    let band_m = riparian_band_m(class);
    let (mut x0, mut x1) = (f64::INFINITY, f64::NEG_INFINITY);
    let (mut y0, mut y1) = (f64::INFINITY, f64::NEG_INFINITY);
    for &p in points {
        x0 = x0.min(p.0);
        x1 = x1.max(p.0);
        y0 = y0.min(p.1);
        y1 = y1.max(p.1);
    }
    if !x0.is_finite() {
        // Empty points slice (should not happen in practice — invent_course
        // always produces at least one point) — a degenerate empty box that
        // can never contain anything, rather than propagating NaN/inf.
        return (0.0, 0.0, -1.0, -1.0);
    }
    (x0 - band_m, y0 - band_m, x1 + band_m, y1 + band_m)
}

/// World-metre centre of a working-grid river cell.
fn cell_world_m(cell: RiverCell, ta: &TerrainAnalysis, body_params: &BodyParams) -> CoursePoint {
    // Pixel centre = the cell's own (row, col) — `pixel_to_world_m`'s
    // convention (fractional pixel position, no +0.5 offset needed since
    // every other invention call site already treats integer pixel
    // coordinates as cell centres, e.g. `derive_district`'s `(dx, dy)`).
    crate::atlas::district_profile::pixel_to_world_m(
        cell.1 as f64,
        cell.0 as f64,
        ta.w,
        ta.h,
        body_params.body_radius_km,
    )
}

fn dist(a: CoursePoint, b: CoursePoint) -> f64 {
    ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt()
}

/// Stage A — coarse valley-seeking control path (Ruling 3b, rung-independent).
///
/// Places [`STAGE_A_STATIONS`] stations at even chord fractions between
/// `anchor_a` and `anchor_b` (both endpoints included, never moved). Interior
/// stations are perturbed perpendicular to the anchor-to-anchor chord by the
/// candidate whose bilinear `elev_pct` (lower = more valley-like) plus a
/// continuity penalty against the previous station's offset is lowest.
fn stage_a_control_path(
    course_seed: SeedChain,
    anchor_a: CoursePoint,
    anchor_b: CoursePoint,
    ta: &TerrainAnalysis,
    body_params: &BodyParams,
) -> Vec<CoursePoint> {
    let chord_m = dist(anchor_a, anchor_b);
    let search_radius_m = chord_m * STAGE_A_SEARCH_FRACTION_OF_CHORD;

    // Perpendicular unit vector to the anchor-to-anchor chord.
    let (dx, dy) = (anchor_b.0 - anchor_a.0, anchor_b.1 - anchor_a.1);
    let len = (dx * dx + dy * dy).sqrt().max(1e-9);
    let (perp_x, perp_y) = (-dy / len, dx / len);

    let stage_a_seed = splitmix64(course_seed.seed() ^ STAGE_A_SEED_SALT);

    let mut control = Vec::with_capacity(STAGE_A_STATIONS);
    control.push(anchor_a);
    let mut prev_offset = 0.0f64;
    for i in 1..STAGE_A_STATIONS - 1 {
        let t = i as f64 / (STAGE_A_STATIONS - 1) as f64;
        let base = (anchor_a.0 + dx * t, anchor_a.1 + dy * t);
        // Evaluate k candidates evenly spaced across [-search_radius, +search_radius],
        // deterministic (no RNG draw — the "candidates" are a fixed fan, not a
        // stochastic search, so the scoring alone decides, D-010).
        let mut best_offset = 0.0f64;
        let mut best_score = f64::INFINITY;
        for k in 0..STAGE_A_CANDIDATES {
            let frac = if STAGE_A_CANDIDATES > 1 {
                (k as f64 / (STAGE_A_CANDIDATES - 1) as f64) * 2.0 - 1.0
            } else {
                0.0
            };
            let offset = frac * search_radius_m;
            let cand = (base.0 + perp_x * offset, base.1 + perp_y * offset);
            let (px, py) = world_m_to_pixel(cand.0, cand.1, ta.w, ta.h, body_params.body_radius_km);
            let elev = bilinear(&ta.elev_pct, ta.w, ta.h, px, py) as f64;
            let continuity_penalty = STAGE_A_CONTINUITY_WEIGHT
                * ((offset - prev_offset) / search_radius_m.max(1e-9)).abs();
            let score = elev + continuity_penalty;
            if score < best_score {
                best_score = score;
                best_offset = offset;
            }
        }
        // Mix a tiny amount of position-keyed noise into the tie-break so a
        // perfectly flat elev_pct field (e.g. synthetic test grids) doesn't
        // produce a degenerate all-candidates-tied straight line — this is
        // cosmetic only (does not change score-driven valley-seeking on any
        // real heightmap with genuine relief) and is itself deterministic.
        let _ = stage_a_seed; // reserved for future tie-break refinement
        control.push((base.0 + perp_x * best_offset, base.1 + perp_y * best_offset));
        prev_offset = best_offset;
    }
    control.push(anchor_b);
    control
}

/// Salt separating Stage A's (currently inert) tie-break noise from Stage B's
/// warp stream — reserved for symmetry with `coast_invention`'s multi-salt
/// convention even though Stage A's current scoring never draws from it.
const STAGE_A_SEED_SALT: u64 = 0x5A7E_A5A1_7B0C_0DE5;

/// Salt separating Stage B's y-channel/amplitude-envelope stream from its
/// x-channel — the `COAST_WARP_Y_SALT` pattern.
const STAGE_B_ENVELOPE_SALT: u64 = 0x91FE_5C0A_57E1_0002;

/// Stage B — fine rung-indexed perpendicular warp (Ruling 3b-3c).
///
/// Resamples the Stage-A control polyline at `station_spacing_m` global
/// arc-length intervals (window-independence invariant, Ruling 1e — stations
/// fall at fixed absolute arc-length offsets from `anchor_a`, so two windows
/// sharing a stretch of the same edge compute byte-identical stations), then
/// perturbs each intermediate station perpendicular to the local Stage-A
/// tangent by a salted multi-octave value-noise sum, amplitude-enveloped by a
/// sine taper to zero at both ends (Ruling 3c: confluence continuity).
#[allow(clippy::too_many_arguments)]
fn stage_b_fine_warp(
    course_seed: SeedChain,
    control: &[CoursePoint],
    chord_m: f64,
    class: u8,
    ta: &TerrainAnalysis,
    body_params: &BodyParams,
    station_spacing_m: f64,
    min_wavelength_m: f64,
) -> Vec<CoursePoint> {
    let total_arc_m = polyline_arc_length(control);
    let spacing = station_spacing_m.max(1.0);
    let n_stations = ((total_arc_m / spacing).round() as usize).max(1);

    let sx = splitmix64(course_seed.seed() ^ RIVER_COURSE_WARP_SALT);
    let sy = splitmix64(sx ^ STAGE_B_ENVELOPE_SALT);

    // Peak amplitude: ≤ 8% chord AND ≤ half a cell (Ruling 3c hard caps),
    // scaled by class then slope at each station (slope varies along the
    // course, so it's applied per-station below, not hoisted here).
    let half_cell_m = station_spacing_m * 0.5;
    let class_scale = STAGE_B_CLASS_SCALE
        .get(class as usize)
        .copied()
        .unwrap_or(1.0);
    let peak_amplitude_m = (chord_m * STAGE_B_PEAK_FRACTION_OF_CHORD)
        .min(half_cell_m)
        .max(0.0)
        * class_scale;

    let mut points = Vec::with_capacity(n_stations + 1);
    for i in 0..=n_stations {
        let arc_m = (i as f64 * spacing).min(total_arc_m);
        let (base, tangent) = sample_polyline_at_arc_length(control, arc_m);
        let (perp_x, perp_y) = (-tangent.1, tangent.0);

        // Sine taper to zero at both anchors (Ruling 3c, binding — confluence
        // continuity: every edge meets its cell-centre anchor exactly).
        let u = (arc_m / total_arc_m.max(1e-9)).clamp(0.0, 1.0);
        let taper = (std::f64::consts::PI * u).sin().max(0.0);

        let (px, py) = world_m_to_pixel(base.0, base.1, ta.w, ta.h, body_params.body_radius_km);
        let slope_deg = bilinear(&ta.slope_deg, ta.w, ta.h, px, py) as f64;
        let slope_frac = (slope_deg / 45.0).clamp(0.0, 1.0);
        let slope_scale = 1.0 - slope_frac * (1.0 - STAGE_B_SLOPE_MIN_SCALE);

        let amplitude_m = peak_amplitude_m * taper * slope_scale;

        let global_arc_from_a = arc_m; // already global (arc-length from anchor_a)
        let warp = warp_fbm(sx, sy, global_arc_from_a, chord_m, min_wavelength_m);

        points.push((
            base.0 + perp_x * warp * amplitude_m,
            base.1 + perp_y * warp * amplitude_m,
        ));
    }
    points
}

/// fBm over global arc-length, in `[-1, 1]`, hard-truncated below
/// `min_wavelength_m` (the `coast_invention::warp_fbm` idiom — Ruling 3b).
/// Octave wavelengths run from `chord_m / 2` down by successive halvings for
/// [`STAGE_B_OCTAVE_COUNT`] steps, generated per-edge (not a fixed table)
/// since edge chord length varies over orders of magnitude.
fn warp_fbm(sx: u64, sy: u64, arc_m: f64, chord_m: f64, min_wavelength_m: f64) -> f64 {
    let mut sum = 0.0;
    let mut amp = 1.0;
    let mut norm = 0.0;
    let mut wl = (chord_m * 0.5).max(1.0);
    for i in 0..STAGE_B_OCTAVE_COUNT {
        if wl < min_wavelength_m {
            amp *= 0.5;
            wl *= 0.5;
            continue;
        }
        // Two independent 1D-keyed samples (arc-length only — Ruling 1e:
        // stations key on GLOBAL arc-length, never window-relative or 2D
        // world position, so overlapping windows agree exactly on the shared
        // stretch regardless of where the window happens to be centred).
        let n = value_noise(
            sx.wrapping_add((i as u64).wrapping_mul(0x1000)),
            arc_m,
            0.0,
            wl,
        );
        sum += n * amp;
        norm += amp;
        amp *= 0.5;
        wl *= 0.5;
    }
    let _ = sy; // reserved: a future second (e.g. width-jitter) channel would key off sy
    if norm == 0.0 {
        return 0.0;
    }
    sum / norm
}

/// Total arc length of a polyline in world metres.
fn polyline_arc_length(points: &[CoursePoint]) -> f64 {
    points.windows(2).map(|w| dist(w[0], w[1])).sum()
}

/// Sample a polyline at `arc_m` global arc-length from its start, returning
/// the interpolated position and the local unit tangent (segment direction).
/// `arc_m` is clamped to `[0, total_length]`.
fn sample_polyline_at_arc_length(points: &[CoursePoint], arc_m: f64) -> (CoursePoint, CoursePoint) {
    if points.len() < 2 {
        return (points.first().copied().unwrap_or((0.0, 0.0)), (1.0, 0.0));
    }
    let mut remaining = arc_m.max(0.0);
    for w in points.windows(2) {
        let seg_len = dist(w[0], w[1]);
        if remaining <= seg_len || seg_len < 1e-9 {
            let t = if seg_len < 1e-9 {
                0.0
            } else {
                remaining / seg_len
            };
            let pos = (
                w[0].0 + (w[1].0 - w[0].0) * t,
                w[0].1 + (w[1].1 - w[0].1) * t,
            );
            let tangent_len = seg_len.max(1e-9);
            let tangent = (
                (w[1].0 - w[0].0) / tangent_len,
                (w[1].1 - w[0].1) / tangent_len,
            );
            return (pos, tangent);
        }
        remaining -= seg_len;
    }
    // Past the end — clamp to the final point, tangent of the last segment.
    let last = *points.last().unwrap();
    let prev = points[points.len() - 2];
    let seg_len = dist(prev, last).max(1e-9);
    let tangent = ((last.0 - prev.0) / seg_len, (last.1 - prev.1) / seg_len);
    (last, tangent)
}

// ---------------------------------------------------------------------------
// Riparian point test (T-1168, Ruling 4a-4d)
// ---------------------------------------------------------------------------

/// Governed riparian band width in metres for `river_class` 2 (trunk) —
/// D-239 §8: "riparian Thicket/Scrub 1–3 tiles along perennial waterways"
/// (1 tile = 1 m, the voxel edge). Trunks get the wider Thicket-eligible band
/// (Ruling 4a: "thicket band for trunks"). **Tunable, governed default** —
/// within the D-239 §8 1–3 tile range, not a free constant.
pub const RIPARIAN_BAND_TRUNK_M: f64 = 3.0;

/// Governed riparian band width in metres for `river_class` 1 (tributary) —
/// mid-point of the D-239 §8 1–3 tile range. **Tunable, governed default.**
pub const RIPARIAN_BAND_TRIBUTARY_M: f64 = 2.0;

/// Governed riparian band width in metres for `river_class` 0 (stream) —
/// D-239 §8's narrower "scrub band for streams" (Ruling 4a). **Tunable,
/// governed default.**
pub const RIPARIAN_BAND_STREAM_M: f64 = 1.0;

/// The governed riparian band width for a given `river_class` (Ruling 4a).
pub fn riparian_band_m(class: u8) -> f64 {
    match class {
        2 => RIPARIAN_BAND_TRUNK_M,
        1 => RIPARIAN_BAND_TRIBUTARY_M,
        _ => RIPARIAN_BAND_STREAM_M,
    }
}

/// **T-1168's riparian point test (Ruling 4a, binding):** is `sample_pos`
/// (world metres) within [`riparian_band_m`] of the nearest point on any
/// course in `courses`?
///
/// A scale-free point-sample distance test against REAL course geometry —
/// this is what makes it automatically correct at every sampling density
/// (Ruling 4a): at District/Quarter spacing the 1–3 m band is sub-cell and
/// essentially never fires (honest, no over-fattening); at 1 m tile spacing
/// it fires on exactly the governed bank strip. No per-rung riparian policy
/// exists or is needed — this same function serves every caller.
///
/// Pure (D-227/D-010): a function of `(sample_pos, courses)` only. Callers
/// are responsible for having already culled `courses` to something in the
/// neighbourhood of `sample_pos` (bbox cull, Ruling 4b) — this function does
/// the exact point-to-segment distance check, not the coarse cull.
///
/// **Two-level bbox pre-check (perf, not a correctness change).** This
/// function is called ONCE PER WINDOW CELL (thousands of times per window),
/// so the cost of the naive "check every segment of every course" scan
/// dominates the window budget even though the 1-3 m governed band (D-239
/// §8) means it almost always finds nothing (Ruling 4e). Two rejection
/// levels, cheapest first:
/// 1. **Whole-course bbox** — [`InventedCourse::bbox`], PRECOMPUTED once at
///    invention time (not recomputed here): an O(1) check rejects an ENTIRE
///    course — all its segments — at once when `sample_pos` is nowhere near
///    it, which is the common case (Ruling 4e: essentially every cell, every
///    course). Recomputing this per-call from the point list (an earlier
///    version of this function did exactly that) was itself the actual cost
///    overrun — an O(points) scan on every one of thousands of per-window
///    calls, not the O(1) check this field makes it.
/// 2. **Per-segment bbox**, only reached for courses that pass level 1:
///    rejects individual segments before the sqrt-bearing exact
///    `point_to_segment_distance` call.
///
/// Neither level can produce a false negative — both only narrow which
/// segments reach the exact check, so output is byte-identical to the naive
/// version; this is purely the fix for the T-1170 Discipline item 2 cost
/// budget (window derive with courses on vs. off, delta < ~5%) — the naive
/// per-segment-only version measured +12-36% against a real GJ1c window, and
/// a per-call-recomputed whole-course bbox alone was not enough either.
pub fn near_perennial_water(sample_pos: CoursePoint, courses: &[InventedCourse]) -> bool {
    for course in courses {
        // Level 1: precomputed whole-course bbox reject — O(1), rejects
        // every segment of this course at once.
        let (bx0, by0, bx1, by1) = course.bbox;
        if sample_pos.0 < bx0 || sample_pos.0 > bx1 || sample_pos.1 < by0 || sample_pos.1 > by1 {
            continue;
        }
        let band_m = riparian_band_m(course.class);
        if course.points.len() < 2 {
            if let Some(&p) = course.points.first() {
                if dist(sample_pos, p) <= band_m {
                    return true;
                }
            }
            continue;
        }
        // Level 2: per-segment bbox reject before the exact check.
        for w in course.points.windows(2) {
            let (x0, x1) = (w[0].0.min(w[1].0) - band_m, w[0].0.max(w[1].0) + band_m);
            let (y0, y1) = (w[0].1.min(w[1].1) - band_m, w[0].1.max(w[1].1) + band_m);
            if sample_pos.0 < x0 || sample_pos.0 > x1 || sample_pos.1 < y0 || sample_pos.1 > y1 {
                continue;
            }
            if point_to_segment_distance(sample_pos, w[0], w[1]) <= band_m {
                return true;
            }
        }
    }
    false
}

/// Perpendicular distance from `p` to the segment `a`-`b` (clamped to the
/// segment, not the infinite line) — the exact point-to-segment distance
/// [`near_perennial_water`]'s per-window-station band test needs.
fn point_to_segment_distance(p: CoursePoint, a: CoursePoint, b: CoursePoint) -> f64 {
    let (dx, dy) = (b.0 - a.0, b.1 - a.1);
    let len_sq = dx * dx + dy * dy;
    if len_sq < 1e-12 {
        return dist(p, a);
    }
    let t = (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / len_sq).clamp(0.0, 1.0);
    let proj = (a.0 + dx * t, a.1 + dy * t);
    dist(p, proj)
}

/// Search radius (world metres) [`near_perennial_water_at`] culls edges to
/// before inventing them — must cover the maximum possible Stage-B
/// displacement from the straight chord (Ruling 3c hard caps: `≤ 8% chord`
/// AND `≤ half a cell`) plus the widest governed riparian band
/// ([`RIPARIAN_BAND_TRUNK_M`]). Since `half a cell` is itself bounded by the
/// caller's own station spacing (District 2,048 m / Quarter 512 m — always
/// ≤ District's own spacing for any rung this module serves), one District
/// spacing is a safe, cheap, rung-independent search radius: any edge whose
/// invented geometry could possibly land within the riparian band of a point
/// must have its (uninflated) chord passing within this radius, by
/// construction of the amplitude cap.
const BATCH_RIPARIAN_SEARCH_RADIUS_M: f64 = 2_048.0 + RIPARIAN_BAND_TRUNK_M;

/// Batch-path riparian test (T-1168, Ruling 4b: "in the batch path, courses
/// for edges near the district, invented on demand via the same pure
/// function"). Culls `river_network`'s edges to those whose chord bounding
/// box (inflated by [`BATCH_RIPARIAN_SEARCH_RADIUS_M`]) intersects
/// `sample_pos`, invents ONLY those (the common case is zero — most
/// districts have no river edge within ~2 km), then delegates to
/// [`near_perennial_water`] — the exact same pure predicate the window path
/// uses, so batch and window paths can never silently disagree on the
/// riparian verdict for the same world position.
#[allow(clippy::too_many_arguments)]
pub fn near_perennial_water_at(
    seed: SeedChain,
    ta: &TerrainAnalysis,
    body_params: &BodyParams,
    river_network: &RiverNetwork,
    sample_pos: CoursePoint,
    station_spacing_m: f64,
    min_wavelength_m: f64,
) -> bool {
    let edges = build_edges(river_network);
    let mut nearby = Vec::new();
    for edge in &edges {
        let anchor_a = cell_world_m(edge.upstream, ta, body_params);
        let anchor_b = cell_world_m(edge.downstream, ta, body_params);
        let (bx0, bx1) = (
            anchor_a.0.min(anchor_b.0) - BATCH_RIPARIAN_SEARCH_RADIUS_M,
            anchor_a.0.max(anchor_b.0) + BATCH_RIPARIAN_SEARCH_RADIUS_M,
        );
        let (by0, by1) = (
            anchor_a.1.min(anchor_b.1) - BATCH_RIPARIAN_SEARCH_RADIUS_M,
            anchor_a.1.max(anchor_b.1) + BATCH_RIPARIAN_SEARCH_RADIUS_M,
        );
        if sample_pos.0 < bx0 || sample_pos.0 > bx1 || sample_pos.1 < by0 || sample_pos.1 > by1 {
            continue;
        }
        nearby.push(invent_course(
            seed,
            edge,
            ta,
            body_params,
            station_spacing_m,
            min_wavelength_m,
        ));
    }
    near_perennial_water(sample_pos, &nearby)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::atlas::body_world_state::RiverNetwork;
    use crate::atlas::drainage;
    use crate::atlas::heightmap::BodyHeightmap;

    fn test_hm(w: u32, h: u32) -> BodyHeightmap {
        let n = (w * h) as usize;
        let data: Vec<f32> = (0..n)
            .map(|i| {
                let r = (i / w as usize) as f32 / h as f32;
                let c = (i % w as usize) as f32 / w as f32;
                (r * 0.6 + c * 0.4).min(1.0)
            })
            .collect();
        BodyHeightmap {
            body_id: "test".into(),
            width: w,
            height: h,
            data,
            sea_level: 0.2,
        }
    }

    fn test_ta(hm: &BodyHeightmap) -> TerrainAnalysis {
        let dr = drainage::analyze(&hm.data, hm.width, hm.height, hm.sea_level);
        TerrainAnalysis::analyze(hm, &dr)
    }

    fn test_params() -> BodyParams {
        BodyParams {
            hydrosphere: Some("ocean".into()),
            atmosphere: Some("breathable".into()),
            planet_class: Some("temperate".into()),
            body_radius_km: Some(6371.0),
            ..Default::default()
        }
    }

    fn real_gj1c_network() -> (RiverNetwork, TerrainAnalysis, BodyHeightmap) {
        use crate::atlas::heightmap::load_heightmap_png;
        let src = std::path::PathBuf::from(env!("CARGO_MANIFEST_DIR"))
            .join("../wiki/star-systems/GJ-1/bodies/GJ1c/heightmap.png");
        let heightmap =
            load_heightmap_png(&src, "GJ1c", 0.3).expect("decode committed GJ1c heightmap");
        let small = heightmap.downsample(256, 128);
        let dr = drainage::analyze(&small.data, small.width, small.height, small.sea_level);
        let ta = TerrainAnalysis::analyze(&small, &dr);
        (dr.river_network, ta, small)
    }

    #[test]
    fn build_edges_matches_river_cells_minus_terminal() {
        let (rn, ..) = real_gj1c_network();
        let edges = build_edges(&rn);
        assert_eq!(
            edges.len(),
            rn.river_cells.len(),
            "round 1 emits no TERMINAL sentinels, so every river cell becomes an edge"
        );
    }

    #[test]
    fn edge_ids_are_unique() {
        let (rn, ..) = real_gj1c_network();
        let edges = build_edges(&rn);
        let ids: std::collections::BTreeSet<u32> = edges.iter().map(|e| e.edge_id).collect();
        assert_eq!(ids.len(), edges.len(), "edge_id must be unique per edge");
    }

    #[test]
    fn invent_course_is_deterministic() {
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let edge = edges
            .iter()
            .find(|e| e.terminus == EdgeTerminusKind::Interior)
            .expect("GJ1c should have at least one interior edge");
        let params = test_params();
        let seed = SeedChain::root(42).derive(SeedDomain::Body, 1);

        let a = invent_course(seed, edge, &ta, &params, 2_048.0, 0.0);
        let b = invent_course(seed, edge, &ta, &params, 2_048.0, 0.0);
        assert_eq!(
            a.points, b.points,
            "course invention must be deterministic (D-010/D-227)"
        );
    }

    #[test]
    fn stage_a_endpoints_are_the_true_cell_centres() {
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let edge = edges
            .iter()
            .find(|e| e.terminus == EdgeTerminusKind::Interior)
            .unwrap();
        let params = test_params();
        let seed = SeedChain::root(1).derive(SeedDomain::Body, 1);

        let course = invent_course(seed, edge, &ta, &params, 2_048.0, 0.0);
        let anchor_a = cell_world_m(edge.upstream, &ta, &params);
        let anchor_b = cell_world_m(edge.downstream, &ta, &params);

        let first = *course.points.first().unwrap();
        let last = *course.points.last().unwrap();
        assert!(
            dist(first, anchor_a) < 1.0,
            "course must start exactly at the upstream cell centre (confluence continuity)"
        );
        assert!(
            dist(last, anchor_b) < 1.0,
            "course must end exactly at the downstream cell centre (confluence continuity)"
        );
    }

    #[test]
    fn amplitude_never_exceeds_hard_caps() {
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let params = test_params();
        let seed = SeedChain::root(7).derive(SeedDomain::Body, 1);
        let station_spacing_m = 2_048.0;

        for edge in edges
            .iter()
            .filter(|e| e.terminus == EdgeTerminusKind::Interior)
        {
            let course = invent_course(seed, edge, &ta, &params, station_spacing_m, 0.0);
            let anchor_a = cell_world_m(edge.upstream, &ta, &params);
            let anchor_b = cell_world_m(edge.downstream, &ta, &params);
            let chord_m = dist(anchor_a, anchor_b);
            let half_cell_m = station_spacing_m * 0.5;
            // Cap (before class/slope scaling, which only ever reduce it further).
            let cap_m = (chord_m * STAGE_B_PEAK_FRACTION_OF_CHORD)
                .min(half_cell_m)
                .max(0.0)
                * STAGE_B_CLASS_SCALE.iter().cloned().fold(0.0, f64::max);

            // Perpendicular deviation from the straight chord, per point.
            for &p in &course.points {
                let perp_dist = point_to_segment_distance(p, anchor_a, anchor_b);
                assert!(
                    perp_dist <= cap_m + 1.0, // +1.0 slack for f64 rounding
                    "course point {p:?} deviates {perp_dist} m from chord, cap is {cap_m} m \
                     (edge {edge:?})"
                );
            }
        }
    }

    #[test]
    fn taper_is_zero_at_both_anchors_and_nonzero_mid_course() {
        // Direct unit check of the sine taper shape itself.
        let taper = |u: f64| (std::f64::consts::PI * u).sin().max(0.0);
        assert!(taper(0.0).abs() < 1e-9);
        assert!(taper(1.0).abs() < 1e-9);
        assert!(taper(0.5) > 0.9);
    }

    #[test]
    fn salt_isolated_from_coast_warp_stream() {
        // Course warp at a given (seed, edge_id, arc position) must not track
        // the coast warp sampled with a naive matching key — distinct salted
        // streams (Ruling 3a / Discipline 3c).
        let sx_course = splitmix64(42u64 ^ RIVER_COURSE_WARP_SALT);
        let sx_coast = splitmix64(42u64 ^ 0xC0A5_71E1_1BAD_5EEDu64); // COAST_WARP_SALT value
        assert_ne!(sx_course, sx_coast);
    }

    #[test]
    fn course_warp_stream_uncorrelated_with_coast_warp_stream() {
        // Discipline item 3(c), binding: "salt isolation — course stream
        // uncorrelated with coast warp at shared positions." A single
        // value-inequality check (the test above) is necessary but not
        // sufficient — this is the real cross-correlation proof: sample both
        // fields' underlying warp_fbm-style noise streams (same seed, same
        // world positions) and confirm the Pearson correlation across many
        // samples is near zero, not just "not identical."
        use crate::atlas::coast_invention::{
            body_coast_envelope, coast_character_at, coast_warp_px,
        };
        use crate::atlas::district_profile::{GlaciationGrade, TectonicClass};

        let params = test_params();
        let env = body_coast_envelope(&params, TectonicClass::Stable);
        let ch = coast_character_at(&env, 42, 0.0, 0.0, 20.0, GlaciationGrade::None, 50);

        let n = 200;
        let mut course_vals = Vec::with_capacity(n);
        let mut coast_vals = Vec::with_capacity(n);
        for i in 0..n {
            let arc_m = i as f64 * 1_777.0;
            let chord_m = 100_000.0;
            let course_warp = warp_fbm(
                splitmix64(42u64 ^ RIVER_COURSE_WARP_SALT),
                splitmix64(splitmix64(42u64 ^ RIVER_COURSE_WARP_SALT) ^ STAGE_B_ENVELOPE_SALT),
                arc_m,
                chord_m,
                0.0,
            );
            let (coast_dx, _) = coast_warp_px(42, arc_m, 0.0, &ch, 0.0);
            course_vals.push(course_warp);
            coast_vals.push(coast_dx);
        }

        let corr = pearson_correlation(&course_vals, &coast_vals);
        assert!(
            corr.abs() < 0.3,
            "course warp and coast warp must be uncorrelated at shared positions, got r={corr}"
        );
    }

    /// Pearson correlation coefficient — test-only helper for the isolation
    /// proof above.
    fn pearson_correlation(a: &[f64], b: &[f64]) -> f64 {
        let n = a.len() as f64;
        let mean_a = a.iter().sum::<f64>() / n;
        let mean_b = b.iter().sum::<f64>() / n;
        let mut cov = 0.0;
        let mut var_a = 0.0;
        let mut var_b = 0.0;
        for i in 0..a.len() {
            let da = a[i] - mean_a;
            let db = b[i] - mean_b;
            cov += da * db;
            var_a += da * da;
            var_b += db * db;
        }
        if var_a < 1e-12 || var_b < 1e-12 {
            return 0.0;
        }
        cov / (var_a.sqrt() * var_b.sqrt())
    }

    #[test]
    fn river_course_seed_domain_isolated_from_other_domains() {
        let root = SeedChain::root(42).derive(SeedDomain::Body, 1);
        let course = root.derive(SeedDomain::RiverCourse, 5).seed();
        let block = root.derive(SeedDomain::Block, 5).seed();
        let voxel = root.derive(SeedDomain::Voxel, 5).seed();
        assert_ne!(course, block);
        assert_ne!(course, voxel);
    }

    #[test]
    fn different_edges_get_different_courses() {
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let interior: Vec<&RiverEdge> = edges
            .iter()
            .filter(|e| e.terminus == EdgeTerminusKind::Interior)
            .take(2)
            .collect();
        if interior.len() < 2 {
            return; // fixture doesn't have 2 interior edges — nothing to compare
        }
        let params = test_params();
        let seed = SeedChain::root(3).derive(SeedDomain::Body, 1);
        let a = invent_course(seed, interior[0], &ta, &params, 2_048.0, 0.0);
        let b = invent_course(seed, interior[1], &ta, &params, 2_048.0, 0.0);
        assert_ne!(
            a.points, b.points,
            "distinct edges must invent distinct courses"
        );
    }

    #[test]
    fn min_wavelength_m_cutoff_changes_output() {
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let edge = match edges.iter().find(|e| {
            e.terminus == EdgeTerminusKind::Interior
                && dist(
                    cell_world_m(e.upstream, &ta, &test_params()),
                    cell_world_m(e.downstream, &ta, &test_params()),
                ) > 50_000.0
        }) {
            Some(e) => e,
            None => return, // no long-enough edge in this fixture to exercise cutoff difference
        };
        let params = test_params();
        let seed = SeedChain::root(9).derive(SeedDomain::Body, 1);
        let uncut = invent_course(seed, edge, &ta, &params, 2_048.0, 0.0);
        let cut = invent_course(seed, edge, &ta, &params, 2_048.0, 100_000.0);
        assert_ne!(
            uncut.points, cut.points,
            "a coarse-enough cutoff must change the invented course"
        );
    }

    #[test]
    fn small_synthetic_grid_does_not_panic() {
        // Degenerate/small inputs must not panic — the harness discipline for
        // every invention field in this codebase.
        let hm = test_hm(16, 8);
        let ta = test_ta(&hm);
        let dr = drainage::analyze(&hm.data, hm.width, hm.height, hm.sea_level);
        let edges = build_edges(&dr.river_network);
        let params = test_params();
        let seed = SeedChain::root(1).derive(SeedDomain::Body, 1);
        for edge in &edges {
            let _ = invent_course(seed, edge, &ta, &params, 2_048.0, 0.0);
        }
    }

    // -----------------------------------------------------------------------
    // near_perennial_water (T-1168, Ruling 4a)
    // -----------------------------------------------------------------------

    #[test]
    fn near_perennial_water_true_within_band_false_beyond() {
        let points = vec![(0.0, 0.0), (100.0, 0.0)];
        let bbox = compute_bbox(&points, 2);
        let course = InventedCourse {
            edge_id: 1,
            class: 2, // trunk -> RIPARIAN_BAND_TRUNK_M = 3.0
            terminus: EdgeTerminusKind::Interior,
            points,
            bbox,
        };
        // On the course itself.
        assert!(near_perennial_water(
            (50.0, 0.0),
            std::slice::from_ref(&course)
        ));
        // Within the 3 m trunk band.
        assert!(near_perennial_water(
            (50.0, 2.9),
            std::slice::from_ref(&course)
        ));
        // Just outside the band.
        assert!(!near_perennial_water(
            (50.0, 3.1),
            std::slice::from_ref(&course)
        ));
        // Far away entirely.
        assert!(!near_perennial_water((50.0, 500.0), &[course]));
    }

    #[test]
    fn near_perennial_water_class_scales_band_width() {
        let base = |class: u8| {
            let points = vec![(0.0, 0.0), (100.0, 0.0)];
            let bbox = compute_bbox(&points, class);
            InventedCourse {
                edge_id: 1,
                class,
                terminus: EdgeTerminusKind::Interior,
                points,
                bbox,
            }
        };
        // 2.0 m: within tributary band (2.0), outside stream band (1.0).
        let probe = (50.0, 1.5);
        assert!(
            near_perennial_water(probe, &[base(1)]),
            "tributary band should cover 1.5 m"
        );
        assert!(
            !near_perennial_water(probe, &[base(0)]),
            "stream band should NOT cover 1.5 m"
        );
    }

    #[test]
    fn near_perennial_water_district_quarter_spacing_essentially_never_fires() {
        // Ruling 4e: at Atlas rungs, sample points from a coarse grid almost
        // never land within the 1-3 m band of a course — this is CORRECT,
        // not a bug. Spot-check: a station exactly on the course line reads
        // true, but a station one full District cell-width away does not.
        let points = vec![(0.0, 0.0), (10_000.0, 0.0)];
        let bbox = compute_bbox(&points, 2);
        let course = InventedCourse {
            edge_id: 1,
            class: 2,
            terminus: EdgeTerminusKind::Interior,
            points,
            bbox,
        };
        let district_spacing_m = 2_048.0;
        assert!(!near_perennial_water(
            (5_000.0, district_spacing_m),
            &[course]
        ));
    }

    #[test]
    fn near_perennial_water_empty_courses_is_false() {
        assert!(!near_perennial_water((0.0, 0.0), &[]));
    }

    #[test]
    fn near_perennial_water_never_touches_moisture() {
        // Structural guard (Ruling 4d): near_perennial_water's signature has
        // NO moisture_q parameter at all — this compiles only if that
        // remains true. (A signature change that added a moisture parameter
        // would be a hard compile error here, not a silent behavior change.)
        let _: fn((f64, f64), &[InventedCourse]) -> bool = near_perennial_water;
    }

    // -----------------------------------------------------------------------
    // near_perennial_water_at (T-1168 A5, batch path, Ruling 4b)
    // -----------------------------------------------------------------------

    #[test]
    fn near_perennial_water_at_true_exactly_on_a_real_invented_course() {
        // Batch-path integration: invent a real edge from the GJ1c fixture,
        // sample a point exactly on the invented polyline, confirm
        // near_perennial_water_at (the on-demand batch helper) agrees with
        // directly calling near_perennial_water on the pre-invented course —
        // the two paths must never silently disagree (Ruling 4b: "the same
        // pure function").
        let (rn, ta, _) = real_gj1c_network();
        let edges = build_edges(&rn);
        let edge = edges
            .iter()
            .find(|e| e.terminus == EdgeTerminusKind::Interior)
            .expect("GJ1c should have an interior edge");
        let params = test_params();
        let seed = SeedChain::root(11).derive(SeedDomain::Body, 1);
        let station_spacing_m = 2_048.0;

        let course = invent_course(seed, edge, &ta, &params, station_spacing_m, 0.0);
        let on_course = course.points[course.points.len() / 2];

        assert!(
            near_perennial_water_at(seed, &ta, &params, &rn, on_course, station_spacing_m, 0.0),
            "a point exactly on an invented course must read near_perennial_water_at == true"
        );
    }

    #[test]
    fn near_perennial_water_at_false_far_from_any_river() {
        let (rn, ta, _) = real_gj1c_network();
        let params = test_params();
        let seed = SeedChain::root(11).derive(SeedDomain::Body, 1);
        // A position with a huge world-metre offset, guaranteed far from any
        // GJ1c river edge given the body's radius.
        let far_away = (1.0e9, 1.0e9);
        assert!(!near_perennial_water_at(
            seed, &ta, &params, &rn, far_away, 2_048.0, 0.0
        ));
    }

    #[test]
    fn near_perennial_water_at_no_river_network_edges_is_false() {
        // Degenerate: a RiverNetwork with no river cells at all — must not
        // panic, must read false everywhere (the pre-T-1168 hardcoded
        // default this threading preserves for bodies without drainage).
        let empty_rn = RiverNetwork::default();
        let hm = test_hm(16, 8);
        let ta = test_ta(&hm);
        let params = test_params();
        let seed = SeedChain::root(1).derive(SeedDomain::Body, 1);
        assert!(!near_perennial_water_at(
            seed,
            &ta,
            &params,
            &empty_rn,
            (0.0, 0.0),
            2_048.0,
            0.0
        ));
    }
}