//! G2 corner patches: the patch that fills a vertex where filleted edges meet, built to repeat each
//! blend's surface up to the second order instead of the rolling ball's sphere.
//!
//! [`crate::blend_corner`] and [`crate::curvecorner`] fill the corner of a vertex with a piece of
//! the ball that rolls along the three edges (a sphere for equal radii, a rational patch for
//! unequal ones). That patch is tangent to every blend along the sections where they stop, so the
//! corner is G1 — but across a seam the ball keeps its own 1/r where a fillet is flat along its
//! axis, so the two surfaces meet with a jump in curvature.
//!
//! Here the corner is one tensor patch whose boundary is the blends' section curves themselves and
//! whose cross-sections carry each blend's own derivatives: the position, the direction the patch
//! leaves the blend's surface in, and the blend's curvature vector along that direction are the
//! first three control points of every cross-section, one control point per continuity condition,
//! exactly as [`crate::bridge`]'s sections do. The other end of every cross-section gathers at
//! `apex` — a point inside the material the corner caps — where the sections arrive in a common
//! tangent plane with a common curvature, so the patch closes there smoothly like a sphere's pole.
//! Along a seam the patch therefore repeats the blend up to the second order (G2).
//!
//! The boundary is a closed loop (the three sections joined corner to corner), so the direction
//! along it is a *closed* cubic B-spline: an interpolating curve with free ends swings wide at the
//! corner where the loop joins and folds the patch there. Where two blends meet at that corner
//! their sections leave in directions 90° apart and ask for different curvature at that one point
//! (a fillet is flat across its own seam, the other is not), which no single surface can hold; the
//! corner is a still point of the construction and [`seam_error`] is measured clear of it. A patch
//! that folds, or a corner with no room, is refused rather than replaced by the G1 ball.

use solidcraft_geom::Vec3;
use solidcraft_geom::bspline::{HPoint, NurbsSurface};
use truck_modeling::{self as mt, builder};

use crate::blend::unsupported;
use crate::body::{from_p3, p3};
use crate::bridge::Frame;
use crate::{KernelError, Result, guard};

/// A corner's name in an error.
fn fail(m: impl std::fmt::Display) -> KernelError {
    KernelError::Failed(format!("G2 corner: {m}"))
}

/// One of the surface's own derivative vectors as ours.
fn vec_of(v: mt::Vector3) -> Vec3 {
    Vec3::new(v.x, v.y, v.z)
}

/// Conditions each end of a cross-section carries: position, tangent and curvature.
const END_CONDITIONS: usize = 3;
/// Degree across the patch: three conditions at each end, less the control point they share.
const ACROSS: usize = 2 * END_CONDITIONS - 1;
/// Degree along the seam loop.
const ALONG: usize = 3;
/// Fewest stations a section curve may be sampled at (a cubic through them).
const MIN_STATIONS: usize = 4;
/// Most stations one section curve may carry.
const MAX_STATIONS: usize = 512;
/// How much of a station's `leave` has to lie across its seam for the patch to have a direction to
/// leave the blend in at all.
const ACROSS_ENOUGH: f64 = 0.2;

/// The G2 corner patch over three blends' section curves: `sides[i]` is one blend's section curve,
/// sampled in order from one corner of the vertex to the next, each side's last station the next
/// side's first (the three close a loop). `apex` is a point inside the material the corner caps;
/// the patch gathers there. The face's boundary edges are the section curves themselves (one edge
/// per side, cut out of the patch's own boundary curve, so they lie on it exactly); the caller
/// re-wraps the wire where it needs its own orientation.
pub fn corner_g2_patch(sides: &[Vec<Frame>], apex: Vec3) -> Result<Vec<mt::Face>> {
    guard("G2 corner patch", || build(sides, apex))
}

/// The seam's own direction at station `i` of a side, from a four-point central difference through
/// the stations (`crate::bridge::frames` samples a chain evenly, and a corner's own callers do the
/// same): a two-point chord would only be second order in the step and would tilt the corner's
/// tangent plane by more than G2 is judged at. `None` within two stations of a side's end.
fn seam_tangent(side: &[Frame], i: usize) -> Option<Vec3> {
    let at = |j: isize| usize::try_from(j).ok().and_then(|j| side.get(j)).map(|f| f.p);
    let (a, b, c, d) = (at(i as isize - 2)?, at(i as isize - 1)?, at(i as isize + 1)?, at(i as isize + 2)?);
    (b * 8.0 - a - c * 8.0 + d).normalized()
}

