BREP_kernel 0.10.0

use super::*;
use crate::{KernelRefusal, KernelStage, OrRefuse, RefusalClass};

// ====================================================================
// OPEN edges (§6.9 transverse construction)
// ====================================================================

/// March an OPEN edge with a small overshoot past both ends (extended
/// curve/surface evaluation) so the support fits reach the end-face
/// crossings.
///
/// Returns the stations together with the INDICES of the two the march
/// snapped exactly onto the edge's start and end parameters — the sections
/// that sit in each vertex's normal plane.  Those indices are the march's own
/// bookkeeping and cannot be recovered from the stations afterwards (see
/// [`vertex_stations`]).
pub(in crate::blend) fn march_open_stations(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
    radius_at: &dyn Fn(f64) -> f64,
    overshoot_fraction: f64,
) -> Result<(Vec<Station>, [usize; 2]), KernelRefusal> {
    march_open_stations_ending(
        edge,
        first,
        second,
        radius_at,
        overshoot_fraction,
        [edge.t0, edge.t1],
        STATIONS,
        [false, false],
        [false, false],
        None,
    )
    .map(|(stations, vertex_indices, _)| (stations, vertex_indices))
}

/// [`march_open_stations`] with the two snapped stations placed at `ends`
/// instead of at the edge's own vertices.
///
/// The snapped stations are where the rows are broken, so they belong where
/// the blend actually ENDS and the marched data stops being the blend: a
/// flush join across a seam stops a stripe at the ball seated on the seam,
/// which is where its side contact runs off its carrier onto the ruled
/// continuation — the same curvature jump [`fit_open_rows`] breaks at the
/// vertex for, one section along.  `ends` must lie inside the overshot span.
///
/// `poles` marks an end whose section is a POINT ([`pole_end`]): the two mates
/// are tangent at that vertex from the ball's side, so the ball seated there
/// touches both at the vertex itself.  The station Newton is singular there —
/// both offset carriers are tangent, so every column of its Jacobian lies in
/// one plane — and past it the two contacts cross over, so such an end is
/// marched with NO overshoot and its station is the pole itself, set from the
/// edge's own pcurves rather than solved. Its end must be the edge's vertex.
///
/// `bounded_ends` permits stopping at the last regular station before a
/// singular carrier pole. Such rows are only usable after the caller verifies
/// that the real end-face intersections lie strictly inside the fitted span;
/// their clamped end indices are fit boundaries, not vertex sections.
#[allow(clippy::too_many_arguments)]
pub(in crate::blend) fn march_open_stations_ending(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
    radius_at: &dyn Fn(f64) -> f64,
    overshoot_fraction: f64,
    ends: [f64; 2],
    station_count: usize,
    poles: [bool; 2],
    bounded_ends: [bool; 2],
    refine_bar: Option<f64>,
) -> Result<(Vec<Station>, [usize; 2], StationRefinement), KernelRefusal> {
    march_open_stations_core(
        edge, first, second, radius_at, overshoot_fraction, ends, station_count, poles, bounded_ends, refine_bar, 0, true,
    )
    .map(|(stations, vertex_indices, refinement, _)| (stations, vertex_indices, refinement))
}

/// The exact ball centres BETWEEN an open march's final stations: for every
/// interval, at fractions `step / subdivisions` of its edge parameter
/// (`step = 1 .. subdivisions - 1`), each the continuation of the interval's
/// LEFT station under the march's own branch guard (`None` where it did not
/// solve, or at a pole). Empty when not asked for.
#[derive(Clone, Debug, Default)]
pub(in crate::blend) struct OpenDenseCentres {
    pub(in crate::blend) subdivisions: usize,
    pub(in crate::blend) centres: Vec<Vec<Option<Vec3>>>,
    /// The final stations' NATIVE edge parameters: interval `i`'s centre `j`
    /// was solved at `t_i + (t_{i+1} - t_i) · j / subdivisions` (never the
    /// wall's chord parameters).
    pub(in crate::blend) station_ts: Vec<f64>,
}

/// The chord road's seed at `t`: the quadratic through the row's last three
/// centres (linear through two), or `None` with fewer than two or a
/// non-finite result, which leaves the caller's own seed in charge.
fn predict_station_uv(history: &[(f64, [f64; 4])], t: f64) -> Option<[f64; 4]> {
    let n = history.len();
    if n < 2 {
        return None;
    }
    let tail = &history[n.saturating_sub(3)..];
    let mut uv = [0.0f64; 4];
    for (i, &(ti, ui)) in tail.iter().enumerate() {
        let mut weight = 1.0;
        for (j, &(tj, _)) in tail.iter().enumerate() {
            if i != j {
                let denominator = ti - tj;
                if denominator == 0.0 {
                    return None;
                }
                weight *= (t - tj) / denominator;
            }
        }
        for k in 0..4 {
            uv[k] += weight * ui[k];
        }
    }
    uv.iter().all(|value| value.is_finite()).then_some(uv)
}

/// `BREP_BLEND_DENSE_CHAIN=0` seeds every dense centre from its interval's
/// start station, as before 2026-10-17; unset or anything else chains them
/// (see the dense loop). Read once.
fn dense_chain_enabled() -> bool {
    static ENABLED: std::sync::OnceLock<bool> = std::sync::OnceLock::new();
    *ENABLED.get_or_init(|| std::env::var("BREP_BLEND_DENSE_CHAIN").as_deref() != Ok("0"))
}

/// [`march_open_stations_ending`], also solving each final interval's
/// [`OpenDenseCentres`] at `dense` subdivisions (0: none) — the exact centres
/// a wall verdict reads the fitted wall against.
#[allow(clippy::too_many_arguments)]
/// `BREP_BLEND_EXACT_STRIPES=1`/`on` translates the stations of a straight
/// stripe between two affine planes from its middle one (lane R2 row 2b,
/// 2026-10-18), `=0`/`off` (the default) marches them as before;
/// `set_exact_stripes_for_test` forces it on a test's own thread.
const EXACT_STRIPES_DEFAULT: bool = false;

thread_local! {
    static EXACT_STRIPES_FOR_TEST: std::cell::Cell<Option<bool>> = const { std::cell::Cell::new(None) };
}

pub(in crate::blend) fn exact_stripes_on() -> bool {
    if let Some(forced) = EXACT_STRIPES_FOR_TEST.with(|forced| forced.get()) {
        return forced;
    }
    match std::env::var("BREP_BLEND_EXACT_STRIPES").as_deref() {
        Ok("0") | Ok("off") => false,
        Ok("1") | Ok("on") => true,
        _ => EXACT_STRIPES_DEFAULT,
    }
}


/// The translation a straight stripe between two affine planes marches by:
/// the edge's unit direction and, per mate, the rows of the plane's
/// left-inverse Jacobian (`(JᵀJ)⁻¹Jᵀ`, J = [s_u s_v]), which map a shift
/// lying in the plane to its (du, dv).
#[derive(Clone, Copy)]
struct TranslatedStations {
    direction: Vec3,
    inverses: [[Vec3; 2]; 2],
}

/// [`TranslatedStations`] when the march may translate its stations
/// (`exact_stripes_on`): both carriers affine planes
/// (`affine_plane_frame`), the edge straight and both planes invariant along
/// it (`extrusion_direction`), the radius the same bits at every station, no
/// pole, and none of the refinement, dense centres or fold probe asked (the
/// exact lane asks none). Never inside the open-pad composition
/// (`open_pad_network_on`), whose landed bodies keep their march.
#[allow(clippy::too_many_arguments)]
fn translated_stations(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
    radius_at: &dyn Fn(f64) -> f64,
    station_ts: &[f64],
    poles: [bool; 2],
    refine_bar: Option<f64>,
    dense: usize,
    probe_fold: bool,
) -> Option<TranslatedStations> {
    if !exact_stripes_on() || crate::blend::open_pad_network_on() || poles != [false, false] || refine_bar.is_some() || dense > 1 || probe_fold {
        return None;
    }
    let radius = radius_at(*station_ts.first()?).to_bits();
    if station_ts.iter().any(|t| radius_at(*t).to_bits() != radius) {
        return None;
    }
    let direction = extrusion_direction(edge, first, second)?;
    let inverse = |surface: &NurbsSurface| -> Option<[Vec3; 2]> {
        let (_, s_u, s_v) = super::open::affine_plane_frame(surface)?;
        let (a, b, c) = (s_u.dot(s_u), s_u.dot(s_v), s_v.dot(s_v));
        let det = a * c - b * b;
        if !(det.abs() > 1e-300) || !det.is_finite() {
            return None;
        }
        Some([s_u.scale(c / det).sub(s_v.scale(b / det)), s_v.scale(a / det).sub(s_u.scale(b / det))])
    };
    Some(TranslatedStations { direction, inverses: [inverse(&first.face.surface)?, inverse(&second.face.surface)?] })
}

pub(in crate::blend) fn march_open_stations_core(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
    radius_at: &dyn Fn(f64) -> f64,
    overshoot_fraction: f64,
    ends: [f64; 2],
    station_count: usize,
    poles: [bool; 2],
    bounded_ends: [bool; 2],
    refine_bar: Option<f64>,
    dense: usize,
    probe_fold: bool,
) -> Result<(Vec<Station>, [usize; 2], StationRefinement, OpenDenseCentres), KernelRefusal> {
    let span = edge.t1 - edge.t0;
    let overshoot = span * overshoot_fraction;
    let [overshoot_before, overshoot_after] =
        poles.map(|pole| if pole { 0.0 } else { overshoot });
    // The overshot span, written `span + 2·overshoot` when neither end is a
    // pole: `(span + o) + o` can round an ulp away from it, and every station
    // parameter and the section step are read off this, so a march with no
    // pole stays the march it was to the bit.
    let marched = if poles == [false, false] {
        span + 2.0 * overshoot
    } else {
        span + overshoot_before + overshoot_after
    };
    for (slot, pole) in poles.into_iter().enumerate() {
        let vertex_t = if slot == 0 { edge.t0 } else { edge.t1 };
        if pole && ends[slot] != vertex_t {
            return Err(KernelRefusal::internal(
                KernelStage::Collect,
                "pole_end_vertex",
                format!(
                    "blend march: edge {}'s pole end is not at its vertex",
                    edge.id
                ),
            ));
        }
    }
    let radius_extent = [0.0, 0.5, 1.0]
        .into_iter()
        .map(|fraction| radius_at(edge.t0 + span * fraction).abs())
        .fold(0.0, f64::max);
    // The overshot span is deliberately NOT used here: `curve_model_scale`
    // clamps out-of-domain probes, so it would measure the same extent while
    // reading as if the overshoot mattered.
    let scale = march_model_scale(&edge.curve, edge.t0, edge.t1, radius_extent)?;
    crate::report_scale_migration("blend_march_open", scale, || {
        let a = edge.curve.evaluate(edge.t0).unwrap_or_default();
        let b = edge.curve.evaluate(edge.t0 + span * 0.5).unwrap_or_default();
        a.length().max(b.length()).max(1.0)
    });
    let surface1 = &first.face.surface;
    let surface2 = &second.face.surface;
    let signs = [first.rho.signum(), second.rho.signum()];
    let rho_at = |t: f64| -> [f64; 2] {
        let radius = radius_at(t);
        [signs[0] * radius, signs[1] * radius]
    };
    // Seed at the middle (the best-conditioned station), then march
    // outward to both ends so every station is seeded by its neighbour.
    let mid_index = station_count / 2;
    // Uniform parameters across the overshot span, with the two marched
    // parameters nearest the BLEND ENDS snapped to exactly t0 / t1: the
    // end stations are analytically knowable (the end radius's tangent
    // offsets on each wall), and making them interpolation points pins
    // the fitted rows THROUGH them — the end trims and any downstream
    // miter boolean read these row endpoints.  A half-step snap keeps the
    // parameter list strictly increasing.
    let (station_ts, vertex_indices) =
        open_station_parameters(edge, overshoot_fraction, ends, station_count, poles);
    let station_t = |index: usize| station_ts[index];
    // The section plane at `t` — the extended point and the edge's unit
    // tangent, the one-sided limit where the parameterization is stationary —
    // advancing past a stationary open end by the march's own chord
    // (`march_section`), so the overshoot stations seat the ball's contacts
    // beyond the vertex instead of on its section plane.
    let station_step = marched.abs() / station_count as f64;
    let section_frame = |t: f64| -> Result<(Vec3, Vec3), KernelRefusal> {
        march_section(edge, t, station_step, "blend march")
    };
    let solve_centered = |t: f64, seed: [f64; 4]| -> Result<([f64; 4], Vec3), KernelRefusal> {
        let (section_point, section_tangent) = section_frame(t)?;
        solve_station_centered(
            surface1,
            surface2,
            rho_at(t),
            seed,
            section_point,
            section_tangent,
            scale,
        )
    };
    let solve_at =
        |t: f64, seed: [f64; 4]| -> Result<[f64; 4], KernelRefusal> { solve_centered(t, seed).map(|(uv, _)| uv) };
    let mid_seed = {
        let t = station_t(mid_index).clamp(edge.t0, edge.t1);
        let uv1 = edge_uv_on_face(first.coedge, edge, t)?;
        let uv2 = edge_uv_on_face(second.coedge, edge, t)?;
        [uv1[0], uv1[1], uv2[0], uv2[1]]
    };
    // A pole's station is the vertex on both mates, read off the edge's own
    // pcurves: the ball there touches both carriers at the one point.
    let pole_index = |index: usize| {
        (poles[0] && index == 0) || (poles[1] && index == station_count)
    };
    let pole_solution = |t: f64| -> Result<[f64; 4], KernelRefusal> {
        let uv1 = edge_uv_on_face(first.coedge, edge, t)?;
        let uv2 = edge_uv_on_face(second.coedge, edge, t)?;
        Ok([uv1[0], uv1[1], uv2[0], uv2[1]])
    };
    // A pole's centre is read off its pcurve station where it can be; the
    // continuation guard skips its centre bar when it cannot (the tangency
    // system is singular there).
    let pole_station = |index: usize| -> Result<([f64; 4], Option<Vec3>), KernelRefusal> {
        let t = station_t(index);
        let uv = pole_solution(t)?;
        let center = section_frame(t).and_then(|(section_point, section_tangent)| {
            tangency_residual(surface1, surface2, rho_at(t), uv, section_point, section_tangent)
                .map(|(_, _, _, center)| center)
        });
        Ok((uv, center.ok()))
    };
    let solve_timing = crate::blend::BlendTimer::start("timing.march.solve_ms");
    let mut solutions = vec![[0.0f64; 4]; station_count + 1];
    let mut centers: Vec<Option<Vec3>> = vec![None; station_count + 1];
    if pole_index(mid_index) {
        let (uv, center) = pole_station(mid_index)?;
        solutions[mid_index] = uv;
        centers[mid_index] = center;
    } else {
        let context = crate::blend::stations::NewtonContext {
            phase: "open-initial",
            stage: "initial",
            index: Some(mid_index),
            t: Some(station_t(mid_index)),
            native: Some([edge.t0, edge.t1]),
            ..Default::default()
        };
        let (uv, center) = crate::blend::stations::with_newton_context(context, || solve_centered(station_t(mid_index), mid_seed))?;
        solutions[mid_index] = uv;
        centers[mid_index] = Some(center);
    }
    // Every other station is the CONTINUATION of the one it is seeded from:
    // an undamped Newton can converge on the far crossing of the section
    // plane, and the guard retries such a hop from a nearer seed
    // (`stations::continue_station`).  The notched crotch arc of a 20°
    // cylinder crossing hopped here at 64 stations (2026-09-26): support 2
    // moved 0.85 of a period and the centre 9.9 for a 2.1 chord.
    let hop = |left: &Continued, right: &Continued| {
        branch_hop([surface1, surface2], edge, &rho_at, left, right)
    };
    let continue_from = |solutions: &mut [[f64; 4]],
                             centers: &mut [Option<Vec3>],
                             index: usize,
                             from: usize|
     -> Result<(), KernelRefusal> {
        if pole_index(index) {
            let (uv, center) = pole_station(index)?;
            solutions[index] = uv;
            centers[index] = center;
            return Ok(());
        }
        let left = Continued { t: station_t(from), uv: solutions[from], center: centers[from] };
        let context = crate::blend::stations::NewtonContext {
            phase: "open-march",
            index: Some(index),
            native: Some([edge.t0, edge.t1]),
            ..Default::default()
        };
        let node = crate::blend::stations::with_newton_context(context, || continue_station(&solve_centered, &hop, &left, station_t(index), 0))?;
        solutions[index] = node.uv;
        centers[index] = node.center;
        Ok(())
    };
    // Near a collapsed carrier boundary the offset surface need not have a
    // continuation past the pole. A free end may already have crossed its
    // end face before that singularity. Keep the solved interval for that
    // case only; the caller must prove both end-face crossings lie inside it.
    let (mut first_kept, mut last_kept) = (0, station_count);
    // Lane R2 row 2b (`exact_stripes_on`): a straight edge between two AFFINE
    // planes at one radius is one section extruded, so every station is the
    // middle one translated along the edge to its section plane — the contacts
    // move in their planes by that shift, the uv by the planes' constant
    // Jacobians — and nothing is solved past the middle. Only where no probe,
    // refinement or dense centre is asked, and never inside the open-pad
    // composition (its landed bodies keep their march).
    let translated = translated_stations(edge, first, second, radius_at, &station_ts, poles, refine_bar, dense, probe_fold);
    if let Some(frame) = translated {
        let (middle_point, _) = section_frame(station_t(mid_index))?;
        for index in (0..=station_count).filter(|index| *index != mid_index) {
            let (point, _) = section_frame(station_t(index))?;
            let shift = frame.direction.scale(point.sub(middle_point).dot(frame.direction));
            let mut uv = solutions[mid_index];
            for (slot, jacobian) in frame.inverses.iter().enumerate() {
                uv[2 * slot] += jacobian[0].dot(shift);
                uv[2 * slot + 1] += jacobian[1].dot(shift);
            }
            solutions[index] = uv;
            centers[index] = centers[mid_index].map(|center| center.add(shift));
        }
        crate::blend::record_count("count.march.translated_stations", station_count as u64);
    }
    for index in (0..mid_index).rev().filter(|_| translated.is_none()) {
        if let Err(error) = continue_from(&mut solutions, &mut centers, index, index + 1) {
            if !bounded_ends[0] {
                return Err(error);
            }
            first_kept = index + 1;
            break;
        }
    }
    for index in (mid_index + 1..=station_count).filter(|_| translated.is_none()) {
        if let Err(error) = continue_from(&mut solutions, &mut centers, index, index - 1) {
            if !bounded_ends[1] {
                return Err(error);
            }
            last_kept = index - 1;
            break;
        }
    }
    solve_timing.stop();
    if last_kept - first_kept < FIT_DEGREE {
        return Err(KernelRefusal::non_convergence(KernelStage::Refine,
            "bounded_march", "blend: too few regular stations before a carrier pole"));
    }
    let station_ts = station_ts[first_kept..=last_kept].to_vec();
    let vertex_indices = vertex_indices.map(|i| i.clamp(first_kept, last_kept) - first_kept);
    let solutions = solutions[first_kept..=last_kept].to_vec();
    let centers = centers[first_kept..=last_kept].to_vec();
    // LOCAL REFINEMENT (2026-09-27). Uniform stations under-sample a wall whose
    // contacts move fast in the edge parameter — a crotch tip, where the ball
    // centre moves ρ/w faster than the edge — and the rows fitted through them
    // bulge off the ball BETWEEN stations while every station is exact (a 15°
    // crossing's rails stood 1.35e-2 off their carriers at 256 uniform
    // stations). Inside the edge's own span, each interval's exact mid-span
    // station is solved and compared with the cubic its four neighbours
    // predict; where they differ by more than `refine_bar` the station is kept,
    // and the two halves are checked next round. Stations only where the wall
    // needs them: every downstream scan samples per knot.
    let mut refinement = StationRefinement::default();
    let refine_timing = refine_bar.map(|_| crate::blend::BlendTimer::start("timing.march.refine_ms"));
    let (station_ts, vertex_indices, solutions, final_centers) = match refine_bar {
        None => (station_ts, vertex_indices, solutions, centers),
        Some(bar) => {
            let contacts = |t: f64, uv: [f64; 4]| -> Result<(Vec3, Vec3), KernelRefusal> {
                let (section_point, section_tangent) = section_frame(t)?;
                let (_, p1, p2, _) = tangency_residual(
                    surface1,
                    surface2,
                    rho_at(t),
                    uv,
                    section_point,
                    section_tangent,
                )?;
                Ok((p1, p2))
            };
            let mut ts = station_ts.clone();
            let mut sols = solutions;
            let mut cens = centers;
            let mut points = ts
                .iter()
                .zip(&sols)
                .map(|(t, uv)| contacts(*t, *uv))
                .collect::<Result<Vec<_>, _>>()?;
            let end_ts = [ts[vertex_indices[0]], ts[vertex_indices[1]]];
            let (low, high) = (end_ts[0].min(end_ts[1]), end_ts[0].max(end_ts[1]));
            let pole_t = |t: f64| (poles[0] && t == ts_first(&station_ts)) || (poles[1] && t == ts_last(&station_ts));
            // At a STATIONARY end (the edge's derivative vanishes at its
            // vertex) the edge parameter does not measure the wall, so a cubic
            // in `t` mis-reads the end interval as under-sampled and keeps
            // halving it into near-coincident stations whose rows cross
            // themselves (zero_tangent_sites' stationary start: a loop that
            // crosses itself at 4 inserted stations). Like a pole, that end's
            // interval is left as marched.
            let stationary_end = end_ts.map(|t| {
                edge.curve
                    .derivatives(t.clamp(edge.t0.min(edge.t1), edge.t0.max(edge.t1)), 1)
                    .map_or(false, |derivatives| derivatives[1].normalized().is_err())
            });
            let at_stationary_end =
                |t: f64| (stationary_end[0] && t == end_ts[0]) || (stationary_end[1] && t == end_ts[1]);
            let ceiling = REFINE_STATION_FACTOR * (station_count + 1);
            let mut dirty = vec![true; ts.len().saturating_sub(1)];
            let mut previous_worst = f64::INFINITY;
            refinement.exit = StationRefineExit::RoundBudget;
            for _round in 0..REFINE_ROUNDS {
                let mut kept: Vec<(usize, Continued, (Vec3, Vec3))> = Vec::new();
                let mut worst = (0.0f64, 0.0f64);
                for interval in 0..ts.len() - 1 {
                    let (a, b) = (ts[interval], ts[interval + 1]);
                    if !dirty[interval] || b <= low || a >= high || pole_t(a) || pole_t(b) || at_stationary_end(a) || at_stationary_end(b) || ts.len() < 4 {
                        continue;
                    }
                    let middle = 0.5 * (a + b);
                    let left = Continued { t: a, uv: sols[interval], center: cens[interval] };
                    // A mid-span station that will not solve is left out: the
                    // refinement may only add stations, never a refusal.
                    let context = crate::blend::stations::NewtonContext {
                        phase: "open-refine",
                        index: Some(interval),
                        native: Some([edge.t0, edge.t1]),
                        ..Default::default()
                    };
                    let Ok(node) = crate::blend::stations::with_newton_context(context, || continue_station(&solve_centered, &hop, &left, middle, 0)) else {
                        continue;
                    };
                    let Ok(exact) = contacts(middle, node.uv) else {
                        continue;
                    };
                    let first = interval.saturating_sub(1).min(ts.len() - 4);
                    let predicted = |pick: fn(&(Vec3, Vec3)) -> Vec3| {
                        lagrange_cubic(&ts[first..first + 4], &points[first..first + 4], pick, middle)
                    };
                    let miss = predicted(|p| p.0)
                        .sub(exact.0)
                        .length()
                        .max(predicted(|p| p.1).sub(exact.1).length());
                    if miss > worst.0 {
                        worst = (miss, middle);
                    }
                    if miss > bar {
                        kept.push((interval, node, exact));
                    }
                }
                refinement.worst_miss = worst.0;
                refinement.worst_t = worst.1;
                if kept.is_empty() {
                    refinement.exit = StationRefineExit::Converged;
                    refinement.worst_miss = 0.0;
                    break;
                }
                if ts.len() + kept.len() > ceiling {
                    refinement.exit = StationRefineExit::StationCeiling;
                    break;
                }
                if !(worst.0 < 0.9 * previous_worst) {
                    refinement.exit = StationRefineExit::NoProgress;
                    break;
                }
                previous_worst = worst.0;
                let mut next_dirty = vec![false; dirty.len() + kept.len()];
                // Insert back to front so earlier indices stay valid.
                for (interval, node, exact) in kept.iter().rev() {
                    let at = interval + 1;
                    ts.insert(at, node.t);
                    sols.insert(at, node.uv);
                    cens.insert(at, node.center);
                    points.insert(at, *exact);
                    refinement.inserted += 1;
                }
                // The halves of every split interval are checked next round.
                let mut shift = 0usize;
                for (interval, _, _) in &kept {
                    let at = interval + shift;
                    next_dirty[at] = true;
                    next_dirty[at + 1] = true;
                    shift += 1;
                }
                dirty = next_dirty;
            }
            let vertex_indices = end_ts.map(|t| ts.iter().position(|x| *x == t).unwrap_or(0));
            (ts, vertex_indices, sols, cens)
        }
    };
    if let Some(refine_timing) = refine_timing {
        refine_timing.stop();
    }
    // CARRIER-BOUNDARY STATIONS (2026-10-17, the face-ring report). A contact
    // that leaves its carrier's patch in a BOUNDED direction runs on past the
    // boundary on the tangent-plane continuation, whose curvature is not the
    // carrier's: a row fitted through stations on both sides of the boundary
    // misses the crossing itself (on the F6 blend of the reported ring the
    // trimmed rail stood 7.1e-5 off its carrier where it met the boundary
    // rail, and no pcurve on the carrier could reach it). The crossing is a
    // station of its own, solved to the march's precision, and the row fit
    // breaks there (`fit_open_rows_breaking_at`) as it breaks at the vertex
    // stations. A periodic direction has no boundary to cross.
    let (station_ts, vertex_indices, solutions, final_centers, boundary_breaks, boundary_inserted, boundary_inserted_count) = {
        let mut station_ts = station_ts;
        let mut vertex_indices = vertex_indices;
        let mut solutions = solutions;
        let mut final_centers = final_centers;
        let mut breaks: [Option<usize>; 4] = [None; 4];
        let mut found = 0usize;
        let mut inserted_indices = [0usize; 20];
        let mut inserted_count = 0usize;
        let limits = |surface: &NurbsSurface| -> Result<[Option<[f64; 2]>; 2], KernelRefusal> {
            let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
            let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
            let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
            Ok([(!closed_u).then_some([u0, u1]), (!closed_v).then_some([v0, v1])])
        };
        let bounds = [limits(surface1)?, limits(surface2)?];
        // Only the open-pad composition's network takes these stations: a
        // plain network keeps the march it always had (its bodies are landed
        // fixtures; `with_open_pad_network`).
        let scoped_on = crate::blend::open_pad_network_on();
        let mut index = 0usize;
        while scoped_on && index + 1 < station_ts.len() && found < breaks.len() {
            let mut crossing: Option<f64> = None;
            'axes: for mate in 0..2 {
                for axis in 0..2 {
                    let Some([low, high]) = bounds[mate][axis] else { continue };
                    let slot = 2 * mate + axis;
                    let (a, b) = (solutions[index][slot], solutions[index + 1][slot]);
                    for limit in [low, high] {
                        if (a - limit) * (b - limit) < 0.0 {
                            // Bisect the march between the two stations.
                            let (mut lo, mut hi) = (station_ts[index], station_ts[index + 1]);
                            let mut lo_value = a - limit;
                            let mut seed = solutions[index];
                            for _ in 0..60 {
                                let mid = 0.5 * (lo + hi);
                                let Ok(uv) = solve_at(mid, seed) else { break };
                                let value = uv[slot] - limit;
                                if value.abs() <= 1e-13 * (1.0 + limit.abs()) {
                                    lo = mid;
                                    hi = mid;
                                    break;
                                }
                                if (value < 0.0) == (lo_value < 0.0) {
                                    lo = mid;
                                    lo_value = value;
                                } else {
                                    hi = mid;
                                }
                                seed = uv;
                            }
                            crossing = Some(0.5 * (lo + hi));
                            break 'axes;
                        }
                    }
                }
            }
            if let Some(t) = crossing {
                let (left, right) = (station_ts[index], station_ts[index + 1]);
                // Within a hair of an existing station: that station is the break.
                let hair = 1e-9 * (right - left).abs().max(1e-300);
                if std::env::var_os("BREP_DEBUG_NETWORK").is_some() {
                    eprintln!(
                        "boundary station: edge {} crossing at t {t:.9} between stations {index}/{} ({left:.9}, {right:.9}), gap fraction {:.3e}",
                        edge.id, index + 1, ((t - left) / (right - left)).min((right - t) / (right - left))
                    );
                }
                if (t - left).abs() <= hair {
                    breaks[found] = Some(index);
                    found += 1;
                } else if (t - right).abs() <= hair {
                    breaks[found] = Some(index + 1);
                    found += 1;
                } else if let Ok((uv, center)) = solve_centered(t, solutions[index]) {
                    // The crossing station, with two plain stations on each
                    // side at a quarter and a half of the old spacing: the
                    // rows are fitted in pieces broken here, and a clamped
                    // cubic's END derivative carries an O(h^3) error that is
                    // the wall's NORMAL error at the break — 5.2e-4 rad on
                    // the reported ring's F6 wall at 64 stations, enough to
                    // fail a later pole test (bar 3e-4 rad) at the trim's
                    // end. Halving the spacing twice takes it 64 times down.
                    let (left_t, right_t) = (station_ts[index], station_ts[index + 1]);
                    let mut inserted: Vec<(f64, [f64; 4], Vec3)> = vec![(t, uv, center)];
                    let mut seed = solutions[index];
                    for fraction in [0.5, 0.75] {
                        let at = left_t + (t - left_t) * fraction;
                        if let Ok((uv, center)) = solve_centered(at, seed) {
                            seed = uv;
                            inserted.push((at, uv, center));
                        }
                    }
                    seed = uv;
                    for fraction in [0.25, 0.5] {
                        let at = t + (right_t - t) * fraction;
                        if let Ok((uv, center)) = solve_centered(at, seed) {
                            seed = uv;
                            inserted.push((at, uv, center));
                        }
                    }
                    inserted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
                    let crossing_offset = inserted.iter().position(|(at, _, _)| *at == t).unwrap_or(0);
                    for (offset, (at, uv, center)) in inserted.iter().enumerate() {
                        station_ts.insert(index + 1 + offset, *at);
                        solutions.insert(index + 1 + offset, *uv);
                        final_centers.insert(index + 1 + offset, Some(*center));
                        if inserted_count < inserted_indices.len() {
                            inserted_indices[inserted_count] = index + 1 + offset;
                            inserted_count += 1;
                        }
                    }
                    let added = inserted.len();
                    for vertex_index in vertex_indices.iter_mut() {
                        if *vertex_index >= index + 1 {
                            *vertex_index += added;
                        }
                    }
                    breaks[found] = Some(index + 1 + crossing_offset);
                    found += 1;
                    index += added;
                }
            }
            index += 1;
        }
        if found > 0 {
            crate::blend::record_count("count.march.boundary_stations", found as u64);
        }
        (station_ts, vertex_indices, solutions, final_centers, breaks, inserted_indices, inserted_count)
    };
    refinement.boundary_breaks = boundary_breaks;
    refinement.boundary_inserted = boundary_inserted;
    refinement.boundary_inserted_count = boundary_inserted_count;
    let station_count = station_ts.len() - 1;
    let station_t = |index: usize| station_ts[index];
    let pole_index = |index: usize| {
        (poles[0] && index == 0) || (poles[1] && index == station_count)
    };
    open_station_trace(edge, &station_ts, &solutions, |t, uv| {
        let (section_point, section_tangent) = section_frame(t)?;
        tangency_residual(surface1, surface2, rho_at(t), uv, section_point, section_tangent)
            .map(|(_, p1, p2, center)| (p1, p2, center, section_point))
    });
    // Does the wall this march would carry FOLD?  Same question the chain and
    // closed-edge marches ask (`blend/fold.rs`), asked here because this is the
    // third march in the kernel and the open-edge and network lanes are the two
    // that reach it.
    //
    // Only the stations INSIDE the edge's own domain are measured.  This march
    // overshoots both ends on extended evaluation, and the centre curve past a
    // carrier's real extent may turn any way at all — refusing on extrapolated
    // geometry would refuse walls that build.  `vertex_indices` are the two
    // stations snapped exactly onto `t0` and `t1`, so they bound the real span.
    //
    // Not asked (`probe_fold` false) by the network's exact extrusion lane:
    // planar mates have no curvature for the ball to exceed.
    if probe_fold {
        let _fold_timing = crate::blend::BlendTimer::start("timing.march.fold_ms");
        let probe = |t: f64,
                     seed: [f64; 4]|
         -> Result<([f64; 4], Vec3, Vec3, Vec3), KernelRefusal> {
            let uv = solve_at(t, seed)?;
            let (section_point, section_tangent) = section_frame(t)?;
            let (_, p1, p2, center) = tangency_residual(
                surface1,
                surface2,
                rho_at(t),
                uv,
                section_point,
                section_tangent,
            )?;
            Ok((uv, p1, p2, center))
        };
        let [low, high] = vertex_indices;
        let (low, high) = (low.min(high), low.max(high));
        // A pole is not probed: the ball has no seat past it to read a
        // curvature from, and its section is a point, which cannot fold.
        let probed: Vec<(f64, [f64; 4])> = (low..=high)
            .filter(|index| !pole_index(*index))
            .map(|index| (station_t(index), solutions[index]))
            .collect();
        let fold_radius = |t: f64| radius_at(t).abs();
        crate::blend::fold::check_wall_fold(
            &fold_radius,
            span,
            &probed,
            &probe,
            edge,
            [first.face, second.face],
            scale,
        )?;
    }
    // The exact centres between the final stations, when a wall verdict asks.
    let dense_timing = (dense > 1).then(|| crate::blend::BlendTimer::start("timing.march.dense_ms"));
    let mut dense_centres = OpenDenseCentres { subdivisions: dense, centres: Vec::new(), station_ts: station_ts.clone() };
    if dense > 1 {
        for interval in 0..station_ts.len().saturating_sub(1) {
            if pole_index(interval) || pole_index(interval + 1) {
                dense_centres.centres.push(vec![None; dense - 1]);
                continue;
            }
            let (a, b) = (station_ts[interval], station_ts[interval + 1]);
            let left = Continued { t: a, uv: solutions[interval], center: final_centers[interval] };
            // CHAINED SEEDING (`BREP_BLEND_DENSE_CHAIN`, 2026-10-17): each
            // centre is continued from the centre before it rather than from
            // the interval's start station, under the same hop guard. The
            // seed is then an interval/32 away instead of up to a whole
            // interval, and the plain Newton converges in two iterations
            // instead of eight: the shelled-box report's F3 solved 18 000
            // dense centres at 40 tangency evaluations each (734 522 in all,
            // 0.76 s of a 1.95 s feature). A centre that fails leaves the
            // next one seeded from the last centre that converged. The
            // centres are the same roots to the solver's own stop
            // (`1e-11 * (1 + scale)`), read by the wall verdict only; the
            // shipped rows are the stations'.
            //
            // And on the same switch the centres are solved by the CHORD road
            // the closed march's wall check took in lane G row 3
            // (`stations::solve_station_state_chord`): the seed is the
            // extrapolation through the last centres of the row, the Jacobian
            // is carried from centre to centre and rebuilt only when a chord
            // step fails to cut the residual by four, and whatever the chord
            // road cannot finish goes to the plain solver from the previous
            // centre's own (u, v) — so it can refuse nothing the plain road
            // would not. The plain road rebuilt a finite-difference Jacobian
            // at every iteration: five evaluations an iteration, thirty a
            // centre even when seeded from the centre before.
            let chain = dense_chain_enabled();
            let mut previous = left;
            let mut row = Vec::with_capacity(dense - 1);
            let jacobian = std::cell::RefCell::new(None::<[[f64; 4]; 4]>);
            let history = std::cell::RefCell::new(vec![(left.t, left.uv)]);
            let solve_chord = |t: f64, seed: [f64; 4]| -> Result<([f64; 4], Vec3), KernelRefusal> {
                let (section_point, section_tangent) = section_frame(t)?;
                let predicted = predict_station_uv(&history.borrow(), t).unwrap_or(seed);
                let mut carried = jacobian.borrow_mut();
                crate::blend::stations::solve_station_state_chord(
                    surface1, surface2, rho_at(t), predicted, seed, section_point, section_tangent, scale, &mut carried,
                )
                .map(|(uv, state)| (uv, state.center))
            };
            for step in 1..dense {
                let t = a + (b - a) * step as f64 / dense as f64;
                let context = crate::blend::stations::NewtonContext {
                    phase: "open-dense",
                    index: Some(interval),
                    native: Some([edge.t0, edge.t1]),
                    ..Default::default()
                };
                let node = if chain {
                    crate::blend::stations::with_newton_context(context, || continue_station(&solve_chord, &hop, &previous, t, 0)).ok()
                } else {
                    crate::blend::stations::with_newton_context(context, || continue_station(&solve_centered, &hop, &left, t, 0)).ok()
                };
                if chain {
                    if let Some(node) = node {
                        previous = node;
                        history.borrow_mut().push((node.t, node.uv));
                    }
                }
                row.push(node.and_then(|node| node.center));
            }
            dense_centres.centres.push(row);
        }
    }
    if let Some(dense_timing) = dense_timing {
        dense_timing.stop();
    }
    let _assemble_timing = crate::blend::BlendTimer::start("timing.march.assemble_ms");
    let mut stations = Vec::with_capacity(station_count + 1);
    for (index, uv) in solutions.iter().enumerate() {
        let t = station_t(index);
        let (section_point, section_tangent) = section_frame(t)?;
        let (_, p1, p2, center) = tangency_residual(
            surface1,
            surface2,
            rho_at(t),
            *uv,
            section_point,
            section_tangent,
        )?;
        let n1 = raw_normal(surface1, uv[0], uv[1])?;
        let n2 = raw_normal(surface2, uv[2], uv[3])?;
        let cos_alpha = signs[0] * signs[1] * n1.dot(n2);
        let weight = ((1.0 + cos_alpha) * 0.5).max(0.0).sqrt();
        if weight <= 1e-6 {
            return Err(KernelRefusal::unsupported(
                KernelStage::Refine,
                "tangent_station",
                "blend: faces are tangent at a station (α = π)",
            ));
        }
        // At a pole the section is a point, so its tangent lines have no
        // meeting point of their own: the section's middle control point is
        // the pole too (a zero-sweep arc).
        let apex = if pole_index(index) {
            p1
        } else {
            apex_point(p1, n1, p2, n2, center)?
        };
        stations.push(Station {
            uv1: [uv[0], uv[1]],
            uv2: [uv[2], uv[3]],
            p1,
            p2,
            center,
            weight,
            apex,
        });
    }
    Ok((stations, vertex_indices, refinement, dense_centres))
}

