//! Building bodies: extrude and revolve planar regions, primitives.

use mt::builder;
use solvecraft_geom::{Loop2, Plane, Region2, Seg2, Vec2, Vec3};
use truck_modeling as mt;

use crate::body::{Body, Solid, from_p3, p3, v3};
use crate::{KernelError, Result, guard};

const MAX_SEGS: usize = 20_000;

fn finite(x: f64, what: &str) -> Result<f64> {
    if x.is_finite() && x.abs() < 1e7 { Ok(x) } else { Err(KernelError::Invalid(format!("{what} = {x}"))) }
}

/// A closed wire in 3D from a sketch loop.
pub(crate) fn wire(plane: &Plane, lp: &Loop2) -> Result<mt::Wire> {
    wire_sources(plane, lp).map(|(wire, _)| wire)
}

fn wire_sources(plane: &Plane, lp: &Loop2) -> Result<(mt::Wire, Vec<(mt::Edge, Seg2)>)> {
    let n = lp.segs.len();
    if n == 0 || n > MAX_SEGS {
        return Err(KernelError::Invalid(format!("loop with {n} segments")));
    }
    if !lp.is_closed(1e-6) {
        return Err(KernelError::Invalid("profile loop is not closed".into()));
    }
    let verts: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(plane.to_world(s.start())))).collect();
    let mut edges = Vec::with_capacity(n);
    let mut sources = Vec::with_capacity(n);
    for (i, s) in lp.segs.iter().enumerate() {
        let (Some(a), Some(b)) = (verts.get(i), verts.get((i + 1) % n)) else { continue };
        let e: mt::Edge = match *s {
            Seg2::Line { .. } => builder::line(a, b),
            Seg2::Arc { sweep, .. } => {
                if sweep.abs() < 1e-9 {
                    return Err(KernelError::Invalid("zero-length arc".into()));
                }
                // Split big arcs so every piece is well conditioned.
                if sweep.abs() > std::f64::consts::PI * 1.5 {
                    let vm = builder::vertex(p3(plane.to_world(s.mid())));
                    let first = builder::circle_arc(a, &vm, p3(plane.to_world(s.point_at(0.25))));
                    sources.push((first.clone(), *s));
                    edges.push(first);
                    builder::circle_arc(&vm, b, p3(plane.to_world(s.point_at(0.75))))
                } else {
                    builder::circle_arc(a, b, p3(plane.to_world(s.mid())))
                }
            }
            Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => {
                crate::freeform::edge(plane, s, a, b).ok_or_else(|| KernelError::Invalid("profile curve".into()))?
            }
        };
        sources.push((e.clone(), *s));
        edges.push(e);
    }
    Ok((edges.into(), sources))
}

/// A planar face from a region (outer loop with holes).
pub(crate) fn face(plane: &Plane, r: &Region2) -> Result<mt::Face> {
    let mut wires = vec![wire(plane, &r.outer.ccw())?];
    for h in &r.holes {
        wires.push(wire(plane, &h.ccw().reversed())?);
    }
    guard("face", || builder::try_attach_plane(&wires).map_err(|e| KernelError::Invalid(format!("profile: {e}"))))
}

fn region_ok(r: &Region2) -> Result<()> {
    if r.area().abs() < 1e-9 {
        return Err(KernelError::Invalid("profile has no area".into()));
    }
    for s in std::iter::once(&r.outer).chain(&r.holes).flat_map(|l| &l.segs) {
        let ok = match *s {
            Seg2::Line { a, b } => a.is_finite() && b.is_finite(),
            Seg2::Arc { center, radius, start, sweep } => {
                center.is_finite() && radius.is_finite() && radius > 1e-9 && start.is_finite() && sweep.is_finite()
            }
            Seg2::Cubic { p0, p1, p2, p3 } => p0.is_finite() && p1.is_finite() && p2.is_finite() && p3.is_finite(),
            Seg2::Conic { a, apex, b, w } => a.is_finite() && apex.is_finite() && b.is_finite() && w.is_finite() && w > 1e-9,
            Seg2::Bezier { n, p } => (4..=solvecraft_geom::MAX_BEZIER_DEGREE as u8).contains(&n) && p.iter().all(|q| q.is_finite()),
        };
        if !ok {
            return Err(KernelError::Invalid("non-finite profile".into()));
        }
    }
    Ok(())
}

/// Extrude each region from `start` to `end` along the plane normal (offsets in mm, `end > start`).
pub fn extrude(plane: &Plane, regions: &[Region2], start: f64, end: f64) -> Result<Vec<Body>> {
    if !plane.origin.is_finite()
        || !plane.x.is_finite()
        || !plane.y.is_finite()
        || (plane.x.len2() - 1.).abs() > 1e-8
        || (plane.y.len2() - 1.).abs() > 1e-8
        || plane.x.dot(plane.y).abs() > 1e-8
    {
        return Err(KernelError::Invalid("extrude frame is not finite and orthonormal".into()));
    }
    let (start, end) = (finite(start, "start")?, finite(end, "end")?);
    let (lo, hi) = if end >= start { (start, end) } else { (end, start) };
    if hi - lo < 1e-6 {
        return Err(KernelError::Invalid("extrude distance is zero".into()));
    }
    if regions.is_empty() {
        return Err(KernelError::Invalid("no profile selected".into()));
    }
    let n = plane.normal();
    let mut out = Vec::new();
    for r in regions {
        region_ok(r)?;
        let base = plane.offset(lo);
        let (outer, mut sources) = wire_sources(&base, &r.outer.ccw())?;
        let mut wires = vec![outer];
        for hole in &r.holes {
            let (wire, mut hole_sources) = wire_sources(&base, &hole.ccw().reversed())?;
            wires.push(wire);
            sources.append(&mut hole_sources);
        }
        let f = guard("extrude profile", || builder::try_attach_plane(&wires).map_err(|e| KernelError::Invalid(format!("profile: {e}"))))?;
        let solid: Solid = guard("extrude", || Ok(builder::tsweep(&f, v3(n * (hi - lo)))))?;
        let upper = guard("extrude cap identities", || {
            let inputs: std::collections::HashSet<_> = sources.iter().map(|(edge, _)| edge.id()).collect();
            let mut visits = 0usize;
            let mut candidate = None;
            for (index, generated) in solid.face_iter().enumerate() {
                if index >= 50_000 {
                    return Ok(None);
                }
                let mut side = false;
                for edge in generated.absolute_boundaries().iter().flatten() {
                    visits += 1;
                    if visits > 200_000 {
                        return Ok(None);
                    }
                    side |= inputs.contains(&edge.id());
                }
                if generated.id() != f.id() && !side {
                    if candidate.is_some() {
                        return Ok(None);
                    }
                    candidate = Some(generated.id());
                }
            }
            Ok(candidate)
        })?;
        let body = Body::new(solid)?;
        // tsweep retains each actual input boundary edge on its generated side.
        // Match topology identity, never face order, surface samples or centroids.
        let provenance = guard("extrude source identities", || {
            let source_edges: std::collections::HashMap<_, _> = sources.iter().map(|(edge, segment)| (edge.id(), *segment)).collect();
            let mut result = Vec::new();
            for generated in body.solid.face_iter() {
                let matching: Vec<_> = generated.absolute_boundaries().iter().flatten().filter_map(|edge| source_edges.get(&edge.id())).collect();
                let source =
                    if generated.id() == f.id() {
                        upper.and_then(|_| {
                            (r.outer.segs.len() + r.holes.iter().map(|hole| hole.segs.len()).sum::<usize>() <= 4096).then(|| {
                                crate::FaceSource::ExtrudeCap {
                                    high: false,
                                    region: r.clone(),
                                    plane: Plane { origin: base.origin, x: plane.x, y: -plane.y },
                                }
                            })
                        })
                    } else if Some(generated.id()) == upper {
                        (r.outer.segs.len() + r.holes.iter().map(|hole| hole.segs.len()).sum::<usize>() <= 4096)
                            .then(|| crate::FaceSource::ExtrudeCap { high: true, region: r.clone(), plane: plane.offset(hi) })
                    } else if let [segment] = matching.as_slice() {
                        let basis = match **segment {
                            Seg2::Line { a, b } => {
                                let d = b - a;
                                d.normalized().and_then(|direction| {
                                    let foot = a - direction * a.dot(direction);
                                    Plane::new(plane.to_world(foot), plane.dir_to_world(direction), n)
                                })
                            }
                            _ => None,
                        };
                        Some(crate::FaceSource::ExtrudeSide { segment: **segment, plane: basis })
                    } else {
                        None
                    };
                result.push(source);
            }
            Ok(result)
        })?;
        out.push(body.with_face_sources(Some(std::sync::Arc::new(provenance))));
    }
    Ok(out)
}