/// The seam's direction at every station of the closed loop (the same difference, wrapping round
/// its ends), falling back to the chord to a neighbour where the four points are not there.
fn loop_tangents(frames: &[Frame]) -> Vec<Option<Vec3>> {
    let m = frames.len();
    if m == 0 {
        return Vec::new();
    }
    let at = |j: isize| frames.get(j.rem_euclid(m as isize) as usize).map(|f| f.p);
    (0..m)
        .map(|i| {
            let four = || {
                let (a, b, c, d) = (at(i as isize - 2)?, at(i as isize - 1)?, at(i as isize + 1)?, at(i as isize + 2)?);
                (b * 8.0 - a - c * 8.0 + d).normalized()
            };
            four().or_else(|| {
                let prev = at(i as isize - 1)?;
                let next = at(i as isize + 1)?;
                (next - prev).normalized()
            })
        })
        .collect()
}

/// The three section curves as one closed loop of stations, the loop index each side starts at,
/// and the corner's own size (twice the mean distance of a station from their centre).
fn corner_loop(sides: &[Vec<Frame>]) -> Result<(Vec<Frame>, Vec<usize>, f64)> {
    if sides.len() != 3 {
        return Err(unsupported("a G2 corner needs three blends meeting at one vertex (a corner of two, or one, still ends on the ball's patch)"));
    }
    let mut frames: Vec<Frame> = Vec::new();
    let mut starts: Vec<usize> = Vec::new();
    for side in sides {
        if side.len() < MIN_STATIONS || side.len() > MAX_STATIONS {
            return Err(KernelError::Invalid(format!("a corner section curve needs {MIN_STATIONS} to {MAX_STATIONS} sampled stations")));
        }
        for f in side {
            if ![f.p, f.leave, f.bend, f.bend_change].into_iter().all(Vec3::is_finite) || !(f.leave.len() > 0.5) {
                return Err(KernelError::Invalid("a corner station must carry finite data and a unit direction".into()));
            }
        }
        starts.push(frames.len());
        // A side runs from one corner to the next; its last station is the next side's first, so the
        // loop takes it once.
        frames.extend(side.iter().take(side.len() - 1).copied());
    }
    let centre = frames.iter().fold(Vec3::ZERO, |a, f| a + f.p) * (1.0 / frames.len().max(1) as f64);
    let size = frames.iter().map(|f| f.p.dist(centre)).fold(0.0, f64::max) * 2.0;
    if !(size > 0.0 && size.is_finite()) {
        return Err(unsupported("a corner whose section curves have no extent"));
    }
    for (i, side) in sides.iter().enumerate() {
        let Some(prev) = sides.get((i + sides.len() - 1) % sides.len()).and_then(|s| s.last()) else { continue };
        let Some(&here) = side.first() else { continue };
        if prev.p.dist(here.p) > size * 1e-3 {
            return Err(unsupported("the section curves do not meet at a corner"));
        }
    }
    Ok((frames, starts, size))
}

/// Every station's `leave` squared up with the seam it belongs to: the part of it across the seam,
/// which is the direction the cross-section has to leave the blend in. A `leave` along the seam
/// itself leaves the cross-section no room, and is refused.
fn transverse(frames: &[Frame], tangents: &[Option<Vec3>]) -> Result<Vec<Frame>> {
    let mut out = Vec::with_capacity(frames.len());
    for (f, t) in frames.iter().zip(tangents) {
        let across = match t {
            Some(t) => f.leave - *t * f.leave.dot(*t),
            // Stations that offer no direction along the seam: the leave stands as it came.
            None => f.leave,
        };
        if across.len() < ACROSS_ENOUGH * f.leave.len() {
            return Err(unsupported("a section curve meets the corner patch along its own tangent, so the patch has no direction to leave it in"));
        }
        let leave = across.normalized().ok_or_else(|| unsupported("a corner station has no direction across its seam"))?;
        out.push(Frame { leave, ..*f });
    }
    Ok(out)
}