fn ts_first(ts: &[f64]) -> f64 {
    ts[0]
}

fn ts_last(ts: &[f64]) -> f64 {
    ts[ts.len() - 1]
}

/// The cubic through four `(t, point)` samples, evaluated at `t`: the
/// prediction a mid-span station is checked against.
fn lagrange_cubic(
    ts: &[f64],
    samples: &[(Vec3, Vec3)],
    pick: fn(&(Vec3, Vec3)) -> Vec3,
    t: f64,
) -> Vec3 {
    let mut sum = Vec3::default();
    for i in 0..4 {
        let mut weight = 1.0;
        for j in 0..4 {
            if i != j {
                weight *= (t - ts[j]) / (ts[i] - ts[j]);
            }
        }
        sum = sum.add(pick(&samples[i]).scale(weight));
    }
    sum
}

/// The section parameters of an open march, as [`march_open_stations_ending`]
/// lays them out: uniform across the overshot span (no overshoot past a pole
/// end), with the two parameters nearest `ends` snapped exactly onto them.
/// Returns the parameters and the indices of the two snapped stations.  A
/// test reading a march's stations against their own section points uses
/// this rather than re-deriving the layout.
pub(in crate::blend) fn open_station_parameters(
    edge: &EdgeRecord,
    overshoot_fraction: f64,
    ends: [f64; 2],
    station_count: usize,
    poles: [bool; 2],
) -> (Vec<f64>, [usize; 2]) {
    let span = edge.t1 - edge.t0;
    let overshoot = span * overshoot_fraction;
    let [overshoot_before, overshoot_after] =
        poles.map(|pole| if pole { 0.0 } else { overshoot });
    // See `march_open_stations_ending` for why `span + 2·overshoot` is
    // written as one product.
    let marched = if poles == [false, false] {
        span + 2.0 * overshoot
    } else {
        span + overshoot_before + overshoot_after
    };
    // An end past the overshoot — a flush join whose seam ball sits before
    // the edge's own start, as a trim curve's does when the pad it continues
    // crossed the wall's floor rail before the trim began (the 2026-10-17
    // face-ring report) — would be snapped onto the first or last station
    // and refused as a join station outside the rows. Such an end takes the
    // march's overshoot past itself; an end inside the overshoot leaves every
    // parameter exactly as it was.
    // The open-pad composition's network only (`with_open_pad_network`): a
    // plain network snaps the nearest station onto an end outside the
    // overshoot as it always did (its bodies are landed fixtures).
    let extend = crate::blend::open_pad_network_on();
    let start = if extend && ends[0] < edge.t0 - overshoot_before { ends[0] - overshoot_before } else { edge.t0 - overshoot_before };
    let finish = if extend && ends[1] > edge.t1 + overshoot_after { ends[1] + overshoot_after } else { edge.t1 + overshoot_after };
    let mut ts: Vec<f64> = if start == edge.t0 - overshoot_before && finish == edge.t1 + overshoot_after {
        (0..=station_count)
            .map(|index| edge.t0 - overshoot_before + marched * index as f64 / station_count as f64)
            .collect()
    } else {
        (0..=station_count)
            .map(|index| start + (finish - start) * index as f64 / station_count as f64)
            .collect()
    };
    let mut snapped = [0usize; 2];
    for (slot, target) in ends.into_iter().enumerate() {
        let mut nearest = 0usize;
        for (index, t) in ts.iter().enumerate() {
            if (t - target).abs() < (ts[nearest] - target).abs() {
                nearest = index;
            }
        }
        ts[nearest] = target;
        snapped[slot] = nearest;
    }
    (ts, snapped)
}

/// `BREP_BLEND_STATION_TRACE=1` prints every solved open-march station, the
/// same way `chain/march.rs::station_trace` prints a chain's: the two contacts,
/// the ball centre, both parameter feet and the centre's step from the previous
/// station.  The fitted rows cannot tell a station from the interpolant through
/// it, so a branch hop between two stations is only visible here.
fn open_station_trace(
    edge: &EdgeRecord,
    station_ts: &[f64],
    solutions: &[[f64; 4]],
    evaluate: impl Fn(f64, [f64; 4]) -> Result<(Vec3, Vec3, Vec3, Vec3), KernelRefusal>,
) {
    if std::env::var("BREP_BLEND_STATION_TRACE").ok().as_deref() != Some("1") {
        return;
    }
    eprintln!(
        "blend open stations: edge {} {} station(s) over t [{:.9}, {:.9}]",
        edge.id,
        solutions.len(),
        station_ts.first().copied().unwrap_or(f64::NAN),
        station_ts.last().copied().unwrap_or(f64::NAN),
    );
    let mut previous: Option<Vec3> = None;
    for (index, uv) in solutions.iter().enumerate() {
        let t = station_ts[index];
        let Ok((p1, p2, center, section)) = evaluate(t, *uv) else {
            eprintln!("  pos {index:>3} t {t:.9} (station could not be evaluated)");
            continue;
        };
        let step = previous
            .map(|point| center.sub(point).length())
            .unwrap_or(f64::NAN);
        eprintln!(
            "  pos {:>3} t {:.9} p1 ({:.9}, {:.9}, {:.9}) p2 ({:.9}, {:.9}, {:.9}) \
             c ({:.9}, {:.9}, {:.9}) s ({:.9}, {:.9}, {:.9}) uv1 ({:.9}, {:.9}) uv2 ({:.9}, {:.9}) \
             dc {:.9}",
            index, t, p1.x, p1.y, p1.z, p2.x, p2.y, p2.z, center.x, center.y, center.z,
            section.x, section.y, section.z, uv[0], uv[1], uv[2], uv[3], step,
        );
        previous = Some(center);
    }
}

/// Whether the stripe along `edge` runs out into a POLE at its `at_start` end:
/// its two mates are tangent at that vertex from the side the ball sits on, so
/// the ball seated there touches both carriers AT the vertex and the blend's
/// section shrinks to that one point.
///
/// This is the end of an edge whose dihedral closes to π along it — the seam
/// between two rounds that meet where both run tangent into a third face, the
/// crease of two tangent-meeting blends.  Read as the two offset carriers
/// meeting at the vertex: `p + ρ·n` on each mate, at the vertex, is one point
/// within `bar`.  Faces tangent from OPPOSITE sides (a knife edge) put the two
/// offsets `2ρ` apart and are not a pole.
pub(in crate::blend) fn pole_end(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
    at_start: bool,
    bar: f64,
) -> Result<bool, KernelRefusal> {
    let t = if at_start { edge.t0 } else { edge.t1 };
    let uv1 = edge_uv_on_face(first.coedge, edge, t)?;
    let uv2 = edge_uv_on_face(second.coedge, edge, t)?;
    let one = blend_offset(&first.face.surface).at(uv1[0], uv1[1], first.rho).or_refuse(KernelStage::Refine, "at")?;
    let two = blend_offset(&second.face.surface).at(uv2[0], uv2[1], second.rho).or_refuse(KernelStage::Refine, "at")?;
    Ok(one.point.sub(two.point).length() <= bar)
}

/// Clamped open interpolation of the blend rows (no wrap, no seam).
/// The row parameters of the two stations the march snapped onto the edge's
/// vertices, start then end.  `vertex_indices` comes straight from
/// [`march_open_stations`], which placed one station EXACTLY at each vertex
/// parameter; the section there is the only one that lies in the vertex's
/// normal plane, and it is what the corner/flush/cap surgery trims to.
///
/// This used to be re-derived geometrically as "the station whose ball centre
/// is nearest the vertex", on the argument that the centre path is an offset
/// of the edge and so turns its closest approach at the vertex.  That is only
/// true when the section's offset direction is constant along the edge (a
/// straight edge, or planar mates).  On a CURVED support the ball centre
/// swings as it rolls, the distance-to-vertex minimum slides off the vertex
/// station, and a NEIGHBOUR wins by a hair: on the reported box+through-hole
/// document the elliptical rim's two arcs picked stations 0.0018 apart in
/// distance but 9.0e-2 apart in space, and the two halves of one smooth rim
/// then refused each other as sitting on "different balls".
pub(in crate::blend) fn vertex_stations(
    parameters: &[f64],
    vertex_indices: [usize; 2],
) -> Option<[f64; 2]> {
    if parameters.len() < 2 || vertex_indices[0] == vertex_indices[1] {
        return None;
    }
    Some([
        *parameters.get(vertex_indices[0])?,
        *parameters.get(vertex_indices[1])?,
    ])
}

/// Fit the rows of an open march.
///
/// `vertex_indices` are the two stations [`march_open_stations`] snapped onto
/// the edge's vertices, when the caller has them.  The rows are interpolated
/// in independent pieces between them and rejoined C0 at those two stations.
///
/// # Why the fit breaks at the vertex stations
///
/// The march overshoots both ends on `derivatives_extended` and
/// `evaluate_extended`, which continue a curve along its end TANGENT and a
/// surface along its boundary TANGENT PLANE.  So on a curved carrier the
/// stations are a circle (say) up to the vertex station and a straight line
/// past it: the data is G1 there and its curvature jumps.  One cubic
/// interpolant through both rings across that jump, and the ringing lands
/// inside the edge's own span, where the corner closures read the rows.
/// Measured on the 2026-09-15 report (a sketch arc notch at r = 0.5, 64
/// stations): the notch stripe's centre row stood 1.9e-4 off tangency one
/// station inside its vertex against 2e-7 mid-span, a miter seam read against
/// it left the strip 1.4e-5 past the sharp edge, the geometrically symmetric
/// corner was built as an asymmetric one with a 1.4e-5 connector, and the
/// network's result was refused for that connector's fit.  Fitted in pieces
/// the same row reads 6.1e-7 at its worst and the seam exits on the sharp edge
/// to 8e-15.
///
/// The pieces need `FIT_DEGREE + 1` samples each; a march too short for that
/// keeps the single fit.
///
/// `extrusion` is [`extrusion_direction`] of the edge and its two mates, when
/// the caller has them: the direction both carriers are invariant along, which
/// is what makes every section of a constant-radius blend the same section.
/// How far a fitted rail stands off the carrier it is tangent to, read
/// BETWEEN the stations it was fitted through.
///
/// The rows are a cubic through the march's stations, so at a station the rail
/// is on its carrier by construction and the reading has to be taken between
/// them — at 0.25, 0.5 and 0.75 of every span, the fit record's between-station
/// sampling, the same places [`crate::blend::is_marched_fit_off_carriers`] reads a
/// marched seam.
///
/// This is the term that decides a filleted solid's volume, and it is not the
/// one any gate was watching.  Measured on the notched cap 2026-09-16, the
/// lobe's 305-degree coaxial stripe against its cylinder, by THIS function (3
/// points per span, in range): 1.419e-6 at 64 stations, where the solid reads
/// -1.618e-5 against its closed form; 3.534e-7 at 128, -1.723e-6; 3.019e-8 at
/// 256, -1.127e-6.  The five straight stripes in the same solid read 7e-15 at
/// every count, because the exact-extrusion lane makes them exact rather than
/// fitted.
///
/// `exact_extrusion_probe closure` reads the same rail as 1.360e-6 / 1.329e-7 /
/// 6.29e-8, because it samples the TRIMMED edge at 129 evenly spaced points
/// rather than three per station span.  Same rail, two instruments; the
/// ladder's numbers are the ones its bar is compared with.
///
/// The reading keeps falling as the march refines, so it has no floor worth
/// naming at these counts; what caps the climb is cost, not the reading
/// ([`crate::blend::stations::MAX_OPEN_STATIONS`]).
pub(in crate::blend) fn rails_off_carriers(
    rows: &FittedRows,
    parameters: &[f64],
    carriers: [&NurbsSurface; 2],
) -> Result<f64, KernelRefusal> {
    // Only BETWEEN the vertex stations.  The march deliberately overshoots
    // past both vertices onto the carriers' tangent continuation, so out there
    // the rail is off the carrier by construction and by a lot -- the notched
    // cap's seven stripes read 0.14 to 1.4 with the overshoot included, on
    // stripes whose in-range rails are 7e-15.  That span is not the blend.
    let [low, high] = match rows.vertex_stations {
        Some([a, b]) if a < b => [a, b],
        Some([a, b]) => [b, a],
        None => [
            *parameters.first().unwrap_or(&0.0),
            *parameters.last().unwrap_or(&1.0),
        ],
    };
    let mut worst = 0.0f64;
    for (rail, surface) in [(&rows.cr, carriers[0]), (&rows.cs, carriers[1])] {
        for pair in parameters.windows(2) {
            if pair[1] <= low || pair[0] >= high {
                continue;
            }
            for fraction in [0.25, 0.5, 0.75] {
                let u = pair[0] + (pair[1] - pair[0]) * fraction;
                if u < low || u > high {
                    continue;
                }
                let point = rail.evaluate(u).or_refuse(KernelStage::Refine, "evaluate")?;
                let off = crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?.distance;
                if !worst.is_nan() && (off.is_nan() || off > worst) {
                    worst = off;
                }
            }
        }
    }
    Ok(worst)
}

pub(in crate::blend) fn fit_open_rows(
    stations: &[Station],
    parameters: &[f64],
    chamfer: bool,
    vertex_indices: Option<[usize; 2]>,
    extrusion: Option<Vec3>,
) -> Result<FittedRows, KernelRefusal> {
    fit_open_rows_breaking_at(stations, parameters, chamfer, vertex_indices, extrusion, &[])
}

