brepkit_math/

analytic_intersection.rs

//! Closed-form and semi-analytic intersections of analytic surfaces with planes.
//!
//! Provides specialized intersection algorithms for cylinder, cone, sphere,
//! and torus surfaces with planes, as well as a general marching approach
//! for analytic-analytic surface intersections.

use std::f64::consts::{FRAC_PI_2, TAU};

use crate::MathError;
use crate::curves::{Circle3D, Ellipse3D};
use crate::frame::Frame3;
use crate::nurbs::fitting::interpolate;
use crate::nurbs::intersection::{IntersectionCurve, IntersectionPoint};
use crate::surfaces::{ConicalSurface, CylindricalSurface, SphericalSurface, ToroidalSurface};
use crate::tolerance::Tolerance;
use crate::vec::{Point3, Vec3};

/// Exact curve type resulting from plane-analytic surface intersection.
#[derive(Debug, Clone)]
pub enum ExactIntersectionCurve {
    /// A circle (plane perpendicular to axis of cylinder/cone/sphere).
    Circle(Circle3D),
    /// An ellipse (plane oblique to cylinder/cone axis).
    Ellipse(Ellipse3D),
    /// Fallback to sampled point chain (torus, degenerate cases).
    Points(Vec<Point3>),
}

/// Compute exact intersection curves between a plane and an analytic surface.
///
/// Returns exact `Circle3D` or `Ellipse3D` where possible, falling back to
/// sampled points for complex cases (torus).
///
/// The plane is defined by `dot(normal, p) = d`.
///
/// # Errors
///
/// Returns an error if the intersection computation fails.
pub fn exact_plane_analytic(
    surface: AnalyticSurface<'_>,
    plane_normal: Vec3,
    plane_d: f64,
) -> Result<Vec<ExactIntersectionCurve>, MathError> {
    match surface {
        AnalyticSurface::Cylinder(cyl) => exact_plane_cylinder(cyl, plane_normal, plane_d),
        AnalyticSurface::Sphere(sphere) => exact_plane_sphere(sphere, plane_normal, plane_d),
        AnalyticSurface::Cone(cone) => exact_plane_cone(cone, plane_normal, plane_d),
        AnalyticSurface::Torus(torus) => {
            // Torus intersections are degree-4 — fall back to sampling.
            let chains = sample_plane_torus(torus, plane_normal, plane_d)?;
            Ok(chains
                .into_iter()
                .map(ExactIntersectionCurve::Points)
                .collect())
        }
    }
}