/// One station's cross-section: `ACROSS + 1` control points of a Bézier from the seam point to
/// `apex`. The first three hold the blend: its position, the direction the patch leaves it in
/// (one control point along `leave`) and its curvature vector across the seam (one more). The last
/// three hold the apex: the point, a tangent in the apex's own plane, and the curvature of a dome
/// over the corner, so the sections arrive smoothly instead of leaving a point.
fn section(f: &Frame, apex: Vec3, n_apex: Vec3, kappa: f64) -> Result<[Vec3; ACROSS + 1]> {
    let k = ACROSS as f64;
    let chord = f.p.dist(apex);
    if !(chord.is_finite() && chord > 0.0) {
        return Err(unsupported("a blend too small for the corner patch to fit"));
    }
    let speed = chord / 3.0;
    let to_apex = (f.p - apex).normalized().ok_or_else(|| unsupported("a section curve starts on the apex"))?;
    let leave_apex =
        (to_apex - n_apex * to_apex.dot(n_apex)).normalized().ok_or_else(|| unsupported("a section curve reaches the apex along the apex normal"))?;
    let mut c = [Vec3::ZERO; ACROSS + 1];
    c[0] = f.p;
    c[1] = f.p + f.leave * (speed / k);
    c[2] = f.bend * (speed * speed / (k * (k - 1.0))) + c[1] * 2.0 - c[0];
    c[ACROSS] = apex;
    c[ACROSS - 1] = apex + leave_apex * (speed / k);
    c[ACROSS - 2] = (-n_apex * kappa) * (speed * speed / (k * (k - 1.0))) + c[ACROSS - 1] * 2.0 - c[ACROSS];
    if c.iter().any(|p| !p.is_finite()) {
        return Err(fail("a cross-section of the corner is not finite"));
    }
    Ok(c)
}