/// [`fit_open_rows`] with the rows and pcurves also broken at `extra`
/// station indices — the carrier-boundary stations a march inserted
/// (`StationRefinement::boundary_breaks`).
pub(in crate::blend) fn fit_open_rows_breaking_at(
    stations: &[Station],
    parameters: &[f64],
    chamfer: bool,
    vertex_indices: Option<[usize; 2]>,
    extrusion: Option<Vec3>,
    extra: &[usize],
) -> Result<FittedRows, KernelRefusal> {
    let mut samples_cr = Vec::new();
    let mut samples_mid = Vec::new();
    let mut samples_cs = Vec::new();
    let mut samples_uv1 = Vec::new();
    let mut samples_uv2 = Vec::new();
    let mut samples_center = Vec::new();
    for station in stations {
        samples_cr.push(Vec4::from_point(station.p1, 1.0));
        samples_cs.push(Vec4::from_point(station.p2, 1.0));
        samples_center.push(Vec4::from_point(station.center, 1.0));
        samples_mid.push(Vec4 {
            x: station.apex.x * station.weight,
            y: station.apex.y * station.weight,
            z: station.apex.z * station.weight,
            w: station.weight,
        });
        samples_uv1.push(Vec4::from_point(
            Vec3::new(station.uv1[0], station.uv1[1], 0.0),
            1.0,
        ));
        samples_uv2.push(Vec4::from_point(
            Vec3::new(station.uv2[0], station.uv2[1], 0.0),
            1.0,
        ));
    }
    let cr_pcurve = fit::interpolate_homogeneous(&samples_uv1, FIT_DEGREE, parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    let cs_pcurve = fit::interpolate_homogeneous(&samples_uv2, FIT_DEGREE, parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    // Rigidly translated sections admit exact extrusion rows, preserving an
    // analytic carrier instead of approximating it with a cubic fitted surface.
    if let Some(rows) =
        exact_extrusion_rows(stations, parameters, chamfer, extrusion, cr_pcurve.clone(), cs_pcurve.clone())?
    {
        return Ok(rows);
    }
    let breaks = row_breaks(stations.len(), vertex_indices, extra);
    let interpolate = |samples: &[Vec4]| interpolate_in_pieces(samples, parameters, &breaks);
    let (cr_pcurve, cs_pcurve) = if breaks.is_empty() {
        (cr_pcurve, cs_pcurve)
    } else {
        (interpolate(&samples_uv1)?, interpolate(&samples_uv2)?)
    };
    let cr = interpolate(&samples_cr)?;
    let cs = interpolate(&samples_cs)?;
    let mid = if chamfer {
        None
    } else {
        Some(interpolate(&samples_mid)?)
    };
    let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
    let surface = crate::blend::rows::surface_from_rows(FIT_DEGREE, &cr, &cs, mid.as_ref(), false)?;
    let center = interpolate(&samples_center)?;
    Ok(FittedRows {
        surface,
        cr,
        cs,
        cr_pcurve,
        cs_pcurve,
        u_domain,
        center: Some(center),
        exact_extrusion: false,
        vertex_stations: None,
        revolution_axis: None,
    })
}

/// Fit the rows of a march whose stations are broken at EVERY index in
/// `breaks` — the runout wall (`runout.rs`), whose carriers switch at the
/// tri-tangent stations: the contact rail is G1 there by the mirror identity
/// but its curvature jumps from a plane's to a pad cylinder's, the same jump
/// [`fit_open_rows`] breaks at the vertex stations for.  No exact-extrusion
/// shortcut (the carriers change along the wall) and no mate pcurves: a
/// station at a switch has one uv on EACH carrier, so the mate pcurves are
/// fitted per stretch by the caller ([`interpolate_piece`]).  Every piece
/// between breaks must carry `FIT_DEGREE` samples or the fit refuses.
pub(in crate::blend) fn fit_rows_breaking(
    stations: &[Station],
    parameters: &[f64],
    breaks: &[usize],
) -> Result<FittedRows, KernelRefusal> {
    let mut bounds = vec![0usize];
    bounds.extend(breaks.iter().copied());
    bounds.push(stations.len().saturating_sub(1));
    if bounds.windows(2).any(|pair| pair[1] < pair[0] + FIT_DEGREE) {
        return Err(KernelRefusal::internal(
            KernelStage::Refine,
            "row_piece_short",
            format!("blend: a row piece between breaks {breaks:?} is too short for the fit degree ({} stations)", stations.len()),
        ));
    }
    let mut samples_cr = Vec::new();
    let mut samples_mid = Vec::new();
    let mut samples_cs = Vec::new();
    let mut samples_center = Vec::new();
    for station in stations {
        samples_cr.push(Vec4::from_point(station.p1, 1.0));
        samples_cs.push(Vec4::from_point(station.p2, 1.0));
        samples_center.push(Vec4::from_point(station.center, 1.0));
        samples_mid.push(Vec4 {
            x: station.apex.x * station.weight,
            y: station.apex.y * station.weight,
            z: station.apex.z * station.weight,
            w: station.weight,
        });
    }
    let interpolate = |samples: &[Vec4]| interpolate_in_pieces(samples, parameters, breaks);
    let cr = interpolate(&samples_cr)?;
    let cs = interpolate(&samples_cs)?;
    let mid = interpolate(&samples_mid)?;
    let center = interpolate(&samples_center)?;
    let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
    let surface = crate::blend::rows::surface_from_rows(FIT_DEGREE, &cr, &cs, Some(&mid), false)?;
    // The mate pcurves are the caller's (one per carrier stretch); these two
    // only keep the record well formed and are never read.
    let placeholder = fit::interpolate_homogeneous(
        &[Vec4::from_point(Vec3::default(), 1.0), Vec4::from_point(Vec3::new(1.0, 0.0, 0.0), 1.0)],
        1,
        &[0.0, 1.0],
    ).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    Ok(FittedRows {
        surface,
        cr,
        cs,
        cr_pcurve: placeholder.clone(),
        cs_pcurve: placeholder,
        u_domain,
        center: Some(center),
        exact_extrusion: false,
        vertex_stations: None,
        revolution_axis: None,
    })
}

/// A cubic (or cubics broken at the LOCAL sample indices in `breaks`) through
/// `samples` at `parameters` — any increasing parameters, not only `[0, 1]`
/// — parameterised over exactly `[parameters[0], parameters[last]]`, the way
/// the pieces of [`interpolate_in_pieces`] are: a mate pcurve for a rail edge
/// that spans a sub-range of its row.
pub(in crate::blend) fn interpolate_piece(
    samples: &[Vec4],
    parameters: &[f64],
    breaks: &[usize],
) -> Result<NurbsCurve, KernelRefusal> {
    if samples.len() < 2 || samples.len() != parameters.len() {
        return Err(KernelRefusal::internal(KernelStage::Refine, "interpolate_piece", "blend: a pcurve piece needs matching samples and parameters"));
    }
    // With breaks inside the piece the row fit is several cubics joined C0,
    // and the pcurve has to be the SAME cubics (on a plane it is their affine
    // image; one cubic over the whole piece stood 7.5e-4 off the rail between
    // stations): `interpolate_in_pieces` maps every piece onto its own
    // parameter range already.
    if !breaks.is_empty() {
        return interpolate_in_pieces(samples, parameters, breaks);
    }
    let (low, high) = (parameters[0], parameters[parameters.len() - 1]);
    if !(high > low) {
        return Err(KernelRefusal::internal(KernelStage::Refine, "interpolate_piece", "blend: a pcurve piece has no parameter span"));
    }
    let local: Vec<f64> = parameters.iter().map(|parameter| (parameter - low) / (high - low)).collect();
    let fitted = fit::interpolate_homogeneous(samples, FIT_DEGREE, &local).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    let knots: Vec<f64> = fitted.knots.iter().map(|knot| low + knot * (high - low)).collect();
    NurbsCurve::new(fitted.degree, knots, fitted.control_points).or_refuse(KernelStage::Refine, "new")
}

/// The interior sample indices the rows break at — the vertex stations that
/// are not already an end of the march — or none when a piece between them
/// would be too short to carry the fit degree.
fn row_breaks(count: usize, vertex_indices: Option<[usize; 2]>, extra: &[usize]) -> Vec<usize> {
    let mut breaks: Vec<usize> = vertex_indices
        .into_iter()
        .flatten()
        .chain(extra.iter().copied())
        .filter(|&index| index > 0 && index + 1 < count)
        .collect();
    if breaks.is_empty() {
        return Vec::new();
    }
    breaks.sort_unstable();
    breaks.dedup();
    let mut bounds = Vec::with_capacity(breaks.len() + 2);
    bounds.push(0);
    bounds.extend(breaks.iter().copied());
    bounds.push(count.saturating_sub(1));
    if bounds.windows(2).all(|pair| pair[1] - pair[0] >= FIT_DEGREE) {
        breaks
    } else {
        Vec::new()
    }
}

/// Interpolate `samples` at `parameters` as one cubic per piece between the
/// sample indices in `breaks`, each clamped over its own parameter interval,
/// joined into one curve with a knot of multiplicity `FIT_DEGREE` at every
/// break.  Both neighbours interpolate the break sample, so the shared control
/// point is that sample and the curve passes through it.  No breaks is the
/// single interpolant.
fn interpolate_in_pieces(
    samples: &[Vec4],
    parameters: &[f64],
    breaks: &[usize],
) -> Result<NurbsCurve, KernelRefusal> {
    if breaks.is_empty() {
        return fit::interpolate_homogeneous(samples, FIT_DEGREE, parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous");
    }
    let mut bounds = Vec::with_capacity(breaks.len() + 2);
    bounds.push(0);
    bounds.extend(breaks.iter().copied());
    bounds.push(samples.len() - 1);
    let pieces = bounds.len() - 1;
    let mut knots: Vec<f64> = Vec::new();
    let mut control: Vec<Vec4> = Vec::new();
    for (piece, pair) in bounds.windows(2).enumerate() {
        let (from, to) = (pair[0], pair[1]);
        let (low, high) = (parameters[from], parameters[to]);
        // `interpolate_homogeneous` clamps over [0, 1], so each piece is
        // fitted on its own normalized parameters and its knots mapped back.
        let local: Vec<f64> = parameters[from..=to]
            .iter()
            .map(|parameter| (parameter - low) / (high - low))
            .collect();
        let fitted = fit::interpolate_homogeneous(&samples[from..=to], FIT_DEGREE, &local).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
        if fitted.degree != FIT_DEGREE {
            return Err(KernelRefusal::internal(
                KernelStage::Refine,
                "row_piece_short",
                "blend: a row piece is too short for the fit degree",
            ));
        }
        if piece == 0 {
            knots.extend(std::iter::repeat(low).take(FIT_DEGREE + 1));
        }
        let interior = &fitted.knots[FIT_DEGREE + 1..fitted.knots.len() - FIT_DEGREE - 1];
        knots.extend(interior.iter().map(|knot| low + knot * (high - low)));
        let multiplicity = if piece + 1 == pieces { FIT_DEGREE + 1 } else { FIT_DEGREE };
        knots.extend(std::iter::repeat(high).take(multiplicity));
        control.extend(fitted.control_points.iter().skip(usize::from(piece > 0)).copied());
    }
    NurbsCurve::new(FIT_DEGREE, knots, control).or_refuse(KernelStage::Refine, "new")
}

/// The direction a blend along `edge` between these two mates is an
/// EXTRUSION along, or `None` when the geometry does not say it is one.
///
/// It is one when the edge is a straight segment and both carriers are
/// invariant under translation along it: a plane containing it, or a cylinder
/// — a `RuledRevolution` whose two radii agree, or a general `Revolution`
/// (a partial sweep, which is what an extruded sketch arc recognizes as)
/// whose generatrix is a straight segment parallel to the axis — with its
/// axis parallel to the edge.  Every comparison is a lateral distance against
/// the kernel's `model` tolerance, the one `straight_segment` is asked with.
///
/// The question is geometric, never the representation's: an edge rejoined
/// with a redundant pole is still straight, and a quarter cylinder recognized
/// as a general revolution is still a cylinder.
pub(in crate::blend) fn extrusion_direction(
    edge: &EdgeRecord,
    first: &BlendMate,
    second: &BlendMate,
) -> Option<Vec3> {
    let tolerance = crate::KernelTolerances::for_scale(
        march_model_scale(&edge.curve, edge.t0, edge.t1, 0.0).ok()?,
        1e-7,
    )
    .model;
    let (start, end) = edge.curve.straight_segment(tolerance)?;
    let chord = end.sub(start);
    let direction = chord.normalized().ok()?;
    // How far `vector` leans off `axis` (a unit vector) over its own length.
    let lean = |vector: Vec3, axis: Vec3| vector.sub(axis.scale(vector.dot(axis))).length();
    let invariant = |mate: &BlendMate| -> Option<()> {
        match mate.face.surface.analytic()? {
            crate::AnalyticSurface::Plane { u_dir, v_dir, .. } => {
                let normal = u_dir.cross(*v_dir).normalized().ok()?;
                (chord.dot(normal).abs() <= tolerance).then_some(())
            }
            crate::AnalyticSurface::RuledRevolution {
                frame, rho0, rho1, ..
            } => ((rho1 - rho0).abs() <= tolerance && lean(chord, frame.axis) <= tolerance)
                .then_some(()),
            crate::AnalyticSurface::Revolution {
                frame, generatrix, ..
            } => {
                let (from, to) = generatrix.straight_segment(tolerance)?;
                (lean(to.sub(from), frame.axis) <= tolerance && lean(chord, frame.axis) <= tolerance)
                    .then_some(())
            }
            _ => None,
        }
    };
    invariant(first)?;
    invariant(second)?;
    Some(direction)
}

/// The exact rows of a translational march, or `None` when the stations are
/// not rigid translates of the first along one direction.
///
/// Exactness bar: every station's contact points, apex and centre must sit
/// within `1e-9·(1 + scale)` of the first station's translated by the
/// centre's own displacement, and the centre displacements must be collinear
/// to the same bar.  The station solver converges to `1e-11·(1 + scale)`, so
/// a genuinely straight×planar edge passes with two orders of margin, and a
/// curved carrier fails by its sagitta at the second station.
///
/// # The section weight, and the geometric door
///
/// The weight is compared too, to an absolute `1e-12`.  It is not an
/// independent quantity — it is the cosine of the section's half angle, set by
/// the contacts and the centre already held to the bar above — and on a curved
/// carrier the solve does not land it that close: the D-prism's plane ×
/// cylinder edge (`open_heal_tests::d_prism`) reads every position within
/// 1.2e-10 of its translate against a 1.2e-8 bar and the weight 6.0e-12 off, so
/// a straight edge between a plane and a cylinder took the 65-row cubic.
///
/// When the caller's `extrusion` says the carriers ARE an extrusion along the
/// edge, every section is the same section by construction, so the weight is
/// not asked of the stations: the rows are the FIRST section and that section
/// translated along the travel, and the stations remain the witness that the
/// march is one — the positions to the same bar, and the travel along the
/// extrusion.  A varying radius law fails that witness at its second station.
/// Without `extrusion`, and wherever the weights do agree, the rows are the
/// first and last stations as they always were.
fn exact_extrusion_rows(
    stations: &[Station],
    parameters: &[f64],
    chamfer: bool,
    extrusion: Option<Vec3>,
    cr_pcurve: NurbsCurve,
    cs_pcurve: NurbsCurve,
) -> Result<Option<FittedRows>, KernelRefusal> {
    let (Some(first), Some(last)) = (stations.first(), stations.last()) else {
        return Ok(None);
    };
    if stations.len() < 3 {
        return Ok(None);
    }
    let travel = last.center.sub(first.center);
    let length = travel.length();
    if !(length > 1e-12) {
        return Ok(None);
    }
    let direction = travel.scale(1.0 / length);
    let scale = stations
        .iter()
        .fold(0.0f64, |worst, station| worst.max(station.center.length()))
        .max(1.0);
    let bar = 1e-9 * (1.0 + scale);
    let mut weights_agree = true;
    for (station, parameter) in stations.iter().zip(parameters) {
        let shift = station.center.sub(first.center);
        // Collinear with the travel, at the row parameter's fraction of it.
        if shift.sub(direction.scale(shift.dot(direction))).length() > bar
            || (shift.dot(direction) - parameter * length).abs() > bar
        {
            return Ok(None);
        }
        let translate = |point: Vec3| point.add(shift);
        if station.p1.sub(translate(first.p1)).length() > bar
            || station.p2.sub(translate(first.p2)).length() > bar
            || station.apex.sub(translate(first.apex)).length() > bar
        {
            return Ok(None);
        }
        if (station.weight - first.weight).abs() > 1e-12 {
            weights_agree = false;
        }
    }
    // The last row: the last station, or the first section carried along the
    // travel when only the geometry vouches for the section.
    let (end, offset) = if weights_agree {
        (last, Vec3::default())
    } else {
        let Some(along) = extrusion else {
            return Ok(None);
        };
        if travel.sub(along.scale(travel.dot(along))).length() > bar {
            return Ok(None);
        }
        (first, travel)
    };
    // Section along v (z), spine along u: the blend's own convention.
    let section = |z_index: usize, station: &Station, offset: Vec3| -> Vec4 {
        match z_index {
            0 => Vec4::from_point(station.p1.add(offset), 1.0),
            1 if !chamfer => {
                let apex = station.apex.add(offset);
                Vec4 {
                    x: apex.x * station.weight,
                    y: apex.y * station.weight,
                    z: apex.z * station.weight,
                    w: station.weight,
                }
            }
            _ => Vec4::from_point(station.p2.add(offset), 1.0),
        }
    };
    let columns = if chamfer { 2 } else { 3 };
    let control = vec![
        (0..columns).map(|j| section(j, first, Vec3::default())).collect::<Vec<_>>(),
        (0..columns).map(|j| section(j, end, offset)).collect::<Vec<_>>(),
    ];
    let (degree_v, knots_v) = crate::blend::rows::cross_section_basis(chamfer);
    let surface = NurbsSurface::new(1, degree_v, vec![0.0, 0.0, 1.0, 1.0], knots_v, control).or_refuse(KernelStage::Refine, "new")?;
    let cr = crate::make_line(first.p1, end.p1.add(offset)).or_refuse(KernelStage::Refine, "make_line")?;
    let cs = crate::make_line(first.p2, end.p2.add(offset)).or_refuse(KernelStage::Refine, "make_line")?;
    let center = crate::make_line(first.center, end.center.add(offset)).or_refuse(KernelStage::Refine, "make_line")?;
    Ok(Some(FittedRows {
        surface,
        cr,
        cs,
        cr_pcurve,
        cs_pcurve,
        u_domain: [0.0, 1.0],
        center: Some(center),
        exact_extrusion: true,
        vertex_stations: None,
        revolution_axis: None,
    }))
}

/// The exact REVOLUTION lane's pre-march decision (2026-10-13): the axis
/// both mates revolve about ([`closed_revolution_axis`]: a plane normal
/// to it counts) when `edge` is a circle about that axis — nine samples
/// of the edge at one radial distance and one axial height, each against
/// the kernel's `model` tolerance at the edge's scale. `None` otherwise,
/// and the ladder marches the stripe as before.
pub(in crate::blend) fn revolution_lane_axis(edge: &EdgeRecord, first: &BlendMate, second: &BlendMate) -> Option<(Vec3, Vec3)> {
    let (origin, axis) = closed_revolution_axis(first, second)?;
    let axis = axis.normalized().ok()?;
    let tolerance = crate::KernelTolerances::for_scale(
        march_model_scale(&edge.curve, edge.t0, edge.t1, 0.0).ok()?,
        1e-7,
    )
    .model;
    let mut reference: Option<(f64, f64)> = None;
    for step in 0..=8 {
        let t = edge.t0 + (edge.t1 - edge.t0) * step as f64 / 8.0;
        let offset = edge.curve.evaluate(t).ok()?.sub(origin);
        let height = offset.dot(axis);
        let radial = offset.sub(axis.scale(height)).length();
        match reference {
            None => reference = Some((radial, height)),
            Some((known_radial, known_height)) => {
                if (radial - known_radial).abs() > tolerance || (height - known_height).abs() > tolerance {
                    return None;
                }
            }
        }
    }
    (reference?.0 > tolerance).then_some((origin, axis))
}

/// What the exact revolution lane hands the network: the rows, every
/// station's parameter on them, and the stations with the OVERSHOOT ones
/// (outside the vertex window) replaced by the rows' own sections — the
/// march solved those on the carriers' ruled extension, where a partial
/// sweep continues as its tangent strip, not as the circle the rows follow.
/// The window's stations are the march's, bit for bit.
pub(in crate::blend) struct ExactRevolution {
    pub(in crate::blend) rows: FittedRows,
    pub(in crate::blend) parameters: Vec<f64>,
    pub(in crate::blend) stations: Vec<Station>,
}

/// The exact rows of a march whose stations between the vertices are ONE
/// section revolved about `axis` (the revolution mates' axis,
/// [`revolution_lane_axis`]), or the reason they are not (the lane then
/// declines to the ladder by name).
///
/// The spine takes the revolution CARRIER's own angular parameterisation:
/// the contact rail on that carrier is its latitude circle at the contact's
/// held coordinate (`iso_curve_v` / `iso_curve_u`), extended past the
/// carrier's sweep along the circle ([`extend_conic_ends`]) to cover the
/// overshoot and affinely on `[0, 1]`, so the support trim on the carrier is
/// a straight line in (u, v) by construction. Every other row is that circle
/// carried by the affine map that takes the carrier contact of the first
/// window station to the row's own point of that section (a radial scaling
/// about the axis with an axial shift: an affine map commutes with the
/// rational construction, so each row is the exact circle of its point and
/// the tensor product is the exact torus patch of the section, in the
/// carrier's parameterisation). The mate's trim is the rail's image through
/// the plane's affine chart when the mate is a plane
/// (`plane_chart_pcurve`), else the uv line through the window's two end
/// contacts; the surgery reads either once as shipped and refits only what
/// does not stand.
///
/// Witness (`1e-9 · (1 + scale)`, the solver's bar with two orders of
/// margin): every window station's two contacts, apex and centre are the
/// rows' at that station's parameter. Declines by name: no revolution
/// carrier among the mates; a contact drifting across the carrier or not
/// monotone along it; a window shorter than three stations; a section that
/// reaches the axis, or a centre circle inside the ball (the wall would
/// fold, which this lane does not probe); a circle whose extension would
/// carry a non-positive weight.
pub(in crate::blend) fn exact_revolution_rows(
    stations: &[Station],
    vertex_indices: [usize; 2],
    chamfer: bool,
    axis: (Vec3, Vec3),
    radius: f64,
    first: &BlendMate,
    second: &BlendMate,
) -> Result<ExactRevolution, String> {
    let (origin, axis) = axis;
    let (lo, hi) = (vertex_indices[0].min(vertex_indices[1]), vertex_indices[0].max(vertex_indices[1]));
    if stations.len() < 3 || hi >= stations.len() || hi < lo + 2 {
        return Err("fewer than three stations between the vertices".to_string());
    }
    let scale = stations.iter().fold(0.0f64, |worst, station| worst.max(station.center.length())).max(1.0);
    let bar = 1e-9 * (1.0 + scale);
    let is_plane = |mate: &BlendMate| matches!(mate.face.surface.analytic(), Some(crate::AnalyticSurface::Plane { .. }));
    // The carrier whose own parameterisation the rows take: the revolution
    // mate, the second when both are.
    let carrier_is_second = if !is_plane(second) {
        true
    } else if !is_plane(first) {
        false
    } else {
        return Err("no revolution carrier among the mates".to_string());
    };
    let (carrier, mate) = if carrier_is_second { (second, first) } else { (first, second) };
    let carrier_uv = |station: &Station| if carrier_is_second { station.uv2 } else { station.uv1 };
    let mate_uv = |station: &Station| if carrier_is_second { station.uv1 } else { station.uv2 };
    let contact = |station: &Station| if carrier_is_second { station.p2 } else { station.p1 };
    let surface = &carrier.face.surface;
    let [u0, u1] = surface.domain_u()?;
    let [v0, v1] = surface.domain_v()?;
    let window = &stations[lo..=hi];
    let (first_uv, last_uv) = (carrier_uv(&window[0]), carrier_uv(&window[window.len() - 1]));
    let along_u = (last_uv[0] - first_uv[0]).abs() >= (last_uv[1] - first_uv[1]).abs();
    let (along, held, held_span) = if along_u { (0usize, 1usize, v1 - v0) } else { (1, 0, u1 - u0) };
    let held_value = window.iter().map(|station| carrier_uv(station)[held]).sum::<f64>() / window.len() as f64;
    if window.iter().any(|station| (carrier_uv(station)[held] - held_value).abs() > 1e-8 * held_span.abs().max(1e-300)) {
        return Err("the carrier contact drifts across the carrier".to_string());
    }
    let along_of = |station: &Station| carrier_uv(station)[along];
    let sweep = along_of(&window[window.len() - 1]) - along_of(&window[0]);
    if !(sweep.abs() > 0.0) || window.windows(2).any(|pair| (along_of(&pair[1]) - along_of(&pair[0])) * sweep <= 0.0) {
        return Err("the carrier contact is not monotone along the carrier".to_string());
    }
    // The carrier's latitude circle at the held coordinate, run the
    // stations' way.
    let circle = if along_u { surface.iso_curve_v(held_value) } else { surface.iso_curve_u(held_value) }?;
    let [c0, c1] = circle.domain()?;
    let reversed = sweep < 0.0;
    let circle = if reversed { circle.reversed()? } else { circle };
    let on_curve = |along: f64| if reversed { c0 + c1 - along } else { along };
    // Every station's along coordinate: the window's own, the overshoot's at
    // the window's mean spacing past its ends.
    let step = sweep / (hi - lo) as f64;
    let along_at = |index: usize| -> f64 {
        if index < lo {
            along_of(&window[0]) - step * (lo - index) as f64
        } else if index > hi {
            along_of(&window[window.len() - 1]) + step * (index - hi) as f64
        } else {
            along_of(&stations[index])
        }
    };
    let (parameter_lo, parameter_hi) = (on_curve(along_at(0)), on_curve(along_at(stations.len() - 1)));
    if !(parameter_hi > parameter_lo) {
        return Err("the rows' parameter range is empty".to_string());
    }
    let rail = extend_conic_ends(&circle, parameter_lo, parameter_hi)?;
    let [r0, r1] = rail.domain()?;
    let mut rail = rail;
    for knot in rail.knots.iter_mut() {
        *knot = ((*knot - r0) / (r1 - r0)).clamp(0.0, 1.0);
    }
    let degree = rail.degree;
    for knot in rail.knots.iter_mut().take(degree + 1) {
        *knot = 0.0;
    }
    let count = rail.knots.len();
    for knot in rail.knots.iter_mut().skip(count - degree - 1) {
        *knot = 1.0;
    }
    let rail = NurbsCurve::new(degree, rail.knots.clone(), rail.control_points.clone())?;
    let mut parameters: Vec<f64> = (0..stations.len()).map(|index| (on_curve(along_at(index)) - parameter_lo) / (parameter_hi - parameter_lo)).collect();
    parameters[0] = 0.0;
    *parameters.last_mut().unwrap() = 1.0;
    if parameters.windows(2).any(|pair| !(pair[1] > pair[0])) {
        return Err("the rows' parameters are not increasing".to_string());
    }
    // The one section: the first window station, carried along the circle.
    let seat = &window[0];
    let axial = |point: Vec3| point.sub(origin).dot(axis);
    let radial = |point: Vec3| {
        let height = axial(point);
        point.sub(origin).sub(axis.scale(height))
    };
    let seat_contact = contact(seat);
    let contact_radius = radial(seat_contact).length();
    if contact_radius <= bar {
        return Err("the carrier contact sits on the axis".to_string());
    }
    let contact_height = axial(seat_contact);
    let map_row = |target: Vec3| -> Result<NurbsCurve, String> {
        let (target_radius, target_height) = (radial(target).length(), axial(target));
        let factor = target_radius / contact_radius;
        let control = rail
            .control_points
            .iter()
            .map(|control| {
                let weight = control.w;
                let point = Vec3::new(control.x / weight, control.y / weight, control.z / weight);
                let image = origin.add(axis.scale(axial(point) - contact_height + target_height)).add(radial(point).scale(factor));
                Vec4 { x: image.x * weight, y: image.y * weight, z: image.z * weight, w: weight }
            })
            .collect();
        NurbsCurve::new(rail.degree, rail.knots.clone(), control)
    };
    let centre_radius = radial(seat.center).length();
    let section_radius = radial(seat.p1).length().min(radial(seat.p2).length()).min(radial(seat.apex).length());
    if section_radius <= 1e-6 * scale {
        return Err("the section reaches the axis".to_string());
    }
    if !chamfer && centre_radius <= radius.abs() * (1.0 + 1e-6) {
        return Err("the centre circle is inside the ball (the wall would fold)".to_string());
    }
    let (cr, cs) = if carrier_is_second { (map_row(seat.p1)?, rail.clone()) } else { (rail.clone(), map_row(seat.p2)?) };
    let apex = map_row(seat.apex)?;
    let apex_weighted = NurbsCurve::new(
        apex.degree,
        apex.knots.clone(),
        apex.control_points.iter().map(|control| control.scale(seat.weight)).collect(),
    )?;
    let center = map_row(seat.center)?;
    let surface = crate::blend::rows::surface_from_rows(rail.degree, &cr, &cs, (!chamfer).then_some(&apex_weighted), false).map_err(|error| error.message)?;
    // The witness: every window station is the rows' section at its parameter.
    for index in lo..=hi {
        let (station, s) = (&stations[index], parameters[index]);
        let off = |curve: &NurbsCurve, point: Vec3| -> Result<f64, String> { Ok(curve.evaluate(s)?.sub(point).length()) };
        let worst = off(&cr, station.p1)?.max(off(&cs, station.p2)?).max(off(&center, station.center)?).max(off(&apex, station.apex)?);
        if !(worst <= bar) {
            return Err(format!("station {index} stands {worst:.3e} off the revolved section (bar {bar:.1e})"));
        }
    }
    // The support trims: the carrier's uv line, exact by construction; the
    // mate's through the plane's chart, or the line through the window's
    // ends when the mate is itself a revolution.
    let line = |from: [f64; 2], to: [f64; 2]| -> Result<NurbsCurve, String> {
        NurbsCurve::new(
            1,
            vec![0.0, 0.0, 1.0, 1.0],
            vec![Vec4::from_point(Vec3::new(from[0], from[1], 0.0), 1.0), Vec4::from_point(Vec3::new(to[0], to[1], 0.0), 1.0)],
        )
    };
    let carrier_point = |along: f64| -> [f64; 2] { if along_u { [along, held_value] } else { [held_value, along] } };
    let carrier_pcurve = line(carrier_point(along_at(0)), carrier_point(along_at(stations.len() - 1)))?;
    let mate_row = if carrier_is_second { &cr } else { &cs };
    let mate_pcurve = match plane_chart_pcurve(&mate.face.surface, mate_row) {
        Some(pcurve) => pcurve,
        None => {
            let (a, b) = (mate_uv(&window[0]), mate_uv(&window[window.len() - 1]));
            let (sa, sb) = (parameters[lo], parameters[hi]);
            let at = |s: f64| -> [f64; 2] {
                let fraction = (s - sa) / (sb - sa);
                [a[0] + fraction * (b[0] - a[0]), a[1] + fraction * (b[1] - a[1])]
            };
            line(at(0.0), at(1.0))?
        }
    };
    let (cr_pcurve, cs_pcurve) = if carrier_is_second { (mate_pcurve, carrier_pcurve) } else { (carrier_pcurve, mate_pcurve) };
    // The overshoot stations, read off the rows.
    let mut stations = stations.to_vec();
    for index in (0..lo).chain(hi + 1..stations.len()) {
        let s = parameters[index];
        let uv = |pcurve: &NurbsCurve| -> Result<[f64; 2], String> {
            let point = pcurve.evaluate(s)?;
            Ok([point.x, point.y])
        };
        stations[index] = Station {
            uv1: uv(&cr_pcurve)?,
            uv2: uv(&cs_pcurve)?,
            p1: cr.evaluate(s)?,
            p2: cs.evaluate(s)?,
            center: center.evaluate(s)?,
            weight: seat.weight,
            apex: apex.evaluate(s)?,
        };
    }
    Ok(ExactRevolution {
        rows: FittedRows {
            surface,
            cr,
            cs,
            cr_pcurve,
            cs_pcurve,
            u_domain: [0.0, 1.0],
            center: Some(center),
            exact_extrusion: false,
            vertex_stations: None,
            revolution_axis: Some((origin, axis)),
        },
        parameters,
        stations,
    })
}

/// `curve` continued past its ends ALONG ITS OWN CONIC pieces to the domain
/// `[a, b]` (one or both ends past the domain), or cut back to it (an end
/// inside the domain). The end piece is rewritten as the same rational
/// polynomial over the wider range — de Casteljau's polygon at the range's
/// local parameter, past 0 or 1 — so a rational quadratic circle piece
/// continues as that circle; the interior pieces are untouched. Refuses
/// when the rewritten polygon carries a non-positive weight (an extension
/// too far for the piece's parameterisation).
pub(in crate::blend) fn extend_conic_ends(curve: &NurbsCurve, a: f64, b: f64) -> Result<NurbsCurve, String> {
    let [d0, d1] = curve.domain()?;
    if !(b > a) || !a.is_finite() || !b.is_finite() {
        return Err(format!("extend_conic_ends: [{a}, {b}] is not a range"));
    }
    let degree = curve.degree;
    let tolerance = crate::KNOT_IDENTITY_TOL;
    // Every interior knot to full multiplicity: piecewise Bezier.
    let mut bezier = curve.clone();
    let mut interior: Vec<f64> = curve.knots.iter().copied().filter(|knot| *knot > d0 + tolerance && *knot < d1 - tolerance).collect();
    interior.dedup_by(|next, previous| (*next - *previous).abs() <= tolerance);
    for knot in interior {
        bezier = bezier.insert_knot(knot, degree)?;
    }
    // de Casteljau's polygons at local parameter `s` on one piece: the left
    // polygon covers local [0, s], the right one [s, 1] — for s past 0 or 1
    // the same algebra gives the piece over the extended range.
    let polygons = |piece: &[Vec4], s: f64| -> (Vec<Vec4>, Vec<Vec4>) {
        let mut levels: Vec<Vec<Vec4>> = vec![piece.to_vec()];
        while levels.last().map_or(false, |level| level.len() > 1) {
            let previous = levels.last().unwrap();
            let next: Vec<Vec4> = previous.windows(2).map(|pair| pair[0].scale(1.0 - s).add(pair[1].scale(s))).collect();
            levels.push(next);
        }
        let left: Vec<Vec4> = levels.iter().map(|level| level[0]).collect();
        let right: Vec<Vec4> = levels.iter().rev().map(|level| level[level.len() - 1]).collect();
        (left, right)
    };
    let mut knots = bezier.knots.clone();
    let mut control = bezier.control_points.clone();
    let count = knots.len();
    if a < d0 - tolerance {
        let piece_end = knots[degree + 1];
        let s = (a - d0) / (piece_end - d0);
        let (_, right) = polygons(&control[..=degree], s);
        control[..=degree].copy_from_slice(&right);
        for knot in knots.iter_mut().take(degree + 1) {
            *knot = a;
        }
    }
    if b > d1 + tolerance {
        let piece_start = knots[count - degree - 2];
        let s = (b - piece_start) / (d1 - piece_start);
        let last = control.len() - 1;
        let (left, _) = polygons(&control[last - degree..=last], s);
        control[last - degree..=last].copy_from_slice(&left);
        for knot in knots.iter_mut().skip(count - degree - 1) {
            *knot = b;
        }
    }
    // The curve constructor's own weight floor (`NurbsCurve::new`: weights
    // above 1e-12), asked here so the refusal is this function's.
    if control.iter().any(|point| !(point.w > 1e-12)) {
        return Err("extend_conic_ends: the extended piece carries a non-positive weight".to_string());
    }
    let mut extended = NurbsCurve::new(degree, knots, control)?;
    let [e0, e1] = extended.domain()?;
    if a > e0 + tolerance {
        extended = extended.split(a)?.1;
    }
    if b < e1 - tolerance {
        extended = extended.split(b)?.0;
    }
    Ok(extended)
}

/// A rail's pcurve on a PLANE mate: the rail's own control points carried
/// through the plane's affine chart — an exact rational image of the rail
/// in (u, v) — read back at 64 stations against the chart at the kernel's
/// `model` tolerance of the chart's scale. `None` when the mate is not a
/// plane, the chart is degenerate or the read-back fails (the closed
/// exact-revolution lane's rule, `stations.rs`, for an open rail).
pub(in crate::blend) fn plane_chart_pcurve(surface: &NurbsSurface, curve: &NurbsCurve) -> Option<NurbsCurve> {
    if !matches!(surface.analytic(), Some(crate::AnalyticSurface::Plane { .. })) {
        return None;
    }
    let [u0, u1] = surface.domain_u().ok()?;
    let [v0, v1] = surface.domain_v().ok()?;
    let origin = surface.evaluate(u0, v0).ok()?;
    let du = surface.evaluate(u1, v0).ok()?.sub(origin).scale(1.0 / (u1 - u0));
    let dv = surface.evaluate(u0, v1).ok()?.sub(origin).scale(1.0 / (v1 - v0));
    let (a, b, c) = (du.dot(du), du.dot(dv), dv.dot(dv));
    let det = a * c - b * b;
    let chart_scale = origin.length().max(du.length()).max(dv.length()).max(1.0);
    if !(det > 1e-24 * chart_scale.powi(4)) {
        return None;
    }
    let affine = |point: Vec3| -> [f64; 2] {
        let offset = point.sub(origin);
        let (r, s) = (offset.dot(du), offset.dot(dv));
        [u0 + (c * r - b * s) / det, v0 + (a * s - b * r) / det]
    };
    let mut pcurve = curve.clone();
    for control in &mut pcurve.control_points {
        let weight = control.w;
        let [u, v] = affine(Vec3::new(control.x / weight, control.y / weight, control.z / weight));
        control.x = u * weight;
        control.y = v * weight;
        control.z = 0.0;
    }
    let pcurve = NurbsCurve::new(pcurve.degree, pcurve.knots.clone(), pcurve.control_points.clone()).ok()?;
    let [t0, t1] = curve.domain().ok()?;
    let floor = crate::KernelTolerances::for_scale(chart_scale, 1e-7).model;
    for step in 0..=64 {
        let t = t0 + (t1 - t0) * step as f64 / 64.0;
        let uv = pcurve.evaluate(t).ok()?;
        let on_chart = surface.evaluate_extended(uv.x, uv.y).ok()?;
        if on_chart.sub(curve.evaluate(t).ok()?).length() > floor {
            return None;
        }
    }
    Some(pcurve)
}

/// One end of an open blend: the mates' boundary edges meeting the end
/// vertex, the single end face across the corner, and where the support
/// curves cross those boundary edges.
pub(in crate::blend) struct EndSurgery {
    /// Support-curve parameters at the crossings.
    pub(in crate::blend) cr_parameter: f64,
    pub(in crate::blend) cs_parameter: f64,
    /// The boundary edges being trimmed and their crossing parameters.
    pub(in crate::blend) first_edge_id: u64,
    pub(in crate::blend) first_edge_parameter: f64,
    pub(in crate::blend) second_edge_id: u64,
    pub(in crate::blend) second_edge_parameter: f64,
    pub(in crate::blend) end_face_id: u64,
    /// New vertices at the crossings (filled during surgery).
    pub(in crate::blend) transverse_curve: NurbsCurve,
    pub(in crate::blend) transverse_blend_pcurve: NurbsCurve,
    pub(in crate::blend) transverse_end_pcurve: NurbsCurve,
    /// True when a support crossing only resolved on a looser rung of
    /// [`support_crossing_tolerances`] than the historical 1e-7 — everything
    /// built from it is NEW, so the caller gates it through
    /// [`SpokeCrossings::gate`].
    pub(in crate::blend) escalated: bool,
}

/// How one end of a stripe is closed.
///
/// A stripe that runs out on its own terminates on the single face across the
/// corner ([`EndSurgery`], the §6.9 transverse construction).  A stripe that
/// runs into a corner shared with other stripes does NOT: it stops at the
/// corner ball's station, and the closure is the corner patch sewn onto every
/// stripe that meets there.  Both occupy the same slot of the blend face's
/// loop, which is why one enum covers them.
pub(in crate::blend) enum EndPlan {
    Free(EndSurgery),
    Corner(CornerEnd),
    Miter(MiterEnd),
    Cap(CapEnd),
}

/// A stripe end closed by a CAP: the stripe stops at the plane normal to
/// the edge at its vertex. The sliver's cross-section — the blend/plane
/// intersection plus one leg on each mate back to the vertex — is filled
/// as a planar bulkhead. The edge continuing past the vertex is left
/// untouched. This is the termination a fillet gets where its edge runs smoothly into an
/// edge that is not selected (tangent chain, no propagation), and it is the
/// cutter's flush end cap built without the cutter.
///
/// The orchestrator commits the arc, the two leg edges, the leg vertices and
/// the bulkhead; the surgery references the arc (like a corner end) and
/// splices each leg into its mate's loop between the rail and the
/// continuing edge.
pub(in crate::blend) struct CapEnd {
    /// Row parameters where the cap plane cuts each rail.
    pub(in crate::blend) cr_parameter: f64,
    pub(in crate::blend) cs_parameter: f64,
    /// Rim vertices: where the cap boundary meets the first and second mate.
    pub(in crate::blend) first_vertex: u64,
    pub(in crate::blend) second_vertex: u64,
    /// The section curve, committed in the blend loop's walk direction.
    pub(in crate::blend) arc_edge_id: u64,
    pub(in crate::blend) arc_blend_pcurve: NurbsCurve,
    /// The leg on each mate, committed from the rim vertex to the sharp
    /// vertex: `(edge id, pcurve on that mate from rim to vertex)`.
    pub(in crate::blend) first_leg: (u64, NurbsCurve),
    pub(in crate::blend) second_leg: (u64, NurbsCurve),
    /// The sharp vertex the legs meet at.
    pub(in crate::blend) sharp_vertex: u64,
}

/// A stripe end that meets ONE other stripe at a vertex whose third edge
/// stays sharp (§6.9.6): the two blend surfaces run past their corner-ball
/// stations and are trimmed against each other along the seam.
///
/// Unlike a star, the two rails of this stripe stop at DIFFERENT stations —
/// the shared-face rail at the ball's tangency point, the other rail where
/// the seam leaves the strip — and the end may be closed by one edge (the
/// seam) or two (the seam and then the sibling's section on this stripe's
/// other mate, when the seam leaves the sibling first).  The orchestrator
/// commits every edge and vertex before the surgery runs; the surgery only
/// references them.
pub(in crate::blend) struct MiterEnd {
    /// Row parameters where the first and second support rows stop.
    pub(in crate::blend) cr_parameter: f64,
    pub(in crate::blend) cs_parameter: f64,
    /// Rim vertices at those stops.
    pub(in crate::blend) first_vertex: u64,
    pub(in crate::blend) second_vertex: u64,
    /// The edges closing this end, in the order the blend face's loop walks
    /// them from the first rim to the second: `(edge id, stored forward along
    /// that walk, pcurve on the blend in walk direction)`.
    pub(in crate::blend) edges: Vec<(u64, bool, NurbsCurve)>,
}

/// A stripe end that stops at a shared corner ball.
///
/// Every field is decided BEFORE any cutting, from the ball solve
/// ([`crate::blend::solve_corner_ball`]): the rows are trimmed at the station
/// where the stripe's rolling ball IS the corner ball, and the rim vertices
/// are that ball's tangency points — the same two vertices the neighbouring
/// stripes on those faces stop at, so the mate loops close with no connector
/// and nothing is left over to identify and delete.
pub(in crate::blend) struct CornerEnd {
    /// Row parameters of the corner station on the two support rows.
    pub(in crate::blend) cr_parameter: f64,
    pub(in crate::blend) cs_parameter: f64,
    /// The ball's tangency vertices on the first and second mate — created
    /// once per corner and shared by every stripe that ends there.
    pub(in crate::blend) first_vertex: u64,
    pub(in crate::blend) second_vertex: u64,
    /// The end cross-section, already committed as an `EdgeRecord` oriented
    /// the way this blend face's loop walks it (so the blend uses it
    /// `forward = true` and the corner patch uses it `forward = false`).
    pub(in crate::blend) arc_edge_id: u64,
    /// Its pcurve in the blend surface's (t, z) space, in that same walk
    /// direction.
    pub(in crate::blend) arc_blend_pcurve: NurbsCurve,
}

impl EndPlan {
    pub(in crate::blend) fn cr_parameter(&self) -> f64 {
        match self {
            EndPlan::Free(end) => end.cr_parameter,
            EndPlan::Corner(end) => end.cr_parameter,
            EndPlan::Miter(end) => end.cr_parameter,
            EndPlan::Cap(end) => end.cr_parameter,
        }
    }

    pub(in crate::blend) fn cs_parameter(&self) -> f64 {
        match self {
            EndPlan::Free(end) => end.cs_parameter,
            EndPlan::Corner(end) => end.cs_parameter,
            EndPlan::Miter(end) => end.cs_parameter,
            EndPlan::Cap(end) => end.cs_parameter,
        }
    }

    pub(in crate::blend) fn free(&self) -> Option<&EndSurgery> {
        match self {
            EndPlan::Free(end) => Some(end),
            EndPlan::Corner(_) | EndPlan::Miter(_) | EndPlan::Cap(_) => None,
        }
    }

    /// The rim vertices an already-planned (corner, miter or cap) end
    /// supplies, first mate then second.
    fn planned_rims(&self) -> Option<(u64, u64)> {
        match self {
            EndPlan::Free(_) => None,
            EndPlan::Corner(end) => Some((end.first_vertex, end.second_vertex)),
            EndPlan::Miter(end) => Some((end.first_vertex, end.second_vertex)),
            EndPlan::Cap(end) => Some((end.first_vertex, end.second_vertex)),
        }
    }

    pub(in crate::blend) fn planned_rim(&self, side_first: bool) -> Option<u64> {
        self.planned_rims()
            .map(|(first, second)| if side_first { first } else { second })
    }
}

/// What one stripe's surgery committed, for the orchestrator to splice
/// against: the ids of the two rail edges it created (a network needs them to
/// insert a miter connector between a rail and the trimmed sharp edge) and
/// the blend face.
///
/// A mate the blend CONSUMED whole is not there any more: the surgery sewed
/// the blend onto the existing edge its rail retraces and dropped the face
/// (`detect_row_coincidence`).  `consumed` says so per side, first then
/// second, and on that side the rail id is that PRE-EXISTING edge rather
/// than one this stripe created — so no later pass may look the mate up or
/// treat the rail as new.
///
/// A rail that COLLAPSED to a point (its two stops within [`consumed_band`])
/// has no edge: its id is `None`, and its two rim vertices are one vertex,
/// listed in `identified` for the orchestrator to merge once every stripe is
/// in (the other rim may be another stripe's, not yet built).
pub(in crate::blend) struct SewnStripe {
    pub(in crate::blend) cr_edge_id: Option<u64>,
    pub(in crate::blend) cs_edge_id: Option<u64>,
    pub(in crate::blend) blend_face_id: u64,
    pub(in crate::blend) consumed: [bool; 2],
    /// The farthest any rail of this stripe was snapped onto an existing
    /// vertex or edge, for [`check_snap_closure`].
    pub(in crate::blend) snap: f64,
    /// Rim vertex pairs a collapsed rail makes one vertex.
    pub(in crate::blend) identified: Vec<(u64, u64)>,
}

/// How one of the four rim points of a stripe resolves.
pub(in crate::blend) enum RimResolution {
    /// Pristine free end: a fresh vertex at the row endpoint, with the
    /// boundary edge trimmed back to it.
    Fresh,
    /// A prior fillet already consumed the whole boundary edge — reuse its far
    /// vertex and delete the boundary.
    Consumed(u64),
    /// A shared corner ball supplies the vertex.  Nothing on the mate is
    /// trimmed here: the neighbouring stripe replaces the boundary outright.
    Corner(u64),
}

impl RimResolution {
    /// The existing vertex to reuse, or `None` when a fresh one is needed.
    pub(in crate::blend) fn existing(&self) -> Option<u64> {
        match self {
            RimResolution::Fresh => None,
            RimResolution::Consumed(id) | RimResolution::Corner(id) => Some(*id),
        }
    }

    pub(in crate::blend) fn is_consumed(&self) -> bool {
        matches!(self, RimResolution::Consumed(_))
    }

    /// True when this rim trims a boundary edge of the mate face.
    pub(in crate::blend) fn trims_boundary(&self) -> bool {
        matches!(self, RimResolution::Fresh)
    }
}

/// The refusal class of a free end whose boundary continues the blended
/// edge smoothly (lane L, 2026-10-14): read by the network's near-tangent
/// junction retry, which caps or flush-joins such an end.
pub(crate) const END_TANGENT_CONTINUATION_WHAT: &str = "end_tangent_continuation";
/// `1 + t_blended·t_boundary` at or under this is a continuation: the two
/// directions leaving the vertex are within 0.81° of opposite (a transverse
/// corner shallower than that is no corner the §6.9 construction can trim).
const END_CONTINUATION_BAND: f64 = 1e-4;

/// `BREP_BLEND_END_CONTINUATION=1` refuses a free end whose boundary
/// continues the blended edge (`END_TANGENT_CONTINUATION_WHAT`), `=0` keeps
/// the phantom-crossing reading. Promoted 2026-10-14 (lane L's record: case
/// manifest OK with it on, the fillet suites green).
const END_CONTINUATION_DEFAULT: bool = true;

thread_local! {
    static END_CONTINUATION_FOR_TEST: std::cell::Cell<Option<bool>> = const { std::cell::Cell::new(None) };
}

pub(crate) fn end_continuation_on() -> bool {
    if let Some(forced) = END_CONTINUATION_FOR_TEST.with(|forced| forced.get()) {
        return forced;
    }
    match std::env::var("BREP_BLEND_END_CONTINUATION").as_deref() {
        Ok("0") | Ok("off") => false,
        Ok("1") | Ok("on") => true,
        _ => END_CONTINUATION_DEFAULT,
    }
}


/// Resolve one FREE end of a stripe: the §6.9 transverse construction on the
/// single face across the corner.
///
/// Returns the end plan together with whether both support crossings landed
/// strictly inside the fitted-row domain — the caller marches with a growing
/// overshoot until they do.
pub(in crate::blend) fn resolve_free_end(
    solid: &BrepSolid,
    edge_id: u64,
    first: &BlendMate,
    second: &BlendMate,
    rows: &FittedRows,
    vertex: u64,
    at_start: bool,
) -> Result<(EndSurgery, bool), KernelRefusal> {
    let first_face = first.face;
    let second_face = second.face;
    let first_edge_id = boundary_edge_at_vertex(solid, first_face, vertex, edge_id)?;
    let second_edge_id = boundary_edge_at_vertex(solid, second_face, vertex, edge_id)?;
    if first_edge_id == second_edge_id {
        return Err(KernelRefusal::unsupported(
            KernelStage::Classify,
            "shared_end_boundary",
            "blend: mates share their end boundary edge (unsupported)",
        ));
    }
    // A boundary that CONTINUES the blended edge smoothly at the vertex — the
    // two leave the vertex in opposite directions — and that the 3D
    // intersector cannot cross is a tangent chain, not a corner: the rail runs
    // a rail offset beside it along its whole extension.  What the
    // parameter-space fallback found there on the 2026-10-12 star boss (lane
    // L, 2026-10-14) was a PHANTOM on the boundary's backward extension, which
    // the EXTEND polish then put on the end face inside the row's domain — a
    // seam wall ending at 55 % of its seam with a 0.24-long cut across the
    // side face, validating.  Refused by name instead — only where the 3D
    // crossing failed AND the parameter-space landing lies past the
    // boundary's own trim (an EXTEND landing) AND was accepted only past the
    // historical bar (escalated: polished, or read against
    // `pcurve_consistency`; the phantom stood 5.6e-2 off before its polish).
    // An exact extend landing on a straight continuation (2e-15 off on
    // `corner_convex`'s non-orthogonal corner) keeps its crossing, as does a
    // collinear boundary the rail does cross (lane L's first cut refused those too and
    // broke six corner_convex tests, a fold refusal and the slit-order
    // control). The network reads this refusal as its near-tangent junction
    // and caps the end.
    let blended = solid
        .edges
        .iter()
        .find(|candidate| candidate.id == edge_id)
        .ok_or_else(|| KernelRefusal::internal(KernelStage::Collect, "blended_edge", "blend: blended edge missing"))?;
    // Past the boundary's own trim: the EXTEND landing, read on two pcurves
    // extrapolated off their images.
    let landed_past_trim = |boundary: &EdgeRecord, t: f64| -> bool {
        let (low, high) = (boundary.t0.min(boundary.t1), boundary.t0.max(boundary.t1));
        let margin = 1e-9 * (high - low).max(1e-12);
        t < low - margin || t > high + margin
    };
    let continuation = |boundary: &EdgeRecord| -> Result<Option<KernelRefusal>, KernelRefusal> {
        if !end_continuation_on() {
            return Ok(None);
        }
        let away_blended = super::super::network::tangent_away_from(blended, vertex)?;
        let away_boundary = super::super::network::tangent_away_from(boundary, vertex)?;
        let gap = 1.0 + away_blended.dot(away_boundary);
        Ok((gap <= END_CONTINUATION_BAND).then(|| {
            KernelRefusal::unsupported(
                KernelStage::Classify,
                END_TANGENT_CONTINUATION_WHAT,
                format!(
                    "blend: boundary edge {} continues the blended edge {edge_id} smoothly at vertex {vertex} \
                     (the two leave it {gap:.2e} from collinear) and the rail does not cross it: a tangent \
                     chain, not a transverse end",
                    boundary.id
                ),
            )
        }))
    };
    let end_face = end_face_id(
        solid,
        first_face.id,
        second_face.id,
        first_edge_id,
        second_edge_id,
    )?;
    let first_boundary = solid
        .edges
        .iter()
        .find(|candidate| candidate.id == first_edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Collect,
                "end_boundary_edge",
                "blend: end boundary edge missing",
            )
        })?;
    let second_boundary = solid
        .edges
        .iter()
        .find(|candidate| candidate.id == second_edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Collect,
                "end_boundary_edge",
                "blend: end boundary edge missing",
            )
        })?;
    let first_boundary_coedge = first_face
        .loops
        .iter()
        .flat_map(|loop_record| &loop_record.coedges)
        .find(|coedge| coedge.edge_id == first_edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Collect,
                "mate_boundary_coedge",
                "blend: mate does not use its end boundary edge",
            )
        })?;
    let second_boundary_coedge = second_face
        .loops
        .iter()
        .flat_map(|loop_record| &loop_record.coedges)
        .find(|coedge| coedge.edge_id == second_edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Collect,
                "mate_boundary_coedge",
                "blend: mate does not use its end boundary edge",
            )
        })?;
    let end_record = solid
        .shells
        .iter()
        .flat_map(|shell| &shell.faces)
        .find(|face| face.id == end_face)
        .ok_or_else(|| {
            KernelRefusal::internal(KernelStage::Collect, "end_face", "blend: end face missing")
        })?;
    let (cr_parameter, first_edge_parameter, cr_escalated) =
        match support_crossing(solid, &rows.cr, first_boundary, at_start) {
            Ok((s, t, _, escalated)) => (s, t, escalated),
            Err(_) => {
                let landing = support_crossing_uv(
                    solid,
                    &rows.cr,
                    &rows.cr_pcurve,
                    first_boundary,
                    first_boundary_coedge,
                    at_start,
                    end_record,
                )?;
                if landing.2 && landed_past_trim(first_boundary, landing.1) {
                    if let Some(refusal) = continuation(first_boundary)? {
                        return Err(refusal);
                    }
                }
                landing
            }
        };
    let (cs_parameter, second_edge_parameter, cs_escalated) =
        match support_crossing(solid, &rows.cs, second_boundary, at_start) {
            Ok((s, t, _, escalated)) => (s, t, escalated),
            Err(_) => {
                let landing = support_crossing_uv(
                    solid,
                    &rows.cs,
                    &rows.cs_pcurve,
                    second_boundary,
                    second_boundary_coedge,
                    at_start,
                    end_record,
                )?;
                if landing.2 && landed_past_trim(second_boundary, landing.1) {
                    if let Some(refusal) = continuation(second_boundary)? {
                        return Err(refusal);
                    }
                }
                landing
            }
        };
    if std::env::var("BREP_DEBUG_NETWORK").is_ok() {
        // The landing, measured: where each trimmed rail ends and how far that
        // point stands off its mate (a rail end off its carrier cannot carry
        // a pcurve on it).
        for (label, row, parameter, mate, boundary, edge_parameter) in [
            ("cr", &rows.cr, cr_parameter, &first_face.surface, first_boundary, first_edge_parameter),
            ("cs", &rows.cs, cs_parameter, &second_face.surface, second_boundary, second_edge_parameter),
        ] {
            let point = row.evaluate_extended(parameter);
            let standoff = point.as_ref().ok().and_then(|p| crate::project_point_to_surface(mate, *p).ok()).map(|projection| projection.distance);
            let on_boundary = boundary.curve.evaluate_extended(edge_parameter);
            let gap = match (&point, &on_boundary) { (Ok(p), Ok(b)) => Some(p.sub(*b).length()), _ => None };
            eprintln!(
                "free end v{vertex} at_start={at_start} {label}: parameter {parameter:.6} of {:?}, rail point {:?}, {:?} off mate {}, boundary {} t={edge_parameter:.6} (edge t in [{:.4}, {:.4}]), rail-to-boundary gap {:?}",
                rows.u_domain, point.ok(), standoff, if label == "cr" { first_face.id } else { second_face.id }, boundary.id, boundary.t0, boundary.t1, gap
            );
            // The row's standoff from its mate along the way to the landing.
            let mut profile = String::new();
            let [d0, d1] = rows.u_domain;
            for step in 0..=20 {
                let u = d0 + (d1 - d0) * step as f64 / 20.0;
                let off = row.evaluate_extended(u).ok().and_then(|p| crate::project_point_to_surface(mate, p).ok()).map(|p| p.distance).unwrap_or(f64::NAN);
                profile.push_str(&format!(" {u:.2}:{off:.1e}"));
            }
            eprintln!("  {label} row standoff profile:{profile}");
        }
    }
    // A crossing past the boundary's far vertex is not an end on this face.
    // Reject it while choosing the end construction, so a tangent continuation
    // can use its cap instead of committing to a transverse end that surgery
    // will later reject. Exact extrusions keep the existing surgery/chain
    // dispatch: a flat cap there could hide an oversized tangent-chain wall.
    // Keep the same geometric snap band as that surgery.
    for (boundary, parameter) in [
        (first_boundary, first_edge_parameter),
        (second_boundary, second_edge_parameter),
    ] {
        let (near, far) = if boundary.start_vertex_id == vertex {
            (boundary.t0, boundary.t1)
        } else {
            (boundary.t1, boundary.t0)
        };
        if !rows.exact_extrusion && (parameter - far) * (far - near) > 0.0 {
            let distance = boundary.curve.evaluate_extended(parameter)
                .or_refuse(KernelStage::Refine, "evaluate_extended")?
                .sub(boundary.curve.evaluate(far).or_refuse(KernelStage::Refine, "evaluate")?)
                .length();
            if distance > consumed_band(solid) {
                return Err(KernelRefusal::unsupported(KernelStage::Refine,
                    "end_crossing_extent", format!(
                        "blend: end crossing is {distance:.3e} past boundary edge {}'s far vertex",
                        boundary.id)));
            }
        }
    }
    // Crossings must land strictly inside the row domain [0, 1] so the
    // surgery can split the support rows.  The start end sits near 0, the
    // finish end near 1.
    let in_range = if at_start {
        cr_parameter >= RIM_MARGIN && cs_parameter >= RIM_MARGIN
    } else {
        cr_parameter <= 1.0 - RIM_MARGIN && cs_parameter <= 1.0 - RIM_MARGIN
    };
    let clamp01 = |value: f64| value.clamp(0.0, 1.0);
    let (transverse, blend_pcurve, end_pcurve) = match exact_section_transverse(
        rows,
        end_record,
        clamp01(cr_parameter),
        clamp01(cs_parameter),
    )? {
        Some(exact) => exact,
        None => {
            // Lane L (2026-10-15): at an EXTEND landing the boundary stub and
            // the end arc leave their shared vertex in the SAME direction (a
            // zero-angle cusp in the end face's loop), so the arc's end tangent
            // must be the exact one or the fitted arc pokes past the boundary
            // (5e-8 on smFilletTest lCount_1, read as a bowtie), and its feet
            // must be on the end face's EXTENSION, not clamped onto its domain.
            let extends = |boundary: &EdgeRecord, parameter: f64| {
                let (low, high) = (boundary.t0.min(boundary.t1), boundary.t0.max(boundary.t1));
                let slack = 1e-9 * (high - low).abs().max(1.0);
                parameter < low - slack || parameter > high + slack
            };
            // Not on a PLANE.
            let planar_end = matches!(end_record.surface.analytic(), Some(crate::AnalyticSurface::Plane { .. }));
            let exact_ends = if crate::blend::planar_chart::ruled_chart_on() && !planar_end {
                [extends(first_boundary, first_edge_parameter), extends(second_boundary, second_edge_parameter)]
            } else {
                [false, false]
            };
            transverse_curve_with_ends(
                rows,
                end_record,
                clamp01(cr_parameter),
                clamp01(cs_parameter),
                exact_ends,
            )?
        }
    };
    Ok((
        EndSurgery {
            cr_parameter,
            cs_parameter,
            first_edge_id,
            first_edge_parameter,
            second_edge_id,
            second_edge_parameter,
            end_face_id: end_face,
            transverse_curve: transverse,
            transverse_blend_pcurve: blend_pcurve,
            transverse_end_pcurve: end_pcurve,
            escalated: cr_escalated || cs_escalated,
        },
        in_range,
    ))
}