/// Choose the per-profile axis sign whose positive rotation follows the sketch normal.
pub fn canonical_revolve_axis(plane: &Plane, region: &Region2, axis_origin: Vec2, axis_dir: Vec2) -> Result<Vec3> {
    if !axis_origin.is_finite() || !axis_dir.is_finite() || !plane.origin.is_finite() || !plane.x.is_finite() || !plane.y.is_finite() {
        return Err(KernelError::Invalid("non-finite revolve frame".into()));
    }
    if (plane.x.len2() - 1.).abs() > 1e-8 || (plane.y.len2() - 1.).abs() > 1e-8 || plane.x.dot(plane.y).abs() > 1e-8 {
        return Err(KernelError::Invalid("revolve frame is not orthonormal".into()));
    }
    region_ok(region)?;
    let direction = axis_dir.normalized().ok_or_else(|| KernelError::Invalid("revolve axis has no direction".into()))?;
    let axis = plane.dir_to_world(direction).normalized().ok_or_else(|| KernelError::Invalid("revolve frame has no axis".into()))?;
    let origin = plane.to_world(axis_origin);
    let center = plane.to_world(region.centroid());
    let radial = center - origin - axis * (center - origin).dot(axis);
    Ok(if axis.cross(radial).dot(plane.normal()) < 0. { -axis } else { axis })
}

/// Revolve each region about the sketch axis by `angle` radians.
/// This legacy API preserves its historical sweep sign behavior.
pub fn revolve(plane: &Plane, regions: &[Region2], axis_origin: Vec2, axis_dir: Vec2, angle: f64) -> Result<Vec<Body>> {
    let angle = finite(angle, "angle")?;
    if angle.abs() < 1e-6 {
        return Err(KernelError::Invalid("revolve angle is zero".into()));
    }
    let d = axis_dir.normalized().ok_or_else(|| KernelError::Invalid("revolve axis has no direction".into()))?;
    if regions.is_empty() {
        return Err(KernelError::Invalid("no profile selected".into()));
    }
    // The profile must lie on one side of the axis.
    let o = plane.to_world(axis_origin);
    let axis = plane.dir_to_world(d);
    let mut out = Vec::new();
    for r in regions {
        region_ok(r)?;
        let (mut neg, mut pos) = (false, false);
        for p in r.outer.polyline(1e-3) {
            let s = d.cross(p - axis_origin);
            neg |= s < -1e-7;
            pos |= s > 1e-7;
        }
        if neg && pos {
            return Err(KernelError::Invalid("the profile crosses the revolve axis".into()));
        }
        // Points on the axis (within solver tolerance) go exactly onto it, so the sweep does not
        // create hair-thin faces.
        let snap = |p: Vec2| {
            let along = (p - axis_origin).dot(d);
            let q = axis_origin + d * along;
            if p.dist(q) < 1e-7 { q } else { p }
        };
        let snap_loop = |l: &Loop2| Loop2 {
            segs: l
                .segs
                .iter()
                .map(|s| match *s {
                    Seg2::Line { a, b } => Seg2::Line { a: snap(a), b: snap(b) },
                    arc => arc,
                })
                .collect(),
        };
        // Full circles start at an odd angle: their seam points then don't lie in the planes
        // through the profile's centre that other bodies tend to cut along.
        let reseam = |l: Loop2| -> Loop2 {
            let c = match l.segs.first() {
                Some(Seg2::Arc { center, radius, .. }) => Some((*center, *radius)),
                _ => None,
            };
            // (A circle touching the axis keeps its vertex there.)
            let touches = c.is_some_and(|(center, radius)| (d.cross(center - axis_origin).abs() - radius).abs() < 1e-7);
            match c {
                Some((center, radius))
                    if !touches
                        && l.segs.iter().all(
                            |s| matches!(s, Seg2::Arc { center: c2, radius: r2, .. } if c2.dist(center) < 1e-9 && (r2 - radius).abs() < 1e-9),
                        )
                        && (l
                            .segs
                            .iter()
                            .map(|s| match s {
                                Seg2::Arc { sweep, .. } => *sweep,
                                _ => 0.0,
                            })
                            .sum::<f64>()
                            .abs()
                            - std::f64::consts::TAU)
                            .abs()
                            < 1e-6 =>
                {
                    let ccw = l.signed_area() > 0.0;
                    let n = Loop2::circle_from(center, radius, 0.4142);
                    if ccw { n } else { n.reversed() }
                }
                _ => l,
            }
        };
        let r = &Region2 { outer: reseam(snap_loop(&r.outer)), holes: r.holes.iter().map(|h| reseam(snap_loop(h))).collect() };
        // A half disc on the axis revolved fully is a sphere: build it without poles.
        if angle.abs() >= std::f64::consts::TAU - 1e-9
            && r.holes.is_empty()
            && let Some((c, rad)) = half_disc(&r.outer, axis_origin, d)
        {
            out.push(sphere_patches(plane.to_world(c), rad)?);
            continue;
        }
        // A full revolution of a profile with one edge on the axis: sweep the rest of the loop
        // as a cone-like shell (no zero-area faces, which booleans can't handle).
        if angle.abs() >= std::f64::consts::TAU - 1e-9
            && r.holes.is_empty()
            && let Some(b) = revolve_touching_axis(plane, &r.outer, axis_origin, d)?
        {
            out.push(b);
            continue;
        }
        let f = face(plane, r)?;
        // Partial revolves sweep toward the side the sketch normal points to (as features
        // extrude toward it).
        let c = plane.to_world(r.centroid());
        let radial = (c - o) - axis * (c - o).dot(axis);
        let toward = axis.cross(radial).dot(plane.normal());
        let angle = if toward < 0.0 { -angle } else { angle };
        // A profile with an edge on the axis: no zero-area face where that edge turns on itself.
        if r.holes.is_empty()
            && let Some(b) = revolve_partial_on_axis(plane, &r.outer, axis_origin, d, angle)?
        {
            out.push(b);
            continue;
        }
        let solid: Solid = guard("revolve", || Ok(builder::rsweep(&f, p3(o), v3(axis), mt::Rad(angle))))?;
        out.push(Body::new(solid)?);
    }
    Ok(out)
}