/// A closed uniform cubic B-spline through `points`, once round the loop: the curve at each station
/// is `(q[i-1] + 4 q[i] + q[i+1]) / 6`, so the stations fix the net by a cyclic tridiagonal system,
/// solved here with the Sherman–Morrison form of the Thomas algorithm (a corner's loop runs to
/// hundreds of stations). The net comes back in the open form the tensor surface wants — the first
/// three control points repeated past the end — so it still closes, with one uniform span per
/// station and no free end anywhere.
fn closed_layer(points: &[Vec3]) -> Result<Vec<Vec3>> {
    let m = points.len();
    if m < 2 * ALONG + 2 {
        return Err(fail("a corner loop needs at least eight stations"));
    }
    // Sub-diagonal 1, diagonal 4, super-diagonal 1, and 1 in both corners.
    let (sub, diag, sup, corner) = (1.0f64, 4.0f64, 1.0f64, 1.0f64);
    let gamma = -diag;
    let mut bb = vec![diag; m];
    if let Some(first) = bb.first_mut() {
        *first = diag - gamma;
    }
    if let Some(last) = bb.last_mut() {
        *last = diag - corner * corner / gamma;
    }
    // Thomas, one right-hand side at a time (three components solved together).
    let thomas = |rhs: &[Vec3]| -> Result<Vec<Vec3>> {
        let mut c = vec![0.0f64; m];
        let mut d = rhs.to_vec();
        let mut denom = bb.first().copied().unwrap_or(0.0);
        if !(denom.abs() > 1e-300) {
            return Err(fail("the corner's stations are degenerate"));
        }
        let (Some(c0), Some(d0)) = (c.first_mut(), d.first_mut()) else { return Err(fail("stations")) };
        *c0 = sup / denom;
        *d0 = *d0 * (1.0 / denom);
        for i in 1..m {
            let (prev_c, prev_d) = (c.get(i - 1).copied().unwrap_or(0.0), d.get(i - 1).copied().unwrap_or(Vec3::ZERO));
            let b = bb.get(i).copied().unwrap_or(diag);
            denom = b - sub * prev_c;
            if !(denom.abs() > 1e-300) {
                return Err(fail("the corner's stations are degenerate"));
            }
            if let Some(slot) = c.get_mut(i) {
                *slot = if i + 1 < m { sup / denom } else { 0.0 };
            }
            if let Some(slot) = d.get_mut(i) {
                *slot = (*slot - prev_d * sub) * (1.0 / denom);
            }
        }
        let mut x = vec![Vec3::ZERO; m];
        let last = d.last().copied().unwrap_or(Vec3::ZERO);
        if let Some(slot) = x.last_mut() {
            *slot = last;
        }
        for i in (0..m - 1).rev() {
            let (Some(ci), Some(xn)) = (c.get(i).copied(), x.get(i + 1).copied()) else { continue };
            let di = d.get(i).copied().unwrap_or(Vec3::ZERO);
            if let Some(slot) = x.get_mut(i) {
                *slot = di - xn * ci;
            }
        }
        Ok(x)
    };
    let rhs: Vec<Vec3> = points.iter().map(|p| *p * 6.0).collect();
    let x = thomas(&rhs)?;
    // The unit right-hand side the corner terms add: gamma at the first station, the corner at the
    // last.
    let mut unit = vec![Vec3::ZERO; m];
    if let Some(slot) = unit.first_mut() {
        *slot = Vec3::new(gamma, gamma, gamma);
    }
    if let Some(slot) = unit.last_mut() {
        *slot = Vec3::new(corner, corner, corner);
    }
    let z = thomas(&unit)?;
    let (x0, xn) = (x.first().copied().unwrap_or(Vec3::ZERO), x.last().copied().unwrap_or(Vec3::ZERO));
    let (z0, zn) = (z.first().copied().unwrap_or(Vec3::ZERO), z.last().copied().unwrap_or(Vec3::ZERO));
    // z answers a scalar right-hand side, so the correction is one scalar factor per component.
    let factor = |x0: f64, xn: f64, z0: f64, zn: f64| (x0 + corner * xn / gamma) / (1.0 + z0 + corner * zn / gamma);
    let f = Vec3::new(factor(x0.x, xn.x, z0.x, zn.x), factor(x0.y, xn.y, z0.y, zn.y), factor(x0.z, xn.z, z0.z, zn.z));
    let mut out = Vec::with_capacity(m + ALONG);
    for i in 0..m {
        let (xi, zi) = (x.get(i).copied().unwrap_or(Vec3::ZERO), z.get(i).copied().unwrap_or(Vec3::ZERO));
        out.push(xi - Vec3::new(f.x * zi.x, f.y * zi.y, f.z * zi.z));
    }
    if !out.iter().all(|p| p.is_finite()) {
        return Err(fail("the corner's boundary is not finite"));
    }
    for i in 0..ALONG {
        let p = out.get(i).copied().unwrap_or(Vec3::ZERO);
        out.push(p);
    }
    Ok(out)
}

/// The stations of the loop as one tensor control net: across (v) each station's cross-section is
/// the Bézier above, along the seam loop (u) each of its control points is a closed cubic through
/// the stations, so the face passes through every section exactly and its parameter line along the
/// loop closes without an end to swing from. The loop's station count comes back with it: station
/// `i` sits at parameter `i + ALONG`, and the loop closes at `count + ALONG`.
fn tensor(layers: &[[Vec3; ACROSS + 1]]) -> Result<(NurbsSurface, f64)> {
    let m = layers.len();
    let nu = m + ALONG;
    let knots_u: Vec<f64> = (0..=(nu + ALONG) as i64).map(|i| i as f64).collect();
    let mut ctrl: Vec<HPoint> = Vec::with_capacity(nu * (ACROSS + 1));
    for j in 0..=ACROSS {
        let points: Vec<Vec3> = layers.iter().map(|c| c.get(j).copied().unwrap_or(Vec3::ZERO)).collect();
        let net = closed_layer(&points)?;
        for p in &net {
            ctrl.push(HPoint::weighted(*p, 1.0));
        }
    }
    let mut knots_v = vec![0.0; ACROSS + 1];
    knots_v.extend(std::iter::repeat_n(1.0, ACROSS + 1));
    let s = NurbsSurface::from_homogeneous([ALONG, ACROSS], [knots_u, knots_v], ctrl).map_err(|e| fail(format!("the corner patch: {e}")))?;
    Ok((s, m as f64))
}