/// The EXACT transverse curve of an exact-extrusion blend on a planar end face
/// perpendicular to the extrusion: the section itself — the surface's own
/// u-iso-curve, a rational arc — with its end-face pcurve mapped control
/// point by control point (an affine map keeps a rational curve rational).
/// `None` when the configuration is not that, and the fitted intersection
/// applies.  This is what keeps a straight×planar fillet's ends exact to the
/// mass integrator's bar instead of a cubic's chord error.
///
/// The control points go through the plane's own affine chart, UNCLAMPED.  A
/// section's middle control point is where the tangents at its two contacts
/// meet; between two planes that is the sharp corner, inside the end face, but
/// a cylinder's tangent meets the plane's past the corner, and past the end
/// face's parameter domain once the section turns far enough.  A projection
/// clamps it back onto the domain and bends the trim: the D-prism's plane ×
/// cylinder edge at r = 1.5 (a 131.8° section) put the middle control point at
/// y = 4.47 on a cap whose domain stops at y = 3, and the solid read 2.04 too
/// light while its blend face kept its exact area.
fn exact_section_transverse(
    rows: &FittedRows,
    end_face: &FaceRecord,
    cr_parameter: f64,
    cs_parameter: f64,
) -> Result<Option<(NurbsCurve, NurbsCurve, NurbsCurve)>, KernelRefusal> {
    if !rows.exact_extrusion || (cr_parameter - cs_parameter).abs() > 1e-9 {
        return Ok(None);
    }
    let Some(crate::AnalyticSurface::Plane {
        origin,
        u_dir,
        v_dir,
        u_domain,
        v_domain,
    }) = end_face.surface.analytic()
    else {
        return Ok(None);
    };
    let (a, b, c) = (u_dir.dot(*u_dir), u_dir.dot(*v_dir), v_dir.dot(*v_dir));
    let determinant = a * c - b * b;
    if determinant.abs() <= 1e-16 * (a * c).max(1.0) {
        return Ok(None);
    }
    let Ok(normal) = u_dir.cross(*v_dir).normalized() else {
        return Ok(None);
    };
    let [u0, u1] = rows.u_domain;
    let direction = rows
        .cr
        .evaluate(u1).or_refuse(KernelStage::Refine, "evaluate")?
        .sub(rows.cr.evaluate(u0).or_refuse(KernelStage::Refine, "evaluate")?)
        .normalized().or_refuse(KernelStage::Refine, "normalized")?;
    if direction.cross(normal).length() > 1e-9 {
        return Ok(None);
    }
    let section = rows.surface.iso_curve_u(cr_parameter).or_refuse(KernelStage::Refine, "iso_curve_u")?;
    let [v0, v1] = section.domain().or_refuse(KernelStage::Refine, "domain")?;
    let blend_pcurve =
        crate::sweep_topology::parameter_line(cr_parameter, v0, cr_parameter, v1).or_refuse(KernelStage::Refine, "parameter_line")?;
    let mut mapped = Vec::with_capacity(section.control_points.len());
    for control in &section.control_points {
        let offset = control.point().or_refuse(KernelStage::Refine, "point")?.sub(*origin);
        let (fu, fv) = (u_dir.dot(offset), v_dir.dot(offset));
        let u = (fu * c - fv * b) / determinant + u_domain[0];
        let v = (fv * a - fu * b) / determinant + v_domain[0];
        mapped.push(Vec4 {
            x: u * control.w,
            y: v * control.w,
            z: 0.0,
            w: control.w,
        });
    }
    let end_pcurve = NurbsCurve::new(section.degree, section.knots.clone(), mapped).or_refuse(KernelStage::Refine, "new")?;
    Ok(Some((section, blend_pcurve, end_pcurve)))
}