/// Axis-aligned box between two corners.
pub fn box_solid(a: Vec3, b: Vec3) -> Result<Body> {
    let (lo, hi) = (a.min(b), a.max(b));
    let s = hi - lo;
    if !(s.x > 1e-6 && s.y > 1e-6 && s.z > 1e-6 && s.is_finite() && s.len() < 1e7) {
        return Err(KernelError::Invalid("box needs a positive size in x, y and z".into()));
    }
    let r = Region2 {
        outer: Loop2::polygon(&[Vec2::new(lo.x, lo.y), Vec2::new(hi.x, lo.y), Vec2::new(hi.x, hi.y), Vec2::new(lo.x, hi.y)]),
        holes: vec![],
    };
    let mut v = extrude(&Plane::XY, &[r], lo.z, hi.z)?;
    let body = v.pop().ok_or_else(|| KernelError::Failed("box".into()))?;
    let sources = body.face_sources().ok_or_else(|| KernelError::Failed("box construction source".into()))?;
    let mut roles = std::collections::HashSet::new();
    let mut certified = Vec::new();
    for source in sources {
        let role = match source {
            Some(crate::FaceSource::ExtrudeCap { high, .. }) => Some((2u8, *high)),
            Some(crate::FaceSource::ExtrudeSide { segment: Seg2::Line { a, b }, .. }) if a.x == b.x && a.x == lo.x => Some((0, false)),
            Some(crate::FaceSource::ExtrudeSide { segment: Seg2::Line { a, b }, .. }) if a.x == b.x && a.x == hi.x => Some((0, true)),
            Some(crate::FaceSource::ExtrudeSide { segment: Seg2::Line { a, b }, .. }) if a.y == b.y && a.y == lo.y => Some((1, false)),
            Some(crate::FaceSource::ExtrudeSide { segment: Seg2::Line { a, b }, .. }) if a.y == b.y && a.y == hi.y => Some((1, true)),
            _ => None,
        }
        .ok_or_else(|| KernelError::Failed("box role correspondence".into()))?;
        if !roles.insert(role) {
            return Err(KernelError::Failed("duplicate box role".into()));
        }
        let (axis, high) = role;
        let sign = if high { 1. } else { -1. };
        let plane = match axis {
            0 => Plane { origin: Vec3::X * if high { hi.x } else { lo.x }, x: Vec3::Y, y: Vec3::Z * sign },
            1 => Plane { origin: Vec3::Y * if high { hi.y } else { lo.y }, x: Vec3::X, y: -Vec3::Z * sign },
            _ => Plane { origin: Vec3::Z * if high { hi.z } else { lo.z }, x: Vec3::X, y: Vec3::Y * sign },
        };
        certified.push(Some(crate::FaceSource::Box { producer: None, axis, high, plane }));
    }
    if roles.len() != 6 {
        return Err(KernelError::Failed("box role count".into()));
    }
    Ok(body.with_face_sources(Some(std::sync::Arc::new(certified))))
}

/// Where primitives start their circles (radians from the frame's x axis).
pub(crate) const SEAM_ANGLE: f64 = 0.6131;

/// Cylinder with base centre `base`, axis direction `axis`, radius and height.
pub fn cylinder(base: Vec3, axis: Vec3, radius: f64, height: f64) -> Result<Body> {
    let plane = Plane::from_normal(base, axis).ok_or_else(|| KernelError::Invalid("cylinder axis".into()))?;
    let radius = finite(radius, "radius")?;
    if radius <= 1e-6 {
        return Err(KernelError::Invalid("radius must be positive".into()));
    }
    // The seams at an odd angle: on the frame's axes they would lie exactly where faces placed
    // on a grid cut the cylinder (along a seam edge, which booleans don't survive).
    let r = Region2 { outer: Loop2::circle_from(Vec2::ZERO, radius, SEAM_ANGLE), holes: vec![] };
    let mut v = extrude(&plane, &[r], 0.0, finite(height, "height")?)?;
    v.pop().ok_or_else(|| KernelError::Failed("cylinder".into()))
}

/// Sphere: revolve a half disc.
pub fn sphere(center: Vec3, radius: f64) -> Result<Body> {
    sphere_patches(center, radius)
}

/// Torus about the Z axis through `center`.
pub fn torus(center: Vec3, major: f64, minor: f64) -> Result<Body> {
    let (major, minor) = (finite(major, "major radius")?, finite(minor, "minor radius")?);
    if !(minor > 1e-6 && major > minor) {
        return Err(KernelError::Invalid("torus needs major > minor > 0".into()));
    }
    let plane = Plane { origin: center, ..Plane::XZ };
    let r = Region2 { outer: Loop2::circle(Vec2::new(major, 0.0), minor), holes: vec![] };
    let mut v = revolve(&plane, &[r], Vec2::ZERO, Vec2::Y, std::f64::consts::TAU)?;
    v.pop().ok_or_else(|| KernelError::Failed("torus".into()))
}

/// Offset every segment of a loop to its right by `d` (outward for a region's outer loop and
/// holes alike). Loops of lines (mitred corners) and loops of arcs on one circle are supported.
fn offset_loop(lp: &Loop2, d: f64) -> Result<Loop2> {
    if d == 0.0 {
        return Ok(lp.clone());
    }
    let all_arcs = lp.segs.iter().all(|s| matches!(s, Seg2::Arc { .. }));
    if all_arcs {
        let segs = lp
            .segs
            .iter()
            .map(|s| match *s {
                Seg2::Arc { center, radius, start, sweep } => {
                    let r = radius + sweep.signum() * d;
                    if r > 1e-6 {
                        Ok(Seg2::Arc { center, radius: r, start, sweep })
                    } else {
                        Err(KernelError::Invalid("the taper closes the profile".into()))
                    }
                }
                _ => Err(KernelError::Failed("not supported yet: tapered profiles mixing lines and arcs".into())),
            })
            .collect::<Result<Vec<_>>>()?;
        return Ok(Loop2 { segs });
    }
    let n = lp.segs.len();
    let mut lines: Vec<(Vec2, Vec2)> = Vec::with_capacity(n);
    for s in &lp.segs {
        let Seg2::Line { a, b } = *s else { return Err(KernelError::Failed("not supported yet: tapered profiles mixing lines and arcs".into())) };
        let dir = (b - a).normalized().ok_or_else(|| KernelError::Invalid("zero-length edge".into()))?;
        let right = Vec2::new(dir.y, -dir.x);
        lines.push((a + right * d, dir));
    }
    let mut pts = Vec::with_capacity(n);
    for i in 0..n {
        let (Some(&(p0, d0)), Some(&(p1, d1))) = (lines.get((i + n - 1) % n), lines.get(i)) else { continue };
        let den = d0.cross(d1);
        if den.abs() < 1e-12 {
            pts.push(p1);
            continue;
        }
        let t = (p1 - p0).cross(d1) / den;
        pts.push(p0 + d0 * t);
    }
    let out = Loop2::polygon(&pts);
    if out.signed_area().signum() != lp.signed_area().signum() || out.signed_area().abs() < 1e-9 {
        return Err(KernelError::Invalid("the taper closes the profile".into()));
    }
    Ok(out)
}

/// Extrude one region by `length` along `dir_sign` × the plane normal with a taper angle
/// (radians; positive grows the profile).
pub fn extrude_tapered(plane: &Plane, region: &Region2, length: f64, dir_sign: f64, taper: f64) -> Result<Body> {
    let (length, taper) = (finite(length, "distance")?, finite(taper, "taper")?);
    if length <= 1e-6 {
        return Err(KernelError::Invalid("extrude distance is zero".into()));
    }
    if taper.abs() >= std::f64::consts::FRAC_PI_2 - 1e-3 {
        return Err(KernelError::Invalid("taper must be between −90° and 90°".into()));
    }
    region_ok(region)?;
    let d = length * taper.tan();
    let top_plane = plane.offset(length * dir_sign.signum());
    let loops0: Vec<Loop2> = std::iter::once(region.outer.ccw()).chain(region.holes.iter().map(|h| h.ccw().reversed())).collect();
    let loops1: Vec<Loop2> = loops0.iter().map(|l| offset_loop(l, d)).collect::<Result<_>>()?;
    guard("tapered extrude", || {
        let w0: Vec<mt::Wire> = loops0.iter().map(|l| wire(plane, l)).collect::<Result<_>>()?;
        let w1: Vec<mt::Wire> = loops1.iter().map(|l| wire(&top_plane, l)).collect::<Result<_>>()?;
        let mut faces: Vec<mt::Face> = Vec::new();
        for (a, b) in w0.iter().zip(&w1) {
            let sides = builder::try_wire_homotopy(a, b).map_err(|e| KernelError::Failed(format!("taper sides: {e}")))?;
            faces.extend(sides.face_iter().cloned());
        }
        let bottom = builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::<Vec<_>>())
            .map_err(|e| KernelError::Failed(format!("taper bottom: {e}")))?;
        let top = builder::try_attach_plane(&w1).map_err(|e| KernelError::Failed(format!("taper top: {e}")))?;
        let mut tries: Vec<(mt::Face, mt::Face)> = vec![(bottom.clone(), top.clone())];
        tries.push((bottom.inverse(), top.inverse()));
        let mut last = String::new();
        for (b, t) in tries {
            let mut all = faces.clone();
            all.push(b);
            all.push(t);
            let shell: mt::Shell = all.into();
            match Solid::try_new(vec![shell]) {
                Ok(s) => return Body::new(s),
                Err(e) => last = e.to_string(),
            }
        }
        Err(KernelError::Failed(format!("tapered extrude: {last}")))
    })
}