/// The surface across its own domain: the patch must keep facing one way, or a blend would be met
/// by a piece of the patch folded back on itself.
fn no_fold(surface: &mt::Surface, loop_range: (f64, f64)) -> Result<()> {
    use mt::ParametricSurface;
    // The u knots run past the loop at both ends (a closed curve has no clamped ends to name), so
    // the loop's own range is brought in rather than read off the knot vector.
    let (u0, u1) = loop_range;
    let Some((v0, v1)) = surface.try_range_tuple().1 else {
        return Err(fail("the corner patch has no parameter domain"));
    };
    let normal = |u: f64, v: f64| -> Option<Vec3> { vec_of(surface.uder(u, v)).cross(vec_of(surface.vder(u, v))).normalized() };
    let mid = normal((u0 + u1) * 0.5, v0 + (v1 - v0) * 0.5).ok_or_else(|| unsupported("the corner patch has no normal at its middle"))?;
    // Clear of the apex (v1), where the patch's parameter line along the seam collapses to a point
    // and the normal is read from the sections approaching it.
    for i in 0..=8 {
        for j in 0..9 {
            let u = u0 + (u1 - u0) * i as f64 / 8.0;
            let v = v0 + (v1 - v0) * j as f64 / 9.0;
            let n = normal(u, v).ok_or_else(|| unsupported("the corner patch folds back on itself"))?;
            if n.dot(mid) <= 0.0 {
                return Err(unsupported("the corner patch folds back on itself with these section curves"));
            }
        }
    }
    Ok(())
}

/// Build the patch: sections, tensor surface, and the wire of the section curves themselves.
fn build(sides: &[Vec<Frame>], apex: Vec3) -> Result<Vec<mt::Face>> {
    if !apex.is_finite() {
        return Err(KernelError::Invalid("the corner apex must be a finite point".into()));
    }
    let (frames, starts, size) = corner_loop(sides)?;
    let frames = transverse(&frames, &loop_tangents(&frames))?;
    // The apex's own plane: the mean direction from the apex to the boundary it caps.
    let n_apex = frames
        .iter()
        .fold(Vec3::ZERO, |a, f| a + (f.p - apex))
        .normalized()
        .ok_or_else(|| unsupported("the corner's apex sits on its own boundary curves"))?;
    // How hard the patch bends into the apex: as much as its blends do across their seams, and
    // never flatter than a dome of the corner's own size (a flat cap would end in a point).
    let mean_bend = frames.iter().map(|f| f.bend.len()).sum::<f64>() / frames.len().max(1) as f64;
    let reach = (frames.iter().map(|f| f.p.dist(apex)).sum::<f64>() / frames.len().max(1) as f64).max(size * 1e-3);
    let kappa = mean_bend.max(1.0 / (2.0 * reach));
    if !(kappa.is_finite() && kappa > 0.0) {
        return Err(unsupported("a corner the patch cannot span"));
    }
    let layers: Vec<[Vec3; ACROSS + 1]> = frames.iter().map(|f| section(f, apex, n_apex, kappa)).collect::<Result<_>>()?;
    let (nurbs, m) = tensor(&layers)?;
    let surface = crate::spline_conv::surface(&nurbs)?;
    no_fold(&surface, (ALONG as f64, m + ALONG as f64))?;
    // The boundary is the loop of section curves; each side is cut out of the patch's own boundary
    // curve so that it lies on the patch exactly.
    let nu = nurbs.count().first().copied().unwrap_or(0);
    let row: Vec<HPoint> = (0..nu).map(|i| nurbs.point(i, 0)).collect();
    let Some(knots_u) = nurbs.knots().first().cloned() else { return Err(fail("knots")) };
    let loop_curve = crate::spline_conv::nurbs_curve(mt::KnotVec::from(knots_u), &row)?;
    let mut cuts: Vec<f64> = starts.iter().map(|i| *i as f64 + ALONG as f64).collect();
    cuts.push(m + ALONG as f64);
    let mut verts: Vec<mt::Vertex> = Vec::with_capacity(sides.len());
    for side in sides {
        let first = side.first().ok_or_else(|| fail("a section curve has no stations"))?;
        verts.push(builder::vertex(p3(first.p)));
    }
    let mut edges: Vec<mt::Edge> = Vec::with_capacity(sides.len());
    for i in 0..sides.len() {
        let (Some(&a), Some(&b)) = (cuts.get(i), cuts.get(i + 1)) else { return Err(fail("the corner's seams")) };
        let mut piece = loop_curve.clone();
        let _tail = mt::Cut::cut(&mut piece, b);
        let piece = mt::Cut::cut(&mut piece, a);
        let (Some(from), Some(to)) = (verts.get(i), verts.get((i + 1) % verts.len())) else { return Err(fail("the corner's vertices")) };
        edges.push(mt::Edge::try_new(from, to, piece).map_err(|e| fail(format!("a corner seam: {e}")))?);
    }
    let face = mt::Face::try_new(vec![mt::Wire::from(edges)], surface).map_err(|e| fail(format!("the corner patch: {e}")))?;
    Ok(vec![face])
}