/// The corner of an END FACE the blend terminates on: the two boundary
/// coedges that meet at the old corner vertex, plus any DEGENERATE coedges
/// standing between them.
///
/// # Why a pole can stand between them
/// A pole is a legitimate corner.  Where the end face's apex sits ON the
/// blended edge's end vertex — revolve a profile about a line through one of
/// its own points and the cone raised on that point has its apex exactly
/// there — the face's loop reads `… first boundary, POLE, second boundary …`
/// and the pair is NOT adjacent.  The blend cuts that apex away, so the pole
/// coedges are located here and dropped by [`splice_end_corner`], and their
/// edges go with the other boundaries the surgery consumes.  With no pole the
/// walk is the plain `index + 1` this replaced.
pub(in crate::blend) struct EndCorner {
    pub(in crate::blend) loop_index: usize,
    /// The coedge the transverse curve is spliced AFTER, and the one it is
    /// spliced BEFORE, by stable id (indices shift when a sibling end shares
    /// this face).
    pub(in crate::blend) x_coedge_id: u64,
    pub(in crate::blend) y_coedge_id: u64,
    /// Whether `x` is the FIRST mate's boundary (it is `y`'s otherwise).
    pub(in crate::blend) x_is_first: bool,
    /// The pole coedges between them, in loop order, and their edges.
    pub(in crate::blend) pole_coedge_ids: Vec<u64>,
    pub(in crate::blend) pole_edge_ids: Vec<u64>,
}

/// Locate [`EndCorner`] on `face`, or `None` when the two boundary edges do
/// not meet at `corner` in any of its loops.
pub(in crate::blend) fn locate_end_corner(
    solid: &BrepSolid,
    face: &FaceRecord,
    corner: u64,
    first_edge_id: u64,
    second_edge_id: u64,
) -> Option<EndCorner> {
    let edge_of = |coedge: &CoedgeRecord| {
        solid
            .edges
            .iter()
            .find(|candidate| candidate.id == coedge.edge_id)
    };
    let traversal = |coedge: &CoedgeRecord| -> Option<(u64, u64)> {
        let edge = edge_of(coedge)?;
        Some(if coedge.forward {
            (edge.start_vertex_id, edge.end_vertex_id)
        } else {
            (edge.end_vertex_id, edge.start_vertex_id)
        })
    };
    // Only a degenerate edge AT THIS CORNER is stepped over — a pole
    // elsewhere on the loop is a boundary like any other.
    let is_pole = |coedge: &CoedgeRecord| -> bool {
        edge_of(coedge).is_some_and(|edge| {
            edge.degenerate && edge.start_vertex_id == corner && edge.end_vertex_id == corner
        })
    };
    for (loop_index, loop_record) in face.loops.iter().enumerate() {
        let count = loop_record.coedges.len();
        for index in 0..count {
            let this = &loop_record.coedges[index];
            if is_pole(this) {
                continue;
            }
            let mut pole_coedge_ids = Vec::new();
            let mut pole_edge_ids = Vec::new();
            let mut step = 1usize;
            while step < count && is_pole(&loop_record.coedges[(index + step) % count]) {
                let pole = &loop_record.coedges[(index + step) % count];
                pole_coedge_ids.push(pole.id);
                pole_edge_ids.push(pole.edge_id);
                step += 1;
            }
            if step >= count {
                continue;
            }
            let next = &loop_record.coedges[(index + step) % count];
            let this_pair = (this.edge_id, next.edge_id);
            if this_pair != (first_edge_id, second_edge_id)
                && this_pair != (second_edge_id, first_edge_id)
            {
                continue;
            }
            let (Some((_, this_end)), Some((next_start, _))) = (traversal(this), traversal(next))
            else {
                continue;
            };
            if this_end != corner || next_start != corner {
                continue;
            }
            return Some(EndCorner {
                loop_index,
                x_coedge_id: this.id,
                y_coedge_id: next.id,
                x_is_first: this.edge_id == first_edge_id,
                pole_coedge_ids,
                pole_edge_ids,
            });
        }
    }
    None
}

/// Splice the transverse coedge into a located [`EndCorner`]: drop the
/// consumed boundary coedge on either side, drop every pole coedge between
/// them, and keep the rest of the loop in order.
pub(in crate::blend) fn splice_end_corner(
    result: &mut BrepSolid,
    end_face: u64,
    corner: &EndCorner,
    transverse: CoedgeRecord,
    x_consumed: bool,
    y_consumed: bool,
) -> Result<(), KernelRefusal> {
    let face = result
        .shells
        .iter_mut()
        .flat_map(|shell| &mut shell.faces)
        .find(|face| face.id == end_face)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "end_face_lost",
                "blend: end face lost during surgery",
            )
        })?;
    let loop_record = face
        .loops
        .get_mut(corner.loop_index)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "corner_loop_lost",
                "blend: end-face corner loop lost during surgery",
            )
        })?;
    let count = loop_record.coedges.len();
    let seg_index = loop_record
        .coedges
        .iter()
        .position(|coedge| coedge.id == corner.x_coedge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "corner_coedge_lost",
                "blend: end-face corner coedge lost during surgery",
            )
        })?;
    for (offset, pole) in corner.pole_coedge_ids.iter().enumerate() {
        if loop_record.coedges[(seg_index + 1 + offset) % count].id != *pole {
            return Err(KernelRefusal::internal(
                KernelStage::Sew,
                "corner_pole_moved",
                "blend: end-face corner pole no longer between the boundaries",
            ));
        }
    }
    let y_index = (seg_index + 1 + corner.pole_coedge_ids.len()) % count;
    if loop_record.coedges[y_index].id != corner.y_coedge_id {
        return Err(KernelRefusal::internal(
            KernelStage::Sew,
            "corner_pair_adjacent",
            "blend: end-face corner pair no longer adjacent",
        ));
    }
    let mut transverse = Some(transverse);
    let mut rebuilt = Vec::with_capacity(count + 1);
    for index in 0..count {
        let coedge = &loop_record.coedges[index];
        if index == seg_index {
            if !x_consumed {
                rebuilt.push(coedge.clone());
            }
            rebuilt.push(transverse.take().expect("transverse spliced once"));
        } else if index == y_index {
            if !y_consumed {
                rebuilt.push(coedge.clone());
            }
        } else if !corner.pole_coedge_ids.contains(&coedge.id) {
            rebuilt.push(coedge.clone());
        }
    }
    check_loop_branches(&face.surface, &rebuilt)?;
    let loop_record = face
        .loops
        .get_mut(corner.loop_index)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "corner_loop_lost",
                "blend: end-face corner loop lost during surgery",
            )
        })?;
    loop_record.coedges = rebuilt;
    Ok(())
}

/// A rebuilt loop must walk ONE branch of a closed carrier's seam.
///
/// # The defect this is the invariant for
/// Where the blend's sliver runs past the corner into material a THIRD face
/// supplies, the §6.9 end cross-section — the genuine blend × end-carrier
/// intersection, exact to 1e-12 in 3D — meets that carrier on the far side of
/// its seam.  The captured case is a profile revolved about a line through one
/// of its own points: the end face is the cone raised on that point, its APEX
/// is the blend's own end vertex, and the arc passes BEHIND the apex, across
/// the part of the cone that lies inside the other body.  Its pcurve then
/// starts at `u = 1` where the boundary it must join ends at `u = 0` — the
/// same 3D point on the same seam, one period apart in parameter — and the
/// loop encloses the COMPLEMENT of the face: a 332-unit² cone came back as
/// 180, the part lost 158 units³ against a 15.6 closed form, and `validate`
/// saw none of it, because every endpoint still coincides in 3D.
///
/// Only a whole branch is refused (half a period or more).  The ordinary
/// junction mismatch is a fit residual — 4e-5 of the domain on that same part
/// — and the §6.9 EXTEND case, whose transverse deliberately lands past a
/// boundary's trim, keeps its endpoints on the branch it started on.
///
/// Representing what the arc actually wants needs the end face SPLIT along the
/// blend and the fragments re-classified; a NURBS carrier's parameter
/// extension does not wrap, so the loop cannot simply be unwrapped across the
/// seam. Until that lands the surgery refuses by name and the selection routes
/// to the cutter, whose boolean does split crossing faces.
fn check_loop_branches(surface: &NurbsSurface, coedges: &[CoedgeRecord]) -> Result<(), KernelRefusal> {
    let (closed_u, closed_v) = surface.closed_directions().unwrap_or((false, false));
    if !closed_u && !closed_v {
        return Ok(());
    }
    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
    let periods = [
        if closed_u { u1 - u0 } else { 0.0 },
        if closed_v { v1 - v0 } else { 0.0 },
    ];
    for index in 0..coedges.len() {
        let this = &coedges[index];
        let next = &coedges[(index + 1) % coedges.len()];
        let [_, this_end] = this.pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
        let [next_start, _] = next.pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
        let leaving = this.pcurve.evaluate(this_end).or_refuse(KernelStage::Refine, "evaluate")?;
        let arriving = next.pcurve.evaluate(next_start).or_refuse(KernelStage::Refine, "evaluate")?;
        for (axis, gap) in [
            (0usize, (leaving.x - arriving.x).abs()),
            (1, (leaving.y - arriving.y).abs()),
        ] {
            if periods[axis] > 0.0 && gap >= periods[axis] * 0.5 {
                return Err(KernelRefusal::unsupported(
                    KernelStage::Sew,
                    "loop_crosses_seam",
                    format!(
                        "blend: the end cross-section joins the end face's loop across its SEAM \
                         (the junction jumps {gap:.6} in the closed direction, a period of \
                         {:.6}) — the blend runs past the corner into material a third face \
                         supplies, and closing it needs that face SPLIT along the blend, which \
                         the rolling-ball surgery does not do",
                        periods[axis]
                    ),
                ));
            }
        }
    }
    Ok(())
}

/// The refusal for a blend whose rail runs PAST the far edge of a face it
/// must lie on, by more than [`consumed_band`] — read at a crossing past the
/// far vertex (`build_open_surgery`) and along a row that leaves the face
/// ([`detect_row_coincidence`]).  Not terminal: the r = 10.001 pinch past an
/// earlier round reads the same way, and the composition's ridge answers it.
pub(in crate::blend) const BLEND_WIDER_THAN_FACE: &str =
    "blend: the blend is wider than the face it must lie on —";

/// The refusal for a rail that shrinks to a point where the surgery has no
/// corner to identify its two rims through.  TERMINAL in the group lane: the
/// selection is a full-width blend, and the cutter composition's answer to a
/// full-width corner the network does not build was measured wrong (a 5F/8E/5V
/// body 8.0e-4 off its closed form at 1e-5 short of the width).
pub(crate) const RAIL_COLLAPSE_UNSUPPORTED: &str =
    "blend: a full-width rail collapse is not built here —";

/// The slug of that refusal (`UnsupportedGeometry`), minted in `edge/open.rs`
/// and `network.rs`.
pub(crate) const RAIL_COLLAPSE_WHAT: &str = "rail_collapse";

/// Is this a full-width rail collapse refusal ([`RAIL_COLLAPSE_UNSUPPORTED`])?
/// Read off the class and its slug, which the mint site and this check share.
pub(crate) fn is_rail_collapse(refusal: &KernelRefusal) -> bool {
    matches!(&refusal.class, RefusalClass::UnsupportedGeometry { what } if what == RAIL_COLLAPSE_WHAT)
}

/// How far inside the fitted-row domain a support crossing must land for the
/// surgery to be able to split the rows there.
pub(in crate::blend) const RIM_MARGIN: f64 = 2.0e-3;

/// The ONE band that decides whether a blend consumes a support face whole:
/// a rail that comes within it of an existing edge is snapped onto that edge,
/// and a rail that stops further short keeps the strip between them as a face.
///
/// Both sides of that choice have a kernel number, and the band is set from
/// them rather than picked:
///
/// * NOT snapping must not build a face narrower than the solid's own
///   [`crate::KernelTolerances::sliver`] — the width below which a feature is a
///   sliver to be repaired, not a face.  So the band is at least that.
/// * Snapping must not leave the shell over its closure bar: a rail snapped
///   across δ leaves its trims δ off the edge they were sewn onto, and
///   [`crate::shell_vector_areas`] allows each trim the refinement floor.
///   That limit depends on the SHAPE (the snapped length against the whole
///   shell's edge length), so it is not read into the band; it is measured on
///   the result by [`check_snap_closure`], and a snap that breaks it refuses.
///
/// Three sites ask the question and all three read this band, because any gap
/// between them is a shape the kernel builds wrong without refusing.  Measured
/// 2026-09-17 on a 20-cube's top-front edge, when the endpoint classifier in
/// `build_open_surgery` read `1e-6·(1 + boundary length)` and
/// [`detect_row_coincidence`] read `max(1e-7, 1e-6·far-vertex separation)`:
/// at r = 20.00002 the ends were consumed (2e-5 ≤ 2.1e-5) and the row was not
/// (2e-5 > 2e-5), so both boundary edges were deleted AND a fresh rail was
/// built coincident with the far edge — a zero-width face, 6F/10E/6V, that
/// `validate()` passes.  The third site is `check_support_extent`
/// (`fillet/analyze.rs`), which refuses a blend wider than its face: a
/// refusal slack WIDER than this band admits a rail that runs past the face
/// and is neither snapped nor refused (7F/11E/6V at r = 20.0000001 with a
/// 1e-7 snap band and the old slack, and no detector flags it).
///
/// That old `1e-6·extent` band was 2.1e-5 on the 20-cube, fifty times the
/// sliver tolerance, and every snap in the upper part of it built a shell over
/// its bar: 2.8e-4 against 2.0e-5 at 1e-5 short of the width, identical on
/// the base commit through the cutter composition's lane.
pub(crate) fn consumed_band(solid: &BrepSolid) -> f64 {
    crate::KernelTolerances::for_solid(solid, 1e-7).sliver
}

/// The refusal for a snap that breaks the shell's closure — see
/// [`check_snap_closure`].
pub(crate) const CONSUMED_SNAP_UNSOUND: &str =
    "blend: the blend is within a sliver of its face's width without matching it";