/// Exact plane-cylinder intersection.
///
/// - Plane perpendicular to axis → `Circle3D`
/// - Plane oblique to axis → `Ellipse3D`
/// - Plane parallel to axis → `Points` fallback (0 or 2 lines)
fn exact_plane_cylinder(
    cyl: &CylindricalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<ExactIntersectionCurve>, MathError> {
    let axis = cyl.axis();
    let cos_theta = normal.dot(axis).abs();
    let r = cyl.radius();

    if cos_theta < 1e-10 {
        // Plane parallel to cylinder axis → 0 or 2 line segments.
        // Fall back to sampled points.
        let chains = sample_plane_cylinder(cyl, normal, d)?;
        return Ok(chains
            .into_iter()
            .map(ExactIntersectionCurve::Points)
            .collect());
    }

    // Find where axis intersects the plane: axis_point + t*axis, dot(normal, P) = d
    // t = (d - dot(normal, origin)) / dot(normal, axis)
    let n_dot_axis = normal.dot(axis);
    let n_dot_origin = dot_np(normal, cyl.origin());
    let t = (d - n_dot_origin) / n_dot_axis;
    let center_on_axis = Point3::new(
        cyl.origin().x() + t * axis.x(),
        cyl.origin().y() + t * axis.y(),
        cyl.origin().z() + t * axis.z(),
    );

    if cos_theta > 1.0 - 1e-10 {
        // Plane perpendicular to axis → Circle
        let circle = Circle3D::new(center_on_axis, normal, r)?;
        Ok(vec![ExactIntersectionCurve::Circle(circle)])
    } else {
        // Oblique plane → Ellipse
        // Semi-minor = r (the cylinder radius, unchanged)
        // Semi-major = r / cos(θ) where θ = angle between plane normal and axis
        let semi_minor = r;
        let semi_major = r / cos_theta;

        // The major axis direction lies in the intersection of the plane
        // with the plane containing the axis and the plane normal.
        // It's the projection of the axis onto the cutting plane, normalized.
        let axis_proj = Vec3::new(
            axis.x() - n_dot_axis * normal.x(),
            axis.y() - n_dot_axis * normal.y(),
            axis.z() - n_dot_axis * normal.z(),
        );
        let u_axis = axis_proj.normalize()?;
        let v_axis = normal.cross(u_axis);

        let ellipse = Ellipse3D::with_axes(
            center_on_axis,
            normal,
            semi_major,
            semi_minor,
            u_axis,
            v_axis,
        )?;
        Ok(vec![ExactIntersectionCurve::Ellipse(ellipse)])
    }
}

/// Exact plane-sphere intersection.
///
/// Always produces a `Circle3D` (or empty if no intersection).
fn exact_plane_sphere(
    sphere: &SphericalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<ExactIntersectionCurve>, MathError> {
    let h = dot_np(normal, sphere.center()) - d;
    let r = sphere.radius();

    if h.abs() > r - 1e-10 {
        return Ok(vec![]);
    }

    let circle_r = (r.mul_add(r, -(h * h))).sqrt();
    let circle_center = Point3::new(
        h.mul_add(-normal.x(), sphere.center().x()),
        h.mul_add(-normal.y(), sphere.center().y()),
        h.mul_add(-normal.z(), sphere.center().z()),
    );

    let circle = Circle3D::new(circle_center, normal, circle_r)?;
    Ok(vec![ExactIntersectionCurve::Circle(circle)])
}

/// Exact plane-cone intersection.
///
/// The conic type is set by the cone's half-opening angle from the axis
/// (`γ = π/2 − half_angle`) versus the plane-axis angle `ψ`:
/// - Plane perpendicular to axis (`ψ = π/2`) → `Circle3D`
/// - `ψ > γ` (ellipse) → closed-form `Ellipse3D`
/// - `ψ ≤ γ` (parabola/hyperbola) → bounded single-branch `Points` (one chain
///   per branch — a hyperbola's two nappes never share a chain)
fn exact_plane_cone(
    cone: &ConicalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<ExactIntersectionCurve>, MathError> {
    let axis = cone.axis();
    let cos_theta = normal.dot(axis).abs();
    let half_angle = cone.half_angle();

    if cos_theta > 1.0 - 1e-10 {
        // Plane perpendicular to axis → Circle
        // Find where axis meets the plane
        let n_dot_axis = normal.dot(axis);
        let n_dot_apex = dot_np(normal, cone.apex());
        let t = (d - n_dot_apex) / n_dot_axis;

        // t is the signed distance from apex to plane along the axis.
        // The real cone is a single nappe; the perpendicular-plane section is a
        // circle whose radius follows from the axial offset |t|.
        // |t| ≈ 0 means the plane passes through the apex → degenerate point.
        if t.abs() < 1e-10 {
            return Ok(vec![]);
        }

        let center = Point3::new(
            cone.apex().x() + t * axis.x(),
            cone.apex().y() + t * axis.y(),
            cone.apex().z() + t * axis.z(),
        );
        // half_angle is the angle from the radial plane to the surface.
        // Axial distance t = v * sin(half_angle), so v = t / sin(half_angle).
        // Radius at v = v * cos(half_angle) = t * cos(half_angle) / sin(half_angle).
        let circle_r = t.abs() * half_angle.cos() / half_angle.sin();
        if circle_r < 1e-15 {
            return Ok(vec![]);
        }

        let circle = Circle3D::new(center, normal, circle_r)?;
        return Ok(vec![ExactIntersectionCurve::Circle(circle)]);
    }

    // Oblique plane. Classify the conic in the plane-aligned frame.
    //
    // Decompose the (unit) axis as a = c·n + p·e1, where c = n·a, e1 is the unit
    // in-plane projection of the axis, and p = |projection| = sqrt(1−c²). Write a
    // point Q on the plane as Q = apex + e·n + s·e1 + t·e2 (e = d − n·apex,
    // e2 = n×e1). The cone equation (w·a)² = cos²γ·(w·w) with k = cos²γ =
    // sin²(half_angle) reduces to (no s·t cross term, since e1/e2 align with the
    // conic axes):
    //     (p²−k)·s² + 2ecp·s + e²(c²−k) = k·t²
    // The s² coefficient A = p²−k = sin²θ − sin²(half_angle) sets the type:
    // A < 0 → ellipse, A = 0 → parabola, A > 0 → hyperbola.
    let c = normal.dot(axis);
    let p2 = (1.0 - c * c).max(0.0);
    let p = p2.sqrt();
    let k = half_angle.sin().powi(2);
    let a_coeff = p2 - k;

    // Build the plane-aligned frame e1 (in-plane axis projection), e2 = n×e1.
    let m = Vec3::new(
        axis.x() - c * normal.x(),
        axis.y() - c * normal.y(),
        axis.z() - c * normal.z(),
    );
    let m_len = m.length();
    if m_len < 1e-12 {
        // Axis parallel to normal — handled by the perpendicular branch above;
        // fall back to sampling for safety.
        let chains = sample_plane_cone(cone, normal, d)?;
        return Ok(chains
            .into_iter()
            .map(ExactIntersectionCurve::Points)
            .collect());
    }
    let e1 = m * (1.0 / m_len);
    let e2 = normal.cross(e1);
    let apex = cone.apex();
    let e = d - dot_np(normal, apex);

    // Ellipse → closed form. A = p²−k < 0 with a margin to keep the
    // near-parabolic regime on the robust sampled path.
    if a_coeff < -1e-9 {
        let abs_a = -a_coeff; // = k − p² > 0
        // Real-nappe guard: in the ellipse regime n·g(u) keeps constant sign(c),
        // so v = e/(n·g) ≥ 0 only when e and c share a sign. When e·c < 0 the
        // plane is offset to the far side of the apex from the cone's opening —
        // the section lies entirely on the phantom nappe, so there is no real
        // curve (RHS below is positive regardless of sign, so it can't catch this).
        if e * c < 0.0 {
            return Ok(vec![]);
        }
        // |A|(s − s_c)² + k·t² = RHS, with s_c = ecp/|A| and
        // RHS = e²·k·(1−k)/|A| (always > 0 for a real ellipse).
        let s_c = e * c * p / abs_a;
        let rhs = e * e * k * (1.0 - k) / abs_a;
        if rhs <= 0.0 {
            return Ok(vec![]);
        }
        let semi_s = (rhs / abs_a).sqrt(); // extent along e1
        let semi_t = (rhs / k).sqrt(); // extent along e2
        if semi_s < 1e-12 || semi_t < 1e-12 {
            return Ok(vec![]);
        }
        let center = apex + normal * e + e1 * s_c;
        let (semi_major, semi_minor, u_axis, v_axis) = if semi_s >= semi_t {
            (semi_s, semi_t, e1, e2)
        } else {
            (semi_t, semi_s, e2, e1)
        };
        let ellipse = Ellipse3D::with_axes(center, normal, semi_major, semi_minor, u_axis, v_axis)?;
        return Ok(vec![ExactIntersectionCurve::Ellipse(ellipse)]);
    }

    // Parabola / hyperbola (and the near-parabolic ellipse margin): the section
    // is unbounded, so emit bounded, branch-separated sample chains.
    let chains = sample_plane_cone(cone, normal, d)?;
    Ok(chains
        .into_iter()
        .map(ExactIntersectionCurve::Points)
        .collect())
}

/// Reference to an analytic surface for intersection dispatch.
#[derive(Clone, Copy)]
pub enum AnalyticSurface<'a> {
    /// Cylindrical surface reference.
    Cylinder(&'a CylindricalSurface),
    /// Conical surface reference.
    Cone(&'a ConicalSurface),
    /// Spherical surface reference.
    Sphere(&'a SphericalSurface),
    /// Toroidal surface reference.
    Torus(&'a ToroidalSurface),
}

/// Compute `n . p` treating a `Point3` as a position vector.
fn dot_np(n: Vec3, p: Point3) -> f64 {
    n.dot(Vec3::new(p.x(), p.y(), p.z()))
}

/// Intersect a plane with an analytic surface.
///
/// The plane is defined by `dot(normal, p) = d`.
///
/// # Errors
///
/// Returns an error if the intersection computation fails.
pub fn intersect_plane_analytic(
    surface: AnalyticSurface<'_>,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, MathError> {
    match surface {
        AnalyticSurface::Cylinder(cyl) => intersect_plane_cylinder(cyl, normal, d),
        AnalyticSurface::Cone(cone) => intersect_plane_cone(cone, normal, d),
        AnalyticSurface::Sphere(sphere) => intersect_plane_sphere(sphere, normal, d),
        AnalyticSurface::Torus(torus) => intersect_plane_torus(torus, normal, d),
    }
}

/// Sample points on the plane-analytic intersection without NURBS curve fitting.
///
/// Returns chains of ordered 3D sample points. Each chain is one connected
/// component of the intersection curve. This is much faster than
/// `intersect_plane_analytic` when only sample points are needed (e.g. for
/// boolean intersection segment generation).
///
/// # Errors
///
/// Returns an error if the intersection computation fails.
pub fn sample_plane_analytic(
    surface: AnalyticSurface<'_>,
    normal: Vec3,
    d: f64,
) -> Result<Vec<Vec<Point3>>, MathError> {
    match surface {
        AnalyticSurface::Cylinder(cyl) => sample_plane_cylinder(cyl, normal, d),
        AnalyticSurface::Cone(cone) => sample_plane_cone(cone, normal, d),
        AnalyticSurface::Sphere(sphere) => sample_plane_sphere(sphere, normal, d),
        AnalyticSurface::Torus(torus) => sample_plane_torus(torus, normal, d),
    }
}

/// Sample the plane-cylinder intersection as ordered 3D points.
#[allow(clippy::cast_precision_loss, clippy::unnecessary_wraps)]
fn sample_plane_cylinder(
    cyl: &CylindricalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<Vec<Point3>>, MathError> {
    let n_samples = 64_usize;
    let mut points = Vec::with_capacity(n_samples + 1);

    for i in 0..=n_samples {
        let u = TAU * (i as f64) / (n_samples as f64);
        let base = cyl.evaluate(u, 0.0);
        let n_dot_axis = normal.dot(cyl.axis());
        let n_dot_base = dot_np(normal, base);

        if n_dot_axis.abs() < 1e-12 {
            if (n_dot_base - d).abs() < 1e-6 {
                points.push(base);
            }
        } else {
            let v = (d - n_dot_base) / n_dot_axis;
            if v.abs() <= 100.0 {
                points.push(cyl.evaluate(u, v));
            }
        }
    }

    if points.len() < 2 {
        Ok(vec![])
    } else {
        Ok(vec![points])
    }
}

/// Sample the plane-sphere intersection as ordered 3D points.
#[allow(clippy::cast_precision_loss)]
fn sample_plane_sphere(
    sphere: &SphericalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<Vec<Point3>>, MathError> {
    let h = dot_np(normal, sphere.center()) - d;
    let r = sphere.radius();

    if h.abs() > r - 1e-10 {
        return Ok(vec![]);
    }

    let circle_r = (r.mul_add(r, -(h * h))).sqrt();
    let circle_center = Point3::new(
        h.mul_add(-normal.x(), sphere.center().x()),
        h.mul_add(-normal.y(), sphere.center().y()),
        h.mul_add(-normal.z(), sphere.center().z()),
    );

    let basis = Frame3::from_normal(circle_center, normal)?;
    let u_dir = basis.x;
    let v_dir = basis.y;

    let n_samples = 64_usize;
    let mut points = Vec::with_capacity(n_samples + 1);

    for i in 0..=n_samples {
        let theta = TAU * (i as f64) / (n_samples as f64);
        let (sin_t, cos_t) = theta.sin_cos();
        points.push(circle_center + u_dir * (circle_r * cos_t) + v_dir * (circle_r * sin_t));
    }

    Ok(vec![points])
}

/// Sample the plane-cone intersection as ordered 3D points.
///
/// The cone is the single real nappe `v >= 0` of `P(u,v) = apex + v·g(u)`.
/// Along each generator `g(u)` the plane `n·P = d` is linear in `v`, so
/// `v = (d − n·apex) / (n·g(u))`. We keep only `v >= 0` (the phantom `v < 0`
/// nappe is geometrically absent) and `v` below a finite bound (near an
/// asymptote `n·g(u) → 0` so `v → ∞` — those points run off the surface and
/// must be excluded). The angular samples that survive form one contiguous arc
/// (ellipse) or two (parabola/hyperbola, one per branch); each contiguous run
/// is returned as a separate ordered chain so the consumer never stitches two
/// disjoint branches into one curve.
#[allow(clippy::cast_precision_loss, clippy::unnecessary_wraps)]
fn sample_plane_cone(
    cone: &ConicalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<Vec<Point3>>, MathError> {
    let apex = cone.apex();
    let n_dot_apex = dot_np(normal, apex);
    let e = d - n_dot_apex;

    // Per-generator solve: along g(u) the plane is linear in v, v = e / (n·g(u)).
    // Sample u densely; keep only the real nappe (v >= 0) and skip near-asymptote
    // generators (n·g(u) ≈ 0 → v → ∞).
    let n_samples = 512_usize;
    let mut vs: Vec<Option<f64>> = Vec::with_capacity(n_samples);
    let mut v_min = f64::INFINITY;
    for i in 0..n_samples {
        let u = TAU * (i as f64) / (n_samples as f64);
        let g = cone.evaluate(u, 1.0) - apex;
        let n_dot_g = normal.dot(Vec3::new(g.x(), g.y(), g.z()));
        if n_dot_g.abs() < 1e-12 {
            vs.push(None);
            continue;
        }
        let v = e / n_dot_g;
        if v >= -1e-12 {
            let v = v.max(0.0);
            v_min = v_min.min(v);
            vs.push(Some(v));
        } else {
            vs.push(None);
        }
    }

    if !v_min.is_finite() {
        return Ok(Vec::new());
    }

    // Bound the arc around the conic vertex (closest approach to the apex, at
    // v_min). An ellipse is naturally bounded; a parabola/hyperbola is not, so
    // cap the cone radius at a generous multiple of the vertex radius. This is
    // scale-invariant and centred on where any finite cone face's overlap lies;
    // the downstream consumer trims the fitted curve to the actual face AABB, so
    // over-coverage is harmless. The floor handles a vertex at the apex (v_min≈0).
    let v_max = (8.0 * v_min).max(v_min + 4.0);

    // Per-sample v within the cap; the raw values stay in `vs` for the
    // boundary solve below.
    let kept: Vec<Option<f64>> = vs.iter().map(|v| v.filter(|&v| v <= v_max)).collect();

    let point_at = |u: f64, v: f64| -> Point3 {
        let g = cone.evaluate(u, 1.0) - apex;
        apex + g * v
    };
    #[allow(clippy::cast_precision_loss)]
    let u_of = |i: usize| TAU * (i as f64) / (n_samples as f64);
    let n_dot_g_at = |u: f64| -> f64 {
        let g = cone.evaluate(u, 1.0) - apex;
        normal.dot(Vec3::new(g.x(), g.y(), g.z()))
    };

    if kept.iter().all(Option::is_some) {
        // Closed loop (ellipse regime): emit all points and repeat the first.
        let mut pts: Vec<Point3> = kept
            .iter()
            .enumerate()
            .filter_map(|(i, v)| v.map(|v| point_at(u_of(i), v)))
            .collect();
        if let Some(&first) = pts.first() {
            pts.push(first);
        }
        return Ok(vec![pts]);
    }

    // A hyperbola/parabola tail diverges as 1/(n·g), so between the last kept
    // sample and its dropped neighbour v can leap far past `v_max` in one
    // uniform-u pitch — and any finite face window inside that leap is lost
    // (a taper cone grazed 0.05 by a prism plane lost its entire 0.5-tall
    // section to exactly this aliasing). Extend each run end to the exact
    // `v_max` boundary: bisect u for `n·g(u) = e/v_max` inside the dropped
    // pitch (n·g is monotone there — its extrema sit at the conic vertex,
    // far from any asymptote), then fill the tail with uniform-u samples.
    let tail = |i_end: usize, forward: bool, kept: &[Option<f64>]| -> Vec<Point3> {
        let Some(v_end) = kept[i_end] else {
            return Vec::new();
        };
        let u_end = u_of(i_end);
        #[allow(clippy::cast_precision_loss)]
        let pitch = TAU / (n_samples as f64);
        let u_next = if forward {
            u_end + pitch
        } else {
            u_end - pitch
        };
        let target = e / v_max;
        let h_end = n_dot_g_at(u_end) - target;
        let h_next = n_dot_g_at(u_next) - target;
        if v_end >= v_max || h_end == 0.0 || h_end.signum() == h_next.signum() {
            return Vec::new();
        }
        let (mut lo, mut hi) = (u_end, u_next);
        for _ in 0..60 {
            let mid = f64::midpoint(lo, hi);
            if (n_dot_g_at(mid) - target).signum() == h_end.signum() {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        let u_star = f64::midpoint(lo, hi);
        let tail_n = 8_usize;
        (1..=tail_n)
            .filter_map(|k| {
                #[allow(clippy::cast_precision_loss)]
                let u = u_end + (u_star - u_end) * (k as f64) / (tail_n as f64);
                let ng = n_dot_g_at(u);
                if ng.abs() < 1e-12 {
                    return None;
                }
                let v = e / ng;
                (v >= -1e-12 && v <= v_max * (1.0 + 1e-9)).then(|| point_at(u, v.max(0.0)))
            })
            .collect()
    };

    // Split into contiguous runs of kept samples, treating the array as
    // circular (rotate past a gap) so a branch straddling u=0 stays whole.
    let gap = kept.iter().position(Option::is_none).unwrap_or(0);
    let mut chains: Vec<Vec<Point3>> = Vec::new();
    let mut run: Vec<usize> = Vec::new();
    let flush = |run: &mut Vec<usize>, chains: &mut Vec<Vec<Point3>>| {
        if run.len() >= 2 {
            let first = run[0];
            let last = run[run.len() - 1];
            let mut pts: Vec<Point3> = tail(first, false, &kept);
            pts.reverse();
            pts.extend(
                run.iter()
                    .filter_map(|&i| kept[i].map(|v| point_at(u_of(i), v))),
            );
            pts.extend(tail(last, true, &kept));
            chains.push(pts);
        }
        run.clear();
    };
    for k in 0..n_samples {
        let idx = (gap + k) % n_samples;
        if kept[idx].is_some() {
            run.push(idx);
        } else {
            flush(&mut run, &mut chains);
        }
    }
    flush(&mut run, &mut chains);
    Ok(chains.into_iter().filter(|c| c.len() >= 2).collect())
}

/// Sample the plane-torus intersection as ordered 3D points.
///
/// Uses the same closed-form crossings and chaining as `intersect_plane_torus`
/// but skips NURBS curve fitting (the callers here only need the points).
#[allow(clippy::unnecessary_wraps)] // sibling match-arms and `?` callers need `Result`
fn sample_plane_torus(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<Vec<Point3>>, MathError> {
    let crossing_pts = plane_torus_crossings(torus, normal, d, 128);
    Ok(chain_torus_crossings(&crossing_pts)
        .into_iter()
        .map(|run| run.into_iter().map(|p| p.point).collect())
        .collect())
}

/// Intersect a plane with a cylindrical surface.
///
/// For each `u` in `[0, 2pi)`, the cylinder point is linear in `v`,
/// so the plane equation `dot(normal, P(u,v)) = d` is linear in `v`
/// and can be solved directly.
///
/// # Errors
///
/// Returns an error if curve fitting fails.
#[allow(clippy::cast_precision_loss)]
pub fn intersect_plane_cylinder(
    cyl: &CylindricalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, MathError> {
    let n_samples = 64_usize;
    let mut points_3d = Vec::new();
    let mut ipoints = Vec::new();

    for i in 0..=n_samples {
        let u = TAU * (i as f64) / (n_samples as f64);
        // P(u, v) = origin + r*(cos(u)*x + sin(u)*y) + v*axis
        // dot(normal, P) = d  =>  dot(normal, base(u)) + v * dot(normal, axis) = d
        let base = cyl.evaluate(u, 0.0);
        let n_dot_axis = normal.dot(cyl.axis());
        let n_dot_base = dot_np(normal, base);

        if n_dot_axis.abs() < 1e-12 {
            // Plane parallel to axis -- check if base is on plane.
            if (n_dot_base - d).abs() < 1e-6 {
                let pt = base;
                points_3d.push(pt);
                ipoints.push(IntersectionPoint {
                    point: pt,
                    param1: (u, 0.0),
                    param2: (0.0, 0.0),
                });
            }
        } else {
            let v = (d - n_dot_base) / n_dot_axis;
            // Only keep points within a reasonable v range.
            if v.abs() <= 100.0 {
                let pt = cyl.evaluate(u, v);
                points_3d.push(pt);
                ipoints.push(IntersectionPoint {
                    point: pt,
                    param1: (u, v),
                    param2: (0.0, 0.0),
                });
            }
        }
    }

    build_curves_from_points(&points_3d, ipoints)
}

/// Intersect a plane with a spherical surface.
///
/// The intersection of a plane with a sphere is a circle (or empty/point).
/// Computes the circle center, radius, and samples points on it.
///
/// # Errors
///
/// Returns an error if curve fitting fails.
#[allow(clippy::cast_precision_loss)]
pub fn intersect_plane_sphere(
    sphere: &SphericalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, MathError> {
    let h = dot_np(normal, sphere.center()) - d;
    let r = sphere.radius();

    // No intersection if plane is too far from center.
    if h.abs() > r - 1e-10 {
        return Ok(vec![]);
    }

    let circle_r = (r.mul_add(r, -(h * h))).sqrt();
    let circle_center = Point3::new(
        h.mul_add(-normal.x(), sphere.center().x()),
        h.mul_add(-normal.y(), sphere.center().y()),
        h.mul_add(-normal.z(), sphere.center().z()),
    );

    // Build a local frame on the plane.
    let basis = Frame3::from_normal(circle_center, normal)?;
    let u_dir = basis.x;
    let v_dir = basis.y;

    let n_samples = 64_usize;
    let mut points_3d = Vec::new();
    let mut ipoints = Vec::new();

    for i in 0..=n_samples {
        let theta = TAU * (i as f64) / (n_samples as f64);
        let (sin_t, cos_t) = theta.sin_cos();
        let pt = circle_center + u_dir * (circle_r * cos_t) + v_dir * (circle_r * sin_t);
        points_3d.push(pt);
        ipoints.push(IntersectionPoint {
            point: pt,
            param1: (theta, 0.0),
            param2: (0.0, 0.0),
        });
    }

    build_curves_from_points(&points_3d, ipoints)
}

/// Intersect a plane with a conical surface.
///
/// Like a cylinder, the cone is linear along each generatrix, so the plane
/// equation is linear in `v` for each fixed `u`.
///
/// # Errors
///
/// Returns an error if curve fitting fails.
#[allow(clippy::cast_precision_loss)]
pub fn intersect_plane_cone(
    cone: &ConicalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, MathError> {
    let n_samples = 64_usize;
    let mut points_3d = Vec::new();
    let mut ipoints = Vec::new();

    for i in 0..n_samples {
        let u = TAU * (i as f64) / (n_samples as f64);
        // P(u, v) = apex + v * dir(u)
        // dot(normal, apex) + v * dot(normal, dir(u)) = d
        let apex = cone.apex();
        let n_dot_apex = dot_np(normal, apex);
        // dir(u) = P(u,1) - apex
        let p1 = cone.evaluate(u, 1.0);
        let dir = p1 - apex;
        let n_dot_dir = normal.dot(dir);

        if n_dot_dir.abs() < 1e-12 {
            continue;
        }

        let v = (d - n_dot_apex) / n_dot_dir;
        // Allow negative v — the cone surface extends in both directions from the apex.
        if v.abs() > 1e-10 && v.abs() < 100.0 {
            let pt = cone.evaluate(u, v);
            points_3d.push(pt);
            ipoints.push(IntersectionPoint {
                point: pt,
                param1: (u, v),
                param2: (0.0, 0.0),
            });
        }
    }

    build_curves_from_points(&points_3d, ipoints)
}

/// Intersect a plane with a toroidal surface.
///
/// The section is a degree-4 curve, but for each `v` the `u` values solve in
/// closed form (see `plane_torus_crossings`), so it is sampled by a v-scan
/// and each connected loop is fitted to a NURBS curve.
///
/// # Errors
///
/// Never returns an error today (curve-fit failures drop the affected loop);
/// the `Result` is kept for signature parity with the other plane-analytic
/// intersectors.
#[allow(clippy::unnecessary_wraps)]
pub fn intersect_plane_torus(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, MathError> {
    // The section satisfies a per-v closed form (see `plane_torus_crossings`),
    // so scan v and solve u directly instead of a 2D sign-change grid with
    // Newton refinement: O(n) rather than O(n²), and every point is exact.
    let crossing_pts = plane_torus_crossings(torus, normal, d, 128);

    let mut curves = Vec::new();
    for ipts in chain_torus_crossings(&crossing_pts) {
        let pts: Vec<Point3> = ipts.iter().map(|p| p.point).collect();
        if let Ok(curve) = interpolate(&pts, 3.min(pts.len() - 1)) {
            curves.push(IntersectionCurve {
                curve,
                points: ipts,
            });
        }
    }

    Ok(curves)
}

/// Greedy nearest-neighbour chaining of torus-plane crossing points into
/// closed section loops. Runs shorter than four points are dropped.
///
/// Plane × full torus is always a set of CLOSED loops, but the greedy walk
/// stops one step short of closing (the first point is already `used`, so it
/// is never re-added and the last point sits ~one step from the start).
/// A loop whose end-to-start gap is within ~2 point-spacings is closed by
/// repeating its first point, so a fitted NURBS closes exactly and downstream
/// consumers see a closed curve. A fragmented chain (greedy walk broke a loop
/// at a near-tangency) ends far from its start and is left open — it must not
/// be force-closed into a wrong loop.
fn chain_torus_crossings(crossing_pts: &[(f64, f64, Point3)]) -> Vec<Vec<IntersectionPoint>> {
    let mut used = vec![false; crossing_pts.len()];
    let mut runs = Vec::new();

    for start in 0..crossing_pts.len() {
        if used[start] {
            continue;
        }
        used[start] = true;
        let mut chain = vec![start];

        loop {
            let last = chain[chain.len() - 1];
            let last_pt = crossing_pts[last].2;
            let mut best_idx = None;
            let mut best_dist = 1.0_f64;

            for (j, &is_used) in used.iter().enumerate() {
                if is_used {
                    continue;
                }
                let dist = (crossing_pts[j].2 - last_pt).length();
                if dist < best_dist {
                    best_dist = dist;
                    best_idx = Some(j);
                }
            }

            if let Some(j) = best_idx {
                used[j] = true;
                chain.push(j);
            } else {
                break;
            }
        }

        if chain.len() < 4 {
            continue;
        }
        let mut ipts: Vec<IntersectionPoint> = chain
            .iter()
            .map(|&i| IntersectionPoint {
                point: crossing_pts[i].2,
                param1: (crossing_pts[i].0, crossing_pts[i].1),
                param2: (0.0, 0.0),
            })
            .collect();

        let closing_gap = (ipts[ipts.len() - 1].point - ipts[0].point).length();
        let median_spacing = {
            let mut spac: Vec<f64> = ipts
                .windows(2)
                .map(|w| (w[1].point - w[0].point).length())
                .collect();
            spac.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
            spac.get(spac.len() / 2).copied().unwrap_or(0.0)
        };
        // Wrap when the closing gap is within ~2 point-spacings (measured
        // ratio ≈ 1.0 for the census ovals) and not already coincident — but
        // only for a SIMPLE loop. A self-touching section (the inner/outer
        // tangent figure-eight) is traced as one chain that folds back through
        // its node and also ends near its start; sealing it would misrepresent
        // a non-manifold singularity as a closed loop, so leave it open.
        if closing_gap > 1e-9
            && median_spacing > 1e-12
            && closing_gap <= 2.0 * median_spacing
            && !chain_self_touches(&ipts, median_spacing)
        {
            ipts.push(ipts[0]);
        }
        runs.push(ipts);
    }

    runs
}

/// Whether a chain folds back on itself in its interior — the signature of a
/// self-touching section (a tangent figure-eight), as opposed to a simple
/// loop whose only near-return is the intended closure at its two ends.
///
/// Checks whether two chain points far apart in index (and both away from the
/// endpoints, so the closure region is excluded) come within ~1.5 spacings of
/// each other. A convex/simple oval never does; a figure-eight does, at its
/// node. Only called when a chain already looks closeable, so the O(m²) scan
/// is rare.
fn chain_self_touches(ipts: &[IntersectionPoint], median_spacing: f64) -> bool {
    let m = ipts.len();
    let k = (m / 4).clamp(1, 6);
    if m < 3 * k || median_spacing <= 0.0 {
        return false;
    }
    let thresh = median_spacing * 1.5;
    for i in k..(m - k) {
        for j in (i + k)..(m - k) {
            if (ipts[i].point - ipts[j].point).length() < thresh {
                return true;
            }
        }
    }
    false
}

/// Closed-form `(u, v, point)` crossings of a plane with a torus.
///
/// In the torus's own frame let `a = n·X`, `b = n·Y`, `c = n·Z`,
/// `s = hypot(a, b)`, `phi = atan2(b, a)`. Substituting the torus
/// parameterization into `n·P = d` gives
///   `(R + r·cos v)·s·cos(u − phi) + r·c·sin v = d − n·center`,
/// so for each `v` the two `u` branches solve directly as
/// `u = phi ± acos((d − n·center − r·c·sin v) / (s·(R + r·cos v)))`.
/// Scanning `v` at `n_v` samples replaces a 2D sign-change grid plus Newton
/// refinement: each point is `torus.evaluate(u, v)` (on the torus by
/// construction) with `u` solved so it lies on the plane to floating-point
/// precision, so no iterative refinement is needed.
///
/// When `s ≈ 0` the plane is perpendicular to the axis and the section is up
/// to two full circles at the `v` values solving `r·c·sin v = d − n·center`;
/// those are sampled by scanning `u`.
#[allow(clippy::cast_precision_loss)]
fn plane_torus_crossings(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
    n_v: usize,
) -> Vec<(f64, f64, Point3)> {
    let big_r = torus.major_radius();
    let small_r = torus.minor_radius();
    let a = normal.dot(torus.x_axis());
    let b = normal.dot(torus.y_axis());
    let c = normal.dot(torus.z_axis());
    let s = a.hypot(b);
    let phi = b.atan2(a);
    let d_local = d - dot_np(normal, torus.center());

    let mut pts: Vec<(f64, f64, Point3)> = Vec::new();

    // Plane perpendicular to the axis: the section is up to two full circles.
    if s < 1e-12 {
        if c.abs() < 1e-12 {
            return pts;
        }
        let sin_v = d_local / (small_r * c);
        if sin_v.abs() > 1.0 + 1e-9 {
            return pts;
        }
        let v0 = sin_v.clamp(-1.0, 1.0).asin();
        let v1 = std::f64::consts::PI - v0;
        let mut vs = vec![v0];
        // Skip the mirror circle when the plane is tangent (v0 == v1).
        if (v1 - v0).abs() > 1e-9 {
            vs.push(v1);
        }
        for v in vs {
            for i in 0..n_v {
                let u = TAU * (i as f64) / (n_v as f64);
                pts.push((u, v, torus.evaluate(u, v)));
            }
        }
        return pts;
    }

    // General plane: scan v, solve the two u branches per v. Offset the scan
    // by half a step so it never lands exactly on a tangency node (e.g. the
    // inner-tangent figure-eight at v = π, where the two u branches collapse
    // to one point) — a coincident node lets greedy chaining thread through
    // and wrongly seal a self-touching section into a closed loop.
    let v_off = TAU / (n_v as f64) * 0.5;
    for i in 0..n_v {
        let v = (i as f64).mul_add(TAU / (n_v as f64), v_off);
        let tube_r = small_r.mul_add(v.cos(), big_r); // R + r·cos v > 0
        let rhs = (d_local - small_r * c * v.sin()) / (s * tube_r);
        if rhs.abs() > 1.0 {
            continue;
        }
        let delta = rhs.clamp(-1.0, 1.0).acos();
        for u in [phi + delta, phi - delta] {
            pts.push((u, v, torus.evaluate(u, v)));
        }
    }
    pts
}

/// Real intersection parameters `t` of the line `origin + t·dir` with a torus.
///
/// A line meets a torus in up to four points (degree-4). Substituting the line
/// into the torus implicit `(a² + b² + c² + R² − r²)² = 4R²(a² + b²)` — where
/// `(a, b, c)` are the line point's coordinates in the torus frame — gives a
/// quartic in `t`, solved here for its real roots (each refined by one Newton
/// step against the implicit). `dir` need not be unit length; `t` is in units of
/// `dir`. Returns the roots sorted ascending (0–4 of them).
///
/// Used by the boolean section trimmer to find where a plane×torus oval exits a
/// box face's straight boundary edge — the exact crossing shared by the two
/// adjacent faces, which is what makes the notch watertight.
#[must_use]
pub fn intersect_line_torus(torus: &ToroidalSurface, origin: Point3, dir: Vec3) -> Vec<f64> {
    let c = torus.center();
    let (xa, ya, za) = (torus.x_axis(), torus.y_axis(), torus.z_axis());
    let big_r = torus.major_radius();
    let small_r = torus.minor_radius();

    // Line point in torus frame: a(t)=a0+a1 t, b(t)=b0+b1 t, c(t)=c0+c1 t.
    let o = Vec3::new(origin.x() - c.x(), origin.y() - c.y(), origin.z() - c.z());
    let (a0, a1) = (xa.dot(o), xa.dot(dir));
    let (b0, b1) = (ya.dot(o), ya.dot(dir));
    let (c0, c1) = (za.dot(o), za.dot(dir));

    // G(t) = a² + b² + c² + R² − r²  (quadratic: g2 t² + g1 t + g0)
    let g2 = a1.mul_add(a1, b1.mul_add(b1, c1 * c1));
    let g1 = 2.0 * a1.mul_add(a0, b1.mul_add(b0, c1 * c0));
    let g0 = a0.mul_add(
        a0,
        b0.mul_add(b0, c0.mul_add(c0, big_r.mul_add(big_r, -small_r * small_r))),
    );

    // H(t) = 4R² (a² + b²)  (quadratic: h2 t² + h1 t + h0)
    let four_rr = 4.0 * big_r * big_r;
    let h2 = four_rr * a1.mul_add(a1, b1 * b1);
    let h1 = four_rr * (2.0 * a1.mul_add(a0, b1 * b0));
    let h0 = four_rr * a0.mul_add(a0, b0 * b0);

    // Quartic G² − H = 0:  e4 t⁴ + e3 t³ + e2 t² + e1 t + e0.
    let e4 = g2 * g2;
    let e3 = 2.0 * g2 * g1;
    let e2 = g1.mul_add(g1, 2.0 * g2 * g0) - h2;
    let e1 = 2.0f64.mul_add(g1 * g0, -h1);
    let e0 = g0.mul_add(g0, -h0);

    let mut roots = real_roots_quartic(e4, e3, e2, e1, e0);
    // One Newton polish against the torus implicit for full precision.
    let impl_f = |t: f64| -> f64 {
        let p = origin + dir * t;
        let q = Vec3::new(p.x() - c.x(), p.y() - c.y(), p.z() - c.z());
        let (a, b, cc) = (xa.dot(q), ya.dot(q), za.dot(q));
        (a.hypot(b) - big_r).hypot(cc) - small_r
    };
    for t in &mut roots {
        let eps = 1e-7;
        let f = impl_f(*t);
        let df = (impl_f(*t + eps) - impl_f(*t - eps)) / (2.0 * eps);
        if df.abs() > 1e-12 {
            *t -= f / df;
        }
    }
    roots.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    roots
}

/// Real roots of `c4 x⁴ + c3 x³ + c2 x² + c1 x + c0` via Durand–Kerner, falling
/// back to the lower-degree solvers when the leading coefficients vanish.
fn real_roots_quartic(c4: f64, c3: f64, c2: f64, c1: f64, c0: f64) -> Vec<f64> {
    // Degenerate leading coefficient → lower degree.
    if c4.abs() < 1e-14 {
        return real_roots_cubic(c3, c2, c1, c0);
    }
    // Monic: x⁴ + a x³ + b x² + c x + d.
    let (a, b, c, d) = (c3 / c4, c2 / c4, c1 / c4, c0 / c4);
    let eval = |z: Complex| -> Complex {
        // Horner.
        let mut acc = Complex::new(1.0, 0.0);
        acc = acc * z + Complex::new(a, 0.0);
        acc = acc * z + Complex::new(b, 0.0);
        acc = acc * z + Complex::new(c, 0.0);
        acc * z + Complex::new(d, 0.0)
    };
    // Durand–Kerner: four roots seeded on a circle, iterated to convergence.
    let seed = Complex::new(0.4, 0.9);
    let mut r = [
        Complex::new(1.0, 0.0),
        seed,
        seed * seed,
        seed * seed * seed,
    ];
    for _ in 0..100 {
        let mut max_step = 0.0_f64;
        for i in 0..4 {
            let mut denom = Complex::new(1.0, 0.0);
            for j in 0..4 {
                if i != j {
                    denom = denom * (r[i] - r[j]);
                }
            }
            if denom.norm() < 1e-300 {
                continue;
            }
            let step = eval(r[i]) / denom;
            r[i] = r[i] - step;
            max_step = max_step.max(step.norm());
        }
        if max_step < 1e-14 {
            break;
        }
    }
    // Keep roots with negligible imaginary part AND a small REAL-polynomial
    // residual — Durand–Kerner stops after a fixed iteration cap whether or not
    // it converged, so a non-converged iterate could otherwise be returned as a
    // spurious root. Evaluate the monic quartic at each candidate (real part) and
    // keep only |p(x)| below a magnitude-scaled tolerance; de-dup near-equal
    // roots (a double root converges to two near-identical iterates).
    let p_real = |x: f64| -> f64 { (((x + a) * x + b) * x + c) * x + d };
    let mut out: Vec<f64> = Vec::new();
    for z in r {
        if z.im.abs() >= 1e-7 {
            continue;
        }
        let x = z.re;
        // Residual tolerance scales with the polynomial's coefficient magnitude
        // and |x|^4 so large-coefficient quartics are not over-rejected.
        let scale = 1.0 + a.abs() + b.abs() + c.abs() + d.abs() + x.abs().powi(4);
        if p_real(x).abs() > 1e-6 * scale {
            continue;
        }
        if out.iter().any(|&y| (y - x).abs() < 1e-9 * (1.0 + x.abs())) {
            continue;
        }
        out.push(x);
    }
    out
}

/// Real roots of `a x³ + b x² + c x + d` (Cardano), with quadratic fallback.
fn real_roots_cubic(a: f64, b: f64, c: f64, d: f64) -> Vec<f64> {
    if a.abs() < 1e-14 {
        return real_roots_quadratic(b, c, d);
    }
    // Depressed cubic t³ + p t + q via x = t − b/(3a).
    let (b, c, d) = (b / a, c / a, d / a);
    let p = c - b * b / 3.0;
    let q = 2.0 * b * b * b / 27.0 - b * c / 3.0 + d;
    let shift = -b / 3.0;
    let disc = q * q / 4.0 + p * p * p / 27.0;
    if disc > 1e-14 {
        let sq = disc.sqrt();
        let u = (-q / 2.0 + sq).cbrt();
        let v = (-q / 2.0 - sq).cbrt();
        vec![u + v + shift]
    } else if disc < -1e-14 {
        // Three real roots (trigonometric).
        let m = 2.0 * (-p / 3.0).sqrt();
        let theta = (3.0 * q / (p * m)).clamp(-1.0, 1.0).acos() / 3.0;
        (0..3)
            .map(|k| {
                m.mul_add(
                    (theta - 2.0 * std::f64::consts::PI * f64::from(k) / 3.0).cos(),
                    shift,
                )
            })
            .collect()
    } else {
        // Repeated roots.
        let u = (-q / 2.0).cbrt();
        vec![2.0 * u + shift, -u + shift]
    }
}

/// Real roots of `a x² + b x + c`, with linear fallback.
fn real_roots_quadratic(a: f64, b: f64, c: f64) -> Vec<f64> {
    if a.abs() < 1e-14 {
        if b.abs() < 1e-14 {
            return Vec::new();
        }
        return vec![-c / b];
    }
    let disc = b * b - 4.0 * a * c;
    if disc < 0.0 {
        Vec::new()
    } else {
        let sq = disc.sqrt();
        vec![(-b - sq) / (2.0 * a), (-b + sq) / (2.0 * a)]
    }
}

/// Minimal complex number for the quartic root finder.
#[derive(Clone, Copy)]
struct Complex {
    re: f64,
    im: f64,
}

impl Complex {
    const fn new(re: f64, im: f64) -> Self {
        Self { re, im }
    }
    fn norm(self) -> f64 {
        self.re.hypot(self.im)
    }
}

impl std::ops::Add for Complex {
    type Output = Self;
    fn add(self, o: Self) -> Self {
        Self::new(self.re + o.re, self.im + o.im)
    }
}

impl std::ops::Sub for Complex {
    type Output = Self;
    fn sub(self, o: Self) -> Self {
        Self::new(self.re - o.re, self.im - o.im)
    }
}

impl std::ops::Mul for Complex {
    type Output = Self;
    fn mul(self, o: Self) -> Self {
        Self::new(
            self.re.mul_add(o.re, -(self.im * o.im)),
            self.re.mul_add(o.im, self.im * o.re),
        )
    }
}

impl std::ops::Div for Complex {
    type Output = Self;
    fn div(self, o: Self) -> Self {
        let den = o.re.mul_add(o.re, o.im * o.im);
        Self::new(
            self.re.mul_add(o.re, self.im * o.im) / den,
            self.im.mul_add(o.re, -(self.re * o.im)) / den,
        )
    }
}

/// Build intersection curves from a collection of ordered 3D points.
///
/// If there are enough points, fits a NURBS curve through them.
fn build_curves_from_points(
    points_3d: &[Point3],
    ipoints: Vec<IntersectionPoint>,
) -> Result<Vec<IntersectionCurve>, MathError> {
    if points_3d.len() < 2 {
        return Ok(vec![]);
    }

    let degree = 3.min(points_3d.len() - 1);
    let curve = interpolate(points_3d, degree)?;
    Ok(vec![IntersectionCurve {
        curve,
        points: ipoints,
    }])
}

// -- Analytic-Analytic Intersection -------------------------------------------

/// Intersect two analytic surfaces using a general marching approach.
///
/// Seeds intersection points by sampling both parameter spaces on a grid,
/// then marches along the intersection curve using the cross product of
/// the two surface normals as the tangent direction.
///
/// # Errors
///
/// Returns an error if curve fitting fails.
#[allow(
    clippy::cast_precision_loss,
    clippy::too_many_lines,
    clippy::similar_names,
    clippy::unnecessary_wraps,
    clippy::type_complexity
)]
pub fn intersect_analytic_analytic(
    a: AnalyticSurface<'_>,
    b: AnalyticSurface<'_>,
    grid_res: usize,
) -> Result<Vec<IntersectionCurve>, MathError> {
    intersect_analytic_analytic_bounded(a, b, grid_res, None, None)
}

/// Intersect two analytic surfaces with optional v-range overrides.
///
/// When `v_range_hint_a` or `v_range_hint_b` is `Some((min, max))`, the
/// marching algorithm searches that v-range instead of the hardcoded default.
/// This is essential for cylinders and cones whose default v-range is small
/// (-1..1 or 0.01..2) but whose actual face may extend much further.
///
/// # Errors
///
/// Returns `MathError` if algebraic intersection fails or marching diverges.
pub fn intersect_analytic_analytic_bounded(
    a: AnalyticSurface<'_>,
    b: AnalyticSurface<'_>,
    grid_res: usize,
    v_range_hint_a: Option<(f64, f64)>,
    v_range_hint_b: Option<(f64, f64)>,
) -> Result<Vec<IntersectionCurve>, MathError> {
    // Try algebraic specialization for known surface pairs before falling
    // back to the general marching approach.
    if let Some(result) = try_algebraic_intersection(&a, &b, v_range_hint_a, v_range_hint_b)? {
        return Ok(result);
    }

    let (surf_a, norm_a, u_range_a, default_v_a) = surface_closures(&a);
    let (surf_b, norm_b, u_range_b, default_v_b) = surface_closures(&b);
    let v_range_a = v_range_hint_a.unwrap_or(default_v_a);
    let v_range_b = v_range_hint_b.unwrap_or(default_v_b);

    // Compute characteristic surface dimensions for adaptive parameters.
    let diag_a = {
        let p00 = surf_a(u_range_a.0, v_range_a.0);
        let p11 = surf_a(u_range_a.1, v_range_a.1);
        (p00 - p11).length()
    };
    let diag_b = {
        let p00 = surf_b(u_range_b.0, v_range_b.0);
        let p11 = surf_b(u_range_b.1, v_range_b.1);
        (p00 - p11).length()
    };
    let char_size = diag_a.min(diag_b).max(0.1);

    // Sample surface A on a grid. For each grid point, project it
    // analytically onto surface B to find the closest point, then check
    // if the distance is below threshold (indicating near-intersection).
    #[allow(clippy::type_complexity)]
    let mut seeds: Vec<(Point3, (f64, f64), (f64, f64))> = Vec::new();
    // Coarse threshold scales with the surface size — the distance from
    // a grid point on A to its projection on B can be large even near
    // the intersection (e.g., sphere R=2 and cylinder R=1 → gap ≈ 1).
    let seed_threshold = diag_a.max(diag_b).max(1.0) * 0.5;
    let mut min_dist = f64::INFINITY;

    #[allow(clippy::cast_precision_loss)]
    for ia in 0..grid_res {
        for ja in 0..grid_res {
            let ua =
                u_range_a.0 + (u_range_a.1 - u_range_a.0) * (ia as f64 + 0.5) / (grid_res as f64);
            let va =
                v_range_a.0 + (v_range_a.1 - v_range_a.0) * (ja as f64 + 0.5) / (grid_res as f64);

            let pa = surf_a(ua, va);

            // Analytically project onto surface B.
            let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
            let pb = surf_b(ub, vb);
            let dist = (pa - pb).length();
            min_dist = min_dist.min(dist);

            if dist < seed_threshold {
                // Use the coarse seed directly. The marching algorithm
                // corrects positions at each step via projection, so seeds
                // don't need to be on the exact intersection — they just
                // need to be close enough for the marcher to converge.
                let mid = Point3::new(
                    (pa.x() + pb.x()) * 0.5,
                    (pa.y() + pb.y()) * 0.5,
                    (pa.z() + pb.z()) * 0.5,
                );
                seeds.push((mid, (ua, va), (ub, vb)));
            }
        }
    }

    // Cheap rejection: the grid samples surface A; the closest sample's
    // distance to B lower-bounds how near the two bounded patches come. A
    // transversal crossing puts a sample within ~one grid cell of it
    // (distance on the order of a cell), so if even the nearest sample is
    // several cells away the patches cannot cross — skip the expensive
    // marching and return empty. Result-preserving: non-crossing pairs
    // already march to nothing, just slowly (this is the gridfinity lip's
    // ~80 inner-wall × outer-wall pairs that dominate pavefiller time).
    let reject_dist = (char_size / grid_res as f64) * 3.0;
    if min_dist > reject_dist {
        return Ok(vec![]);
    }

    if seeds.is_empty() {
        return Ok(vec![]);
    }

    // Aggressively deduplicate seeds — we only need 1-2 per intersection
    // branch. Scale dedup radius to ~2% of characteristic surface size
    // (at least 10× the march step size) to avoid redundant marches.
    let march_step = (char_size * 0.02).clamp(0.005, 0.5);
    let dedup_radius = march_step * 10.0;
    let mut unique_seeds = Vec::new();
    for seed in &seeds {
        let dominated = unique_seeds
            .iter()
            .any(|s: &(Point3, (f64, f64), (f64, f64))| (s.0 - seed.0).length() < dedup_radius);
        if !dominated {
            unique_seeds.push(*seed);
        }
    }

    // March from each seed.
    let mut curves = Vec::new();
    let mut used_seeds = vec![false; unique_seeds.len()];

    for si in 0..unique_seeds.len() {
        if used_seeds[si] {
            continue;
        }
        used_seeds[si] = true;

        let march_result = march_analytic_intersection(
            &a,
            &b,
            surf_a.as_ref(),
            norm_a.as_ref(),
            surf_b.as_ref(),
            norm_b.as_ref(),
            unique_seeds[si].0,
            u_range_a,
            v_range_a,
            u_range_b,
            v_range_b,
            march_step,
            is_u_periodic(&a),
            is_u_periodic(&b),
        );

        if march_result.len() >= 2 {
            for (sj, other) in unique_seeds.iter().enumerate() {
                if !used_seeds[sj]
                    && march_result
                        .iter()
                        .any(|p| (*p - other.0).length() < dedup_radius)
                {
                    used_seeds[sj] = true;
                }
            }

            let ipts: Vec<IntersectionPoint> = march_result
                .iter()
                .map(|&pt| IntersectionPoint {
                    point: pt,
                    param1: (0.0, 0.0),
                    param2: (0.0, 0.0),
                })
                .collect();

            let degree = 3.min(march_result.len() - 1);
            if let Ok(curve) = interpolate(&march_result, degree) {
                curves.push(IntersectionCurve {
                    curve,
                    points: ipts,
                });
            }
        }
    }

    Ok(curves)
}

/// Try algebraic (closed-form or semi-algebraic) intersection for known
/// surface pairs before falling back to general marching.
///
/// Returns `Some(curves)` if a specialized method exists, `None` otherwise.
///
/// Currently handles:
/// - **Sphere-sphere**: intersection is a circle (plane through the two centers)
/// - **Coaxial cylinders**: same axis → circle(s) or empty
/// - **Sphere-cylinder**: reduce to quadratic in one parameter
#[allow(clippy::too_many_lines)]
fn try_algebraic_intersection(
    a: &AnalyticSurface<'_>,
    b: &AnalyticSurface<'_>,
    v_range_a: Option<(f64, f64)>,
    v_range_b: Option<(f64, f64)>,
) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
    match (a, b) {
        (AnalyticSurface::Cone(cone), AnalyticSurface::Cylinder(cyl)) => {
            algebraic_parallel_cone_cylinder(cone, cyl, v_range_a, v_range_b)
        }
        (AnalyticSurface::Cylinder(cyl), AnalyticSurface::Cone(cone)) => {
            algebraic_parallel_cone_cylinder(cone, cyl, v_range_b, v_range_a)
        }
        (AnalyticSurface::Sphere(s1), AnalyticSurface::Sphere(s2)) => {
            algebraic_sphere_sphere(s1, s2).map(Some)
        }
        (AnalyticSurface::Cylinder(c1), AnalyticSurface::Cylinder(c2)) => {
            let axis_dot = c1.axis().dot(c2.axis()).abs();
            if axis_dot > 1.0 - 1e-10 {
                // Axes are parallel — check if they're the same line.
                let delta = c2.origin() - c1.origin();
                let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
                let along = delta_vec.dot(c1.axis());
                let perp = (delta_vec - c1.axis() * along).length();
                if perp < 1e-8 {
                    // Coaxial: same axis, different radii → no intersection
                    // (unless equal radius → degenerate overlap, skip)
                    if (c1.radius() - c2.radius()).abs() < 1e-8 {
                        return Ok(None); // Overlapping — let marcher handle
                    }
                    return Ok(Some(vec![])); // Coaxial, different radii
                }
            }
            // Non-coaxial: algebraic quadratic in v.
            algebraic_cylinder_cylinder(c1, c2)
        }
        // Sphere-cylinder (both orderings).
        (AnalyticSurface::Sphere(s), AnalyticSurface::Cylinder(c)) => {
            algebraic_sphere_cylinder(s, c, true)
        }
        (AnalyticSurface::Cylinder(c), AnalyticSurface::Sphere(s)) => {
            algebraic_sphere_cylinder(s, c, false)
        }
        (AnalyticSurface::Cone(c1), AnalyticSurface::Cone(c2)) => algebraic_cone_cone(c1, c2),
        (AnalyticSurface::Torus(t), AnalyticSurface::Cylinder(c)) => {
            Ok(parallel_axis_torus_cylinder(t, c, true))
        }
        (AnalyticSurface::Cylinder(c), AnalyticSurface::Torus(t)) => {
            Ok(parallel_axis_torus_cylinder(t, c, false))
        }
        _ => Ok(None),
    }
}

/// A torus and a cylinder whose axes are parallel but distinct (a drill
/// through a ring parallel to its axis), traced along the cylinder's
/// rulings. A ruling stays at one distance `ρ` from the torus axis, so it
/// meets the tube where `(ρ − R)² + z² = r²`: a quadratic in its axial
/// parameter. `None` for any other pair (tilted or coaxial axes) and when
/// the sampling misses a window narrower than itself.
fn parallel_axis_torus_cylinder(
    torus: &ToroidalSurface,
    cyl: &CylindricalSurface,
    torus_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    let axis = torus.z_axis();
    let along = cyl.axis().dot(axis);
    if along.abs() < 1.0 - 1e-10 {
        return None;
    }
    let offset = cyl.origin() - torus.center();
    if (offset - axis * offset.dot(axis)).length() < Tolerance::new().linear {
        return None;
    }
    let (major, minor) = (torus.major_radius(), torus.minor_radius());
    let roots = |u: f64| {
        let q = cyl.evaluate(u, 0.0) - torus.center();
        let height = q.dot(axis);
        let rho = (q - axis * height).length();
        let reach = minor * minor - (rho - major) * (rho - major);
        ruling_quadratic(1.0, 2.0 * along.signum() * height, height * height - reach)
    };
    let samples = ruling_samples(cyl, &roots);
    let loops = if samples.iter().all(Option::is_some) {
        closed_ruling_loops(&samples)
    } else {
        partial_ruling_loops(cyl, &roots, &samples)
    };
    if loops.is_empty() {
        return None;
    }
    Some(fit_ruling_loops(&loops, |p| {
        in_order(torus.project_point(p), cyl.project_point(p), torus_first)
    }))
}

/// Exact coaxial cone-cone intersection: returns the shared circle.
///
/// Two cones that share an axis are concentric circles at every axial
/// station, so they meet only where their radii are equal. Each cone's
/// radius is linear in the axial coordinate `t` (measured along the shared
/// axis from cone 1's apex): `r1 = m1·t` and `r2 = m2·σ·(t − d2)`, where
/// `m_i = cot(half_angle_i)`, `σ = sign(axis2·axis1)`, and `d2` is cone 2's
/// apex position in that coordinate. Equating gives a single crossing `t*`
/// → one circle (the shared rim). The general marcher mishandles this case:
/// at the radii-crossing the surfaces are nearly tangent, so a grid-seeded
/// march fragments the clean circle into dozens of degenerate micro-curves.
///
/// Returns `Some(vec![circle])` for a genuine crossing, `Some(vec![])` when
/// the cones do not meet (parallel radius lines or a crossing on the wrong
/// nappe), and `None` for the identical-cone overlap or a degenerate
/// (near-flat) cone — both of which fall through to the general path.
/// Parallel-but-offset axes with equal half-angle tangents reduce to a
/// radical-plane conic (`offset_parallel_cone_cone`); other offset
/// configurations defer to the marcher with `None`.
///
/// # Errors
///
/// Returns [`MathError`] if the shared-rim `Circle3D` cannot be constructed
/// (e.g. a non-finite center or radius from a malformed cone).
pub fn exact_cone_cone(
    c1: &ConicalSurface,
    c2: &ConicalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let axis = c1.axis();
    let axis2 = c2.axis();

    // Coaxial check: parallel axes and the second apex lies on the first axis.
    if axis.dot(axis2).abs() < 1.0 - 1e-10 {
        return Ok(None); // Non-coaxial: quartic curve, let the marcher handle.
    }
    let apex1 = c1.apex();
    let apex2 = c2.apex();
    let delta = apex2 - apex1;
    let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
    let along = delta_v.dot(axis);
    if (delta_v - axis * along).length() > 1e-8 {
        return offset_parallel_cone_cone(c1, c2);
    }

    let (s1, s2) = (c1.half_angle().sin(), c2.half_angle().sin());
    if s1.abs() < 1e-12 || s2.abs() < 1e-12 {
        return Ok(None); // Degenerate (near-flat) cone.
    }
    let m1 = c1.half_angle().cos() / s1;
    let m2 = c2.half_angle().cos() / s2;
    let sigma = if axis.dot(axis2) >= 0.0 { 1.0 } else { -1.0 };
    let d2 = along; // apex2 position along `axis`, measured from apex1.

    let denom = m1 - m2 * sigma;
    if denom.abs() < 1e-12 {
        // Parallel radius lines: identical cones (coincident apex, same opening)
        // overlap — defer to the general/same-domain path; otherwise no meeting.
        if sigma > 0.0 && d2.abs() < 1e-9 {
            return Ok(None);
        }
        return Ok(Some(vec![]));
    }

    let t_star = (-m2 * sigma * d2) / denom;
    let radius = m1 * t_star;
    if radius < 1e-12 {
        return Ok(Some(vec![])); // Crossing on the wrong nappe / no real circle.
    }

    let center = Point3::new(
        apex1.x() + axis.x() * t_star,
        apex1.y() + axis.y() * t_star,
        apex1.z() + axis.z() * t_star,
    );
    let circle = Circle3D::new(center, axis, radius)?;
    Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
}

/// Parallel-axis (or anti-parallel), offset-apex cones with equal half-angle
/// tangents: subtracting the two quadric equations cancels both the radial
/// and the axial quadratic terms (their coefficients depend only on
/// `tan²(half_angle)`), so every intersection point lies on a plane — the
/// degenerate member of the quadric pencil — and plane ∩ cone is an exact
/// conic. The gridfinity spacer lip fuse hits this exactly: opposed 45°
/// corner cones offset 0.25mm, which the marcher shreds into ~64 closed
/// micro-loops per pair (#1570). Unequal angles keep a genuine quadratic
/// term, and an unbounded section (hyperbola/parabola) has no closed-form
/// win over the marcher — both defer with `None`.
fn offset_parallel_cone_cone(
    c1: &ConicalSurface,
    c2: &ConicalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    if c1.half_angle().sin().abs() < 1e-12 || c2.half_angle().sin().abs() < 1e-12 {
        return Ok(None); // Degenerate (near-flat) cone, as in the coaxial path.
    }
    let t1 = c1.half_angle().tan();
    let t2 = c2.half_angle().tan();
    if !t1.is_finite() || !t2.is_finite() {
        return Ok(None);
    }
    if (t1 - t2).abs() > 1e-9 * (1.0 + t1.abs().max(t2.abs())) {
        return Ok(None);
    }

    let w = c1.axis();
    let apex1 = c1.apex();
    let apex2 = c2.apex();
    let delta = apex2 - apex1;
    let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
    let s = delta_v.dot(w);
    let tm = 0.5 * (t1 + t2);
    let k = 1.0 + tm * tm;

    // In the apex1 frame each cone is |P|² − k(P·w)² = 0 (shifted by δ for
    // cone 2; the axis SIGN drops out since only (P·w)² appears). Their
    // difference: P·(2δ − 2ksw) = |δ|² − ks².
    let n = (delta_v - w * (k * s)) * 2.0;
    let n_len = n.length();
    if n_len < 1e-12 {
        return Ok(None);
    }
    let n_hat = n * (1.0 / n_len);
    let d = (dot_np(n, apex1) + delta_v.dot(delta_v) - k * s * s) / n_len;

    // `exact_plane_cone` already rejects sections on cone 1's phantom nappe;
    // cone 2's nappe must be checked here. A conic on the shared quadric
    // pencil cannot cross between nappes except exactly through apex 2, so
    // sampled quarter-points either all pass or all fail; a mixed verdict
    // means an apex-touching degeneracy — defer to the marcher.
    let axis2 = c2.axis();
    let scale = 1.0 + delta_v.length();
    let mut out = Vec::new();
    for curve in exact_plane_cone(c1, n_hat, d)? {
        let samples: Vec<Point3> = match &curve {
            ExactIntersectionCurve::Circle(c) => (0..4)
                .map(|i| crate::traits::ParametricCurve::evaluate(c, TAU * f64::from(i) / 4.0))
                .collect(),
            ExactIntersectionCurve::Ellipse(e) => (0..4)
                .map(|i| crate::traits::ParametricCurve::evaluate(e, TAU * f64::from(i) / 4.0))
                .collect(),
            ExactIntersectionCurve::Points(_) => return Ok(None),
        };
        let on_real_nappe = |p: &Point3| {
            let rel = *p - apex2;
            Vec3::new(rel.x(), rel.y(), rel.z()).dot(axis2) >= -1e-9 * scale
        };
        let hits = samples.iter().filter(|p| on_real_nappe(p)).count();
        match hits {
            0 => {}
            4 => out.push(curve),
            _ => return Ok(None),
        }
    }
    Ok(Some(out))
}

/// Exact coaxial cone-cylinder intersection: returns the shared circle.
///
/// A cone and a cylinder sharing an axis are concentric circles at every
/// axial station, so they meet only where the cone's radius equals the
/// cylinder's. The cone radius is linear in the axial coordinate `t` from its
/// apex (`r = m·t`, `m = cot(half_angle)`), the cylinder radius is the
/// constant `R`, so `m·t = R` gives a single crossing `t*` → one circle. This
/// is the gridfinity lip's top knife edge (inner tapered corner = cone, outer
/// corner = cylinder, concentric, radii matching at `Z_PEAK`); the general
/// marcher fragments that near-tangent contact into dozens of degenerate
/// micro-curves.
///
/// Returns `Some(vec![circle])` for a genuine crossing, `Some(vec![])` when
/// the crossing degenerates to the apex, and `None` (defer to the marcher)
/// when the surfaces are not coaxial or the cone is near-flat / near-axial.
///
/// # Errors
///
/// Returns [`MathError`] if the shared `Circle3D` cannot be constructed.
pub fn exact_cone_cylinder(
    cone: &ConicalSurface,
    cyl: &CylindricalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let axis = cone.axis();
    let cyl_axis = cyl.axis();

    // Coaxial check: parallel axes and the cone apex on the cylinder's axis.
    if axis.dot(cyl_axis).abs() < 1.0 - 1e-10 {
        return Ok(None);
    }
    let apex = cone.apex();
    let delta = apex - cyl.origin();
    let delta_v = Vec3::new(delta.x(), delta.y(), delta.z());
    let along = delta_v.dot(cyl_axis);
    if (delta_v - cyl_axis * along).length() > 1e-8 {
        return Ok(None);
    }

    let s = cone.half_angle().sin();
    if s.abs() < 1e-12 {
        return Ok(None); // near-flat cone.
    }
    let m = cone.half_angle().cos() / s; // dr/dt along the cone axis.
    if m.abs() < 1e-12 {
        return Ok(None); // near-axial cone: radius ~constant.
    }

    let t_star = cyl.radius() / m; // where the cone radius m·t equals R.
    if t_star.abs() < 1e-12 {
        return Ok(Some(vec![])); // crossing at the apex — no real circle.
    }
    let center = Point3::new(
        apex.x() + axis.x() * t_star,
        apex.y() + axis.y() * t_star,
        apex.z() + axis.z() * t_star,
    );
    let circle = Circle3D::new(center, axis, cyl.radius())?;
    Ok(Some(vec![ExactIntersectionCurve::Circle(circle)]))
}

/// Algebraic cone-cone intersection (NURBS form for the general bounded
/// path). Delegates to [`exact_cone_cone`] and samples each exact conic
/// (coaxial circle or offset-parallel radical-plane ellipse) into an
/// interpolated NURBS `IntersectionCurve`, mirroring the
/// sphere-cylinder algebraic path. phase FF prefers the exact circle form
/// directly (so the section edge links to the coincident boundary), but a
/// caller of `intersect_analytic_analytic_bounded` still gets one clean
/// curve instead of the marcher's fragments.
fn algebraic_cone_cone(
    c1: &ConicalSurface,
    c2: &ConicalSurface,
) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
    let Some(exacts) = exact_cone_cone(c1, c2)? else {
        return Ok(None);
    };
    let mut curves = Vec::new();
    for exact in exacts {
        let n_samples = 33;
        let mut positions = Vec::with_capacity(n_samples);
        let mut points = Vec::with_capacity(n_samples);
        #[allow(clippy::cast_precision_loss)]
        for i in 0..n_samples {
            let theta = TAU * i as f64 / (n_samples - 1) as f64;
            let pt = match &exact {
                ExactIntersectionCurve::Circle(circle) => {
                    crate::traits::ParametricCurve::evaluate(circle, theta)
                }
                ExactIntersectionCurve::Ellipse(ellipse) => {
                    crate::traits::ParametricCurve::evaluate(ellipse, theta)
                }
                ExactIntersectionCurve::Points(_) => break,
            };
            positions.push(pt);
            points.push(IntersectionPoint {
                point: pt,
                param1: (0.0, 0.0),
                param2: (0.0, 0.0),
            });
        }
        if positions.is_empty() {
            continue;
        }
        let degree = 3.min(positions.len() - 1);
        let curve = interpolate(&positions, degree)?;
        curves.push(IntersectionCurve { curve, points });
    }
    Ok(Some(curves))
}

/// Exact coaxial sphere-cylinder intersection: returns the shared circle(s).
///
/// A sphere of radius `R` centered at `C` and a cylinder of radius `r` whose
/// axis passes through `C` meet in concentric circles of radius `r` at the
/// axial stations where `sqrt(R² − z²) = r`, i.e. `z = ±sqrt(R² − r²)`
/// measured from `C` along the axis. A proper crossing yields two circles; a
/// tangent contact (`r = R`) yields one; a cylinder wider than the sphere, or
/// a non-coaxial configuration (quartic curve), yields none/defers.
///
/// Mirrors [`exact_cone_cylinder`] so phase FF can emit the section as an
/// exact `Circle3D` (which the closed-circle split + seam adoption recognise)
/// rather than the marcher's NURBS fragments.
///
/// Returns `Some(vec![..])` (0, 1, or 2 circles) for the coaxial case, and
/// `None` (defer to the general marcher) when the axes are not coaxial.
///
/// # Errors
///
/// Returns [`MathError`] if a shared `Circle3D` cannot be constructed.
pub fn exact_sphere_cylinder(
    sphere: &SphericalSurface,
    cyl: &CylindricalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let sc = sphere.center();
    let r_sphere = sphere.radius();
    let co = cyl.origin();
    let axis = cyl.axis();
    let r_cyl = cyl.radius();

    // Project sphere center onto the cylinder axis.
    let delta = sc - co;
    let delta_vec = Vec3::new(delta.x(), delta.y(), delta.z());
    let along = delta_vec.dot(axis);
    let perp_vec = delta_vec - axis * along;
    let d_perp = perp_vec.length();

    // Non-coaxial sphere-cylinder intersections produce quartic curves;
    // defer those to the general marcher.
    if d_perp > 1e-7 {
        return Ok(None);
    }

    // Coaxial: the sphere center lies on the cylinder axis. No real circle
    // when the cylinder is wider than the sphere or they are tangent-internal.
    if r_cyl > r_sphere + 1e-10 {
        return Ok(Some(vec![]));
    }
    let z_sq = r_sphere * r_sphere - r_cyl * r_cyl;
    if z_sq < 0.0 {
        return Ok(Some(vec![]));
    }
    let z = z_sq.sqrt();

    // The sphere center projected onto the axis is the midpoint of the two
    // section circles, each offset by ±z along the axis with radius `r_cyl`.
    let center_axis_pt = Point3::new(
        co.x() + axis.x() * along,
        co.y() + axis.y() * along,
        co.z() + axis.z() * along,
    );

    let mut circles = Vec::new();
    let offsets: &[f64] = if z < 1e-10 { &[0.0] } else { &[z, -z] };
    for &z_offset in offsets {
        let center = Point3::new(
            center_axis_pt.x() + axis.x() * z_offset,
            center_axis_pt.y() + axis.y() * z_offset,
            center_axis_pt.z() + axis.z() * z_offset,
        );
        let circle = Circle3D::new(center, axis, r_cyl)?;
        circles.push(ExactIntersectionCurve::Circle(circle));
    }
    Ok(Some(circles))
}

/// Algebraic sphere-cylinder intersection (NURBS form for the general bounded
/// path). A coaxial pair delegates to [`exact_sphere_cylinder`] and samples
/// each exact circle into an interpolated NURBS `IntersectionCurve`. phase FF
/// prefers the exact circle form directly (so the section edge links to the
/// coincident boundary and the closed-circle splitter can carve the spherical
/// band), but a caller of `intersect_analytic_analytic_bounded` still gets
/// clean curves instead of the marcher's fragments. Any other pair is traced
/// along the cylinder's rulings ([`off_axis_sphere_cylinder`]).
fn algebraic_sphere_cylinder(
    sphere: &SphericalSurface,
    cyl: &CylindricalSurface,
    sphere_first: bool,
) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
    let Some(exacts) = exact_sphere_cylinder(sphere, cyl)? else {
        return Ok(off_axis_sphere_cylinder(sphere, cyl, sphere_first));
    };

    let mut curves = Vec::new();
    for exact in exacts {
        let ExactIntersectionCurve::Circle(circle) = exact else {
            continue;
        };
        let n_samples = 33;
        let mut points = Vec::with_capacity(n_samples);
        let mut positions = Vec::with_capacity(n_samples);
        #[allow(clippy::cast_precision_loss)]
        for i in 0..n_samples {
            let theta = TAU * i as f64 / (n_samples - 1) as f64;
            let pt = crate::traits::ParametricCurve::evaluate(&circle, theta);
            positions.push(pt);
            let (param1, param2) = in_order(
                sphere.project_point(pt),
                cyl.project_point(pt),
                sphere_first,
            );
            points.push(IntersectionPoint {
                point: pt,
                param1,
                param2,
            });
        }
        let degree = 3.min(positions.len() - 1);
        let curve = interpolate(&positions, degree)?;
        curves.push(IntersectionCurve { curve, points });
    }

    Ok(Some(curves))
}

/// A sphere and a cylinder whose axis misses the sphere's centre (a drill
/// entering a ball off its axis), traced along the cylinder's rulings: the
/// ruling `c(u) + v·a` meets the sphere where `v² + 2(q·a)·v + |q|² − R² = 0`,
/// with `q = c(u) − C`. When every ruling meets the sphere (the cylinder
/// passes wholly through it) the roots trace an entry and an exit loop;
/// otherwise each window of meeting rulings carries one loop. `Some(empty)`
/// when the two cannot meet, `None` when the sampling misses a window
/// narrower than itself.
fn off_axis_sphere_cylinder(
    sphere: &SphericalSurface,
    cyl: &CylindricalSurface,
    sphere_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    let (centre, radius) = (sphere.center(), sphere.radius());
    let axis = cyl.axis();
    let offset = centre - cyl.origin();
    let axis_distance = (offset - axis * offset.dot(axis)).length();
    let lin_tol = Tolerance::new().linear;
    if axis_distance > radius + cyl.radius() + lin_tol
        || axis_distance + radius < cyl.radius() - lin_tol
    {
        return Some(Vec::new());
    }
    let roots = |u: f64| {
        let q = cyl.evaluate(u, 0.0) - centre;
        ruling_quadratic(1.0, 2.0 * q.dot(axis), q.dot(q) - radius * radius)
    };
    let samples = ruling_samples(cyl, &roots);
    let loops = if samples.iter().all(Option::is_some) {
        closed_ruling_loops(&samples)
    } else {
        partial_ruling_loops(cyl, &roots, &samples)
    };
    if loops.is_empty() {
        return None;
    }
    Some(fit_ruling_loops(&loops, |p| {
        in_order(sphere.project_point(p), cyl.project_point(p), sphere_first)
    }))
}

/// Parameters on the pair's first and second surfaces, from those on `a`
/// and `b` and whether `a` came first.
const fn in_order(a: (f64, f64), b: (f64, f64), a_first: bool) -> ((f64, f64), (f64, f64)) {
    if a_first { (a, b) } else { (b, a) }
}

/// Algebraic cylinder-cylinder intersection for non-coaxial cylinders.
///
/// For two cylinders with axes that are NOT parallel, the intersection
/// consists of up to two closed space curves. These are found by
/// parameterizing one cylinder's angular coordinate `u ∈ [0, 2π]` and
/// solving a quadratic in the axial parameter `v` to find where each
/// "ring" of cylinder A sits on cylinder B.
///
/// The quadratic is:
///   `v²·(1 - α²) + 2v·(q·a₁ - α·q·a₂) + (|q|² - (q·a₂)² - r₂²) = 0`
/// where `α = a₁·a₂`, `q(u)` is the radial point on cylinder 1 minus
/// cylinder 2's origin, `a₁`/`a₂` are the cylinder axes, and `r₂` is
/// cylinder 2's radius.
#[allow(clippy::too_many_lines, clippy::unnecessary_wraps)]
fn algebraic_cylinder_cylinder(
    c1: &CylindricalSurface,
    c2: &CylindricalSurface,
) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
    let alpha = c1.axis().dot(c2.axis());
    let a_coeff = 1.0 - alpha * alpha;

    // Should only be called for non-parallel axes.
    if a_coeff.abs() < 1e-12 {
        return Ok(None);
    }

    let r1 = c1.radius();
    let r2 = c2.radius();
    let o1 = c1.origin();
    let o2 = c2.origin();
    let a1 = c1.axis();
    let a2 = c2.axis();

    // Separation check: distance between axes vs sum of radii.
    // Closest approach of two skew lines:
    let delta = Vec3::new(o1.x() - o2.x(), o1.y() - o2.y(), o1.z() - o2.z());
    let cross = a1.cross(a2);
    let cross_len = cross.length();
    if cross_len > 1e-12 {
        let axis_dist = delta.dot(cross).abs() / cross_len;
        if axis_dist > r1 + r2 + Tolerance::new().linear {
            return Ok(Some(vec![])); // No intersection
        }
    }

    // Solve every ruling of one cylinder against the other. When EVERY
    // ruling of the swept cylinder meets the other (the thinner of two
    // crossing tubes), the two roots trace the curve's two closed loops.
    // Swept the other way only a window of rulings meets, and each root
    // traces an open arc of one loop.
    let roots = |sweep: &CylindricalSurface, other: &CylindricalSurface| {
        let (o, a, radius) = (other.origin(), other.axis(), other.radius());
        let alpha = sweep.axis().dot(a);
        let quad = 1.0 - alpha * alpha;
        let (axis, sweep) = (sweep.axis(), sweep.clone());
        move |u: f64| {
            let q = sweep.evaluate(u, 0.0) - o;
            let (q_a1, q_a2) = (q.dot(axis), q.dot(a));
            let b = 2.0 * (q_a1 - alpha * q_a2);
            let c = q.dot(q) - q_a2 * q_a2 - radius * radius;
            ruling_quadratic(quad, b, c)
        }
    };
    let (roots1, roots2) = (roots(c1, c2), roots(c2, c1));
    let samples1 = ruling_samples(c1, &roots1);
    let loops = if samples1.iter().all(Option::is_some) {
        closed_ruling_loops(&samples1)
    } else {
        let samples2 = ruling_samples(c2, &roots2);
        if samples2.iter().all(Option::is_some) {
            closed_ruling_loops(&samples2)
        } else if samples1.iter().any(Option::is_some) {
            partial_ruling_loops(c1, &roots1, &samples1)
        } else {
            partial_ruling_loops(c2, &roots2, &samples2)
        }
    };
    if loops.is_empty() {
        return Ok(None);
    }
    Ok(Some(fit_ruling_loops(&loops, |p| {
        (c1.project_point(p), c2.project_point(p))
    })))
}

/// Rulings sampled around a swept cylinder, half a step off u = 0 so the
/// branches of a self-touching curve (equal crossing cylinders) do not share
/// a sample.
const RULING_SAMPLES: usize = 128;

#[allow(clippy::cast_precision_loss)]
fn ruling_u(i: usize) -> f64 {
    TAU * (i as f64 + 0.5) / RULING_SAMPLES as f64
}

/// The discriminant and roots of `quad·v² + b·v + c = 0`.
fn ruling_quadratic(quad: f64, b: f64, c: f64) -> (f64, f64, f64) {
    let disc = b * b - 4.0 * quad * c;
    let root = disc.max(0.0).sqrt();
    (disc, (-b + root) / (2.0 * quad), (-b - root) / (2.0 * quad))
}

/// The two points where each sampled ruling of `sweep` meets the other
/// surface, from `roots(u)` (the discriminant and the two axial parameters),
/// or `None` for a ruling that misses it.
fn ruling_samples(
    sweep: &CylindricalSurface,
    roots: &impl Fn(f64) -> (f64, f64, f64),
) -> Vec<Option<(Point3, Point3)>> {
    let lin_tol = Tolerance::new().linear;
    (0..RULING_SAMPLES)
        .map(|i| {
            let u = ruling_u(i);
            let (disc, vp, vm) = roots(u);
            (disc >= -lin_tol).then(|| (sweep.evaluate(u, vp), sweep.evaluate(u, vm)))
        })
        .collect()
}

/// Every ruling meets the other surface: each root traces a closed loop.
fn closed_ruling_loops(samples: &[Option<(Point3, Point3)>]) -> Vec<Vec<Point3>> {
    let mut plus: Vec<Point3> = samples.iter().flatten().map(|s| s.0).collect();
    let mut minus: Vec<Point3> = samples.iter().flatten().map(|s| s.1).collect();
    plus.push(plus[0]);
    minus.push(minus[0]);
    vec![plus, minus]
}

/// Only windows of rulings meet the other surface: each cyclic window
/// carries one loop, out along one root and back along the other, the two
/// joined where the discriminant vanishes. Empty when no sample meets it (a
/// window narrower than the sampling).
fn partial_ruling_loops(
    sweep: &CylindricalSurface,
    roots: &impl Fn(f64) -> (f64, f64, f64),
    samples: &[Option<(Point3, Point3)>],
) -> Vec<Vec<Point3>> {
    let branch_point = |inside: usize, outside: usize| -> Point3 {
        let (mut lo, mut hi) = (ruling_u(inside), ruling_u(outside));
        if (hi - lo).abs() > std::f64::consts::PI {
            hi += if hi < lo { TAU } else { -TAU };
        }
        for _ in 0..60 {
            let mid = 0.5 * (lo + hi);
            if roots(mid).0 >= 0.0 {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        let (_, vp, vm) = roots(lo);
        sweep.evaluate(lo, 0.5 * (vp + vm))
    };
    let Some(first_gap) = samples.iter().position(Option::is_none) else {
        return Vec::new();
    };
    let mut loops = Vec::new();
    let mut k = 0;
    while k < RULING_SAMPLES {
        let i = (first_gap + k) % RULING_SAMPLES;
        if samples[i].is_none() {
            k += 1;
            continue;
        }
        let start = i;
        let mut run = Vec::new();
        while k < RULING_SAMPLES {
            let j = (first_gap + k) % RULING_SAMPLES;
            let Some(pair) = samples[j] else { break };
            run.push(pair);
            k += 1;
        }
        let end = (start + run.len() - 1) % RULING_SAMPLES;
        let head = branch_point(start, (start + RULING_SAMPLES - 1) % RULING_SAMPLES);
        let tail = branch_point(end, (end + 1) % RULING_SAMPLES);
        let mut pts = vec![head];
        pts.extend(run.iter().map(|p| p.0));
        pts.push(tail);
        pts.extend(run.iter().rev().map(|p| p.1));
        pts.push(head);
        loops.push(pts);
    }
    loops
}

/// Cubic interpolants through the swept loops, with `params(p)` giving each
/// point's parameters on the two surfaces.
fn fit_ruling_loops(
    loops: &[Vec<Point3>],
    params: impl Fn(Point3) -> ((f64, f64), (f64, f64)),
) -> Vec<IntersectionCurve> {
    let mut curves = Vec::new();
    for pts in loops {
        if pts.len() < 4 {
            continue;
        }
        let ipts: Vec<IntersectionPoint> = pts
            .iter()
            .map(|&p| {
                let (param1, param2) = params(p);
                IntersectionPoint {
                    point: p,
                    param1,
                    param2,
                }
            })
            .collect();
        let degree = 3.min(pts.len() - 1);
        if let Ok(curve) = interpolate(pts, degree) {
            curves.push(IntersectionCurve {
                curve,
                points: ipts,
            });
        }
    }
    curves
}

/// Algebraic cone-cylinder intersection for PARALLEL (or antiparallel) axes.
///
/// When the axes are parallel, every plane perpendicular to them cuts the cone
/// in a circle of radius `rho = v * cos(half_angle)` about a FIXED centre and
/// the cylinder in a circle of radius `R` about a second FIXED centre, so the
/// axis separation `d` is constant in `v`. Two coplanar circles meet at
/// `u = phi0 +/- acos((d^2 + rho^2 - R^2) / (2*d*rho))`, giving two branches
/// parameterised exactly by the cone's own `v`. The branches exist only where
/// `rho` lies in `[|d - R|, d + R]`, which bounds the curve naturally.
///
/// This replaces the general grid-seeded marcher for the configuration, which
/// mis-handles it badly: seeds are accepted anywhere within half the surface
/// diagonal of the partner, the march-result dedup only consumes seeds the
/// traced polyline passes near, and the survivors are dozens of overlapping
/// partial traces of the same curve. Those fragments carry no usable in-face
/// span, so a cone corner-round crossed by a boss cylinder never splits (a
/// counterbore/countersink meeting a pad — the gridfinity lightweight base).
///
/// Returns `None` (defer to the caller's other paths) when the axes are not
/// parallel, or when they are coaxial — a coaxial pair degenerates to shared
/// circles, which [`exact_cone_cylinder`] emits exactly and phase FF calls
/// directly. Note that `intersect_analytic_analytic_bounded` does NOT consult
/// `exact_cone_cylinder`, so a coaxial pair reaching this path through that
/// caller falls through to the marcher; only the FF path gets the exact circles.
// Result-wrapped to match the other `try_algebraic_intersection` arms' shape.
#[allow(clippy::unnecessary_wraps)]
fn algebraic_parallel_cone_cylinder(
    cone: &ConicalSurface,
    cyl: &CylindricalSurface,
    v_range_cone: Option<(f64, f64)>,
    v_range_cyl: Option<(f64, f64)>,
) -> Result<Option<Vec<IntersectionCurve>>, MathError> {
    let axis = cone.axis();
    if axis.dot(cyl.axis()).abs() < 1.0 - 1e-10 {
        return Ok(None); // Skew/oblique — general marcher.
    }

    let apex = cone.apex();
    let delta = cyl.origin() - apex;
    let along = delta.dot(axis);
    let perp = delta - axis * along;
    let d = perp.length();
    if d < 1e-9 {
        return Ok(None); // Coaxial — `exact_cone_cylinder` owns this.
    }

    let (e1, e2) = (cone.x_axis(), cone.y_axis());
    let phi0 = perp.dot(e2).atan2(perp.dot(e1));

    let (sin_t, cos_t) = cone.half_angle().sin_cos();
    if cos_t < 1e-12 || sin_t < 1e-12 {
        return Ok(None);
    }
    let r = cyl.radius();

    // Branch existence: |d - R| <= rho <= d + R, with rho = v * cos(half_angle).
    let mut v_min = (d - r).abs() / cos_t;
    let mut v_max = (d + r) / cos_t;
    if v_max <= v_min {
        return Ok(Some(vec![]));
    }

    // Narrow the sampled span to the faces' own extents so the fixed sample
    // budget resolves the in-face part of the curve rather than spreading over
    // a loop that mostly lies off both patches. A face's crossing can be a
    // fraction of a degree of the cone's sweep (the corner-round case above),
    // and an unnarrowed sampling puts fewer than one sample across it.
    let mut lo = v_min;
    let mut hi = v_max;
    // Clip EXACTLY to the hints, not to a padded window: an endpoint that lands
    // exactly on the face's own v-limit lies ON that boundary rim, so the
    // downstream pave machinery anchors it to the rim edge instead of leaving
    // the section dangling just past the face.
    if let Some((a, b)) = v_range_cone {
        let (a, b) = if a <= b { (a, b) } else { (b, a) };
        lo = lo.max(a);
        hi = hi.min(b);
    }
    if let Some((a, b)) = v_range_cyl {
        // The cylinder's v is a signed distance along its axis from its origin;
        // convert both ends to the cone's v via the shared axial direction.
        let flip = cyl.axis().dot(axis);
        let to_cone_v = |cv: f64| (along + cv * flip) / sin_t;
        let (a, b) = (to_cone_v(a), to_cone_v(b));
        let (a, b) = if a <= b { (a, b) } else { (b, a) };
        lo = lo.max(a);
        hi = hi.min(b);
    }
    v_min = lo.max(v_min);
    v_max = hi.min(v_max);
    if v_max - v_min <= 1e-12 {
        return Ok(Some(vec![]));
    }

    let n_samples = 128;
    let mut plus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
    let mut minus: Vec<Point3> = Vec::with_capacity(n_samples + 1);
    #[allow(clippy::cast_precision_loss)]
    for i in 0..=n_samples {
        let v = v_min + (v_max - v_min) * (i as f64) / (n_samples as f64);
        let rho = v * cos_t;
        if rho < 1e-12 {
            // The apex. `cos_alpha` has rho in its denominator, so it is only
            // meaningful in the limit: it tends to 0 (alpha -> pi/2) when the
            // cylinder passes exactly through the apex (d == R), and diverges
            // otherwise — where the clamp would manufacture a spurious alpha of
            // 0 or pi. So keep the apex only in the d == R case, where it is a
            // genuine point of the intersection and the shared endpoint at
            // which the two branches meet.
            if (d - r).abs() < 1e-12 {
                let apex = cone.evaluate(phi0, v);
                plus.push(apex);
                minus.push(apex);
            }
            continue;
        }
        let cos_alpha = ((d * d + rho * rho - r * r) / (2.0 * d * rho)).clamp(-1.0, 1.0);
        let alpha = cos_alpha.acos();
        plus.push(cone.evaluate(phi0 + alpha, v));
        minus.push(cone.evaluate(phi0 - alpha, v));
    }

    let mut curves = Vec::new();
    for pts in [&plus, &minus] {
        // Fewer than four samples in range means this branch does not cross the
        // bounded region at all (the other branch may still).
        if pts.len() < 4 {
            continue;
        }
        let ipts: Vec<IntersectionPoint> = pts
            .iter()
            .map(|&p| IntersectionPoint {
                point: p,
                param1: cone.project_point(p),
                param2: cyl.project_point(p),
            })
            .collect();
        let degree = 3.min(pts.len() - 1);
        match interpolate(pts, degree) {
            Ok(curve) => curves.push(IntersectionCurve {
                curve,
                points: ipts,
            }),
            // Emitting only the branch that happened to fit would starve the
            // section chain of exactly the piece this path exists to supply —
            // the same silent half-answer the marcher's fragments produced.
            // Defer the whole pair to the caller's other paths instead.
            Err(_) => return Ok(None),
        }
    }

    Ok(Some(curves))
}

/// Algebraic sphere-sphere intersection.
///
/// Two spheres intersect in a circle lying in the radical plane.
/// The radical plane is perpendicular to the line connecting the centers,
/// at a distance d1 from center1 where:
///   d1 = (D² + R1² - R2²) / (2D)
/// and D is the distance between centers.
fn algebraic_sphere_sphere(
    s1: &SphericalSurface,
    s2: &SphericalSurface,
) -> Result<Vec<IntersectionCurve>, MathError> {
    let c1 = s1.center();
    let c2 = s2.center();
    let r1 = s1.radius();
    let r2 = s2.radius();

    let delta = c2 - c1;
    let d_sq = delta.x() * delta.x() + delta.y() * delta.y() + delta.z() * delta.z();
    let d = d_sq.sqrt();

    if d < 1e-12 {
        // Concentric spheres: no intersection (unless same radius → degenerate).
        return Ok(vec![]);
    }

    // Check separation conditions.
    if d > r1 + r2 + 1e-10 {
        return Ok(vec![]); // Too far apart
    }
    if d + r2.min(r1) + 1e-10 < r1.max(r2) {
        return Ok(vec![]); // One inside the other
    }

    // Distance from c1 to the radical plane along the center line.
    let d1 = (d_sq + r1 * r1 - r2 * r2) / (2.0 * d);

    // Radius of the intersection circle.
    let r_circle_sq = r1 * r1 - d1 * d1;
    if r_circle_sq < 0.0 {
        // Tangent or no intersection (numerical noise).
        if r_circle_sq > -1e-10 {
            // Tangent: single point.
            let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
            let tangent_pt = Point3::new(
                c1.x() + axis.x() * d1,
                c1.y() + axis.y() * d1,
                c1.z() + axis.z() * d1,
            );
            let ipt = IntersectionPoint {
                point: tangent_pt,
                param1: (0.0, 0.0),
                param2: (0.0, 0.0),
            };
            // Single-point "curve" — not very useful but correct.
            return Ok(vec![IntersectionCurve {
                curve: interpolate(&[tangent_pt, tangent_pt], 1)?,
                points: vec![ipt],
            }]);
        }
        return Ok(vec![]);
    }

    let r_circle = r_circle_sq.sqrt();
    let axis = Vec3::new(delta.x() / d, delta.y() / d, delta.z() / d);
    let center = Point3::new(
        c1.x() + axis.x() * d1,
        c1.y() + axis.y() * d1,
        c1.z() + axis.z() * d1,
    );

    // Build a reference frame for the circle.
    let basis = Frame3::from_normal(center, axis)?;
    let u_dir = basis.x;
    let v_dir = basis.y;

    // Sample the circle for the IntersectionCurve representation.
    let n_samples = 33; // Odd for symmetry
    let mut points = Vec::with_capacity(n_samples);
    let mut positions = Vec::with_capacity(n_samples);
    #[allow(clippy::cast_precision_loss)]
    for i in 0..n_samples {
        let theta = TAU * i as f64 / (n_samples - 1) as f64;
        let (sin_t, cos_t) = theta.sin_cos();
        let pt = Point3::new(
            center.x() + (u_dir.x() * cos_t + v_dir.x() * sin_t) * r_circle,
            center.y() + (u_dir.y() * cos_t + v_dir.y() * sin_t) * r_circle,
            center.z() + (u_dir.z() * cos_t + v_dir.z() * sin_t) * r_circle,
        );
        positions.push(pt);
        points.push(IntersectionPoint {
            point: pt,
            param1: (0.0, 0.0),
            param2: (0.0, 0.0),
        });
    }

    let degree = 3.min(positions.len() - 1);
    let curve = interpolate(&positions, degree)?;

    Ok(vec![IntersectionCurve { curve, points }])
}

/// Newton correction: project a point back onto the intersection curve
/// of two analytic surfaces. Solves the 3×3 system:
///   δ · na = -da  (eliminate distance to surface A)
///   δ · nb = -db  (eliminate distance to surface B)
///   δ · t  = 0    (minimal correction, perpendicular to tangent)
#[allow(clippy::too_many_arguments)]
fn correct_to_intersection(
    a: &AnalyticSurface<'_>,
    b: &AnalyticSurface<'_>,
    surf_a: &dyn Fn(f64, f64) -> Point3,
    norm_a: &dyn Fn(f64, f64) -> Vec3,
    surf_b: &dyn Fn(f64, f64) -> Point3,
    norm_b: &dyn Fn(f64, f64) -> Vec3,
    point: Point3,
    u_range_a: (f64, f64),
    v_range_a: (f64, f64),
    u_range_b: (f64, f64),
    v_range_b: (f64, f64),
    max_iters: usize,
) -> Point3 {
    let mut p = point;
    for _ in 0..max_iters {
        let (ua, va) = project_analytic(a, p, u_range_a, v_range_a);
        let (ub, vb) = project_analytic(b, p, u_range_b, v_range_b);
        let pa = surf_a(ua, va);
        let pb = surf_b(ub, vb);
        let na = norm_a(ua, va);
        let nb = norm_b(ub, vb);
        let pv = Vec3::new(p.x(), p.y(), p.z());

        let da = (pv - Vec3::new(pa.x(), pa.y(), pa.z())).dot(na);
        let db = (pv - Vec3::new(pb.x(), pb.y(), pb.z())).dot(nb);

        if da.abs() < 1e-7 && db.abs() < 1e-7 {
            break;
        }

        let t = na.cross(nb);
        let t_len = t.length();
        if t_len < 1e-10 {
            // Surfaces are tangent — fall back to midpoint.
            return Point3::new(
                (pa.x() + pb.x()) * 0.5,
                (pa.y() + pb.y()) * 0.5,
                (pa.z() + pb.z()) * 0.5,
            );
        }
        let t_hat = t * (1.0 / t_len);

        // Solve [na; nb; t_hat] · δ = [-da, -db, 0] via Cramer's rule.
        let det = na.x() * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
            - na.y() * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
            + na.z() * (nb.x() * t_hat.y() - nb.y() * t_hat.x());
        if det.abs() < 1e-15 {
            return Point3::new(
                (pa.x() + pb.x()) * 0.5,
                (pa.y() + pb.y()) * 0.5,
                (pa.z() + pb.z()) * 0.5,
            );
        }
        let inv = 1.0 / det;
        // Cramer's rule: replace each column of A with rhs = (-da, -db, 0).
        let dx = inv
            * (-da * (nb.y() * t_hat.z() - nb.z() * t_hat.y())
                + db * (na.y() * t_hat.z() - na.z() * t_hat.y()));
        let dy = inv
            * (da * (nb.x() * t_hat.z() - nb.z() * t_hat.x())
                - db * (na.x() * t_hat.z() - na.z() * t_hat.x()));
        let dz = inv
            * (-da * (nb.x() * t_hat.y() - nb.y() * t_hat.x())
                + db * (na.x() * t_hat.y() - na.y() * t_hat.x()));
        let candidate = Point3::new(p.x() + dx, p.y() + dy, p.z() + dz);

        // Divergence guard: if the correction moves farther from both
        // surfaces, abandon Newton and return the best point so far.
        let (uc, vc) = project_analytic(a, candidate, u_range_a, v_range_a);
        let (ud, vd) = project_analytic(b, candidate, u_range_b, v_range_b);
        let pc_a = surf_a(uc, vc);
        let pc_b = surf_b(ud, vd);
        let cv = Vec3::new(candidate.x(), candidate.y(), candidate.z());
        let da_new = (cv - Vec3::new(pc_a.x(), pc_a.y(), pc_a.z()))
            .dot(norm_a(uc, vc))
            .abs();
        let db_new = (cv - Vec3::new(pc_b.x(), pc_b.y(), pc_b.z()))
            .dot(norm_b(ud, vd))
            .abs();
        if da_new > da.abs() && db_new > db.abs() {
            return p;
        }

        p = candidate;
    }
    p
}

/// March along the intersection of two surfaces from a seed point.
///
/// Uses the cross product of surface normals as the tangent direction
/// and projects back onto both surfaces using analytical projection
/// (for cylinders/spheres) or grid search (fallback).
#[allow(clippy::too_many_arguments)]
fn march_analytic_intersection(
    a: &AnalyticSurface<'_>,
    b: &AnalyticSurface<'_>,
    surf_a: &dyn Fn(f64, f64) -> Point3,
    norm_a: &dyn Fn(f64, f64) -> Vec3,
    surf_b: &dyn Fn(f64, f64) -> Point3,
    norm_b: &dyn Fn(f64, f64) -> Vec3,
    seed: Point3,
    u_range_a: (f64, f64),
    v_range_a: (f64, f64),
    u_range_b: (f64, f64),
    v_range_b: (f64, f64),
    initial_step: f64,
    u_periodic_a: bool,
    u_periodic_b: bool,
) -> Vec<Point3> {
    let max_steps = 500;
    let h_min = 1e-6;
    let h_max = initial_step * 4.0;
    // Fixed closure threshold: the adaptive step `h` varies with curvature
    // and can shrink below the actual miss distance at the seed re-approach.
    // Use `initial_step * 5` to robustly detect closure on the first pass.
    let closure_dist = initial_step * 5.0;
    // Angular thresholds for curvature-adaptive stepping.
    let max_angle = 10.0_f64.to_radians();
    let min_angle = 2.0_f64.to_radians();

    // March forward from seed, collecting points.
    let mut forward = Vec::new();
    // March backward from seed, collecting points (reversed at end).
    let mut backward = Vec::new();

    for (direction, points) in [(1.0_f64, &mut forward), (-1.0_f64, &mut backward)] {
        let mut current = seed;
        let mut h = initial_step;
        let mut prev_tangent: Option<Vec3> = None;

        for _ in 0..max_steps {
            let (ua, va) = project_analytic(a, current, u_range_a, v_range_a);
            let (ub, vb) = project_analytic(b, current, u_range_b, v_range_b);

            let na = norm_a(ua, va);
            let nb = norm_b(ub, vb);

            let tangent = na.cross(nb);
            let t_len = tangent.length();
            if t_len < 1e-10 {
                break;
            }
            let t_dir = tangent * (direction / t_len);

            // Curvature-adaptive step: check angular deviation from previous tangent.
            if let Some(prev_t) = prev_tangent {
                let cos_angle = prev_t.dot(t_dir).clamp(-1.0, 1.0);
                let angle = cos_angle.acos();
                if angle > max_angle && h > h_min {
                    h = (h * 0.5).max(h_min);
                } else if angle < min_angle {
                    h = (h * 2.0).min(h_max);
                }
            }
            prev_tangent = Some(t_dir);

            let next = Point3::new(
                h.mul_add(t_dir.x(), current.x()),
                h.mul_add(t_dir.y(), current.y()),
                h.mul_add(t_dir.z(), current.z()),
            );

            let (ua2, va2) = project_analytic(a, next, u_range_a, v_range_a);
            let (ub2, vb2) = project_analytic(b, next, u_range_b, v_range_b);

            let pa = surf_a(ua2, va2);
            let pb = surf_b(ub2, vb2);
            let mid = Point3::new(
                (pa.x() + pb.x()) * 0.5,
                (pa.y() + pb.y()) * 0.5,
                (pa.z() + pb.z()) * 0.5,
            );
            let out_a = (!u_periodic_a && (ua2 <= u_range_a.0 || ua2 >= u_range_a.1))
                || va2 <= v_range_a.0
                || va2 >= v_range_a.1;
            let out_b = (!u_periodic_b && (ub2 <= u_range_b.0 || ub2 >= u_range_b.1))
                || vb2 <= v_range_b.0
                || vb2 >= v_range_b.1;

            if out_a || out_b {
                break;
            }

            // Check for loop closure — if we've collected enough points and
            // the current point is close to the seed, the curve is closed.
            // Require ≥10 steps to avoid premature closure near the seed.
            let dist_to_seed = (mid - seed).length();
            if points.len() > 10 && dist_to_seed < closure_dist {
                points.push(seed);
                break;
            }

            points.push(mid);
            current = mid;
        }
    }

    // Assemble result: backward (reversed) + seed + forward
    backward.reverse();
    let mut result = backward;
    result.push(seed);
    result.append(&mut forward);

    // Refine all points onto the intersection curve via Newton correction.
    for pt in &mut result {
        *pt = correct_to_intersection(
            a, b, surf_a, norm_a, surf_b, norm_b, *pt, u_range_a, v_range_a, u_range_b, v_range_b,
            5,
        );
    }

    result
}

/// Project a 3D point onto an analytic surface using the surface's
/// analytical projection method. Falls back to grid search for surface
/// types without analytical projection.
fn project_analytic(
    surface: &AnalyticSurface<'_>,
    point: Point3,
    u_range: (f64, f64),
    v_range: (f64, f64),
) -> (f64, f64) {
    match surface {
        AnalyticSurface::Cylinder(cyl) => {
            let (u, v) = cyl.project_point(point);
            (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
        }
        AnalyticSurface::Sphere(sphere) => {
            let (u, v) = sphere.project_point(point);
            (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
        }
        AnalyticSurface::Cone(cone) => {
            let (u, v) = cone.project_point(point);
            (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
        }
        AnalyticSurface::Torus(torus) => {
            let (u, v) = torus.project_point(point);
            (u.clamp(u_range.0, u_range.1), v.clamp(v_range.0, v_range.1))
        }
    }
}

/// Returns `true` if the surface's u-parameter is periodic (wraps around 2π).
/// All current `AnalyticSurface` variants have periodic u — this is trivially
/// true today but exists as a guard for future non-periodic analytic types.
fn is_u_periodic(surface: &AnalyticSurface<'_>) -> bool {
    matches!(
        surface,
        AnalyticSurface::Cylinder(_)
            | AnalyticSurface::Cone(_)
            | AnalyticSurface::Sphere(_)
            | AnalyticSurface::Torus(_)
    )
}

/// Extract closures and parameter ranges for an analytic surface.
#[allow(clippy::type_complexity)]
fn surface_closures<'a>(
    surface: &'a AnalyticSurface<'a>,
) -> (
    Box<dyn Fn(f64, f64) -> Point3 + 'a>,
    Box<dyn Fn(f64, f64) -> Vec3 + 'a>,
    (f64, f64),
    (f64, f64),
) {
    match surface {
        AnalyticSurface::Cylinder(cyl) => (
            Box::new(|u, v| cyl.evaluate(u, v)),
            Box::new(|u, v| cyl.normal(u, v)),
            (0.0, TAU),
            (-1.0, 1.0),
        ),
        AnalyticSurface::Cone(cone) => (
            Box::new(|u, v| cone.evaluate(u, v)),
            Box::new(|u, v| cone.normal(u, v)),
            (0.0, TAU),
            (0.01, 2.0),
        ),
        AnalyticSurface::Sphere(sphere) => (
            Box::new(|u, v| sphere.evaluate(u, v)),
            Box::new(|u, v| sphere.normal(u, v)),
            (0.0, TAU),
            (-FRAC_PI_2, FRAC_PI_2),
        ),
        AnalyticSurface::Torus(torus) => (
            Box::new(|u, v| torus.evaluate(u, v)),
            Box::new(|u, v| torus.normal(u, v)),
            (0.0, TAU),
            (0.0, TAU),
        ),
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used)]
mod tests {
    use super::*;
    use crate::tolerance::Tolerance;

    #[test]
    fn plane_cylinder_perpendicular() {
        let cyl =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
                .unwrap();

        // Horizontal plane at z=3 -- produces a circle at height 3.
        let curves = intersect_plane_cylinder(&cyl, Vec3::new(0.0, 0.0, 1.0), 3.0).unwrap();
        assert!(!curves.is_empty(), "should find intersection curve");
        assert!(
            curves[0].points.len() > 10,
            "should have many sample points"
        );

        let tol = Tolerance::loose();
        for pt in &curves[0].points {
            assert!(
                tol.approx_eq(pt.point.z(), 3.0),
                "z should be ~3.0, got {}",
                pt.point.z()
            );
            let r = pt.point.x().hypot(pt.point.y());
            assert!(tol.approx_eq(r, 2.0), "radius should be ~2.0, got {r}");
        }
    }

    #[test]
    fn plane_sphere_equator() {
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();

        let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
        assert!(!curves.is_empty());

        let tol = Tolerance::loose();
        for pt in &curves[0].points {
            assert!(
                tol.approx_eq(pt.point.z(), 0.0),
                "z should be ~0, got {}",
                pt.point.z()
            );
            let r = pt.point.x().hypot(pt.point.y());
            assert!(tol.approx_eq(r, 3.0), "radius should be ~3.0, got {r}");
        }
    }

    #[test]
    fn plane_sphere_no_intersection() {
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0).unwrap();

        let curves = intersect_plane_sphere(&sphere, Vec3::new(0.0, 0.0, 1.0), 5.0).unwrap();
        assert!(curves.is_empty());
    }

    #[test]
    fn plane_cone_cross_section() {
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();

        let curves = intersect_plane_cone(&cone, Vec3::new(0.0, 0.0, 1.0), 1.0).unwrap();
        assert!(!curves.is_empty(), "should find intersection with cone");
    }

    /// The 1u gridfinity spacer lip fuse corner (#1570): the body's lip
    /// recess cone (45 deg, opening downward) meets the tool's lip cone
    /// (45 deg, opening upward) with axes offset 0.25mm in x and y. Equal
    /// half-angle tangents put the whole intersection on the radical plane,
    /// so the section is one exact ellipse; the marcher shredded this into
    /// ~64 closed micro-loops per pair.
    #[test]
    fn offset_parallel_equal_angle_cones_give_one_exact_ellipse() {
        let c1 = ConicalSurface::new(
            Point3::new(
                -16.999_999_999_999_975,
                -16.999_999_999_999_975,
                5.849_999_999_999_951,
            ),
            Vec3::new(0.0, 0.0, -1.0),
            0.785_398_163_397_433_5,
        )
        .unwrap();
        let c2 = ConicalSurface::new(
            Point3::new(
                -16.750_000_000_000_036,
                -16.750_000_000_000_018,
                0.749_999_999_999_881,
            ),
            Vec3::new(0.0, 0.0, 1.0),
            0.785_398_163_397_467_6,
        )
        .unwrap();

        let curves = exact_cone_cone(&c1, &c2)
            .unwrap()
            .expect("offset parallel equal-angle cones must take the radical-plane path");
        assert_eq!(curves.len(), 1, "expected exactly one section conic");
        assert!(
            matches!(curves[0], ExactIntersectionCurve::Ellipse(_)),
            "expected an ellipse section, got {:?}",
            curves[0]
        );
        let ExactIntersectionCurve::Ellipse(ellipse) = &curves[0] else {
            return;
        };

        // Every sample must lie on BOTH cones: distance to the axis equals
        // tan(half_angle) times the axial distance from the apex, on the
        // real nappe of each.
        for i in 0..16 {
            let p = crate::traits::ParametricCurve::evaluate(ellipse, TAU * f64::from(i) / 16.0);
            for (cone, label) in [(&c1, "c1"), (&c2, "c2")] {
                let rel = p - cone.apex();
                let rel_v = Vec3::new(rel.x(), rel.y(), rel.z());
                let axial = rel_v.dot(cone.axis());
                let radial = (rel_v - cone.axis() * axial).length();
                assert!(
                    axial > 0.0,
                    "{label}: sample on phantom nappe (axial {axial})"
                );
                let expect = cone.half_angle().tan() * axial;
                assert!(
                    (radial - expect).abs() < 1e-9,
                    "{label}: sample off surface by {}",
                    (radial - expect).abs()
                );
            }
        }
    }

    /// Opposed cones whose real nappes occupy disjoint half-spaces share a
    /// radical-plane conic only on the phantom nappe — the exact path must
    /// report a definitive empty intersection, not defer to the marcher.
    #[test]
    fn offset_parallel_cones_opening_apart_have_no_real_intersection() {
        let c1 = ConicalSurface::new(
            Point3::new(0.0, 0.0, 5.0),
            Vec3::new(0.0, 0.0, -1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let c2 = ConicalSurface::new(
            Point3::new(0.25, 0.25, 20.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let curves = exact_cone_cone(&c1, &c2)
            .unwrap()
            .expect("radical-plane path");
        assert!(curves.is_empty(), "disjoint nappes must yield no curves");
    }

    /// Unequal half-angles keep a quadratic term in the pencil — no plane
    /// reduction exists, so the exact path must defer to the marcher.
    #[test]
    fn offset_parallel_cones_with_unequal_angles_defer() {
        let c1 = ConicalSurface::new(
            Point3::new(0.0, 0.0, 5.0),
            Vec3::new(0.0, 0.0, -1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let c2 = ConicalSurface::new(Point3::new(0.25, 0.25, 0.5), Vec3::new(0.0, 0.0, 1.0), 0.6)
            .unwrap();
        assert!(exact_cone_cone(&c1, &c2).unwrap().is_none());
    }

    #[test]
    fn coaxial_cones_cross_at_single_circle() {
        // Two coaxial truncated cones (outer base r10->top r8, inner r9->r8
        // over height 10) cross where their radii match: z=10, r=8. The
        // intersection must be ONE clean circle, not the dozens of degenerate
        // micro-curves the general marcher produces at near-tangency.
        let outer = ConicalSurface::new(
            Point3::new(0.0, 0.0, 50.0),
            Vec3::new(0.0, 0.0, -1.0),
            5.0_f64.atan(),
        )
        .unwrap();
        let inner = ConicalSurface::new(
            Point3::new(0.0, 0.0, 90.0),
            Vec3::new(0.0, 0.0, -1.0),
            10.0_f64.atan(),
        )
        .unwrap();

        let curves = intersect_analytic_analytic_bounded(
            AnalyticSurface::Cone(&outer),
            AnalyticSurface::Cone(&inner),
            32,
            None,
            None,
        )
        .unwrap();

        assert_eq!(
            curves.len(),
            1,
            "coaxial cones crossing at one circle must yield exactly one curve, got {}",
            curves.len()
        );
        for p in &curves[0].points {
            let r = p.point.x().hypot(p.point.y());
            assert!(
                (p.point.z() - 10.0).abs() < 1e-6 && (r - 8.0).abs() < 1e-6,
                "intersection point off the expected z=10,r=8 circle: {:?}",
                p.point
            );
        }
    }

    #[test]
    fn plane_torus_cross_section() {
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 5.0, 1.0).unwrap();

        let curves = intersect_plane_torus(&torus, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
        assert!(
            !curves.is_empty(),
            "should find intersection curves with torus"
        );
    }

    /// Signed distance of a point to a z-axis torus centred at the origin:
    /// `sqrt((sqrt(x^2+y^2) - R)^2 + z^2) - r`.
    fn torus_implicit(p: Point3, major: f64, minor: f64) -> f64 {
        let rho = p.x().hypot(p.y());
        ((rho - major).hypot(p.z())) - minor
    }

    /// The gridfinity lightweight base's failing corner, reduced: a cavity
    /// corner-round cone (apex below the floor, 45 deg, axis +z) crossed by a
    /// parallel-axis boss cylinder. The general marcher returned ~49 overlapping
    /// partial traces of one curve here; the algebraic path must return exactly
    /// the two branches, each ON both surfaces and inside the cone's v-hint.
    #[test]
    fn parallel_cone_cylinder_gives_two_exact_branches() {
        use crate::traits::ParametricCurve;
        let cone = ConicalSurface::new(
            Point3::new(-5.45, -36.55, -4.85),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let cyl = CylindricalSurface::new(
            Point3::new(-8.0, -34.0, -5.0),
            Vec3::new(0.0, 0.0, 1.0),
            4.45,
        )
        .unwrap();
        // The cone face spans z in [-3.8, -3.0]; v = (z - apex_z) / sin(45 deg).
        let v_hint = (1.484_924_240_492_058, 2.616_295_090_390_43);
        let curves = intersect_analytic_analytic_bounded(
            AnalyticSurface::Cone(&cone),
            AnalyticSurface::Cylinder(&cyl),
            32,
            Some(v_hint),
            Some((0.0, 2.5)),
        )
        .unwrap();

        assert_eq!(curves.len(), 2, "expected exactly the two branches");
        for c in &curves {
            let (t0, t1) = c.curve.domain();
            for k in 0..=32 {
                let t = (t1 - t0).mul_add(f64::from(k) / 32.0, t0);
                let p = ParametricCurve::evaluate(&c.curve, t);
                // On the cylinder: radial distance from its axis is the radius.
                let radial = ((p.x() + 8.0).powi(2) + (p.y() + 34.0).powi(2)).sqrt();
                assert!((radial - 4.45).abs() < 1e-6, "off cylinder: {radial}");
                // On the cone: radial distance from its axis is z - apex_z.
                let cone_r = ((p.x() + 5.45).powi(2) + (p.y() + 36.55).powi(2)).sqrt();
                assert!((cone_r - (p.z() + 4.85)).abs() < 1e-6, "off cone at {p:?}");
                // Inside the cone face's own v-window (the hint is respected).
                assert!(p.z() >= -3.8 - 1e-9 && p.z() <= -3.0 + 1e-9, "z={}", p.z());
            }
        }
    }

    /// A coaxial pair has no radical line; the algebraic path must defer rather
    /// than divide by a zero axis separation.
    #[test]
    fn coaxial_cone_cylinder_defers_to_other_paths() {
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let cyl =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 2.0)
                .unwrap();
        assert!(
            algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
                .unwrap()
                .is_none()
        );
    }

    #[test]
    fn oblique_cone_cylinder_defers_to_other_paths() {
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let cyl =
            CylindricalSurface::new(Point3::new(3.0, 0.0, 1.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
                .unwrap();
        assert!(
            algebraic_parallel_cone_cylinder(&cone, &cyl, None, None)
                .unwrap()
                .is_none()
        );
    }

    #[test]
    fn plane_torus_lobe_closes_and_stays_on_surface() {
        use crate::traits::ParametricCurve;
        let (major, minor) = (10.0, 3.0);
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();

        // The census cutting planes (y=-4, x=6) each cut the +x and -x tube lobes
        // in a CLOSED oval. The greedy marcher stops one grid step short of
        // closing; the wrap-close must make every fitted lobe close exactly.
        for (n, d) in [
            (Vec3::new(0.0, -1.0, 0.0), 4.0),  // y = -4
            (Vec3::new(-1.0, 0.0, 0.0), -6.0), // x = 6
            (Vec3::new(0.0, 0.0, 1.0), 0.0),   // z = 0 -> two concentric circles
        ] {
            let curves = intersect_plane_torus(&torus, n, d).unwrap();
            assert!(!curves.is_empty(), "plane n={n:?} d={d} found no curves");
            for c in &curves {
                let p0 = ParametricCurve::evaluate(&c.curve, 0.0);
                let p1 = ParametricCurve::evaluate(&c.curve, 1.0);
                assert!(
                    (p0 - p1).length() < 1e-7,
                    "lobe not closed: gap={} (n={n:?} d={d})",
                    (p0 - p1).length()
                );
                // Every fitted sample stays on the torus (shape-preserving).
                for k in 0..=64 {
                    let t = f64::from(k) / 64.0;
                    let p = ParametricCurve::evaluate(&c.curve, t);
                    assert!(
                        torus_implicit(p, major, minor).abs() < 1e-2,
                        "off-surface point {p:?} implicit={}",
                        torus_implicit(p, major, minor)
                    );
                }
            }
        }
    }

    #[test]
    fn plane_torus_inner_tangent_figure_eight_stays_open() {
        use crate::traits::ParametricCurve;
        let (major, minor) = (10.0, 3.0);
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();

        // A plane tangent to the inner equator (x = major - minor = 7) cuts a
        // self-touching figure-eight. The marcher traces it as a single chain
        // whose end lands on the opposite lobe — FAR from its start (gap is many
        // point-spacings). The wrap-close must NOT force-close this into a wrong
        // loop; it must stay OPEN so a self-touching curve is never sealed.
        let curves =
            intersect_plane_torus(&torus, Vec3::new(-1.0, 0.0, 0.0), -(major - minor)).unwrap();
        assert!(!curves.is_empty(), "inner-tangent plane found no curves");
        let max_gap = curves
            .iter()
            .map(|c| {
                let p0 = ParametricCurve::evaluate(&c.curve, 0.0);
                let p1 = ParametricCurve::evaluate(&c.curve, 1.0);
                (p0 - p1).length()
            })
            .fold(0.0_f64, f64::max);
        assert!(
            max_gap > 1e-2,
            "figure-eight chain was wrongly force-closed (max end-gap={max_gap})"
        );
    }

    #[test]
    fn line_torus_box_edge_crossing_is_exact() {
        // The census box edge x=6, y=-4 (z varying) crosses the torus (R=10,r=3)
        // at z = ±sqrt(r² − (rho−R)²), rho = hypot(6,4) ≈ 7.2111 → z ≈ ±1.1055.
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
        let ts = intersect_line_torus(
            &torus,
            Point3::new(6.0, -4.0, -5.0),
            Vec3::new(0.0, 0.0, 1.0),
        );
        // Vertical line through (6,-4) meets the tube twice.
        assert_eq!(ts.len(), 2, "expected 2 crossings, got {ts:?}");
        let zs: Vec<f64> = ts.iter().map(|t| -5.0 + t).collect();
        let rho = 6.0_f64.hypot(4.0);
        let z_exp = (9.0 - (rho - 10.0).powi(2)).sqrt();
        assert!(
            (zs[0] - (-z_exp)).abs() < 1e-9,
            "z0={} exp={}",
            zs[0],
            -z_exp
        );
        assert!((zs[1] - z_exp).abs() < 1e-9, "z1={} exp={}", zs[1], z_exp);
        // Each crossing lies on the torus.
        for &t in &ts {
            let p = Point3::new(6.0, -4.0, -5.0 + t);
            let rho = p.x().hypot(p.y());
            let impl_v = (rho - 10.0).hypot(p.z()) - 3.0;
            assert!(impl_v.abs() < 1e-9, "off-torus impl={impl_v}");
        }
    }

    #[test]
    fn line_torus_miss_and_tangent() {
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 10.0, 3.0).unwrap();
        // A vertical line at rho beyond the outer rim (x=20) misses entirely.
        let miss = intersect_line_torus(
            &torus,
            Point3::new(20.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
        );
        assert!(miss.is_empty(), "expected no crossings, got {miss:?}");
        // The z-axis (rho=0) passes through the hole — no intersection.
        let axis =
            intersect_line_torus(&torus, Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
        assert!(axis.is_empty(), "z-axis should miss the tube, got {axis:?}");
    }

    #[test]
    fn dispatch_via_analytic_surface() {
        let cyl =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
                .unwrap();
        let curves = intersect_plane_analytic(
            AnalyticSurface::Cylinder(&cyl),
            Vec3::new(0.0, 0.0, 1.0),
            0.0,
        )
        .unwrap();
        assert!(!curves.is_empty());
    }

    #[test]
    fn perpendicular_cylinders_intersect() {
        let cyl_z =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
                .unwrap();
        let cyl_x =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
                .unwrap();

        let curves = intersect_analytic_analytic(
            AnalyticSurface::Cylinder(&cyl_z),
            AnalyticSurface::Cylinder(&cyl_x),
            16,
        )
        .unwrap();

        assert!(
            !curves.is_empty(),
            "perpendicular cylinders should intersect"
        );

        for c in &curves {
            assert!(
                c.points.len() >= 2,
                "intersection curve should have >= 2 points, got {}",
                c.points.len()
            );
        }
    }

    /// Neither cylinder's rulings all meet the other: the curve is one loop
    /// joined at its two branch points.
    #[test]
    fn partially_overlapping_cylinders_meet_in_one_closed_loop() {
        let cyl_z =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
                .unwrap();
        let cyl_x =
            CylindricalSurface::new(Point3::new(0.0, 1.2, 0.0), Vec3::new(1.0, 0.0, 0.0), 1.0)
                .unwrap();
        let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
            .unwrap()
            .unwrap();
        assert_eq!(curves.len(), 1);
        let curve = &curves[0].curve;
        let (t0, t1) = curve.domain();
        assert!((curve.evaluate(t0) - curve.evaluate(t1)).length() < 1e-9);
        let off = |p: Point3| {
            let on_z = (p.x().hypot(p.y()) - 1.0).abs();
            let on_x = ((p.y() - 1.2).hypot(p.z()) - 1.0).abs();
            on_z.max(on_x)
        };
        let worst = (0..=400)
            .map(|k| off(curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0)))
            .fold(0.0, f64::max);
        assert!(worst < 2e-4, "curve leaves the cylinders by {worst}");
    }

    /// Near tangency the thick cylinder's window of rulings (0.02 either side
    /// of a quarter turn) falls between its samples; the thin one's sweep
    /// finds the loop.
    #[test]
    fn near_tangent_cylinders_find_their_loop_on_the_thinner_sweep() {
        let cyl_z =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
                .unwrap();
        let cyl_x =
            CylindricalSurface::new(Point3::new(0.0, 1.1998, 0.0), Vec3::new(1.0, 0.0, 0.0), 0.2)
                .unwrap();
        let curves = algebraic_cylinder_cylinder(&cyl_z, &cyl_x)
            .unwrap()
            .expect("the thin cylinder's sweep finds the loop");
        assert_eq!(curves.len(), 1);
    }

    #[test]
    fn sphere_cylinder_intersect() {
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
        let cyl =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 1.0)
                .unwrap();

        let curves = intersect_analytic_analytic(
            AnalyticSurface::Sphere(&sphere),
            AnalyticSurface::Cylinder(&cyl),
            16,
        )
        .unwrap();

        // A sphere of radius 2 and a cylinder of radius 1, both centered
        // at the origin, should intersect (the cylinder passes through
        // the sphere).
        assert!(!curves.is_empty(), "sphere and cylinder should intersect");
    }

    #[test]
    fn exact_sphere_cylinder_coaxial_two_circles() {
        // Sphere r=6 at origin, coaxial cylinder r=3 along z: two latitude
        // circles at z = ±sqrt(36-9) = ±sqrt(27), each of radius 3.
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
        let cyl =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
                .unwrap();
        let circles = exact_sphere_cylinder(&sphere, &cyl)
            .unwrap()
            .expect("coaxial case returns Some");
        assert_eq!(circles.len(), 2, "through-bore meets the sphere twice");
        let mut zs: Vec<f64> = circles
            .iter()
            .filter_map(|c| match c {
                ExactIntersectionCurve::Circle(circle) => {
                    assert!(
                        (circle.radius() - 3.0).abs() < 1e-9,
                        "rim radius == cyl radius"
                    );
                    Some(circle.center().z())
                }
                _ => None,
            })
            .collect();
        assert_eq!(zs.len(), 2, "both sections must be exact circles");
        zs.sort_by(f64::total_cmp);
        let z = 27.0_f64.sqrt();
        assert!((zs[0] + z).abs() < 1e-9 && (zs[1] - z).abs() < 1e-9);
    }

    #[test]
    fn exact_sphere_cylinder_non_coaxial_defers() {
        // Cylinder axis offset from the sphere center → quartic curve, deferred.
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 6.0).unwrap();
        let cyl =
            CylindricalSurface::new(Point3::new(2.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 3.0)
                .unwrap();
        assert!(
            exact_sphere_cylinder(&sphere, &cyl).unwrap().is_none(),
            "non-coaxial sphere/cylinder defers to the marcher"
        );
    }

    /// Loops of an off-axis sphere-cylinder pair: `(count, worst distance
    /// from either surface)`.
    fn off_axis_loops(cylinder_origin: Point3, cylinder_radius: f64) -> (usize, f64) {
        let sphere = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 2.0).unwrap();
        let cyl =
            CylindricalSurface::new(cylinder_origin, Vec3::new(0.0, 0.0, 1.0), cylinder_radius)
                .unwrap();
        let curves = algebraic_sphere_cylinder(&sphere, &cyl, true)
            .unwrap()
            .unwrap();
        let mut worst: f64 = 0.0;
        for c in &curves {
            for ip in &c.points {
                let on_sphere = sphere.evaluate(ip.param1.0, ip.param1.1);
                let on_cylinder = cyl.evaluate(ip.param2.0, ip.param2.1);
                worst = worst
                    .max((on_sphere - ip.point).length())
                    .max((on_cylinder - ip.point).length());
            }
            let (t0, t1) = c.curve.domain();
            assert!((c.curve.evaluate(t0) - c.curve.evaluate(t1)).length() < 1e-9);
            for k in 0..=400 {
                let p = c.curve.evaluate(t0 + (t1 - t0) * f64::from(k) / 400.0);
                let on_sphere = ((p - Point3::new(0.0, 0.0, 0.0)).length() - 2.0).abs();
                let on_cylinder = ((p.x() - cylinder_origin.x())
                    .hypot(p.y() - cylinder_origin.y())
                    - cylinder_radius)
                    .abs();
                worst = worst.max(on_sphere).max(on_cylinder);
            }
        }
        (curves.len(), worst)
    }

    /// A drill off the ball's axis passes through it: an entry and an exit
    /// loop.
    #[test]
    fn off_axis_drill_through_a_sphere_meets_it_in_two_loops() {
        let (count, worst) = off_axis_loops(Point3::new(0.5, 0.0, 0.0), 0.2);
        assert_eq!(count, 2);
        assert!(worst < 1e-5, "loops leave the surfaces by {worst}");
    }

    /// A cylinder over the ball's side: one loop joined at its branch points.
    #[test]
    fn cylinder_over_a_spheres_side_meets_it_in_one_loop() {
        let (count, worst) = off_axis_loops(Point3::new(1.8, 0.0, 0.0), 0.5);
        assert_eq!(count, 1);
        assert!(worst < 5e-4, "loop leaves the surfaces by {worst}");
    }

    #[test]
    fn disjoint_cylinders_no_intersection() {
        let cyl_a =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
                .unwrap();
        let cyl_b =
            CylindricalSurface::new(Point3::new(5.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5)
                .unwrap();

        let curves = intersect_analytic_analytic(
            AnalyticSurface::Cylinder(&cyl_a),
            AnalyticSurface::Cylinder(&cyl_b),
            16,
        )
        .unwrap();

        assert!(curves.is_empty(), "disjoint cylinders should not intersect");
    }

    // ── Oblique plane × cone conic (ellipse / parabola / hyperbola) ──────

    /// Collect 3D points from a returned exact curve, sampling analytic forms.
    fn collect_points(curve: &ExactIntersectionCurve) -> Vec<Point3> {
        use crate::traits::ParametricCurve;
        match curve {
            ExactIntersectionCurve::Circle(c) => (0..=64)
                .map(|i| ParametricCurve::evaluate(c, TAU * f64::from(i) / 64.0))
                .collect(),
            ExactIntersectionCurve::Ellipse(e) => (0..=64)
                .map(|i| ParametricCurve::evaluate(e, TAU * f64::from(i) / 64.0))
                .collect(),
            ExactIntersectionCurve::Points(pts) => pts.clone(),
        }
    }

    /// Assert every returned point lies on the plane and the cone surface, on
    /// the real (`v >= 0`) nappe, and within a sane axial bound.
    fn assert_on_plane_and_cone(
        curves: &[ExactIntersectionCurve],
        cone: &ConicalSurface,
        n: Vec3,
        d: f64,
        z_bound: (f64, f64),
    ) {
        assert!(!curves.is_empty(), "expected at least one section curve");
        let mut total = 0;
        for curve in curves {
            for p in collect_points(curve) {
                total += 1;
                let plane_err = (n.x() * p.x() + n.y() * p.y() + n.z() * p.z() - d).abs();
                assert!(
                    plane_err < 1e-9,
                    "point off plane by {plane_err:.2e}: {p:?}"
                );
                let (u, v) = cone.project_point(p);
                let q = cone.evaluate(u, v);
                let cone_err =
                    ((p.x() - q.x()).powi(2) + (p.y() - q.y()).powi(2) + (p.z() - q.z()).powi(2))
                        .sqrt();
                assert!(cone_err < 1e-7, "point off cone by {cone_err:.2e}: {p:?}");
                assert!(v >= -1e-9, "point on phantom nappe (v={v:.4}): {p:?}");
                assert!(
                    p.z() >= z_bound.0 - 1e-6 && p.z() <= z_bound.1 + 1e-6,
                    "point z={:.4} outside sane bound {z_bound:?}: {p:?}",
                    p.z()
                );
            }
        }
        assert!(total >= 8, "too few section points ({total})");
    }

    #[test]
    fn oblique_plane_cone_ellipse_is_exact_and_on_both() {
        // 45°-half-angle cone (axis +z). A plane tilted only ~16.7° off horizontal
        // has plane-axis angle ≈ 73° > 45° (the cone's half-opening from axis) →
        // ellipse. Must come back as an exact Ellipse, fully on both surfaces.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
        // Plane through (0,0,5): d = n·(0,0,5).
        let d = n.z() * 5.0;
        let curves = exact_plane_cone(&cone, n, d).unwrap();
        assert!(
            curves
                .iter()
                .any(|c| matches!(c, ExactIntersectionCurve::Ellipse(_))),
            "oblique steep plane × cone must yield an exact Ellipse"
        );
        // The ellipse straddles z=5; with the 0.3 tilt the z-extent stays modest.
        assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 12.0));
    }

    #[test]
    fn oblique_plane_cone_wrong_nappe_is_empty() {
        // Same ellipse-regime plane as above, but offset to the FAR side of the
        // apex (z=-5). The +z cone's real (v≥0) nappe is not met — only the
        // phantom v<0 nappe — so the result must be EMPTY, not a phantom ellipse.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let n = Vec3::new(0.3, 0.0, 1.0).normalize().unwrap();
        let d = n.z() * -5.0;
        let curves = exact_plane_cone(&cone, n, d).unwrap();
        assert!(
            curves.is_empty(),
            "plane on the phantom-nappe side must yield no real curve, got {}",
            curves.len()
        );
    }

    #[test]
    fn oblique_plane_cone_parabola_on_both_single_branch() {
        // Plane normal at exactly 45° to the axis (= the cone half-opening) → the
        // plane is parallel to a generator → parabola. One unbounded branch.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let n = Vec3::new(1.0, 0.0, 1.0).normalize().unwrap();
        let d = n.x() * 3.0 + n.z() * 3.0; // through (3,0,3)
        let curves = exact_plane_cone(&cone, n, d).unwrap();
        assert_eq!(
            curves.len(),
            1,
            "a parabola is a single branch, got {}",
            curves.len()
        );
        // Bounded by r_max = 32·|e|; |e| here is O(few), so allow a wide z window.
        assert_on_plane_and_cone(&curves, &cone, n, d, (0.0, 400.0));
    }

    #[test]
    fn oblique_plane_cone_hyperbola_real_nappe_only() {
        // Faithful scooplabel lip-foot geometry: a 45° cone with axis −z and
        // apex at (−59,−59,15.85) (a bin corner), cut by the upper ramp tread
        // plane n=(0,0.99518,0.09802), d=−58.36056. The plane is nearly parallel
        // to the axis (cos≈0.098) → plane-axis angle ≈ 5.6° < 45° → hyperbola.
        // The downward real nappe is hit by exactly one branch; the phantom
        // upward nappe (and the asymptote runaway) must NOT appear, and the arc
        // must stay near the apex (the plane is ~1.2 mm from it).
        let cone = ConicalSurface::new(
            Point3::new(-59.0, -59.0, 15.85),
            Vec3::new(0.0, 0.0, -1.0),
            std::f64::consts::FRAC_PI_4,
        )
        .unwrap();
        let n = Vec3::new(0.0, 0.995_18, 0.098_02).normalize().unwrap();
        let d = -58.360_56;
        let cos_theta = n.dot(cone.axis()).abs();
        assert!(cos_theta < 0.2, "expected a shallow (hyperbola) plane");
        let curves = exact_plane_cone(&cone, n, d).unwrap();
        // Real downward nappe only: never above the apex (z=15.85). The vertex is
        // ~1.2 mm from the apex, so the bounded arc stays within a few mm of it.
        assert_on_plane_and_cone(&curves, &cone, n, d, (5.0, 15.85));
        // Every returned curve is sampled Points (no false Circle/Ellipse).
        for c in &curves {
            assert!(
                matches!(c, ExactIntersectionCurve::Points(_)),
                "hyperbola must be sampled Points, not a closed conic"
            );
        }
    }
}