/// What the patch achieves along the seams of `sides`: the largest distance from a seam station to
/// the patch, the largest angle (radians) between the patch's normal there and the blend's, and the
/// largest difference between the curvature across the seam and the blend's own, relative to
/// `reference` — the curvature a rolling-ball corner shows across the same seam, which is what G2
/// removes. Measured on the built surface alone, by the same finite differences of its own samples
/// that `crate::bridge`'s tests read a cross-section's shape with, so it is evidence about the
/// geometry and not about the construction.
pub fn seam_error(surface: &mt::Surface, sides: &[Vec<Frame>], reference: f64) -> Result<(f64, f64, f64)> {
    use mt::{ParametricSurface, SearchNearestParameter};
    if !(reference.is_finite() && reference > 0.0) {
        return Err(KernelError::Invalid("a reference curvature must be finite and positive".into()));
    }
    let (Some((u0, u1)), Some((v0, v1))) = surface.try_range_tuple() else {
        return Err(fail("the corner patch has no parameter domain"));
    };
    // Samples across the section for a five-point difference, as `bridge`'s tests use.
    const INWARDS: usize = 5;
    let hv = (v1 - v0).abs() * 1e-3;
    let hu = (u1 - u0).abs() * 1e-3;
    if !(hv > 0.0 && hu > 0.0) {
        return Err(fail("the corner patch has no extent to measure"));
    }
    let (mut gap, mut angle, mut curvature) = (0.0f64, 0.0f64, 0.0f64);
    for side in sides {
        if side.len() < 5 {
            return Err(KernelError::Invalid("a seam needs at least five stations to measure its direction".into()));
        }
        let mut hint: Option<(f64, f64)> = None;
        for (i, f) in side.iter().enumerate() {
            let Some(t) = seam_tangent(side, i) else { continue };
            let found = surface
                .search_nearest_parameter(p3(f.p), hint, 100)
                .or_else(|| surface.search_nearest_parameter(p3(f.p), None, 200))
                .ok_or_else(|| fail("a seam station is not on the corner patch"))?;
            hint = Some(found);
            let (u, v) = found;
            let at = from_p3(surface.subs(u, v));
            gap = gap.max(at.dist(f.p));
            // The cross-section, sampled inwards from the seam (or from the other end of the patch,
            // with the odd differences turned round).
            let out = v <= (v0 + v1) * 0.5;
            let step = if out { hv } else { -hv };
            let q: Vec<Vec3> = (0..INWARDS).map(|k| from_p3(surface.subs(u, v + step * k as f64))).collect();
            let (Some(&f0), Some(&f1), Some(&f2), Some(&f3), Some(&f4)) = (q.first(), q.get(1), q.get(2), q.get(3), q.get(4)) else {
                return Err(fail("the corner patch has no room for a sample"));
            };
            let sign = if out { 1.0 } else { -1.0 };
            let d1 = (f0 * -25.0 + f1 * 48.0 + f2 * -36.0 + f3 * 16.0 + f4 * -3.0) * (sign / (12.0 * hv));
            let d2 = (f0 * 35.0 + f1 * -104.0 + f2 * 114.0 + f3 * -56.0 + f4 * 11.0) * (1.0 / (12.0 * hv * hv));
            // Along the seam, to stand the tangent plane up.
            let du = if u - hu >= u0 && u + hu <= u1 {
                (from_p3(surface.subs(u + hu, v)) - from_p3(surface.subs(u - hu, v))) * (0.5 / hu)
            } else {
                let (a, b) = if u + hu <= u1 { (u, u + hu) } else { (u, u - hu) };
                (from_p3(surface.subs(b, v)) - from_p3(surface.subs(a, v))) * (1.0 / hu)
            };
            let n = du.cross(d1).normalized().ok_or_else(|| fail("the corner patch has no normal at a seam"))?;
            let blend_n = t.cross(f.leave).normalized().ok_or_else(|| fail("a blend has no normal at a seam"))?;
            angle = angle.max(n.dot(blend_n).abs().clamp(0.0, 1.0).acos());
            // The curvature vector of the cross-seam curve: what the blend's own surface does along
            // the direction the patch leaves it in.
            let speed2 = d1.len2().max(1e-300);
            let bend = (d2 - d1 * (d2.dot(d1) / speed2)) * (1.0 / speed2);
            curvature = curvature.max((bend - f.bend).len() / reference.max(f.bend.len()).max(1e-30));
        }
    }
    Ok((gap, angle, curvature))
}