/// The slug of that refusal (`IllPosed`).
pub(crate) const CONSUMED_SNAP_WHAT: &str = "consumed_snap";

/// Is this the sliver-snap refusal ([`CONSUMED_SNAP_UNSOUND`])?
/// Read off the class and its slug, which the mint site and this check share.
pub(crate) fn is_consumed_snap(refusal: &KernelRefusal) -> bool {
    matches!(&refusal.class, RefusalClass::IllPosed { what } if what == CONSUMED_SNAP_WHAT)
}

/// Accept a body whose rails were snapped across `snap` only while its shells
/// still close: a snap is a trim moved off the edge it claims, and a shell
/// whose trims no longer enclose zero vector area is a wrong body that
/// `validate()` passes.  A snap inside the refinement floor is exactly what
/// every trim is allowed, so it is not measured.
///
/// Refusing here is a statement about the WINDOW between the two rules of
/// [`consumed_band`]: this rail is too close to the far edge to keep the strip
/// as a face (narrower than the sliver tolerance) and too far to be snapped
/// onto it soundly on this shape.  The input is measured only when the result
/// fails, so a body that arrived already open is not blamed on the snap.
pub(in crate::blend) fn check_snap_closure(
    input: &BrepSolid,
    result: &BrepSolid,
    snap: f64,
) -> Result<(), KernelRefusal> {
    if !(snap > crate::pcurve::PCURVE_REFINEMENT_TOLERANCE) {
        return Ok(());
    }
    let report = crate::shell_vector_areas(result);
    if !report.is_flagged() && report.unreadable.is_empty() {
        return Ok(());
    }
    if crate::shell_vector_areas(input).is_flagged() {
        return Ok(());
    }
    let (residual, bar) = report
        .worst_shell()
        .map_or((f64::NAN, f64::NAN), |shell| (shell.residual(), shell.bar));
    Err(KernelRefusal::ill_posed(
        KernelStage::Validate,
        CONSUMED_SNAP_WHAT,
        format!(
            "{CONSUMED_SNAP_UNSOUND} — its rail ends {snap:.3e} from the face's far edge, inside \
             the {:.3e} sliver band, so the strip between them is no face; and snapped onto that \
             edge, the shell's trims close to {residual:.3e} against a bar of {bar:.3e}{}. Make \
             the size the width, or move it away from the width by more than the band",
            consumed_band(input),
            if report.unreadable.is_empty() {
                String::new()
            } else {
                format!(" ({} faces unreadable)", report.unreadable.len())
            }
        ),
    ))
}

/// A detected full-span coincidence of a blend support row with an EXISTING
/// mate-face boundary edge (OCCT PR #1449 lesson 3, the r = face-width
/// class).  The surgery sews the blend face straight onto `edge_id` instead
/// of creating a coincident fresh support edge, and drops the fully-consumed
/// mate face whose loop collapsed between the two coincident curves.
pub(in crate::blend) struct RowSew {
    /// The existing boundary edge the blend row retraces.
    pub(in crate::blend) edge_id: u64,
    /// Traversal sense the blend face's coedge on that edge must carry
    /// (equal to the dropped face's old sense, so manifold pairing with the
    /// surviving neighbour face is preserved).
    pub(in crate::blend) forward: bool,
    /// Lesson-7 leftovers of the collapsing loop: zero-span edges wholly
    /// owned by the dropped face, collapsed (deleted) with it.
    pub(in crate::blend) collapsed_edges: Vec<u64>,
    /// The largest distance from the row to the edge it was sewn onto: how
    /// far the snap moved the rail.
    pub(in crate::blend) deviation: f64,
}

/// Detect that the trimmed support `row` coincides with an existing boundary
/// edge of the mate face over that edge's FULL span, so the mate face is
/// fully consumed (OCCT PR #1449 lesson 3 — the curve-level analog of the
/// endpoint-level `classify` rule in open.rs).
///
/// Fail-safe contract: ANY ambiguity returns `Ok(None)` and the caller keeps
/// the existing surgery path (the sliver face survives, downstream validation
/// still guards it).  `Some` requires ALL of:
///
/// * both end crossings already endpoint-consumed onto DISTINCT far vertices;
/// * the mate face is single-loop (dropping a face with holes would orphan
///   them) and its loop consists of the blended edge, the two consumed
///   boundary edges, exactly ONE existing edge joining the far vertices, and
///   (lesson 7) only zero-span collapse leftovers wholly owned by this face;
/// * the row retraces that edge under the AFFINE parameter map the sewn
///   coedge will carry — sampled in 3D, the exact metric `validate`'s
///   `adaptive_coedge_error` measures (the pr1449 doc sketches this check in
///   the mate's UV space; 3D avoids the anisotropic UV-metric conversion and
///   validates the true downstream contract) — within `band`, the
///   [`consumed_band`] the caller classified both ends against, with the
///   OPPOSITE map refused;
/// * relative orientation (the parameter-delta sign of the matching map)
///   agrees with the vertex correspondence, with the dropped face's own use
///   of the edge, and with exactly one opposing use by the surviving
///   neighbour (manifold pairing).
#[allow(clippy::too_many_arguments)]
pub(in crate::blend) fn detect_row_coincidence(
    solid: &BrepSolid,
    mate: &BlendMate,
    blended_edge_id: u64,
    boundary_start_id: u64,
    boundary_finish_id: u64,
    start_consumed: bool,
    start_far_vertex: u64,
    finish_consumed: bool,
    finish_far_vertex: u64,
    row: &NurbsCurve,
    row_traversed_from_start: bool,
    band: f64,
) -> Result<Option<RowSew>, KernelRefusal> {
    if !(start_consumed && finish_consumed) {
        return Ok(None);
    }
    if start_far_vertex == finish_far_vertex {
        return Ok(None);
    }
    if mate.face.loops.len() != 1 || mate.loop_index != 0 {
        return Ok(None);
    }
    let edge_by_id = |id: u64| solid.edges.iter().find(|candidate| candidate.id == id);
    let vertex_point = |id: u64| {
        solid
            .vertices
            .iter()
            .find(|candidate| candidate.id == id)
            .map(|vertex| vertex.point)
    };
    // The ONE existing edge joining the two far vertices, plus any leftovers.
    let mut sewn_coedge: Option<&CoedgeRecord> = None;
    let mut leftovers: Vec<u64> = Vec::new();
    for coedge in &mate.face.loops[0].coedges {
        if coedge.edge_id == blended_edge_id
            || coedge.edge_id == boundary_start_id
            || coedge.edge_id == boundary_finish_id
        {
            continue;
        }
        let Some(candidate) = edge_by_id(coedge.edge_id) else {
            return Err(KernelRefusal::internal(
                KernelStage::Classify,
                "mate_loop_edge",
                "blend: mate loop references a missing edge",
            ));
        };
        let joins = (candidate.start_vertex_id == start_far_vertex
            && candidate.end_vertex_id == finish_far_vertex)
            || (candidate.start_vertex_id == finish_far_vertex
                && candidate.end_vertex_id == start_far_vertex);
        if joins {
            if sewn_coedge.is_some() {
                return Ok(None);
            }
            sewn_coedge = Some(coedge);
        } else if !leftovers.contains(&candidate.id) {
            leftovers.push(candidate.id);
        }
    }
    let Some(sewn_coedge) = sewn_coedge else {
        return Ok(None);
    };
    let Some(target) = edge_by_id(sewn_coedge.edge_id) else {
        return Ok(None);
    };
    let (Some(start_point), Some(finish_point)) = (
        vertex_point(start_far_vertex),
        vertex_point(finish_far_vertex),
    ) else {
        return Ok(None);
    };
    let extent = finish_point.sub(start_point).length();
    if extent <= band * 4.0 {
        // Far vertices too close to discriminate orientation — ambiguous.
        return Ok(None);
    }
    // Orientation by vertex correspondence; confirmed below by the matching
    // affine map (the parameter-delta sign) and the manifold senses.
    let aligned = target.start_vertex_id == start_far_vertex;
    let [row_t0, row_t1] = row.domain().or_refuse(KernelStage::Refine, "domain")?;
    // Which way is INTO the face from the far edge: toward the blended edge.
    // A row that leaves the band on the other side runs outside its face —
    // both ends read consumed only because each end's crossing is clamped to
    // the far vertex, and the pristine path would build a lens face between
    // the rail and the edge where the part has no material at all
    // (7F/11E/6V at r = 20.0000008 on a 20-cube, closure clean).
    let blended_middle = edge_by_id(blended_edge_id)
        .map(|blended| blended.curve.evaluate(0.5 * (blended.t0 + blended.t1)))
        .transpose()
        .or_refuse(KernelStage::Classify, "evaluate")?;
    const ROW_COINCIDENCE_SAMPLES: usize = 16;
    let mut deviation = 0.0f64;
    for sample in 0..=ROW_COINCIDENCE_SAMPLES {
        let fraction = sample as f64 / ROW_COINCIDENCE_SAMPLES as f64;
        let on_row = row.evaluate(row_t0 + (row_t1 - row_t0) * fraction).or_refuse(KernelStage::Refine, "evaluate")?;
        let edge_fraction = if aligned { fraction } else { 1.0 - fraction };
        let on_edge = target
            .curve
            .evaluate(target.t0 + (target.t1 - target.t0) * edge_fraction).or_refuse(KernelStage::Refine, "evaluate")?;
        let gap = on_row.sub(on_edge).length();
        if gap > band {
            let outside = blended_middle
                .is_some_and(|middle| on_row.sub(on_edge).dot(middle.sub(on_edge)) < 0.0);
            if outside {
                return Err(KernelRefusal::unsupported(
                    KernelStage::Classify,
                    "blend_wider_than_face",
                    format!(
                        "{BLEND_WIDER_THAN_FACE} its rail runs {gap:.3e} past face {}'s far edge \
                         {} (more than the {band:.3e} it can be snapped across), outside the face",
                        mate.face.id, target.id
                    ),
                ));
            }
            return Ok(None);
        }
        deviation = deviation.max(gap);
    }
    // The OPPOSITE map must NOT also fit — a both-ways match means the edge
    // is degenerate at this band and orientation is ambiguous.
    let quarter_row = row.evaluate(row_t0 + (row_t1 - row_t0) * 0.25).or_refuse(KernelStage::Refine, "evaluate")?;
    let opposite_fraction = if aligned { 0.75 } else { 0.25 };
    let quarter_opposite = target
        .curve
        .evaluate(target.t0 + (target.t1 - target.t0) * opposite_fraction).or_refuse(KernelStage::Refine, "evaluate")?;
    if quarter_row.sub(quarter_opposite).length() <= band {
        return Ok(None);
    }
    // Lesson 7: leftover loop edges are tolerated only as ZERO-SPAN collapse
    // edges (both ends on one vertex, no geometric extent) wholly owned by
    // the dropped face — collapsed with it, never a generic deletion.
    for id in &leftovers {
        let Some(leftover) = edge_by_id(*id) else {
            return Ok(None);
        };
        if leftover.start_vertex_id != leftover.end_vertex_id {
            return Ok(None);
        }
        let a = leftover.curve.evaluate(leftover.t0).or_refuse(KernelStage::Refine, "evaluate")?;
        let mid = leftover
            .curve
            .evaluate(leftover.t0 + (leftover.t1 - leftover.t0) * 0.5).or_refuse(KernelStage::Refine, "evaluate")?;
        let b = leftover.curve.evaluate(leftover.t1).or_refuse(KernelStage::Refine, "evaluate")?;
        if a.sub(b).length() > band || a.sub(mid).length() > band {
            return Ok(None);
        }
        let used_elsewhere = solid.shells.iter().flat_map(|shell| &shell.faces).any(|face| {
            face.id != mate.face.id
                && face
                    .loops
                    .iter()
                    .flat_map(|loop_record| &loop_record.coedges)
                    .any(|coedge| coedge.edge_id == *id)
        });
        if used_elsewhere {
            return Ok(None);
        }
    }
    // Manifold pairing: the blend's traversal of the row determines the sense
    // its coedge on the existing edge must carry; it must equal the dropped
    // face's old sense (whose partner is the surviving neighbour) and be
    // opposed by exactly one use elsewhere.
    let traversal_start = if row_traversed_from_start {
        start_far_vertex
    } else {
        finish_far_vertex
    };
    let forward = target.start_vertex_id == traversal_start;
    if sewn_coedge.forward != forward {
        return Ok(None);
    }
    let mut opposing_uses = 0usize;
    for face in solid.shells.iter().flat_map(|shell| &shell.faces) {
        if face.id == mate.face.id {
            continue;
        }
        for coedge in face.loops.iter().flat_map(|loop_record| &loop_record.coedges) {
            if coedge.edge_id != target.id {
                continue;
            }
            if coedge.forward == forward {
                return Ok(None);
            }
            opposing_uses += 1;
        }
    }
    if opposing_uses != 1 {
        return Ok(None);
    }
    Ok(Some(RowSew {
        edge_id: target.id,
        forward,
        collapsed_edges: leftovers,
        deviation,
    }))
}

/// The boundary edge of `face` that meets the blended `exclude` edge at
/// `vertex`.  Uses loop adjacency: the answer is the `exclude` coedge's
/// loop-neighbour that shares `vertex`.  This DISAMBIGUATES the case where a
/// prior fillet has left more than one edge touching `vertex` (the original
/// plus a fresh fillet/support edge) — a plain incidence scan would see
/// several and could not tell which continues this face's border.
pub(in crate::blend) fn boundary_edge_at_vertex(
    solid: &BrepSolid,
    face: &FaceRecord,
    vertex: u64,
    exclude: u64,
) -> Result<u64, KernelRefusal> {
    let edge_at = |edge_id: u64| -> Option<&EdgeRecord> {
        solid.edges.iter().find(|edge| edge.id == edge_id)
    };
    // A collapsed pole coedge carries a chart boundary but no transverse
    // curve. Skip it when locating the physical boundary of the fillet.
    // Primary: the loop-neighbour of the blended edge that shares the vertex.
    for loop_record in &face.loops {
        let count = loop_record.coedges.len();
        if count == 0 {
            continue;
        }
        for index in 0..count {
            if loop_record.coedges[index].edge_id != exclude {
                continue;
            }
            // Prefer the neighbour on the side of `vertex`, then the other.
            for step in [count - 1, 1] {
                let neighbour = &loop_record.coedges[(index + step) % count];
                if neighbour.edge_id == exclude {
                    continue;
                }
                if let Some(edge) = edge_at(neighbour.edge_id) {
                    if !edge.degenerate
                        && (edge.start_vertex_id == vertex || edge.end_vertex_id == vertex)
                    {
                        return Ok(edge.id);
                    }
                }
            }
        }
    }
    // Fallback: a plain incidence scan (used when the blended edge is not in
    // this face's loops, e.g. degenerate inputs) — accept a single unique hit.
    let mut found = None;
    for loop_record in &face.loops {
        for coedge in &loop_record.coedges {
            if coedge.edge_id == exclude {
                continue;
            }
            let edge = edge_at(coedge.edge_id).ok_or_else(|| {
                KernelRefusal::internal(
                    KernelStage::Collect,
                    "loop_edge",
                    "blend: loop references missing edge",
                )
            })?;
            if !edge.degenerate
                && (edge.start_vertex_id == vertex || edge.end_vertex_id == vertex)
            {
                if found.is_some() && found != Some(edge.id) {
                    return Err(KernelRefusal::ill_posed(
                        KernelStage::Collect,
                        "boundary_edge_ambiguous",
                        "blend: multiple boundary edges at the end vertex",
                    ));
                }
                found = Some(edge.id);
            }
        }
    }
    found.ok_or_else(|| {
        KernelRefusal::internal(
            KernelStage::Collect,
            "boundary_edge_none",
            "blend: no boundary edge at the end vertex",
        )
    })
}

pub(in crate::blend) fn end_face_id(
    solid: &BrepSolid,
    first_face: u64,
    second_face: u64,
    first_edge: u64,
    second_edge: u64,
) -> Result<u64, KernelRefusal> {
    let mut found = None;
    for shell in &solid.shells {
        for face in &shell.faces {
            if face.id == first_face || face.id == second_face {
                continue;
            }
            let uses_first = face
                .loops
                .iter()
                .flat_map(|l| &l.coedges)
                .any(|coedge| coedge.edge_id == first_edge);
            let uses_second = face
                .loops
                .iter()
                .flat_map(|l| &l.coedges)
                .any(|coedge| coedge.edge_id == second_edge);
            if uses_first && uses_second {
                if found.is_some() {
                    return Err(KernelRefusal::ill_posed(
                        KernelStage::Collect,
                        "end_face_ambiguous",
                        "blend: ambiguous end face",
                    ));
                }
                found = Some(face.id);
            }
        }
    }
    found.ok_or_else(|| {
        KernelRefusal::unsupported(
            KernelStage::Collect,
            "end_face_none",
            "blend: no single end face across the corner",
        )
    })
}

/// Crossing of the fitted support pcurve with a boundary edge's pcurve in
/// the MATE FACE's parameter space — robust for free-form boundaries, and
/// capable of landing BEYOND the boundary's trim (the §6.9 "trim or
/// EXTEND" case).  Returns (support parameter, boundary EDGE parameter).
pub(in crate::blend) fn support_crossing_uv(
    solid: &BrepSolid,
    support: &NurbsCurve,
    support_pcurve: &NurbsCurve,
    boundary: &EdgeRecord,
    boundary_coedge: &CoedgeRecord,
    at_start: bool,
    end_face: &FaceRecord,
) -> Result<(f64, f64, bool), KernelRefusal> {
    let [s0, s1] = support_pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
    let [q0, q1] = boundary_coedge.pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
    let seed_s = if at_start { s0 } else { s1 };
    let mut solved = None;
    for seed_fraction in [0.0f64, 1.0, 0.5] {
        let mut x = [seed_s, q0 + (q1 - q0) * seed_fraction];
        let mut converged = false;
        for _ in 0..NEWTON_ITERATIONS {
            let a = support_pcurve.evaluate_extended(x[0]).or_refuse(KernelStage::Refine, "evaluate_extended")?;
            let b = boundary_coedge.pcurve.evaluate_extended(x[1]).or_refuse(KernelStage::Refine, "evaluate_extended")?;
            let residual = [a.x - b.x, a.y - b.y];
            if residual[0].abs().max(residual[1].abs()) <= 1e-12 {
                converged = true;
                break;
            }
            let step = 1e-8;
            let a_probe = support_pcurve.evaluate_extended(x[0] + step).or_refuse(KernelStage::Refine, "evaluate_extended")?;
            let b_probe = boundary_coedge.pcurve.evaluate_extended(x[1] + step).or_refuse(KernelStage::Refine, "evaluate_extended")?;
            let jacobian = [
                [(a_probe.x - a.x) / step, -(b_probe.x - b.x) / step],
                [(a_probe.y - a.y) / step, -(b_probe.y - b.y) / step],
            ];
            let determinant = jacobian[0][0] * jacobian[1][1] - jacobian[0][1] * jacobian[1][0];
            if determinant.abs() <= 1e-16 {
                break;
            }
            x[0] -= (residual[0] * jacobian[1][1] - residual[1] * jacobian[0][1]) / determinant;
            x[1] -= (residual[1] * jacobian[0][0] - residual[0] * jacobian[1][0]) / determinant;
        }
        if converged {
            // Prefer the solution nearest the support's queried end.
            let rank = (x[0] - seed_s).abs();
            if solved
                .map(|(existing_s, _): (f64, f64)| (existing_s - seed_s).abs() > rank)
                .unwrap_or(true)
            {
                solved = Some((x[0], x[1]));
            }
        }
    }
    let (s, q) = solved.ok_or_else(|| {
        KernelRefusal::unsupported(
            KernelStage::Intersect,
            "support_crossing_uv_none",
            "blend: support curve does not reach the end boundary edge",
        )
    })?;
    // The row was marched with an overshoot past both vertices, so a genuine
    // crossing always lies within its domain; a solution outside it is the
    // pcurve extrapolated off its surface image.
    let span = (s1 - s0).abs();
    if s < s0 - 1e-6 * span || s > s1 + 1e-6 * span {
        return Err(KernelRefusal::unsupported(
            KernelStage::Intersect,
            "support_crossing_uv_outside",
            format!(
                "blend: support curve does not reach the end boundary edge (the parameter-space \
                 crossing lands at {s:.4}, outside the row's [{s0:.4}, {s1:.4}])"
            ),
        ));
    }
    // Map the pcurve parameter to the boundary EDGE parameter through the
    // affine traversal contract.
    let fraction = (q - q0) / (q1 - q0);
    let fraction = if boundary_coedge.forward {
        fraction
    } else {
        1.0 - fraction
    };
    let t = boundary.t0 + (boundary.t1 - boundary.t0) * fraction;
    // A crossing that lands PAST the boundary's trim (the §6.9 EXTEND
    // landing) was solved on two pcurves EXTRAPOLATED off their surface
    // images: the boundary's past its own end, the support's past the march.
    // That extrapolation is the pcurve fit's polynomial continuation, not the
    // boundary's carrier intersection, and its error grows with the square
    // of the reach — on the 2026-10-12 report (a free-form body's short cap
    // edge continued past the vertex a rolling ball of r = 0.5 … 1.0 ends
    // beyond) the GENUINE crossing read 4.1e-3 … 1.6e-2 off in 3D and was
    // refused below as a phantom, while r ≤ 0.4 passed on the same path.
    // The boundary's 3D curve is no arbiter there either: it is the fit of
    // the carriers' intersection over the edge and a short margin, and
    // `evaluate_extended` continues it LINEARLY past that, so the closest
    // approach of the two 3D curves still read 1.4e-3 … 3.6e-3 and the rim
    // vertex it placed sat off the end face — the cap's loop then crossed
    // itself at r = 1.  What the landing IS: the point where the support row
    // (on the mate) meets the END FACE's surface, since the boundary is that
    // face's intersection with the mate.  The row parameter is polished onto
    // the end surface (three equations in `s` and the end face's `uv`, read
    // through the extensions: the landing lies outside the end face's trim),
    // and the boundary parameter is the nearest point of its extended curve
    // to that landing — the surgery uses it only to extend the edge's range,
    // and re-marches the carriers to the rim vertex itself.  A boundary that
    // CONTINUES the blended edge keeps its rail a rail offset away along its
    // whole extension, so the polish either does not converge or leaves the
    // row's domain, and the historical reading (the phantom test below)
    // stands.  Inside the trim the pcurves are on their images and nothing
    // here changes.
    let extended = fraction < -1e-9 || fraction > 1.0 + 1e-9;
    // Model scale from the boundary's own extent — never from the evaluated
    // points, which an extrapolation can throw arbitrarily far and thereby
    // make any gap look small.
    let scale = boundary
        .curve
        .evaluate(boundary.t0).or_refuse(KernelStage::Refine, "evaluate")?
        .length()
        .max(boundary.curve.evaluate(boundary.t1).or_refuse(KernelStage::Refine, "evaluate")?.length())
        .max(1.0);
    let historical = 1e-6 * scale;
    let band = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
    // The reading as it always was: the two 3D curves at the pcurve landing.
    let raw_gap = support
        .evaluate_extended(s).or_refuse(KernelStage::Refine, "evaluate_extended")?
        .sub(boundary.curve.evaluate_extended(t).or_refuse(KernelStage::Refine, "evaluate_extended")?)
        .length();
    // Only an extended landing the historical bar would not pass as it
    // stands is polished; one inside that bar ships exactly as it always has.
    let mut polished = false;
    let (s, t) = if extended && raw_gap > historical {
        let vertex_t = if fraction < 0.0 { boundary.t0 } else { boundary.t1 };
        match polish_crossing_onto_end_face(support, boundary, end_face, s, t, vertex_t)? {
            Some((s_polished, t_polished))
                if s_polished >= s0 - 1e-6 * span && s_polished <= s1 + 1e-6 * span =>
            {
                polished = true;
                (s_polished, t_polished)
            }
            _ => (s, t),
        }
    } else {
        (s, t)
    };
    // A UV solution is only a crossing if the two 3D curves meet there.  A
    // rational pcurve evaluated past its domain sweeps away from its surface
    // image, so a boundary the rail never reaches — one that CONTINUES the
    // blended edge tangentially, the rail running parallel to it — can still
    // hand Newton a phantom (s, q).  The §6.9 EXTEND case, a short straight
    // boundary met just past its end, keeps its 3D point on the boundary's
    // line and passes.
    let gap = if polished {
        // The polish converged on the end face to the Newton bar; the
        // distance to the boundary's LINEAR continuation is not a reading of
        // the landing (see above) and is left out of the phantom test.
        0.0
    } else {
        raw_gap
    };
    if std::env::var("BREP_DEBUG_BLEND_MARCH").is_ok() {
        let on_support = support.evaluate_extended(s).or_refuse(KernelStage::Refine, "evaluate_extended")?;
        let on_boundary = boundary.curve.evaluate_extended(t).or_refuse(KernelStage::Refine, "evaluate_extended")?;
        let domain = boundary.curve.domain().unwrap_or([f64::NAN, f64::NAN]);
        eprintln!(
            "  support_crossing_uv: boundary {} edge t=[{:.4},{:.4}] curve domain [{:.4},{:.4}] fraction {fraction:.4} extended={extended} polished={polished} s={s:.6} of [{s0:.4},{s1:.4}] t={t:.6} raw_gap={raw_gap:.3e} gap={gap:.3e} on_support=({:.6},{:.6},{:.6}) on_boundary=({:.6},{:.6},{:.6})",
            boundary.id, boundary.t0, boundary.t1, domain[0], domain[1],
            on_support.x, on_support.y, on_support.z, on_boundary.x, on_boundary.y, on_boundary.z
        );
    }
    // Both sides of that subtraction are APPROXIMATIONS of the same rail: the
    // support is a cubic fit through the marched stations, and the crossing
    // parameter came out of the pcurve fit through the SAME stations — two
    // interpolants that agree at every station and drift between them by the
    // fit's own accuracy, which the surface's parameterisation amplifies (a
    // revolution's u is projective, not proportional to angle).  On a 208°
    // arc of a r≈15 revolve that drift is 2e-4, so `1e-6·scale` was demanding
    // more precision of the row than the row carries and a genuinely
    // transverse 43° corner was refused.  The kernel already names the size
    // of that disagreement — see [`support_crossing_tolerances`], which makes
    // exactly this argument for the 3D intersector — so the phantom test is
    // judged against `pcurve_consistency`, the same contract, and never
    // against a literal of its own.  A crossing accepted only past the
    // historical bar is ESCALATED: the caller must validate what it builds.
    if gap > band.max(historical) {
        return Err(KernelRefusal::unsupported(
            KernelStage::Intersect,
            "support_crossing_uv_phantom",
            format!(
                "blend: support curve does not reach the end boundary edge (the parameter-space \
                 crossing is {gap:.3e} off in 3D — the boundary continues the blended edge rather \
                 than crossing its rail)"
            ),
        ));
    }
    Ok((s, t, polished || gap > historical))
}