/// Revolve 360° a loop that has exactly one line segment on the axis, without degenerate faces.
/// A loop made of a half circle and its diameter on the axis: (centre, radius).
fn half_disc(lp: &Loop2, axis_origin: Vec2, d: Vec2) -> Option<(Vec2, f64)> {
    let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-7 * (1.0 + p.len());
    let mut arcs = Vec::new();
    for s in &lp.segs {
        match *s {
            Seg2::Line { a, b } if on_axis(a) && on_axis(b) => {}
            Seg2::Arc { center, radius, sweep, .. } if on_axis(center) => arcs.push((center, radius, sweep)),
            _ => return None,
        }
    }
    let (c, r, _) = *arcs.first()?;
    let total: f64 = arcs.iter().map(|a| a.2.abs()).sum();
    let same = arcs.iter().all(|(c2, r2, _)| c2.dist(c) < 1e-7 * (1.0 + r) && (r2 - r).abs() < 1e-7 * (1.0 + r));
    (same && (total - std::f64::consts::PI).abs() < 1e-6).then_some((c, r))
}

/// A partial revolve of a loop with one edge on the axis: the other edges sweep their faces,
/// the two caps (the profile at the start and at the end) share the axis edge. `None` when the
/// loop isn't like that.
fn revolve_partial_on_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2, angle: f64) -> Result<Option<Body>> {
    use mt::{BoundedCurve, Invertible, ParametricCurve, ParametricSurface3D, SearchNearestParameter, Transformed};
    let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-9;
    let n = lp.segs.len();
    let axis_segs: Vec<usize> =
        (0..n).filter(|i| lp.segs.get(*i).is_some_and(|s| matches!(s, Seg2::Line { a, b } if on_axis(*a) && on_axis(*b)))).collect();
    let [k] = axis_segs[..] else { return Ok(None) };
    // The open chain after the axis segment, back around to it: only its ends on the axis.
    let chain: Vec<Seg2> = (1..n).filter_map(|j| lp.segs.get((k + j) % n).copied()).collect();
    if chain.is_empty() || chain.iter().skip(1).any(|s| on_axis(s.start())) || chain.iter().any(|s| on_axis(s.mid())) {
        return Ok(None);
    }
    // Built with the axis through the world origin, then moved into place.
    let shift = plane.to_world(axis_origin);
    let local = Plane { origin: plane.origin - shift, ..*plane };
    let plane = &local;
    let axis = plane.dir_to_world(d);
    let rot = mt::Matrix4::from_axis_angle(v3(axis), mt::Rad(angle));
    let half = mt::Matrix4::from_axis_angle(v3(axis), mt::Rad(angle * 0.5));
    let turn = |m: &mt::Matrix4, p: Vec3| from_p3(mt::EuclideanSpace::from_vec((m * p3(p).to_homogeneous()).truncate()));
    guard("revolve", || {
        let fail = |m: &str| KernelError::Failed(format!("revolve: {m}"));
        let pts: Vec<Vec3> = chain.iter().map(|s| plane.to_world(s.start())).chain(chain.last().map(|s| plane.to_world(s.end()))).collect();
        let last = pts.len() - 1;
        let pv: Vec<mt::Vertex> = pts.iter().map(|p| builder::vertex(p3(*p))).collect();
        // The ends stay where they are (on the axis).
        let qv: Vec<mt::Vertex> = pts
            .iter()
            .enumerate()
            .map(|(i, p)| {
                if i == 0 || i == last { pv.get(i).cloned().unwrap_or_else(|| builder::vertex(p3(*p))) } else { builder::vertex(p3(turn(&rot, *p))) }
            })
            .collect();
        let mut e: Vec<mt::Edge> = Vec::new();
        let mut er: Vec<mt::Edge> = Vec::new();
        for (i, sg) in chain.iter().enumerate() {
            let (Some(a), Some(b), Some(qa), Some(qb)) = (pv.get(i), pv.get(i + 1), qv.get(i), qv.get(i + 1)) else { return Err(fail("chain")) };
            let edge = match *sg {
                Seg2::Line { .. } => builder::line(a, b),
                Seg2::Arc { .. } => builder::circle_arc(a, b, p3(plane.to_world(sg.mid()))),
                Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => {
                    crate::freeform::edge(plane, sg, a, b).unwrap_or_else(|| builder::line(a, b))
                }
            };
            let curve = edge.oriented_curve().transformed(rot);
            er.push(mt::Edge::new(qa, qb, curve));
            e.push(edge);
        }
        // Arcs swept by the vertices off the axis.
        let mut arcs: Vec<Option<mt::Edge>> = vec![None; pts.len()];
        for i in 1..last {
            let (Some(a), Some(b), Some(p)) = (pv.get(i), qv.get(i), pts.get(i)) else { continue };
            if let Some(slot) = arcs.get_mut(i) {
                *slot = Some(builder::circle_arc(a, b, p3(turn(&half, *p))));
            }
        }
        // Side faces: each edge's swept surface, oriented to agree with its boundary.
        let mut faces: Vec<mt::Face> = Vec::new();
        for (i, sg) in chain.iter().enumerate() {
            let (Some(ei), Some(eri)) = (e.get(i), er.get(i)) else { return Err(fail("edges")) };
            let mut ws: Vec<mt::Edge> = vec![ei.clone()];
            if let Some(Some(c)) = arcs.get(i + 1) {
                ws.push(c.clone());
            }
            ws.push(eri.inverse());
            if let Some(Some(c)) = arcs.get(i) {
                ws.push(c.inverse());
            }
            let wire: mt::Wire = ws.into();
            let radial = matches!(*sg, Seg2::Line { a, b } if (b - a).normalized().is_some_and(|u| u.dot(d).abs() < 1e-9));
            let f = if radial {
                builder::try_attach_plane(&[wire]).map_err(|er| fail(&er.to_string()))?
            } else {
                let template = builder::rsweep(ei, mt::Point3::new(0.0, 0.0, 0.0), v3(axis), mt::Rad(angle));
                let mut surface = template.face_iter().next().map(|f| f.oriented_surface()).ok_or_else(|| fail("swept face"))?;
                // The loop's turning (Newell) against the surface normal halfway round.
                let mut nw = Vec3::ZERO;
                let poly: Vec<Vec3> = wire
                    .edge_iter()
                    .flat_map(|ed| {
                        let c = ed.oriented_curve();
                        let (t0, t1) = c.range_tuple();
                        (0..8).map(move |j| from_p3(c.subs(t0 + (t1 - t0) * j as f64 / 8.0))).collect::<Vec<_>>()
                    })
                    .collect();
                for (j, p) in poly.iter().enumerate() {
                    if let Some(q) = poly.get((j + 1) % poly.len()) {
                        nw += p.cross(*q);
                    }
                }
                let probe = turn(&half, plane.to_world(sg.mid()));
                let sn = surface.search_nearest_parameter(p3(probe), None, 100).map(|(u, v)| surface.normal(u, v));
                if sn.is_some_and(|m| nw.dot(Vec3::new(m.x, m.y, m.z)) < 0.0) {
                    surface = surface.inverse();
                }
                mt::Face::try_new(vec![wire], surface).map_err(|er| fail(&er.to_string()))?
            };
            faces.push(f);
        }
        // Caps: the profile at the start and at the end, sharing the axis edge.
        let (Some(a0), Some(an)) = (pv.first(), pv.get(last)) else { return Err(fail("ends")) };
        let axis_edge = builder::line(an, a0);
        let start: mt::Wire = e.iter().rev().map(|x| x.inverse()).chain(std::iter::once(axis_edge.inverse())).collect::<Vec<_>>().into();
        let end: mt::Wire = er.iter().cloned().chain(std::iter::once(axis_edge.clone())).collect::<Vec<_>>().into();
        let cs = builder::try_attach_plane(&[start]).map_err(|er| fail(&er.to_string()))?;
        let ce = builder::try_attach_plane(&[end]).map_err(|er| fail(&er.to_string()))?;
        let mut last_err = String::new();
        for flip in 0..4 {
            let mut all = faces.clone();
            all.push(if flip & 1 == 0 { cs.clone() } else { cs.inverse() });
            all.push(if flip & 2 == 0 { ce.clone() } else { ce.inverse() });
            match Solid::try_new(vec![all.into()]) {
                Ok(solid) => {
                    let solid = if shift.len() > 0.0 { builder::translated(&solid, v3(shift)) } else { solid };
                    return Body::new(solid).map(Some);
                }
                Err(er) => last_err = er.to_string(),
            }
        }
        Err(fail(&last_err))
    })
}

fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2) -> Result<Option<Body>> {
    let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-9;
    let n = lp.segs.len();
    let axis_segs: Vec<usize> =
        (0..n).filter(|i| lp.segs.get(*i).is_some_and(|s| matches!(s, Seg2::Line { a, b } if on_axis(*a) && on_axis(*b)))).collect();
    let [k] = axis_segs[..] else { return Ok(None) };
    // The open chain after the axis segment, back around to it.
    let mut chain: Vec<Seg2> = (1..n).filter_map(|j| lp.segs.get((k + j) % n).copied()).collect();
    // A single arc (a sphere) is split off-centre: a seam on the equator would lie in the very
    // planes other bodies are often cut with.
    if let [Seg2::Arc { center, radius, start, sweep }] = chain[..] {
        let f = 0.37;
        chain = vec![
            Seg2::Arc { center, radius, start, sweep: sweep * f },
            Seg2::Arc { center, radius, start: start + sweep * f, sweep: sweep * (1.0 - f) },
        ];
    }
    if chain.is_empty() || chain.iter().skip(1).take(chain.len().saturating_sub(2)).any(|s| on_axis(s.start()) && on_axis(s.end())) {
        return Ok(None);
    }
    // Built with the axis through the world origin (truck tessellates such cones reliably), then
    // moved into place.
    let shift = plane.to_world(axis_origin);
    let local = Plane { origin: plane.origin - shift, ..*plane };
    let plane = &local;
    // Radial lines at the ends become planar caps (no degenerate disc centres).
    let radial = |sg: &Seg2| matches!(*sg, Seg2::Line { a, b } if (b - a).normalized().is_some_and(|u| u.dot(d).abs() < 1e-9));
    let mut mid: &[Seg2] = &chain;
    if mid.first().is_some_and(|f| radial(f) && on_axis(f.start())) {
        mid = mid.get(1..).unwrap_or(&[]);
    }
    if mid.last().is_some_and(|l| radial(l) && on_axis(l.end())) {
        mid = mid.get(..mid.len().saturating_sub(1)).unwrap_or(&[]);
    }
    if mid.is_empty() {
        return Ok(None);
    }
    let apex_start = mid.first().is_some_and(|f| on_axis(f.start()));
    let apex_end = mid.last().is_some_and(|l| on_axis(l.end()));
    guard("revolve", || {
        let verts: Vec<mt::Vertex> = mid
            .iter()
            .map(|s| builder::vertex(p3(plane.to_world(s.start()))))
            .chain(mid.last().map(|s| builder::vertex(p3(plane.to_world(s.end())))))
            .collect();
        let mut edges: Vec<mt::Edge> = Vec::new();
        for (i, s) in mid.iter().enumerate() {
            let (Some(a), Some(b)) = (verts.get(i), verts.get(i + 1)) else { continue };
            edges.push(match *s {
                Seg2::Line { .. } => builder::line(a, b),
                Seg2::Arc { .. } => builder::circle_arc(a, b, p3(plane.to_world(s.mid()))),
                Seg2::Cubic { .. } | Seg2::Conic { .. } | Seg2::Bezier { .. } => match crate::freeform::edge(plane, s, a, b) {
                    Some(e) => e,
                    None => builder::line(a, b),
                },
            });
        }
        let mut w: mt::Wire = edges.into();
        let axis = plane.dir_to_world(d);
        let mut shell = if apex_start || apex_end {
            if !apex_start {
                // builder::cone wants the apex first.
                w = w.inverse();
            }
            builder::cone(&w, v3(axis), mt::Rad(std::f64::consts::TAU))
        } else {
            builder::rsweep(&w, p3(plane.to_world(axis_origin)), v3(axis), mt::Rad(std::f64::consts::TAU))
        };
        let sides = shell.len();
        for b in shell.extract_boundaries() {
            let cap = builder::try_attach_plane(&[b]).map_err(|e| KernelError::Failed(format!("revolve cap: {e}")))?;
            shell.push(cap);
        }
        let solid = Solid::try_new(vec![shell.clone()]).or_else(|_| {
            // Caps may come out facing inwards: flip them.
            let mut s2 = shell.clone();
            for f in s2.face_iter_mut().skip(sides) {
                f.invert();
            }
            Solid::try_new(vec![s2])
        });
        let solid = solid.map_err(|e| KernelError::Failed(format!("revolve: {e}")))?;
        let solid = if shift.len() > 0.0 { builder::translated(&solid, v3(shift)) } else { solid };
        Body::new(solid).map(Some)
    })
}

/// Close a ruled shell (side faces between matching wires) with planar caps into a body.
pub(crate) fn cap_and_close(mut faces: Vec<mt::Face>, w0: &[mt::Wire], w1: &[mt::Wire]) -> Result<Body> {
    let bottom =
        builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::<Vec<_>>()).map_err(|e| KernelError::Failed(format!("cap: {e}")))?;
    let top = builder::try_attach_plane(w1).map_err(|e| KernelError::Failed(format!("cap: {e}")))?;
    let mut last = String::new();
    for (b, t) in [(bottom.clone(), top.clone()), (bottom.inverse(), top.inverse()), (bottom.clone(), top.inverse()), (bottom.inverse(), top)] {
        let mut all = std::mem::take(&mut faces);
        let keep = all.clone();
        all.push(b);
        all.push(t);
        let shell: mt::Shell = all.into();
        match Solid::try_new(vec![shell]) {
            Ok(s) => return Body::new(s),
            Err(e) => last = e.to_string(),
        }
        faces = keep;
    }
    Err(KernelError::Failed(format!("closing the shell: {last}")))
}