#[cfg(test)]
#[path = "corner_g2_tests.rs"]
mod tests;

/// The G2 corner through the fillet path: a box's three edges at one corner rounded, rolling ball
/// against G2 patch, measured on the built bodies with the continuity analysis.
#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod fillet_corner_tests {
    use super::*;
    use crate::{ContinuityTolerances, SeamReport, measure, seam_diagnostics};

    /// A 20 mm box with its three edges at one corner rounded by 4: the rolling ball fills that
    /// corner with a piece of its sphere, the G2 route with a patch on the blends' own sections.
    fn cornered(g2: bool) -> crate::Body {
        let b = crate::box_solid(Vec3::ZERO, Vec3::new(20.0, 20.0, 20.0)).unwrap();
        let edges = [Vec3::new(10.0, 0.0, 0.0), Vec3::new(0.0, 10.0, 0.0), Vec3::new(0.0, 0.0, 10.0)];
        if g2 { crate::fillet_g2_corners(&b, &edges, 4.0).unwrap() } else { crate::fillet(&b, &edges, 4.0).unwrap() }
    }

    #[test]
    fn the_ball_corner_is_measurably_g1_and_a_g2_corner_is_refused_until_it_sews() {
        // The rolling ball: tangent to every blend, keeping its own 1/r across each corner seam.
        let body = cornered(false);
        let volume = measure(&body).unwrap().volume;
        assert!(volume > 7000.0 && volume < 8000.0, "the box is still a solid: {volume}");
        let report = seam_diagnostics(&body, &[], 33, ContinuityTolerances { gap_mm: 1e-3, angle_deg: 1.0, curvature_per_mm: 1e-6 }).unwrap();
        // Every seam of a filleted box is tangent; the corner's three are the short ones (the
        // blend-to-face seams run the edge's own 12 mm and more).
        let tangent: Vec<&SeamReport> = report.seams.iter().filter(|s| s.normal_angle_deg.value < 0.01).collect();
        let corner: Vec<&&SeamReport> = tangent.iter().filter(|s| s.length < 8.0).collect();
        assert_eq!(
            corner.len(),
            3,
            "three corner seams among {} tangent ones: {:?}",
            tangent.len(),
            tangent.iter().map(|s| s.length).collect::<Vec<_>>()
        );
        let worst = corner.iter().map(|s| s.curvature_jump_per_mm.value).fold(0.0f64, f64::max);
        eprintln!("rolling ball: corner seams jump {worst:.6} per mm (volume {volume:.3})");
        assert!((worst - 0.25).abs() < 1e-3, "the ball keeps its own 1/r (r = 4) across the seam: {worst}");

        // The G2 corner patch is built on the blends' own sections, but sewing it into the shell
        // is not done yet: the fillet with `corner: g2` refuses with a clear message rather than
        // leaving the shell open, and leaves the ball the corner.
        let b = crate::box_solid(Vec3::ZERO, Vec3::new(20.0, 20.0, 20.0)).unwrap();
        let edges = [Vec3::new(10.0, 0.0, 0.0), Vec3::new(0.0, 10.0, 0.0), Vec3::new(0.0, 0.0, 10.0)];
        let err = crate::fillet_g2_corners(&b, &edges, 4.0).unwrap_err().to_string();
        eprintln!("G2 corner: {err}");
        assert!(err.contains("not supported yet:") && err.contains("G2 corner"), "{err}");
    }
}