/// An EXTEND landing polished onto the end face: the row parameter `s` and
/// the end face's `(u, v)` at which `support(s) = end_surface(u, v)`, by
/// Newton from the pcurve landing's `s` and the end face's own uv of the
/// boundary's vertex (the end nearest the landing, `vertex_t`), both read
/// through their extensions.  Returns the polished `(s, t)`, `t` the nearest
/// parameter of the boundary's extended curve to the landing, or `None` when
/// the system does not converge, the landing moves further than the reach it
/// polishes (another crossing of the end surface, not this one), or the end
/// face does not use the boundary edge.
fn polish_crossing_onto_end_face(
    support: &NurbsCurve,
    boundary: &EdgeRecord,
    end_face: &FaceRecord,
    seed_s: f64,
    seed_t: f64,
    vertex_t: f64,
) -> Result<Option<(f64, f64)>, KernelRefusal> {
    let Some(end_coedge) = end_face
        .loops
        .iter()
        .flat_map(|loop_record| &loop_record.coedges)
        .find(|coedge| coedge.edge_id == boundary.id)
    else {
        return Ok(None);
    };
    let seed_uv = edge_uv_on_face(end_coedge, boundary, vertex_t)?;
    let vertex_point = boundary.curve.evaluate(vertex_t).or_refuse(KernelStage::Refine, "evaluate")?;
    let landing = support.evaluate_extended(seed_s).or_refuse(KernelStage::Refine, "evaluate_extended")?;
    let reach = landing.sub(vertex_point).length();
    let scale = landing.length().max(vertex_point.length()).max(1.0);
    let tolerance = 1e-12 * scale;
    let surface = &end_face.surface;
    let mut x = [seed_s, seed_uv[0], seed_uv[1]];
    let mut converged = false;
    for _ in 0..NEWTON_ITERATIONS {
        let a = support.derivatives_extended(x[0], 1).or_refuse(KernelStage::Refine, "derivatives_extended")?;
        let d = surface.derivatives_extended(x[1], x[2], 1).or_refuse(KernelStage::Refine, "derivatives_extended")?;
        let residual = a[0].sub(d[0][0]);
        if residual.length() <= tolerance {
            converged = true;
            break;
        }
        let matrix = [
            [a[1].x, -d[1][0].x, -d[0][1].x],
            [a[1].y, -d[1][0].y, -d[0][1].y],
            [a[1].z, -d[1][0].z, -d[0][1].z],
        ];
        let Ok(delta) = fit::solve_small::<3>(matrix, [residual.x, residual.y, residual.z], 3) else {
            return Ok(None);
        };
        if delta.iter().any(|value| !value.is_finite()) {
            return Ok(None);
        }
        for (value, correction) in x.iter_mut().zip(delta) {
            *value -= correction;
        }
    }
    if !converged {
        let a = support.evaluate_extended(x[0]).or_refuse(KernelStage::Refine, "evaluate_extended")?;
        let d = surface.evaluate_extended(x[1], x[2]).or_refuse(KernelStage::Refine, "evaluate_extended")?;
        if a.sub(d).length() > 1e-9 * scale {
            return Ok(None);
        }
    }
    let polished = support.evaluate_extended(x[0]).or_refuse(KernelStage::Refine, "evaluate_extended")?;
    // The polish corrects the pcurves' extrapolation, which is a fraction of
    // the extension itself; a landing that moves further than its whole
    // reach has found another crossing of the end surface.
    if polished.sub(landing).length() > reach.max(1e-9 * scale) {
        return Ok(None);
    }
    // The boundary parameter: the nearest point of the extended curve.
    let mut t = seed_t;
    for _ in 0..NEWTON_ITERATIONS {
        let b = boundary.curve.derivatives_extended(t, 2).or_refuse(KernelStage::Refine, "derivatives_extended")?;
        let offset = b[0].sub(polished);
        let g = offset.dot(b[1]);
        let dg = b[1].dot(b[1]) + offset.dot(b[2]);
        if !(dg.abs() > 0.0) || !dg.is_finite() {
            break;
        }
        let step = g / dg;
        if !step.is_finite() {
            break;
        }
        t -= step;
        if step.abs() <= 1e-14 * (1.0 + t.abs()) {
            break;
        }
    }
    if !t.is_finite() {
        return Ok(None);
    }
    Ok(Some((x[0], t)))
}

/// Where a fitted support curve crosses a boundary edge (nearest genuine
/// crossing to the given end of the support domain).
pub(in crate::blend) fn support_crossing(
    solid: &BrepSolid,
    support: &NurbsCurve,
    boundary: &EdgeRecord,
    at_start: bool,
) -> Result<(f64, f64, Vec3, bool), KernelRefusal> {
    for tolerance in support_crossing_tolerances(solid) {
        if let Some((s, t, point)) = support_crossing_at(support, boundary, at_start, tolerance)? {
            return Ok((s, t, point, tolerance > 1e-7));
        }
    }
    Err(KernelRefusal::unsupported(
        KernelStage::Intersect,
        "support_crossing_none",
        "blend: support curve does not reach the end boundary edge",
    ))
}

/// One rung of [`support_crossing`]'s ladder.
fn support_crossing_at(
    support: &NurbsCurve,
    boundary: &EdgeRecord,
    at_start: bool,
    tolerance: f64,
) -> Result<Option<(f64, f64, Vec3)>, KernelRefusal> {
    let hits = crate::intersect_curves(support, &boundary.curve, tolerance).or_refuse(KernelStage::Refine, "intersect_curves")?;
    let [s0, s1] = support.domain().or_refuse(KernelStage::Refine, "domain")?;
    let mut best: Option<(f64, f64, Vec3)> = None;
    for hit in hits {
        if hit.t < boundary.t0 - 1e-9 || hit.t > boundary.t1 + 1e-9 {
            continue;
        }
        // The support is a fitted row over [s0, s1]; a hit reported outside
        // it is the intersector extrapolating the row's carrier (a tangent
        // continuation lets a circle's extension meet a line far away) and
        // is not a crossing of the rail.
        let margin = 1e-9 * (1.0 + (s1 - s0).abs());
        if hit.s < s0 - margin || hit.s > s1 + margin {
            continue;
        }
        let rank = if at_start { hit.s - s0 } else { s1 - hit.s };
        if best
            .map(|(existing, ..)| {
                let existing_rank = if at_start {
                    existing - s0
                } else {
                    s1 - existing
                };
                rank < existing_rank
            })
            .unwrap_or(true)
        {
            best = Some((hit.s, hit.t, hit.point));
        }
    }
    Ok(best)
}

/// The tolerance ladder a fitted support row is intersected against, tightest
/// first.
///
/// # Why a ladder and not one number
/// A support row is a FIT through marched stations; the spoke/boundary edge it
/// crosses is the carrier's own EXACT curve.  Both lie on the same carrier by
/// construction, so the crossing exists geometrically — but the fitted row only
/// reproduces the carrier to the fit's accuracy, so in 3D the two curves pass
/// each other with a residual gap of that size instead of meeting exactly.
/// [`crate::intersect_curves`] discards any Newton answer whose 3D gap exceeds
/// the tolerance it was given, so demanding 1e-7 of a fitted row asks for more
/// precision than the row carries and the crossing is reported as ABSENT.
///
/// The kernel already names the size of that disagreement, and names it twice:
/// [`crate::KernelTolerances::intersection_fit`] ("maximum geometric error
/// accepted while fitting SSI curves") and
/// [`crate::KernelTolerances::pcurve_consistency`] ("maximum edge/pcurve/
/// carrier disagreement in a valid in-memory BREP", whose own comment records
/// that *existing fitted NURBS intersections can deviate by a few microns*).
/// Those are the contract; 1e-7 was tighter than any of them.
///
/// The ladder starts at the historical 1e-7 so every crossing that resolves
/// today resolves identically — the looser rungs are reached ONLY where the
/// tight one finds nothing, i.e. only where the caller used to hard-refuse —
/// and it stops at the FIRST rung that produces a crossing, so a crossing is
/// never accepted at a looser tolerance than it actually needs.
pub(in crate::blend) fn support_crossing_tolerances(solid: &BrepSolid) -> [f64; 3] {
    let policy = crate::KernelTolerances::for_solid(solid, 1e-7);
    [1e-7, policy.intersection_fit, policy.pcurve_consistency]
}

/// Where a fitted support row crosses a chain SPOKE edge.
///
/// Returns the crossing nearest the chain vertex as `(support parameter, spoke
/// edge parameter, escalated)`, where `escalated` is true when the crossing
/// only resolved on a looser rung than the historical 1e-7 — the caller MUST
/// then validate what the surgery built (see [`SpokeCrossings`]).  `None` means
/// the row genuinely does not cross the spoke on any rung; the caller renders
/// that refusal, because what it means differs between the closed and open
/// chain lanes.
pub(in crate::blend) fn spoke_crossing(
    solid: &BrepSolid,
    support: &NurbsCurve,
    spoke: &EdgeRecord,
    vertex_point: Vec3,
) -> Result<Option<(f64, f64, bool)>, KernelRefusal> {
    for tolerance in support_crossing_tolerances(solid) {
        let hits = crate::intersect_curves(support, &spoke.curve, tolerance).or_refuse(KernelStage::Refine, "intersect_curves")?;
        let mut best: Option<(f64, f64, f64)> = None;
        for hit in hits {
            if hit.t < spoke.t0 - 1e-9 || hit.t > spoke.t1 + 1e-9 {
                continue;
            }
            let distance = hit.point.sub(vertex_point).length();
            if best.map(|(.., known)| distance < known).unwrap_or(true) {
                best = Some((hit.s, hit.t, distance));
            }
        }
        if let Some((crossing_support, crossing_spoke, _)) = best {
            return Ok(Some((crossing_support, crossing_spoke, tolerance > 1e-7)));
        }
    }
    Ok(None)
}

/// Whether any spoke crossing in one chain surgery needed a looser rung than
/// the historical 1e-7, and the gate that fact imposes on the result.
///
/// # Why the escalation is gated and the tight rung is not
/// The tight rung reproduces master exactly, INCLUDING whatever it already
/// builds — so it is never re-judged here; re-judging it would turn an existing
/// output into a refusal, which is a regression whatever the output's quality.
/// The looser rungs, by contrast, only ever run where master hard-refused, so
/// everything they produce is NEW.  A new output is worth having only if it is
/// a valid solid: a refusal is recoverable (the caller's maximal-valid-subset
/// search, or the user's own edit), a silently invalid solid is not.  So the
/// escalated path must earn its answer, and a radius that marches past what the
/// supports can carry goes back to being refused — with a message that says so
/// — instead of shipping a torn blend.
#[derive(Clone, Copy, Default)]
pub(in crate::blend) struct SpokeCrossings {
    escalated: bool,
}

impl SpokeCrossings {
    pub(in crate::blend) fn note(&mut self, escalated: bool) {
        self.escalated |= escalated;
    }

    /// Accept `result` unless it was built on an escalated crossing AND does
    /// not validate.
    pub(in crate::blend) fn gate(&self, result: BrepSolid) -> Result<BrepSolid, KernelRefusal> {
        if !self.escalated {
            return Ok(result);
        }
        let issues = result.validate();
        if issues.is_empty() {
            return Ok(result);
        }
        Err(KernelRefusal::new(
            RefusalClass::InvalidResultTopology {
                issues: issues.len() as u32,
            },
            KernelStage::Validate,
            format!(
                "blend: chain surgery across a fit-tolerance spoke crossing did not \
                 produce a valid solid ({})",
                issues
                    .first()
                    .map(|issue| issue.message.clone())
                    .unwrap_or_default()
            ),
        ))
    }
}

/// The transverse curve where the blend surface crosses the end face:
/// per-z Newton in t on the signed distance to the end carrier, endpoints
/// pinned to the rim crossings.
pub(in crate::blend) fn transverse_curve(
    rows: &FittedRows,
    end_face: &FaceRecord,
    cr_parameter: f64,
    cs_parameter: f64,
) -> Result<(NurbsCurve, NurbsCurve, NurbsCurve), KernelRefusal> {
    transverse_curve_with_ends(rows, end_face, cr_parameter, cs_parameter, [false, false])
}

/// [`transverse_curve`], with the z = 0 (cr) and z = 1 (cs) ends EXACT where
/// `exact_ends` says so: the 3D curve's end derivative along the wall ∩ end
/// face direction (`n_wall × n_end`, which is the mate ∩ end-face boundary's
/// direction — the wall is G1 with its mate along the rail), the end face's
/// feet on its EXTENSION rather than its clamped domain, and the end pcurve's
/// end derivative that direction pulled back through the extension's
/// Jacobian.  An EXTEND landing needs both (see the caller).
pub(in crate::blend) fn transverse_curve_with_ends(
    rows: &FittedRows,
    end_face: &FaceRecord,
    cr_parameter: f64,
    cs_parameter: f64,
    exact_ends: [bool; 2],
) -> Result<(NurbsCurve, NurbsCurve, NurbsCurve), KernelRefusal> {
    // A LADDER on the section count (lane R2, 2026-10-18). The edge and its
    // two pcurves are fitted through the same samples, each on its own: they
    // agree AT the samples and part company between them where the section's
    // parameter on the wall bends hard (smFilletTest lCount_224's free end at
    // v6, whose rail meets its boundary on that boundary's extension: the
    // blend pcurve's image stood 1.449e-3 off the edge between its first two
    // samples, and the shell's closure read 0.0377 against a bar of 0.0036).
    // The 16-section build is kept, bit for bit, whenever both pcurves'
    // images stand within the pcurve floor of the edge at four points per
    // section; otherwise the sections double, at most to 256, and the closest
    // build is kept.
    const FIRST_SECTIONS: usize = 16;
    const MAX_SECTIONS: usize = 256;
    if !transverse_ladder_on() {
        return transverse_curve_sections(rows, end_face, cr_parameter, cs_parameter, exact_ends, FIRST_SECTIONS);
    }
    let mut sections = FIRST_SECTIONS;
    let mut best: Option<((NurbsCurve, NurbsCurve, NurbsCurve), f64)> = None;
    loop {
        let built = transverse_curve_sections(rows, end_face, cr_parameter, cs_parameter, exact_ends, sections)?;
        let miss = transverse_pcurve_miss(rows, &built, 4 * sections)?;
        if std::env::var("BREP_DEBUG_TRANSVERSE_FIT").is_ok() {
            let (sliver, length) = transverse_sliver(rows, &built, 4 * sections);
            eprintln!(
                "transverse ladder: end face {:?}, exact ends {exact_ends:?}, {sections} sections, pcurve images {miss:.3e} off the edge; \
                 sliver vector area {sliver:.3e} over length {length:.3e} ({:.3e} per length)",
                end_face.name,
                sliver / length
            );
        }
        if sections > FIRST_SECTIONS {
            crate::blend::record_count("count.transverse.refined_rungs", 1);
        }
        let done = miss <= crate::pcurve::PCURVE_REFINEMENT_TOLERANCE || sections >= MAX_SECTIONS;
        if best.as_ref().map_or(true, |(_, known)| miss < *known) {
            best = Some((built, miss));
        }
        if done {
            break;
        }
        sections *= 2;
    }
    Ok(best.expect("one build at least").0)
}

/// `BREP_BLEND_TRANSVERSE_LADDER=1`/`on` climbs [`transverse_curve_with_ends`]'s
/// section ladder, `=0`/`off` (the default) builds 16 sections as before;
/// `set_transverse_ladder_for_test` forces it on a test's own thread.
const TRANSVERSE_LADDER_DEFAULT: bool = false;

thread_local! {
    static TRANSVERSE_LADDER_FOR_TEST: std::cell::Cell<Option<bool>> = const { std::cell::Cell::new(None) };
}

pub(in crate::blend) fn transverse_ladder_on() -> bool {
    if let Some(forced) = TRANSVERSE_LADDER_FOR_TEST.with(|forced| forced.get()) {
        return forced;
    }
    match std::env::var("BREP_BLEND_TRANSVERSE_LADDER").as_deref() {
        Ok("0") | Ok("off") => false,
        Ok("1") | Ok("on") => true,
        _ => TRANSVERSE_LADDER_DEFAULT,
    }
}


/// How far the transverse section's BLEND pcurve's image stands off its
/// edge: the worst `|C(z) − S_wall(b(z))|` over `samples` midpoints of
/// `z ∈ [0, 1]`, the wall read on its extension (the overshoot is where the
/// section lives). The end face's pcurve is not read here: off an exact
/// landing its samples are the clamped feet, which the planar-chart pass
/// rebuilds on the widened chart (`fit_planar_charts_to_trims`). An
/// unreadable point reads infinitely far.
fn transverse_pcurve_miss(
    rows: &FittedRows,
    built: &(NurbsCurve, NurbsCurve, NurbsCurve),
    samples: usize,
) -> Result<f64, KernelRefusal> {
    let (curve, blend_pcurve, _) = built;
    let mut worst = 0.0f64;
    for k in 0..samples {
        let z = (k as f64 + 0.5) / samples as f64;
        let read = (|| -> Result<f64, String> {
            let point = curve.evaluate(z)?;
            let on_wall = blend_pcurve.evaluate(z)?;
            Ok(point.sub(rows.surface.evaluate_extended(on_wall.x, on_wall.y)?).length())
        })();
        match read {
            Ok(miss) if miss.is_finite() => worst = worst.max(miss),
            _ => return Ok(f64::INFINITY),
        }
    }
    Ok(worst)
}

/// BREP_DEBUG_TRANSVERSE_FIT only: the vector area of the sliver between the
/// section edge and its blend pcurve's image, `|∫ (image − edge) × dedge|`,
/// and the edge's length — the shell closure's per-edge share, which its bar
/// reads per length.
fn transverse_sliver(rows: &FittedRows, built: &(NurbsCurve, NurbsCurve, NurbsCurve), samples: usize) -> (f64, f64) {
    let (curve, blend_pcurve, _) = built;
    let (mut area, mut length) = (Vec3::new(0.0, 0.0, 0.0), 0.0);
    let at = |z: f64| -> Option<(Vec3, Vec3)> {
        let on_wall = blend_pcurve.evaluate(z).ok()?;
        Some((curve.evaluate(z).ok()?, rows.surface.evaluate_extended(on_wall.x, on_wall.y).ok()?))
    };
    for k in 0..samples {
        let (z0, z1) = (k as f64 / samples as f64, (k + 1) as f64 / samples as f64);
        let (Some((p0, i0)), Some((p1, i1))) = (at(z0), at(z1)) else { return (f64::NAN, f64::NAN) };
        let gap = i0.sub(p0).add(i1.sub(p1)).scale(0.5);
        area = area.add(gap.cross(p1.sub(p0)));
        length += p1.sub(p0).length();
    }
    (area.length(), length)
}