/// Re-split a loop into exactly `n` segments: a full circle gets breaks at `angles` (or evenly),
/// other loops split their longest segments in half.
fn resample(lp: &Loop2, n: usize, angles: &[f64]) -> Loop2 {
    let circle = lp.segs.iter().all(|s| matches!(s, Seg2::Arc { .. }))
        && lp.segs.windows(2).all(|w| match (w[0], w[1]) {
            (Seg2::Arc { center: c1, radius: r1, .. }, Seg2::Arc { center: c2, radius: r2, .. }) => c1.dist(c2) < 1e-9 && (r1 - r2).abs() < 1e-9,
            _ => false,
        });
    if circle && let Some(Seg2::Arc { center, radius, .. }) = lp.segs.first().copied() {
        let mut a: Vec<f64> = if angles.len() == n { angles.to_vec() } else { (0..n).map(|i| std::f64::consts::TAU * i as f64 / n as f64).collect() };
        // Keep the given order but make them increase.
        for i in 1..a.len() {
            while a[i] <= a[i - 1] {
                a[i] += std::f64::consts::TAU;
            }
        }
        let segs = (0..n)
            .map(|i| {
                let s = a[i];
                let e = if i + 1 < n { a[i + 1] } else { a[0] + std::f64::consts::TAU };
                Seg2::Arc { center, radius, start: s, sweep: e - s }
            })
            .collect();
        return Loop2 { segs };
    }
    let mut segs = lp.segs.clone();
    while segs.len() < n {
        let Some((k, _)) = segs.iter().enumerate().max_by(|a, b| a.1.length().total_cmp(&b.1.length())) else { break };
        let s = segs.remove(k);
        let (p, q) = s.split_half();
        segs.insert(k, q);
        segs.insert(k, p);
    }
    Loop2 { segs }
}

/// Ruled loft through two or more planar sections (outer loops only). The sections and the
/// same sections in reverse order give the same body.
pub fn loft(sections: &[(Plane, Loop2)]) -> Result<Body> {
    if sections.len() < 2 || sections.len() > 50 {
        return Err(KernelError::Invalid("a loft needs 2…50 sections".into()));
    }
    let n = sections.iter().map(|(_, l)| l.segs.len()).max().unwrap_or(0);
    if n == 0 || n > MAX_SEGS {
        return Err(KernelError::Invalid("empty section".into()));
    }
    // A canonical order, so that reversing the sections changes nothing: the one whose
    // first-to-last step has a positive largest component. The side faces' parameterisation
    // follows the order, and the boolean worked on one and failed on the other (a cut whose
    // loft touches the body along an edge: an internal panic, or "shell not oriented and
    // closed").
    let centre = |(p, l): &(Plane, Loop2)| p.to_world(l.centroid());
    let (Some(first), Some(last)) = (sections.first(), sections.last()) else { return Err(KernelError::Invalid("no sections".into())) };
    let step = centre(last) - centre(first);
    let big = [step.x, step.y, step.z].into_iter().max_by(|a, b| a.abs().total_cmp(&b.abs())).unwrap_or(0.0);
    let ordered: Vec<&(Plane, Loop2)> = if big < 0.0 { sections.iter().rev().collect() } else { sections.iter().collect() };
    // Orient every section counter-clockwise about the first section's normal (`Body::new`
    // turns the solid outward), then match segment counts and starting points to the section
    // with the most segments.
    let n0 = ordered.first().map(|(p, _)| p.normal()).unwrap_or(Vec3::Z);
    let oriented: Vec<(Plane, Loop2)> = ordered
        .iter()
        .map(|(p, l)| {
            let l = if p.normal().dot(n0) >= 0.0 { l.ccw() } else { l.ccw().reversed() };
            (*p, l)
        })
        .collect();
    let lead = oriented.iter().max_by_key(|(_, l)| l.segs.len()).cloned().ok_or_else(|| KernelError::Invalid("no sections".into()))?;
    let lead_c = lead.0.to_world(lead.1.centroid());
    let lead_dirs: Vec<Vec3> = lead.1.segs.iter().map(|s| lead.0.to_world(s.start()) - lead_c).collect();
    let mut wires: Vec<mt::Wire> = Vec::new();
    for (p, l) in &oriented {
        let c = l.centroid();
        // Break angles: the lead section's vertex directions seen in this plane.
        let angles: Vec<f64> = lead_dirs.iter().map(|d| Vec2::new(d.dot(p.x), d.dot(p.y)).angle()).collect();
        let mut r = resample(l, n, &angles);
        // Start at the vertex closest in direction to the lead's first vertex.
        if let Some(d0) = lead_dirs.first() {
            let target = Vec2::new(d0.dot(p.x), d0.dot(p.y)).angle();
            let best = (0..r.segs.len()).min_by(|a, b| {
                let ang = |k: usize| {
                    let v = r.segs[k].start() - c;
                    let d = (v.angle() - target).rem_euclid(std::f64::consts::TAU);
                    d.min(std::f64::consts::TAU - d)
                };
                ang(*a).total_cmp(&ang(*b))
            });
            if let Some(k) = best {
                r.segs.rotate_left(k);
            }
        }
        wires.push(uniform_wire(p, &r)?);
    }
    // Three sections or more: one smooth surface through them all for each segment.
    if wires.len() >= 3
        && let Ok(b) = guard("loft", || smooth_loft(&wires))
    {
        return Ok(b);
    }
    guard("loft", || {
        let mut faces: Vec<mt::Face> = Vec::new();
        for w in wires.windows(2) {
            let sides = builder::try_wire_homotopy(&w[0], &w[1]).map_err(|e| KernelError::Failed(format!("loft: {e}")))?;
            faces.extend(sides.face_iter().cloned());
        }
        let (Some(first), Some(last)) = (wires.first(), wires.last()) else { return Err(KernelError::Failed("loft".into())) };
        cap_and_close(faces, std::slice::from_ref(first), std::slice::from_ref(last))
    })
}

/// Points across the sections (one per section, in order) as a smooth B-spline: cubic through
/// four or more, quadratic through three.
fn across(pts: &[Vec3]) -> Option<mt::BSplineCurve<mt::Point3>> {
    match pts {
        [a, b, c] => {
            let m = *b * 2.0 - (*a + *c) * 0.5;
            mt::BSplineCurve::try_new(mt::KnotVec::bezier_knot(2), vec![p3(*a), p3(m), p3(*c)]).ok()
        }
        _ => interpolate_cubic(pts),
    }
}

/// A loft whose side faces are B-spline surfaces through every section (smooth across the
/// middle sections, which leave no edges). Every section wire comes from `uniform_wire`, so the
/// k-th edges of all sections share one knot vector and their control points pair up.
fn smooth_loft(wires: &[mt::Wire]) -> Result<Body> {
    let fail = |m: &str| KernelError::Failed(format!("smooth loft: {m}"));
    let edges: Vec<Vec<mt::Edge>> = wires.iter().map(|w| w.edge_iter().cloned().collect()).collect();
    let n = edges.first().map(Vec::len).ok_or_else(|| fail("no sections"))?;
    if edges.iter().any(|e| e.len() != n) {
        return Err(fail("sections with different edge counts"));
    }
    let curve_of = |e: &mt::Edge| -> Option<mt::BSplineCurve<mt::Point3>> {
        match e.curve() {
            mt::Curve::BSplineCurve(c) if e.orientation() => Some(c),
            _ => None,
        }
    };
    let (Some(first), Some(last)) = (edges.first(), edges.last()) else { return Err(fail("sections")) };
    // Seams: through the k-th vertex of every section.
    let mut seams: Vec<mt::Edge> = Vec::with_capacity(n);
    for k in 0..n {
        let pts: Vec<Vec3> = edges.iter().filter_map(|es| es.get(k).map(|e| from_p3(e.front().point()))).collect();
        let c = across(&pts).ok_or_else(|| fail("seam"))?;
        let (Some(a), Some(b)) = (first.get(k), last.get(k)) else { return Err(fail("seam")) };
        seams.push(mt::Edge::new(a.front(), b.front(), mt::Curve::BSplineCurve(c)));
    }
    let mut faces = Vec::with_capacity(n);
    for k in 0..n {
        let curves: Vec<mt::BSplineCurve<mt::Point3>> =
            edges.iter().map(|es| es.get(k).and_then(curve_of)).collect::<Option<_>>().ok_or_else(|| fail("section curves"))?;
        let c0 = curves.first().ok_or_else(|| fail("section curves"))?;
        let uk = c0.knot_vec().clone();
        let m = c0.control_points().len();
        if curves.iter().any(|c| c.control_points().len() != m || c.knot_vec() != &uk) {
            return Err(fail("sections parameterised differently"));
        }
        let mut ctrl: Vec<Vec<mt::Point3>> = Vec::with_capacity(m);
        let mut vk = None;
        for j in 0..m {
            let pts: Vec<Vec3> = curves.iter().filter_map(|c| c.control_points().get(j).map(|p| from_p3(*p))).collect();
            let c = across(&pts).ok_or_else(|| fail("surface"))?;
            vk = Some(c.knot_vec().clone());
            ctrl.push(c.control_points().clone());
        }
        let vk = vk.ok_or_else(|| fail("surface"))?;
        let surface = mt::BSplineSurface::try_new((uk, vk), ctrl).map_err(|e| fail(&e.to_string()))?;
        let (Some(e0), Some(e1), Some(s0), Some(s1)) = (first.get(k), last.get(k), seams.get(k), seams.get((k + 1) % n)) else {
            return Err(fail("face"));
        };
        // u along the section edge, v across the sections: counter-clockwise in (u, v).
        let w: mt::Wire = vec![e0.clone(), s1.clone(), e1.inverse(), s0.inverse()].into();
        faces.push(mt::Face::try_new(vec![w], mt::Surface::BSplineSurface(surface)).map_err(|e| fail(&e.to_string()))?);
    }
    let (Some(w0), Some(w1)) = (wires.first(), wires.last()) else { return Err(fail("sections")) };
    cap_and_close(faces, std::slice::from_ref(w0), std::slice::from_ref(w1))
}