/// One build of [`transverse_curve_with_ends`] on `sections` sections.
fn transverse_curve_sections(
    rows: &FittedRows,
    end_face: &FaceRecord,
    cr_parameter: f64,
    cs_parameter: f64,
    exact_ends: [bool; 2],
    sections: usize,
) -> Result<(NurbsCurve, NurbsCurve, NurbsCurve), KernelRefusal> {
    let exact = exact_ends[0] || exact_ends[1];
    let mut points3 = Vec::with_capacity(sections + 1);
    let mut blend_uv = Vec::with_capacity(sections + 1);
    let mut end_uv: Vec<[f64; 2]> = Vec::with_capacity(sections + 1);
    let mut parameters = Vec::with_capacity(sections + 1);
    // The curve is ONE branch of the end carrier: after the first section every
    // foot continues from the last one found — the previous section's, the
    // previous iterate's, and for the difference quotient the iterate's own —
    // so the level the Newton differences is the distance to that sheet, and no
    // sample of the end pcurve is the nearest point of another.
    let project = |point: Vec3, seed: Option<[f64; 2]>| match seed {
        None => crate::project_point_to_surface(&end_face.surface, point)
            .or_refuse(KernelStage::Refine, "project_point_to_surface"),
        Some(seed) => crate::projection::project_point_to_surface_from_seed(&end_face.surface, point, seed)
            .or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed"),
    };
    for section in 0..=sections {
        let z = section as f64 / sections as f64;
        let mut t = cr_parameter + (cs_parameter - cr_parameter) * z;
        let mut foot = end_uv.last().copied();
        if section > 0 && section < sections {
            for _ in 0..NEWTON_ITERATIONS {
                let point = rows.surface.evaluate_extended(t, z).or_refuse(KernelStage::Refine, "evaluate_extended")?;
                let projection = project(point, foot)?;
                foot = Some([projection.u, projection.v]);
                let normal = raw_normal(&end_face.surface, projection.u, projection.v)?;
                let distance = point.sub(projection.point).dot(normal);
                if distance.abs() <= 1e-11 {
                    break;
                }
                let step = 1e-7;
                let probe = rows.surface.evaluate_extended(t + step, z).or_refuse(KernelStage::Refine, "evaluate_extended")?;
                let probe_projection = project(probe, foot)?;
                let probe_distance = probe.sub(probe_projection.point).dot(raw_normal(
                    &end_face.surface,
                    probe_projection.u,
                    probe_projection.v,
                )?);
                let derivative = (probe_distance - distance) / step;
                if derivative.abs() <= 1e-14 {
                    return Err(KernelRefusal::non_convergence(
                        KernelStage::Refine,
                        "transverse_newton",
                        "blend: transverse Newton stalled",
                    ));
                }
                t -= distance / derivative;
            }
        }
        let point = rows.surface.evaluate_extended(t, z).or_refuse(KernelStage::Refine, "evaluate_extended")?;
        let projection = project(point, foot)?;
        points3.push(Vec4::from_point(point, 1.0));
        blend_uv.push(Vec4::from_point(Vec3::new(t, z, 0.0), 1.0));
        if exact {
            // From the branch-continuing foot: only an excursion past the
            // domain moves it.
            let (u, v, _) = crate::blend::planar_chart::extended_foot_from(&end_face.surface, point, [projection.u, projection.v])?;
            if std::env::var("BREP_DEBUG_TRANSVERSE_FIT").is_ok() {
                let du = end_face.surface.domain_u().unwrap_or([f64::NAN; 2]);
                let dv = end_face.surface.domain_v().unwrap_or([f64::NAN; 2]);
                eprintln!(
                    "transverse foot z {z:.4}: seeded ({:.6}, {:.6}) {:.3e} off; extension ({u:.6}, {v:.6}); end face {:?} kind {} domain u {du:?} v {dv:?}",
                    projection.u, projection.v, projection.distance, end_face.name,
                    end_face.surface.analytic().map(|_| "analytic").unwrap_or("nurbs")
                );
            }
            end_uv.push([u, v]);
        } else {
            end_uv.push([projection.u, projection.v]);
        }
        parameters.push(z);
    }
    unwrap_seam_branch(&end_face.surface, &mut end_uv)?;
    if exact {
        // The exact ends: the end derivative's DIRECTION is the wall ∩ end
        // face direction (n_wall × n_end), its MAGNITUDE the free fit's end
        // derivative along it; a non-exact end keeps the free fit's.  The
        // end pcurve takes the same derivative pulled back through the
        // extension's Jacobian.  Interpolated through every sample.
        let free3 = fit::interpolate_homogeneous(&points3, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
        let end_uv4: Vec<Vec4> = end_uv.iter().map(|uv| Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)).collect();
        let free_uv = fit::interpolate_homogeneous(&end_uv4, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
        let free_blend = fit::interpolate_homogeneous(&blend_uv, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
        let mut derivatives3 = [Vec3::default(); 2];
        let mut derivatives_uv = [Vec3::default(); 2];
        for (side, index) in [(0usize, 0usize), (1, sections)] {
            let z = parameters[index];
            let free_d3 = free3.derivatives(z, 1).or_refuse(KernelStage::Refine, "derivatives")?[1];
            let free_duv = free_uv.derivatives(z, 1).or_refuse(KernelStage::Refine, "derivatives")?[1];
            if !exact_ends[side] {
                derivatives3[side] = free_d3;
                derivatives_uv[side] = free_duv;
                continue;
            }
            let wall = raw_normal(&rows.surface, blend_uv[index].x, z)?;
            let [u, v] = end_uv[index];
            let d = end_face.surface.derivatives_extended(u, v, 1).or_refuse(KernelStage::Refine, "derivatives_extended")?;
            let (su, sv) = (d[1][0], d[0][1]);
            let end_normal = su.cross(sv).normalized().or_refuse(KernelStage::Refine, "normalized")?;
            let mut direction = wall.cross(end_normal).normalized().or_refuse(KernelStage::Refine, "normalized")?;
            if direction.dot(free_d3) < 0.0 {
                direction = direction.scale(-1.0);
            }
            let tangent = direction.scale(free_d3.dot(direction));
            let (a11, a12, a22) = (su.dot(su), su.dot(sv), sv.dot(sv));
            let (b1, b2) = (su.dot(tangent), sv.dot(tangent));
            let det = a11 * a22 - a12 * a12;
            if !(det.abs() > 1e-300) {
                return Err(KernelRefusal::non_convergence(KernelStage::Refine, "transverse_exact_end", "blend: the end face's extension is singular at a transverse end"));
            }
            derivatives3[side] = tangent;
            derivatives_uv[side] = Vec3::new((b1 * a22 - b2 * a12) / det, (a11 * b2 - a12 * b1) / det, 0.0);
            if std::env::var("BREP_DEBUG_RULED_CHART").is_ok() {
                let tilt = free_d3.normalized().map(|free| (1.0 - free.dot(direction)).max(0.0)).unwrap_or(f64::NAN);
                eprintln!("transverse exact end {side}: the free fit's end tangent stood 1 - cos = {tilt:.3e} off the wall ∩ end-face direction");
            }
        }
        let points: Vec<Vec3> = points3.iter().map(|p| p.point()).collect::<Result<_, _>>().or_refuse(KernelStage::Refine, "point")?;
        let uv_points: Vec<Vec3> = end_uv.iter().map(|uv| Vec3::new(uv[0], uv[1], 0.0)).collect();
        // All three on ONE knot vector, as the free fits are: the edge and
        // its two pcurves agree between the samples only if they share it
        // (with the blend pcurve left on the free fit's knots).
        let blend_points: Vec<Vec3> = blend_uv.iter().map(|p| p.point()).collect::<Result<_, _>>().or_refuse(KernelStage::Refine, "point")?;
        let blend_derivatives = [
            free_blend.derivatives(parameters[0], 1).or_refuse(KernelStage::Refine, "derivatives")?[1],
            free_blend.derivatives(parameters[sections], 1).or_refuse(KernelStage::Refine, "derivatives")?[1],
        ];
        if let Ok(path) = std::env::var("BREP_DUMP_TRANSVERSE_ENDS") {
            // The samples and end derivatives as JSON lines: the regression
            // case for `fit::interpolate_curve_with_end_tangents`'s knot fix.
            let row = |v: &Vec3| format!("[{:.17e},{:.17e},{:.17e}]", v.x, v.y, v.z);
            let line = format!(
                "{{\"points\":[{}],\"parameters\":[{}],\"start_derivative\":{},\"end_derivative\":{}}}\n",
                points.iter().map(row).collect::<Vec<_>>().join(","),
                parameters.iter().map(|t| format!("{t:.17e}")).collect::<Vec<_>>().join(","),
                row(&derivatives3[0]),
                row(&derivatives3[1]),
            );
            use std::io::Write as _;
            if let Ok(mut file) = std::fs::OpenOptions::new().create(true).append(true).open(path) {
                let _ = file.write_all(line.as_bytes());
            }
        }
        // Lane Y's knot fix (2026-10-15) made the shared fitter the Piegl–Tiller
        // end-derivative interpolant this lane had duplicated; lane Y's
        // differential matched every captured input exactly (knots, controls,
        // positions, derivatives).
        let refuse = |message: String| KernelRefusal::non_convergence(KernelStage::Refine, "transverse_exact_end", format!("blend: {message}"));
        let shared = |points: &[Vec3], derivatives: [Vec3; 2]| {
            crate::interpolate_curve_with_end_tangents(points, &parameters, derivatives[0], derivatives[1]).map_err(refuse)
        };
        let curve = shared(&points, derivatives3)?;
        let blend_pcurve = shared(&blend_points, blend_derivatives)?;
        let end_pcurve = shared(&uv_points, derivatives_uv)?;
        if std::env::var("BREP_DEBUG_TRANSVERSE_FIT").is_ok() {
            // Between the samples, how far each fit stands off the two
            // surfaces it is the intersection of (the end face on its
            // extension, the wall by projection).
            let off = |fitted: &NurbsCurve| -> Result<(f64, f64), KernelRefusal> {
                let (mut end_off, mut wall_off) = (0.0_f64, 0.0_f64);
                for k in 0..sections * 8 {
                    let z = (k as f64 + 0.5) / (sections * 8) as f64;
                    let point = fitted.evaluate(z).or_refuse(KernelStage::Refine, "evaluate")?;
                    end_off = end_off.max(crate::blend::planar_chart::extended_foot(&end_face.surface, point)?.2);
                    wall_off = wall_off.max(crate::project_point_to_surface(&rows.surface, point).or_refuse(KernelStage::Refine, "project")?.distance);
                }
                Ok((end_off, wall_off))
            };
            let (free_end, free_wall) = off(&free3)?;
            let (exact_end, exact_wall) = off(&curve)?;
            eprintln!("transverse fit {exact_ends:?}: free fit off end face {free_end:.3e}, off wall {free_wall:.3e}; exact-end fit off end face {exact_end:.3e}, off wall {exact_wall:.3e}");
        }
        return Ok((curve, blend_pcurve, end_pcurve));
    }
    let end_uv: Vec<Vec4> = end_uv
        .into_iter()
        .map(|uv| Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0))
        .collect();
    let curve = fit::interpolate_homogeneous(&points3, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    let blend_pcurve = fit::interpolate_homogeneous(&blend_uv, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    let end_pcurve = fit::interpolate_homogeneous(&end_uv, FIT_DEGREE, &parameters).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
    Ok((curve, blend_pcurve, end_pcurve))
}

/// Make a sampled parameter sequence CONTINUOUS across a closed carrier's
/// seam, in place.
///
/// # Why a projection alone cannot decide this
/// On a surface closed in u, the parameters `u0` and `u1` name the SAME
/// meridian, so a point sitting on it projects to either, and the projector's
/// answer is arbitrary — measured on the captured revolve-about-its-own-edge
/// part, the rim point where the blend's row meets its mate came back at
/// `u = 1.000000` while every neighbouring sample answered near `u = 0`.  That
/// happens at every free end whose mate boundary IS the end face's seam, which
/// a revolve about a line through the profile produces by construction (the
/// shared generatrix is both).  One sample left on the far branch makes the
/// fitted pcurve walk the long way round the carrier.
///
/// The seam is a parameter artefact, so the fix is the parametric one: choose,
/// for each sample, the branch nearest its neighbour.  The walk starts at the
/// MIDDLE sample and runs outward in both directions — an end sample is the
/// one that lands on a seam, the middle is the one that does not.  Whether the
/// branch the WHOLE curve lands on is the one the end face's loop can join is
/// a separate question, and [`check_loop_branches`] is where it is asked.
fn unwrap_seam_branch(surface: &NurbsSurface, samples: &mut [[f64; 2]]) -> Result<(), KernelRefusal> {
    let (closed_u, closed_v) = surface.closed_directions().unwrap_or((false, false));
    if (!closed_u && !closed_v) || samples.len() < 2 {
        return Ok(());
    }
    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
    let periods = [
        if closed_u { u1 - u0 } else { 0.0 },
        if closed_v { v1 - v0 } else { 0.0 },
    ];
    let snap = |value: f64, reference: f64, period: f64| -> f64 {
        if !(period > 0.0) {
            return value;
        }
        let shifted = value - period * ((value - reference) / period).round();
        if (shifted - reference).abs() < (value - reference).abs() {
            shifted
        } else {
            value
        }
    };
    let middle = samples.len() / 2;
    for index in (0..middle).rev() {
        for axis in 0..2 {
            samples[index][axis] = snap(samples[index][axis], samples[index + 1][axis], periods[axis]);
        }
    }
    for index in middle + 1..samples.len() {
        for axis in 0..2 {
            samples[index][axis] = snap(samples[index][axis], samples[index - 1][axis], periods[axis]);
        }
    }
    Ok(())
}

/// Locate a trim's physical point on its own carrier. Reconstructed boundary
/// curves may use a native carrier parameter while the inherited pcurve keeps
/// its original parameterization; equal traversal fractions are then seeds,
/// not a correspondence. Restrict the stored pcurve, without refitting it.
pub(in crate::blend) fn trim_parameter_at_point(
    surface: &crate::NurbsSurface,
    pcurve: &NurbsCurve,
    seed: f64,
    target: Vec3,
    band: f64,
) -> Result<f64, KernelRefusal> {
    let [lo, hi] = pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
    let at = |q: f64| -> Result<(Vec3, Vec3), KernelRefusal> {
        let d = pcurve.derivatives(q, 1).or_refuse(KernelStage::Refine, "derivatives")?;
        let (point, su, sv) = surface.deriv1_extended(d[0].x, d[0].y)
            .or_refuse(KernelStage::Refine, "deriv1_extended")?;
        Ok((point, su.scale(d[1].x).add(sv.scale(d[1].y))))
    };
    if !band.is_finite() || band <= 0.0 || !seed.is_finite() {
        return Err(KernelRefusal::unsupported(KernelStage::Sew, "trim_pcurve_crossing",
            "blend: boundary trim requires a finite seed and positive finite consistency bar"));
    }
    // Scale stopping precision by local travel, never by the world origin.
    // Every success still passes the existing acceptance bar below.
    let (_, initial_tangent) = at(seed)?;
    let precision = (1e-11 * (1.0 + initial_tangent.length() * (hi - lo))).min(band * 0.01);
    let mut seeds = vec![seed];
    let mut breaks = vec![lo, hi];
    breaks.extend(pcurve.knots.iter().copied().filter(|q| *q > lo && *q < hi));
    breaks.sort_by(f64::total_cmp);
    breaks.dedup();
    for span in breaks.windows(2) {
        for i in 0..=4 { seeds.push(span[0] + (span[1] - span[0]) * i as f64 / 4.0); }
    }
    let mut best = (f64::INFINITY, seed);
    let mut visits = Vec::new();
    for (index, mut q) in seeds.into_iter().enumerate() {
        let mut local_best = (f64::INFINITY, q);
        for _ in 0..32 {
            let (point, tangent) = at(q)?;
            let difference = point.sub(target);
            let error = difference.length();
            if error < local_best.0 { local_best = (error, q); }
            if error <= precision { break; }
            let speed2 = tangent.dot(tangent);
            if speed2 <= 1e-300 { break; }
            let step = difference.dot(tangent) / speed2;
            let mut next = q;
            for j in 0..16 {
                let probe = (q - step * 0.5_f64.powi(j)).clamp(lo, hi);
                if at(probe)?.0.sub(target).length() < error { next = probe; break; }
            }
            if next == q { break; }
            q = next;
        }
        // Only the traversal-seeded solve can select the split parameter.
        // Other seeds discover ambiguity; they never substitute another visit.
        if index == 0 { best = local_best; }
        if local_best.0.is_finite() && local_best.0 <= band {
            visits.push(local_best.1);
        }
    }
    // No global fallback may replace a failed traversal-seeded crossing.
    if !best.0.is_finite() || best.0 > band {
        return Err(KernelRefusal::unsupported(KernelStage::Sew, "trim_pcurve_crossing",
            format!("blend: boundary trim cannot reach its new vertex on its carrier (miss {:.3e}, bar {band:.3e})", best.0)));
    }
    // A lower residual cannot choose between separate visits to the same
    // physical point. Reject discovered overlap instead of cutting a different
    // traversal interval. Nearby iterates of one crossing stay inside the bar.
    for q in visits {
        let a = q.min(best.1);
        let b = q.max(best.1);
        let mut intervals = vec![a];
        intervals.extend(breaks.iter().copied().filter(|t| *t > a && *t < b));
        intervals.push(b);
        for span in intervals.windows(2) {
            for i in 1..4 {
                let between = span[0] + (span[1] - span[0]) * i as f64 / 4.0;
                let miss = at(between)?.0.sub(target).length();
                if !miss.is_finite() || miss > band {
                    return Err(KernelRefusal::unsupported(KernelStage::Sew, "trim_pcurve_ambiguous",
                        "blend: boundary trim has separate physical crossing visits; traversal correspondence is ambiguous"));
                }
            }
        }
    }
    Ok(best.1)
}

/// Trim `edge_id` at `parameter`, moving the `old_vertex` endpoint to
/// `new_vertex`, and split every referencing coedge pcurve to the
/// same physical point on its carrier. The affine traversal fraction is
/// only a seed when reconstruction has changed the edge parameterization.
pub(in crate::blend) fn trim_edge_at(
    solid: &mut BrepSolid,
    edge_id: u64,
    parameter: f64,
    old_vertex: u64,
    new_vertex: u64,
) -> Result<(), KernelRefusal> {
    let (old_t0, old_t1, trims_start) = {
        let edge = solid
            .edges
            .iter()
            .find(|edge| edge.id == edge_id)
            .ok_or_else(|| {
                KernelRefusal::internal(
                    KernelStage::Sew,
                    "trim_edge",
                    "blend: edge to trim is missing",
                )
            })?;
        (edge.t0, edge.t1, edge.start_vertex_id == old_vertex)
    };
    let target = solid.vertices.iter().find(|vertex| vertex.id == new_vertex)
        .ok_or_else(|| KernelRefusal::internal(KernelStage::Sew, "trim_vertex", "blend: new trim vertex is missing"))?.point;
    let pcurve_band = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
    let fraction = (parameter - old_t0) / (old_t1 - old_t0);
    let interior = (1e-9..=1.0 - 1e-9).contains(&fraction);
    let extends = if trims_start {
        fraction < 1e-9
    } else {
        fraction > 1.0 - 1e-9
    };
    if !interior && !extends {
        return Err(KernelRefusal::unsupported(
            KernelStage::Sew,
            "trim_degenerate",
            "blend: end trim degenerates the boundary edge",
        ));
    }
    if extends {
        let edge = solid.edges.iter().find(|edge| edge.id == edge_id).unwrap();
        let linear = edge.curve.degree == 1
            && edge.curve.control_points.len() == 2
            && (edge.curve.control_points[0].w - edge.curve.control_points[1].w).abs() <= 1e-12;
        let curved = !linear || solid
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .flat_map(|face| &face.loops)
            .flat_map(|loop_record| &loop_record.coedges)
            .any(|coedge| coedge.edge_id == edge_id && coedge.pcurve.degree != 1);
        if curved {
            return extend_edge_on_carriers(solid, edge_id, parameter, old_vertex, new_vertex);
        }
    }
    for shell in &mut solid.shells {
        for face in &mut shell.faces {
            for loop_record in &mut face.loops {
                for coedge in &mut loop_record.coedges {
                    if coedge.edge_id != edge_id {
                        continue;
                    }
                    let [q0, q1] = coedge.pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
                    // Pcurve fraction runs with the traversal: forward
                    // coedges map fraction f -> q0 + (q1-q0)f, reversed
                    // ones map f -> q1 - (q1-q0)f.
                    let split_q = if coedge.forward {
                        q0 + (q1 - q0) * fraction
                    } else {
                        q1 - (q1 - q0) * fraction
                    };
                    if interior {
                        let split_q = trim_parameter_at_point(&face.surface, &coedge.pcurve, split_q, target, pcurve_band)?;
                        let keep_upper = trims_start == coedge.forward;
                        let (low, high) = coedge.pcurve.split(split_q).or_refuse(KernelStage::Refine, "split")?;
                        coedge.pcurve = if keep_upper { high } else { low };
                    } else {
                        // EXTEND (§6.9): the boundary continues past its old
                        // end to meet the blend.  Exact for linear pcurves;
                        // curved pcurves would need re-fitting.
                        if coedge.pcurve.degree != 1 {
                            return Err(KernelRefusal::unsupported(
                                KernelStage::Sew,
                                "extend_curved_pcurve",
                                "blend: cannot extend a curved boundary pcurve",
                            ));
                        }
                        let a = coedge.pcurve.evaluate(q0).or_refuse(KernelStage::Refine, "evaluate")?;
                        let b = coedge.pcurve.evaluate(q1).or_refuse(KernelStage::Refine, "evaluate")?;
                        let direction = b.sub(a);
                        let extended_fraction = (split_q - q0) / (q1 - q0);
                        let target = a.add(direction.scale(extended_fraction));
                        let controls = &mut coedge.pcurve.control_points;
                        let end_index = if (trims_start == coedge.forward) == false {
                            controls.len() - 1
                        } else {
                            0
                        };
                        let w = controls[end_index].w;
                        controls[end_index] = crate::Vec4::from_point(target, w);
                    }
                }
            }
        }
    }
    let edge = solid
        .edges
        .iter_mut()
        .find(|edge| edge.id == edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(KernelStage::Sew, "trim_edge", "blend: edge to trim is missing")
        })?;
    if extends {
        // Updating the trim alone cannot extend a bounded NURBS line:
        // ordinary evaluation clamps at its original domain. Rebuild the
        // line over the new range, preserving the affine edge parameter.
        let t0 = if trims_start { parameter } else { old_t0 };
        let t1 = if trims_start { old_t1 } else { parameter };
        let a = edge.curve.evaluate_extended(t0).or_refuse(KernelStage::Refine, "extend_line")?;
        let b = edge.curve.evaluate_extended(t1).or_refuse(KernelStage::Refine, "extend_line")?;
        edge.curve = NurbsCurve::new(1, vec![t0, t0, t1, t1],
            vec![Vec4::from_point(a, 1.0), Vec4::from_point(b, 1.0)])
            .or_refuse(KernelStage::Refine, "extend_line")?;
    }
    if trims_start {
        edge.t0 = parameter;
        edge.start_vertex_id = new_vertex;
    } else {
        edge.t1 = parameter;
        edge.end_vertex_id = new_vertex;
    }
    Ok(())
}

/// Stations per original edge span, and per extension, of the refit in
/// [`extend_edge_on_carriers`].
const EXTENSION_REFIT_SAMPLES: usize = 96;
const EXTENSION_SAMPLES: usize = 8;

/// The §6.9 EXTEND on a CURVED boundary: the edge continues past `old_vertex`
/// to `new_vertex` along the intersection of the two faces it bounds, and its
/// curve and both pcurves are refitted over the extended range.
///
/// # Why the intersection, and not the curve's own continuation
/// A boundary is where its two carriers meet, so that is the only curve its
/// continuation can be: a cone's rim runs on round the cone, not along its end
/// tangent (whose departure is `½·d²·κ` — 4e-4 over the 0.08 a concave r = 0.1
/// blend's rail runs past a radius-7.4 rim), and a pcurve extrapolated in
/// parameter space leaves its surface image the same way. So every extension
/// point is SOLVED on both carriers, Newton on the two faces' (u, v) with the
/// point held in a plane across the chord, and the pcurves take the solved
/// parameters; the original span is resampled from the edge as it stands.
///
/// A curved 3D boundary or curved pcurve comes here. An affine line with
/// affine pcurves extends exactly by rebuilding the line's domain and moving
/// the pcurve end controls.
fn extend_edge_on_carriers(
    solid: &mut BrepSolid,
    edge_id: u64,
    parameter: f64,
    old_vertex: u64,
    new_vertex: u64,
) -> Result<(), KernelRefusal> {
    let edge = solid
        .edges
        .iter()
        .find(|edge| edge.id == edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "extend_edge",
                "blend: edge to extend is missing",
            )
        })?
        .clone();
    let trims_start = edge.start_vertex_id == old_vertex;
    let target = solid
        .vertices
        .iter()
        .find(|vertex| vertex.id == new_vertex)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "extend_vertex",
                "blend: the extended boundary's new vertex is missing",
            )
        })?
        .point;
    // The two faces the boundary bounds, with their coedges' pcurves.
    let mut uses: Vec<(u64, CoedgeRecord)> = Vec::new();
    for face in solid.shells.iter().flat_map(|shell| &shell.faces) {
        for coedge in face.loops.iter().flat_map(|loop_record| &loop_record.coedges) {
            if coedge.edge_id == edge_id {
                uses.push((face.id, coedge.clone()));
            }
        }
    }
    let [(face_a, coedge_a), (face_b, coedge_b)] = <[(u64, CoedgeRecord); 2]>::try_from(uses)
        .map_err(|uses| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "extend_uses",
                format!(
                    "blend: cannot extend boundary edge {edge_id} along its carriers — it is used \
                     {} times, not twice",
                    uses.len()
                ),
            )
        })?;
    let surface_of = |face_id: u64| {
        solid
            .shells
            .iter()
            .flat_map(|shell| &shell.faces)
            .find(|face| face.id == face_id)
            .map(|face| face.surface.clone())
            .ok_or_else(|| {
                KernelRefusal::internal(
                    KernelStage::Sew,
                    "extend_face",
                    "blend: a face of the extended boundary is missing",
                )
            })
    };
    let (surface_a, surface_b) = (surface_of(face_a)?, surface_of(face_b)?);
    let (t0, t1) = (edge.t0, edge.t1);
    // A pcurve read at an EDGE parameter of the original range, honouring the
    // affine traversal contract.
    let uv_at = |coedge: &CoedgeRecord, t: f64| -> Result<[f64; 2], KernelRefusal> {
        let point = edge_uv_on_face(coedge, &edge, t)?;
        Ok(point)
    };
    let end_t = if trims_start { t0 } else { t1 };
    let end_point = edge.curve.evaluate(end_t).or_refuse(KernelStage::Refine, "evaluate")?;
    let chord = target.sub(end_point);
    let reach = chord.length();
    let direction = chord
        .normalized()
        .map_err(|_| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "extend_length",
                "blend: the boundary extension has no length",
            )
        })?;
    let scale = target.length().max(end_point.length()).max(1.0);
    let tolerance = 1e-11 * (1.0 + scale);
    // March the carriers' intersection from the old end to the new vertex.
    let mut seed = {
        let a = uv_at(&coedge_a, end_t)?;
        let b = uv_at(&coedge_b, end_t)?;
        [a[0], a[1], b[0], b[1]]
    };
    let mut extension: Vec<(f64, Vec3, [f64; 4])> = Vec::with_capacity(EXTENSION_SAMPLES);
    for step in 1..=EXTENSION_SAMPLES {
        let fraction = step as f64 / EXTENSION_SAMPLES as f64;
        let level = reach * fraction;
        let mut x = seed;
        let mut solved = None;
        for _ in 0..NEWTON_ITERATIONS {
            let da = surface_a.derivatives_extended(x[0], x[1], 1).or_refuse(KernelStage::Refine, "derivatives_extended")?;
            let db = surface_b.derivatives_extended(x[2], x[3], 1).or_refuse(KernelStage::Refine, "derivatives_extended")?;
            let mismatch = da[0][0].sub(db[0][0]);
            let plane = da[0][0].sub(end_point).dot(direction) - level;
            if mismatch.length() <= tolerance && plane.abs() <= tolerance {
                solved = Some(da[0][0]);
                break;
            }
            let matrix = [
                [da[1][0].x, da[0][1].x, -db[1][0].x, -db[0][1].x],
                [da[1][0].y, da[0][1].y, -db[1][0].y, -db[0][1].y],
                [da[1][0].z, da[0][1].z, -db[1][0].z, -db[0][1].z],
                [da[1][0].dot(direction), da[0][1].dot(direction), 0.0, 0.0],
            ];
            let delta = fit::solve_small::<4>(
                matrix,
                [mismatch.x, mismatch.y, mismatch.z, plane],
                4,
            )
            .map_err(|error| {
                KernelRefusal::ill_posed(
                    KernelStage::Refine,
                    "extend_tangent_carriers",
                    format!(
                        "blend: the carriers of boundary edge {edge_id} do not cross transversally \
                         where it must extend ({error})"
                    ),
                )
            })?;
            for (value, correction) in x.iter_mut().zip(delta) {
                *value -= correction;
            }
        }
        let point = solved.ok_or_else(|| {
            KernelRefusal::non_convergence(
                KernelStage::Refine,
                "extend_march",
                format!(
                    "blend: the carriers' intersection could not be followed along boundary edge \
                     {edge_id} to the blend"
                ),
            )
        })?;
        let t = end_t + (parameter - end_t) * fraction;
        extension.push((t, point, x));
        seed = x;
    }
    // The last solved point is the carriers' own crossing of the plane through
    // the new vertex, which sits on that crossing to the rail's fit.
    let arrival = extension
        .last()
        .map(|(_, point, _)| point.sub(target).length())
        .unwrap_or(f64::INFINITY);
    let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
    if arrival > bar {
        return Err(KernelRefusal::internal(
            KernelStage::Sew,
            "extend_arrival",
            format!(
                "blend: boundary edge {edge_id} extended along its carriers arrives {arrival:.3e} \
                 from the blend's rim"
            ),
        ));
    }
    // Every sample over the new range, in increasing edge parameter.
    let mut samples: Vec<(f64, Vec3, [f64; 2], [f64; 2])> =
        Vec::with_capacity(EXTENSION_REFIT_SAMPLES + EXTENSION_SAMPLES + 1);
    for index in 0..=EXTENSION_REFIT_SAMPLES {
        let t = t0 + (t1 - t0) * index as f64 / EXTENSION_REFIT_SAMPLES as f64;
        samples.push((
            t,
            edge.curve.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?,
            uv_at(&coedge_a, t)?,
            uv_at(&coedge_b, t)?,
        ));
    }
    for (t, point, x) in extension {
        samples.push((t, point, [x[0], x[1]], [x[2], x[3]]));
    }
    samples.sort_by(|a, b| a.0.total_cmp(&b.0));
    let (new_t0, new_t1) = if trims_start {
        (parameter, t1)
    } else {
        (t0, parameter)
    };
    let parameters: Vec<f64> = samples.iter().map(|sample| sample.0).collect();
    if parameters.windows(2).any(|pair| pair[1] <= pair[0]) {
        return Err(KernelRefusal::internal(
            KernelStage::Sew,
            "extend_direction",
            format!(
                "blend: boundary edge {edge_id} extends in the wrong direction"
            ),
        ));
    }
    let curve = fit::interpolate_curve(
        &samples.iter().map(|sample| sample.1).collect::<Vec<_>>(),
        3,
        &parameters,
    ).or_refuse(KernelStage::Refine, "interpolate_curve")?;
    // Each pcurve keeps its own domain and the affine contract over the NEW
    // edge range.
    let refit = |coedge: &CoedgeRecord, side: usize| -> Result<NurbsCurve, KernelRefusal> {
        let [q0, q1] = coedge.pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
        let mut pairs: Vec<(f64, Vec3)> = samples
            .iter()
            .map(|sample| {
                let fraction = (sample.0 - new_t0) / (new_t1 - new_t0);
                let q = if coedge.forward {
                    q0 + (q1 - q0) * fraction
                } else {
                    q1 - (q1 - q0) * fraction
                };
                let uv = if side == 0 { sample.2 } else { sample.3 };
                (q, Vec3::new(uv[0], uv[1], 0.0))
            })
            .collect();
        pairs.sort_by(|a, b| a.0.total_cmp(&b.0));
        fit::interpolate_curve(
            &pairs.iter().map(|pair| pair.1).collect::<Vec<_>>(),
            3,
            &pairs.iter().map(|pair| pair.0).collect::<Vec<_>>(),
        )
        .or_refuse(KernelStage::Sew, "interpolate_curve")
    };
    let pcurve_a = refit(&coedge_a, 0)?;
    let pcurve_b = refit(&coedge_b, 1)?;
    for shell in &mut solid.shells {
        for face in &mut shell.faces {
            for loop_record in &mut face.loops {
                for coedge in &mut loop_record.coedges {
                    if coedge.id == coedge_a.id {
                        coedge.pcurve = pcurve_a.clone();
                    } else if coedge.id == coedge_b.id {
                        coedge.pcurve = pcurve_b.clone();
                    }
                }
            }
        }
    }
    let record = solid
        .edges
        .iter_mut()
        .find(|edge| edge.id == edge_id)
        .ok_or_else(|| {
            KernelRefusal::internal(
                KernelStage::Sew,
                "extend_edge",
                "blend: edge to extend is missing",
            )
        })?;
    record.curve = curve;
    record.t0 = new_t0;
    record.t1 = new_t1;
    if trims_start {
        record.start_vertex_id = new_vertex;
    } else {
        record.end_vertex_id = new_vertex;
    }
    Ok(())
}