/// A 3D path segment for sweeps.
#[derive(Clone, Copy, Debug)]
pub enum PathSeg {
    Line {
        a: Vec3,
        b: Vec3,
    },
    /// Arc around `center`, rotating by `angle` about `axis` (right hand) from `a`.
    Arc {
        a: Vec3,
        center: Vec3,
        axis: Vec3,
        angle: f64,
    },
}

impl PathSeg {
    /// Unit directions of the path where the segment starts and ends.
    fn tangents(&self) -> Option<(Vec3, Vec3)> {
        match *self {
            PathSeg::Line { a, b } => {
                let t = (b - a).normalized()?;
                Some((t, t))
            }
            PathSeg::Arc { a, center, axis, angle } => {
                let k = axis.normalized()?;
                let end = rotate_about(a, center, k, angle);
                Some((k.cross(a - center).normalized()?, k.cross(end - center).normalized()?))
            }
        }
    }

    /// Points along the segment (an arc within `tol` of its chords), without its start.
    fn points_after_start(&self, tol: f64) -> Vec<Vec3> {
        match *self {
            PathSeg::Line { b, .. } => vec![b],
            PathSeg::Arc { a, center, axis, angle } => {
                let Some(k) = axis.normalized() else { return vec![a] };
                let r = (a - center).len();
                let step = if r > tol { 2.0 * (1.0 - tol / r).acos() } else { angle.abs() };
                let n = ((angle.abs() / step.max(1e-6)).ceil() as usize).clamp(1, 512);
                (1..=n).map(|i| rotate_about(a, center, k, angle * i as f64 / n as f64)).collect()
            }
        }
    }
}

/// `p` turned by `angle` about the unit axis `k` through `center` (right hand).
fn rotate_about(p: Vec3, center: Vec3, k: Vec3, angle: f64) -> Vec3 {
    let v = p - center;
    let (s, c) = angle.sin_cos();
    center + v * c + k.cross(v) * s + k * (k.dot(v) * (1.0 - c))
}

/// Sweep a planar region (outer loop) along a chain of path segments, keeping the profile
/// perpendicular to the path (lines translate it, arcs rotate it). A path with sharp corners
/// is swept as a polyline, mitred at the corners.
pub fn sweep(plane: &Plane, region: &Region2, path: &[PathSeg]) -> Result<Body> {
    region_ok(region)?;
    if path.is_empty() || path.len() > 1000 {
        return Err(KernelError::Invalid("the path needs 1…1000 segments".into()));
    }
    // Segment by segment the profile only follows a path whose direction runs on across each
    // joint: at a corner the next line would slide the profile within its own plane.
    let ends: Vec<(Vec3, Vec3)> =
        path.iter().map(PathSeg::tangents).collect::<Option<_>>().ok_or_else(|| KernelError::Invalid("a path segment has no direction".into()))?;
    let corner = ends.windows(2).any(|w| w[0].1.dot(w[1].0) < (1e-3f64).cos());
    if corner {
        let first = match path.first() {
            Some(PathSeg::Line { a, .. } | PathSeg::Arc { a, .. }) => *a,
            None => return Err(KernelError::Invalid("the path needs 1…1000 segments".into())),
        };
        let mut pts = vec![first];
        for seg in path {
            pts.extend(seg.points_after_start(1e-2));
        }
        return crate::sweep_path::sweep_path(plane, region, &pts);
    }
    // The outer loop and each hole, swept along the path; the ends are capped with holes.
    let mut starts = vec![wire(plane, &region.outer.ccw())?];
    for h in &region.holes {
        starts.push(wire(plane, &h.ccw().reversed())?);
    }
    guard("sweep", || {
        let mut shell = mt::Shell::new();
        let mut ends = Vec::new();
        for w0 in &starts {
            let mut cur = w0.clone();
            for seg in path {
                let mut part: mt::Shell = match *seg {
                    PathSeg::Line { a, b } => builder::tsweep(&cur, v3(b - a)),
                    PathSeg::Arc { center, axis, angle, .. } => {
                        let ax = axis.normalized().ok_or_else(|| KernelError::Invalid("arc axis".into()))?;
                        builder::rsweep(&cur, p3(center), v3(ax), mt::Rad(angle))
                    }
                };
                let next = part
                    .extract_boundaries()
                    .into_iter()
                    .find(|w| w.edge_iter().all(|e| !cur.edge_iter().any(|c| c.id() == e.id())))
                    .ok_or_else(|| KernelError::Failed("sweep: lost the profile".into()))?;
                shell.append(&mut part);
                cur = next.inverse();
            }
            ends.push(cur);
        }
        let faces: Vec<mt::Face> = shell.face_iter().cloned().collect();
        cap_and_close(faces, &starts, &ends)
    })
}

/// Cox–de Boor basis values N_{j,3}(u) for all j.
fn cubic_basis(knots: &[f64], n_ctrl: usize, u: f64) -> Vec<f64> {
    let p = 3;
    let m = knots.len();
    let mut n = vec![0.0; m.saturating_sub(1)];
    let last = n_ctrl.saturating_sub(1);
    for (j, v) in n.iter_mut().enumerate() {
        let (Some(&a), Some(&b)) = (knots.get(j), knots.get(j + 1)) else { continue };
        // Half-open spans, with the last span closed at the end.
        if (a <= u && u < b) || (u >= b && (b - a) > 0.0 && j == last && (u - b).abs() < 1e-12) {
            *v = 1.0;
        }
    }
    for k in 1..=p {
        for j in 0..m.saturating_sub(1 + k) {
            let (Some(&tj), Some(&tjk), Some(&tj1), Some(&tjk1)) = (knots.get(j), knots.get(j + k), knots.get(j + 1), knots.get(j + k + 1)) else {
                continue;
            };
            let left = if tjk - tj > 0.0 { (u - tj) / (tjk - tj) * n.get(j).copied().unwrap_or(0.0) } else { 0.0 };
            let right = if tjk1 - tj1 > 0.0 { (tjk1 - u) / (tjk1 - tj1) * n.get(j + 1).copied().unwrap_or(0.0) } else { 0.0 };
            if let Some(v) = n.get_mut(j) {
                *v = left + right;
            }
        }
    }
    n.truncate(n_ctrl);
    n
}

/// Cubic B-spline through points at uniform parameters (clamped, knots by averaging).
pub(crate) fn interpolate_cubic(pts: &[Vec3]) -> Option<mt::BSplineCurve<mt::Point3>> {
    let m = pts.len();
    if m < 4 {
        return None;
    }
    let u: Vec<f64> = (0..m).map(|i| i as f64 / (m - 1) as f64).collect();
    let mut knots = vec![0.0; 4];
    for j in 1..(m - 3) {
        knots.push((j..j + 3).filter_map(|i| u.get(i)).sum::<f64>() / 3.0);
    }
    knots.extend([1.0; 4]);
    // Dense solve of N · P = Q for each coordinate.
    let mut a: Vec<Vec<f64>> = u.iter().map(|ui| cubic_basis(&knots, m, *ui)).collect();
    let mut q: Vec<[f64; 3]> = pts.iter().map(|p| [p.x, p.y, p.z]).collect();
    for col in 0..m {
        let piv = (col..m).max_by(|x, y| a[*x][col].abs().total_cmp(&a[*y][col].abs()))?;
        if a[piv][col].abs() < 1e-14 {
            return None;
        }
        a.swap(col, piv);
        q.swap(col, piv);
        for r in 0..m {
            if r == col {
                continue;
            }
            let f = a[r][col] / a[col][col];
            if f == 0.0 {
                continue;
            }
            let prow = a[col].clone();
            for (x, v) in a[r].iter_mut().zip(&prow).skip(col) {
                *x -= f * v;
            }
            let qc = q[col];
            for k in 0..3 {
                q[r][k] -= f * qc[k];
            }
        }
    }
    let ctrl: Vec<mt::Point3> = (0..m).map(|i| mt::Point3::new(q[i][0] / a[i][i], q[i][1] / a[i][i], q[i][2] / a[i][i])).collect();
    mt::BSplineCurve::try_new(mt::KnotVec::from(knots), ctrl).ok()
}

/// A closed wire whose edges are cubic B-splines parameterised (nearly) by arc length, so that
/// ruled surfaces between sections pair points proportionally along each edge.
pub(crate) fn uniform_wire(plane: &Plane, lp: &Loop2) -> Result<mt::Wire> {
    const SAMPLES: usize = 24;
    let n = lp.segs.len();
    let verts: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(plane.to_world(s.start())))).collect();
    let mut edges = Vec::with_capacity(n);
    for (i, s) in lp.segs.iter().enumerate() {
        let (Some(a), Some(b)) = (verts.get(i), verts.get((i + 1) % n)) else { continue };
        let pts: Vec<Vec3> = (0..SAMPLES).map(|k| plane.to_world(s.point_at(k as f64 / (SAMPLES - 1) as f64))).collect();
        let curve = interpolate_cubic(&pts).ok_or_else(|| KernelError::Failed("loft curve".into()))?;
        edges.push(mt::Edge::new(a, b, mt::Curve::BSplineCurve(curve)));
    }
    Ok(edges.into())
}

/// A sphere as six patches (the cube's faces projected onto it) on revolved surfaces whose poles
/// and seams lie outside each patch. No face has a singular point, which the boolean needs; the
/// patches are one analytic sphere, so measured topology counts it as one face.
pub fn sphere_patches(center: Vec3, radius: f64) -> Result<Body> {
    sphere_patches_rotated(center, radius, [Vec3::X, Vec3::Y, Vec3::Z])
}

/// [`sphere_patches`] with the patch layout turned to the frame `rot` (orthonormal columns).
pub fn sphere_patches_rotated(center: Vec3, radius: f64, rot: [Vec3; 3]) -> Result<Body> {
    let radius = finite(radius, "radius")?;
    let turn = |v: Vec3| rot[0] * v.x + rot[1] * v.y + rot[2] * v.z;
    if radius <= 1e-6 || !center.is_finite() {
        return Err(KernelError::Invalid("radius must be positive".into()));
    }
    let k = radius / 3f64.sqrt();
    let corner = |sx: f64, sy: f64, sz: f64| center + turn(Vec3::new(sx * k, sy * k, sz * k));
    // Cube faces: outward direction, rotation axis for the patch's surface (perpendicular to it),
    // and the four corner sign triples counter-clockwise seen from outside.
    let faces: [(Vec3, Vec3, [[f64; 3]; 4]); 6] = [
        (Vec3::X, Vec3::Z, [[1.0, -1.0, -1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0]]),
        (-Vec3::X, Vec3::Z, [[-1.0, 1.0, -1.0], [-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, 1.0]]),
        (Vec3::Y, Vec3::Z, [[1.0, 1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]),
        (-Vec3::Y, Vec3::Z, [[-1.0, -1.0, -1.0], [1.0, -1.0, -1.0], [1.0, -1.0, 1.0], [-1.0, -1.0, 1.0]]),
        (Vec3::Z, Vec3::X, [[-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0]]),
        (-Vec3::Z, Vec3::X, [[-1.0, 1.0, -1.0], [1.0, 1.0, -1.0], [1.0, -1.0, -1.0], [-1.0, -1.0, -1.0]]),
    ];
    guard("sphere", || {
        let key = |s: [f64; 3]| ((s[0] > 0.0) as usize) | (((s[1] > 0.0) as usize) << 1) | (((s[2] > 0.0) as usize) << 2);
        let mut verts: std::collections::HashMap<usize, mt::Vertex> = std::collections::HashMap::new();
        let mut vert = |s: [f64; 3]| verts.entry(key(s)).or_insert_with(|| builder::vertex(p3(corner(s[0], s[1], s[2])))).clone();
        let mut edges: std::collections::HashMap<(usize, usize), mt::Edge> = std::collections::HashMap::new();
        let mut out: Vec<mt::Face> = Vec::new();
        for (dir, axis, cs) in faces {
            let (dir, axis) = (turn(dir), turn(axis));
            let mut wire: Vec<mt::Edge> = Vec::new();
            for i in 0..4 {
                let (a, b) = (cs[i], cs[(i + 1) % 4]);
                let (ka, kb) = (key(a), key(b));
                let e = match edges.get(&(ka.min(kb), ka.max(kb))) {
                    Some(e) => e.clone(),
                    None => {
                        let (lo, hi) = if ka < kb { (a, b) } else { (b, a) };
                        let (va, vb) = (vert(lo), vert(hi));
                        let mid = (corner(lo[0], lo[1], lo[2]) + corner(hi[0], hi[1], hi[2])) * 0.5 - center;
                        let transit = center + mid.normalized().unwrap_or(dir) * radius;
                        let e = builder::circle_arc(&va, &vb, p3(transit));
                        edges.insert((ka.min(kb), ka.max(kb)), e.clone());
                        e
                    }
                };
                wire.push(if ka < kb { e } else { e.inverse() });
            }
            // The meridian runs pole to pole through the side away from the patch.
            let (np, sp, back) = (center + axis * radius, center - axis * radius, center - dir * radius);
            let meridian = builder::circle_arc(&builder::vertex(p3(np)), &builder::vertex(p3(sp)), p3(back));
            let mut surface =
                mt::Surface::RevolutedCurve(mt::Processor::new(mt::RevolutedCurve::by_revolution(meridian.oriented_curve(), p3(center), v3(axis))));
            // Outward normal at the patch centre.
            let probe = p3(center + dir * radius);
            use mt::{ParametricSurface3D, SearchParameter};
            if let Some((u, v)) = surface.search_parameter(probe, mt::SPHint2D::None, 100) {
                let n = surface.normal(u, v);
                if mt::InnerSpace::dot(n, v3(dir)) < 0.0 {
                    mt::Invertible::invert(&mut surface);
                }
            }
            out.push(mt::Face::new(vec![wire.into()], surface));
        }
        let shell: mt::Shell = out.into();
        let solid = Solid::try_new(vec![shell]).map_err(|e| KernelError::Failed(format!("sphere: {e}")))?;
        Body::new(solid)
    })
}