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::aabb::Aabb3;
use crate::curves::{Circle3D, Ellipse3D};
use crate::frame::Frame3;
use crate::nurbs::curve::NurbsCurve;
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> {
    exact_plane_analytic_reaching(surface, plane_normal, plane_d, 0.0)
}

/// [`exact_plane_analytic`] with a cone's sampled hyperbola or parabola
/// carried at least `reach` from the apex, so it spans whatever faces the
/// caller will trim it to.
///
/// # Errors
///
/// Returns an error if the intersection computation fails.
pub fn exact_plane_analytic_reaching(
    surface: AnalyticSurface<'_>,
    plane_normal: Vec3,
    plane_d: f64,
    reach: 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, reach),
        AnalyticSurface::Torus(torus) => {
            if let Some(circles) = exact_plane_torus(torus, plane_normal, plane_d)? {
                return Ok(circles);
            }
            if let Some(loops) = plane_torus_winding_loops(torus, plane_normal, plane_d, 128) {
                return Ok(loops
                    .into_iter()
                    .map(ExactIntersectionCurve::Points)
                    .collect());
            }
            // Other torus sections 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())
        }
    }
}

/// The plane-torus sections that are circles:
///
/// - a plane across the axis at height `h` from the centre, `|h| < r`: the
///   two circles of radius `R ± sqrt(r² − h²)` about the axis;
/// - a plane through the axis: the two tube cross-sections of radius `r`,
///   `R` either side of the axis.
///
/// `Some` of no curves for a plane across the axis that misses the tube;
/// `None` for any other plane, a plane tangent to the tube, or a torus whose
/// tube reaches its axis.
fn exact_plane_torus(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let len = normal.length();
    let n = normal.normalize()?;
    let d = d / len;
    let axis = torus.z_axis();
    let center = torus.center();
    let (big, small) = (torus.major_radius(), torus.minor_radius());
    let height = d - dot_np(n, center);
    let along = n.dot(axis);
    if along.abs() > 1.0 - 1e-10 {
        if height.abs() >= small - 1e-10 * small {
            return Ok(if height.abs() > small + 1e-10 * small {
                Some(Vec::new())
            } else {
                None
            });
        }
        let reach = small.mul_add(small, -(height * height)).sqrt();
        if big - reach <= 1e-10 * big {
            return Ok(None);
        }
        let middle = center + n * height;
        return Ok(Some(vec![
            ExactIntersectionCurve::Circle(Circle3D::new(middle, n, big + reach)?),
            ExactIntersectionCurve::Circle(Circle3D::new(middle, n, big - reach)?),
        ]));
    }
    if along.abs() < 1e-10 && height.abs() < 1e-10 * (big + small) {
        let out = axis.cross(n).normalize()?;
        return Ok(Some(vec![
            ExactIntersectionCurve::Circle(Circle3D::new(center + out * big, n, small)?),
            ExactIntersectionCurve::Circle(Circle3D::new(center - out * big, n, small)?),
        ]));
    }
    Ok(None)
}

/// 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,
    reach: 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, reach)?;
        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, reach)?;
    Ok(chains
        .into_iter()
        .map(ExactIntersectionCurve::Points)
        .collect())
}

/// The exact arc from `from` to `to` of a plane's parabola or hyperbola
/// section of a cone, both ends on one branch of it, as a rational quadratic
/// NURBS.
///
/// A hyperbola `x = a cosh φ, y = b sinh φ` (in the plane frame of
/// `exact_plane_cone`) is cut into pieces of at most one unit of `φ`, each
/// a conic Bézier: its middle point is where the end tangents meet and its
/// middle weight is the cosh of half the piece's span. A parabola is one
/// polynomial quadratic. `None` for an elliptic or circular section, a plane
/// through the apex, or ends off one branch of the section.
///
/// # Errors
///
/// Returns an error if the plane normal is zero or the curve cannot be built.
#[allow(clippy::many_single_char_names)]
pub fn plane_cone_conic_arc(
    cone: &ConicalSurface,
    normal: Vec3,
    d: f64,
    from: Point3,
    to: Point3,
) -> Result<Option<NurbsCurve>, MathError> {
    let len = normal.length();
    if len < 1e-15 {
        return Err(MathError::ZeroVector);
    }
    let (normal, d) = (normal * (1.0 / len), d / len);
    let axis = cone.axis();
    let c = normal.dot(axis);
    let p2 = (1.0 - c * c).max(0.0);
    let p = p2.sqrt();
    let k = cone.half_angle().sin().powi(2);
    let a_coeff = p2 - k;
    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 || a_coeff < -1e-9 {
        return Ok(None);
    }
    let e1 = m * (1.0 / m_len);
    let e2 = normal.cross(e1);
    let apex = cone.apex();
    let e = d - dot_np(normal, apex);
    let origin = apex + normal * e;
    let plane_st = |q: Point3| {
        let w = q - origin;
        (w.dot(e1), w.dot(e2))
    };
    let ((s0, t0), (s1, t1)) = (plane_st(from), plane_st(to));
    let scale = s0.abs().max(t0.abs()).max(s1.abs()).max(t1.abs()).max(1.0);
    if e.abs() < 1e-9 * scale || (from - to).length() <= 1e-9 * scale {
        return Ok(None);
    }
    let point = |s: f64, t: f64| origin + e1 * s + e2 * t;
    let on_curve = |q: Point3, r: Point3| (q - r).length() <= 1e-6 * scale;
    let (control, weights) = if a_coeff.abs() <= 1e-9 {
        // (p² − k) s² vanishes: 2ecp·s + e²(c² − k) = k·t², s = α t² + β.
        let lin = 2.0 * e * c * p;
        if lin.abs() < 1e-12 * scale {
            return Ok(None);
        }
        let (alpha, beta) = (k / lin, -e * e * (c * c - k) / lin);
        if !on_curve(point(alpha * t0 * t0 + beta, t0), from)
            || !on_curve(point(alpha * t1 * t1 + beta, t1), to)
        {
            return Ok(None);
        }
        let mid = point(alpha * t0 * t1 + beta, 0.5 * (t0 + t1));
        (vec![from, mid, to], vec![1.0; 3])
    } else {
        // A (s − s_c)² − k t² = R with R = e² k (1 − k) / A.
        let s_c = -e * c * p / a_coeff;
        let r = e * e * k * (1.0 - k) / a_coeff;
        if r <= 0.0 {
            return Ok(None);
        }
        let (a, b) = ((r / a_coeff).sqrt(), (r / k).sqrt());
        let (x0, x1) = (s0 - s_c, s1 - s_c);
        if x0 * x1 <= 0.0 {
            return Ok(None);
        }
        let side = x0.signum();
        let hyperbola = |phi: f64| point(s_c + side * a * phi.cosh(), b * phi.sinh());
        let (phi0, phi1) = ((t0 / b).asinh(), (t1 / b).asinh());
        if !on_curve(hyperbola(phi0), from) || !on_curve(hyperbola(phi1), to) {
            return Ok(None);
        }
        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
        let pieces = ((phi1 - phi0).abs().ceil() as usize).max(1);
        let mut control = vec![from];
        let mut weights = vec![1.0];
        for i in 0..pieces {
            #[allow(clippy::cast_precision_loss)]
            let (fa, fb) = (i as f64 / pieces as f64, (i + 1) as f64 / pieces as f64);
            let (pa, pb) = (phi0 + (phi1 - phi0) * fa, phi0 + (phi1 - phi0) * fb);
            let (mid, half) = (0.5 * (pa + pb), 0.5 * (pb - pa));
            let w = half.cosh();
            control.push(point(s_c + side * a * mid.cosh() / w, b * mid.sinh() / w));
            weights.push(w);
            control.push(if i + 1 == pieces { to } else { hyperbola(pb) });
            weights.push(1.0);
        }
        (control, weights)
    };
    let pieces = (control.len() - 1) / 2;
    let mut knots = vec![0.0; 3];
    for i in 1..pieces {
        #[allow(clippy::cast_precision_loss)]
        knots.extend([i as f64; 2]);
    }
    #[allow(clippy::cast_precision_loss)]
    knots.extend([pieces as f64; 3]);
    let curve = NurbsCurve::new(2, knots, control, weights)?;
    // The closed forms drop terms that vanish only on the exact conic (a
    // barely elliptic section read as a parabola), so the arc must meet the
    // cone between its ends too: a point at radius ρ and height h off the
    // apex lies |ρ sin α − |h| cos α| from it.
    let (sin_a, cos_a) = cone.half_angle().sin_cos();
    let off_cone = |q: Point3| {
        let w = q - apex;
        let h = w.dot(axis);
        (w - axis * h)
            .length()
            .mul_add(sin_a, -(h.abs() * cos_a))
            .abs()
    };
    for i in 0..pieces {
        for f in [0.25, 0.5, 0.75] {
            #[allow(clippy::cast_precision_loss)]
            if off_cone(curve.evaluate(i as f64 + f)) > 1e-9 * scale {
                return Ok(None);
            }
        }
    }
    Ok(Some(curve))
}

/// 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, 0.0),
        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,
    reach: 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).
    // A caller that knows its faces asks for their reach: an open hyperbola
    // (a vertex close to the axis) crosses a rim far past eight vertex radii.
    let v_max = (8.0 * v_min).max(v_min + 4.0).max(reach);

    // 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 loops 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> {
    Ok(plane_torus_loops(torus, normal, d, 128)
        .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_loops`), so it is sampled by a v-scan
/// and each 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_loops`),
    // 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 mut curves = Vec::new();
    for ipts in plane_torus_loops(torus, normal, d, 128) {
        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)
}

/// The fewest and most steps along each branch of a loop a plane cuts from
/// a torus between its two turns.
const PLANE_TORUS_LOOP_SAMPLES: (f64, f64) = (24.0, 512.0);

/// The loops a plane cuts from a torus, each as points on the torus in
/// order, closed by repeating its first point.
///
/// 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 at each tube angle `v` where `|rhs(v)| <= 1`, with
/// `rhs(v) = (d − n·center − r·c·sin v) / (s·(R + r·cos v))`, the section's
/// two branches are `u = phi ± acos(rhs(v))`. Each run of `v` where the
/// section exists is one loop: out along one branch to where the run ends
/// and the branches meet, and back along the other. A section present at
/// every `v` is two loops winding the tube, one per branch. Every point
/// comes from the closed form, so the loops need no chaining.
///
/// A loop whose branches touch inside its run (a plane tangent to an
/// equator) is a self-touching section, not a simple loop, and is left open.
///
/// 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`,
/// sampled by scanning `u`.
#[allow(clippy::cast_precision_loss, clippy::too_many_lines)]
fn plane_torus_loops(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
    n_v: usize,
) -> Vec<Vec<IntersectionPoint>> {
    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 point = |u: f64, v: f64| IntersectionPoint {
        point: torus.evaluate(u, v),
        param1: (u, v.rem_euclid(TAU)),
        param2: (0.0, 0.0),
    };
    let closed = |mut run: Vec<IntersectionPoint>| {
        run.push(run[0]);
        run
    };

    // Plane perpendicular to the axis: the section is up to two full circles.
    if s < 1e-12 {
        if c.abs() < 1e-12 {
            return Vec::new();
        }
        let sin_v = d_local / (small_r * c);
        if sin_v.abs() > 1.0 + 1e-9 {
            return Vec::new();
        }
        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);
        }
        return vs
            .into_iter()
            .map(|v| {
                closed(
                    (0..n_v)
                        .map(|i| point(TAU * (i as f64) / (n_v as f64), v))
                        .collect(),
                )
            })
            .collect();
    }

    // The scan is offset by half a step so it never lands on a node where
    // the branches touch (the inner-tangent figure-eight at v = π).
    let step = TAU / (n_v as f64);
    let v_off = step * 0.5;
    // R + r·cos v > 0 on a ring torus.
    let rhs_at = |v: f64| (d_local - small_r * c * v.sin()) / (s * small_r.mul_add(v.cos(), big_r));
    let branch = |v: f64, sign: f64| point(sign.mul_add(rhs_at(v).clamp(-1.0, 1.0).acos(), phi), v);
    let inside = |v: f64| rhs_at(v).abs() <= 1.0;
    let scan: Vec<f64> = (0..n_v).map(|i| (i as f64).mul_add(step, v_off)).collect();
    // Whether |rhs| reaches 1 between two scan samples inside the section:
    // the branches touch there.
    let touches = |lo: f64, hi: f64| {
        let golden = 0.5 * (5.0_f64.sqrt() - 1.0);
        let (mut lo, mut hi) = (lo, hi);
        for _ in 0..80 {
            let (m1, m2) = (hi - golden * (hi - lo), lo + golden * (hi - lo));
            if rhs_at(m1).abs() > rhs_at(m2).abs() {
                hi = m2;
            } else {
                lo = m1;
            }
        }
        1.0 - rhs_at(f64::midpoint(lo, hi)).abs() < 1e-12
    };
    // Where the section's run ends between an inside and an outside sample.
    let turn = |v_in: f64, v_out: f64| {
        let (mut lo, mut hi) = (v_in, v_out);
        for _ in 0..60 {
            let mid = f64::midpoint(lo, hi);
            if inside(mid) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        lo
    };
    let in_scan: Vec<bool> = scan.iter().map(|&v| inside(v)).collect();
    if in_scan.iter().all(|&x| x) {
        let touching = scan.iter().any(|&v| touches(v, v + step));
        return [1.0, -1.0]
            .into_iter()
            .map(|sign| {
                let run: Vec<IntersectionPoint> = scan.iter().map(|&v| branch(v, sign)).collect();
                if touching { run } else { closed(run) }
            })
            .collect();
    }
    let Some(first) = (0..n_v).find(|&i| in_scan[i] && !in_scan[(i + n_v - 1) % n_v]) else {
        return Vec::new();
    };
    let mut loops = Vec::new();
    let mut k = 0;
    while k < n_v {
        let i = (first + k) % n_v;
        if !in_scan[i] {
            k += 1;
            continue;
        }
        // The run from scan sample `i`, its `v` unwrapped past a turn.
        let len = (0..n_v - k).take_while(|&j| in_scan[(i + j) % n_v]).count();
        let v_a = scan[i];
        let v_b = ((len - 1) as f64).mul_add(step, v_a);
        let run_v = |j: usize| (j as f64).mul_add(step, v_a);
        let (t_lo, t_hi) = (turn(v_a, v_a - step), turn(v_b, v_b + step));
        let touching = (0..len - 1).any(|j| touches(run_v(j), run_v(j + 1)));
        // `u` moves like the square root of the distance in `v` to a turn, so
        // the loop is sampled at `v = t_lo + (t_hi - t_lo) (1 - cos θ) / 2`
        // for even steps of `θ`: `v` then moves like `θ²` at each turn and the
        // points fall at near-even steps along the curve, turns included. A
        // loop that sweeps far round the axis on a short run of `v` (a plane
        // through the centre, tilted a little) takes a step per half scan
        // step of `u` its branches sweep.
        let (u_lo, u_hi) = (0..len)
            .map(run_v)
            .chain([t_lo, t_hi])
            .map(|v| rhs_at(v).clamp(-1.0, 1.0).acos())
            .fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), u| {
                (lo.min(u), hi.max(u))
            });
        let m = (len as f64)
            .max((n_v as f64) * (u_hi - u_lo) / std::f64::consts::PI)
            .max(PLANE_TORUS_LOOP_SAMPLES.0)
            .min(PLANE_TORUS_LOOP_SAMPLES.1)
            .ceil();
        let at = |k: f64| {
            let f = 0.5 * (1.0 - (std::f64::consts::PI * k / m).cos());
            (t_hi - t_lo).mul_add(f, t_lo)
        };
        let steps = m as usize;
        let mut pts: Vec<IntersectionPoint> =
            (0..=steps).map(|k| branch(at(k as f64), 1.0)).collect();
        pts.extend((1..steps).rev().map(|k| branch(at(k as f64), -1.0)));
        loops.push(if touching { pts } else { closed(pts) });
        k += len;
    }
    loops
}

/// The two sections of a plane that crosses every tube cross-section of a
/// torus twice (one parallel to the axis within `R − r` of it, or tilted a
/// little from that): both branches `u = phi ± acos(rhs(v))` of
/// [`plane_torus_loops`] are then defined for every `v`, so each closes
/// into a loop that winds once around the tube. Sampled at `n_v` steps from
/// `v = 0`, the outer equator, so every such loop on a torus starts on one
/// latitude, as the tube cross-sections of a plane through the axis do.
/// `None` when a branch lapses somewhere or the two come close to meeting.
#[allow(clippy::cast_precision_loss)]
fn plane_torus_winding_loops(
    torus: &ToroidalSurface,
    normal: Vec3,
    d: f64,
    n_v: usize,
) -> Option<Vec<Vec<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);
    if s < 1e-12 * normal.length() || small_r >= big_r {
        return None;
    }
    let phi = b.atan2(a);
    let d_local = d - dot_np(normal, torus.center());
    let rhs = |v: f64| (d_local - small_r * c * v.sin()) / (s * small_r.mul_add(v.cos(), big_r));
    let dense = 8 * n_v;
    if (0..dense).any(|i| rhs(TAU * i as f64 / dense as f64).abs() > 1.0 - 1e-3) {
        return None;
    }
    let mut loops = [Vec::with_capacity(n_v + 1), Vec::with_capacity(n_v + 1)];
    for i in 0..n_v {
        let v = TAU * i as f64 / n_v as f64;
        let delta = rhs(v).acos();
        loops[0].push(torus.evaluate(phi + delta, v));
        loops[1].push(torus.evaluate(phi - delta, v));
    }
    Some(
        loops
            .into_iter()
            .map(|mut run| {
                run.push(run[0]);
                run
            })
            .collect(),
    )
}

/// 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> {
    intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, None)
}

/// [`intersect_analytic_analytic_bounded`] with its general marcher confined
/// to `region`.
///
/// The marcher starts only from seeds that converge into the region (with a
/// margin of about two grid cells) and stops where a curve leaves it. Pairs
/// solved in closed form are returned whole, as by the unconfined call. Two
/// bounded faces meet only inside both their boxes, so their overlap is a
/// sound region for a face pair.
///
/// # Errors
///
/// Returns `MathError` if algebraic intersection fails or marching diverges.
pub fn intersect_analytic_analytic_in_region(
    a: AnalyticSurface<'_>,
    b: AnalyticSurface<'_>,
    grid_res: usize,
    v_range_hint_a: Option<(f64, f64)>,
    v_range_hint_b: Option<(f64, f64)>,
    region: Aabb3,
) -> Result<Vec<IntersectionCurve>, MathError> {
    intersect_analytic_analytic_impl(a, b, grid_res, v_range_hint_a, v_range_hint_b, Some(region))
}

fn intersect_analytic_analytic_impl(
    a: AnalyticSurface<'_>,
    b: AnalyticSurface<'_>,
    grid_res: usize,
    v_range_hint_a: Option<(f64, f64)>,
    v_range_hint_b: Option<(f64, f64)>,
    region: Option<Aabb3>,
) -> 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);
        }
    }

    // A seed is only near the intersection; settle it onto the curve before
    // asking whether it lies in the region, so a coarse grid point beside a
    // short in-region span still starts its march, and march from there, so
    // the march does not begin outside the region and stop at once.
    let region = region.map(|r| r.expanded(2.0 * char_size / grid_res as f64));
    if let Some(r) = region {
        for seed in &mut unique_seeds {
            let mut p = seed.0;
            for _ in 0..8 {
                let (ua, va) = project_analytic(&a, p, u_range_a, v_range_a);
                let pa = surf_a(ua, va);
                let (ub, vb) = project_analytic(&b, pa, u_range_b, v_range_b);
                let pb = surf_b(ub, vb);
                p = Point3::new(
                    (pa.x() + pb.x()) * 0.5,
                    (pa.y() + pb.y()) * 0.5,
                    (pa.z() + pb.z()) * 0.5,
                );
                if (pa - pb).length() < 1e-9 {
                    break;
                }
            }
            seed.0 = p;
        }
        unique_seeds.retain(|seed| r.contains_point(seed.0));
    }

    // 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),
            region,
        );

        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
/// - **Cone-cylinder**: parallel axes in the cone's own `v`, other axes
///   along the cylinder's rulings
/// - **Cone-sphere**: off the cone's axis, along the cone's generators
/// - **Torus-cylinder**: along the cylinder's rulings
#[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)) => Ok(
            algebraic_parallel_cone_cylinder(cone, cyl, v_range_a, v_range_b)?
                .or_else(|| ruling_cone_cylinder(cone, cyl, true)),
        ),
        (AnalyticSurface::Cylinder(cyl), AnalyticSurface::Cone(cone)) => Ok(
            algebraic_parallel_cone_cylinder(cone, cyl, v_range_b, v_range_a)?
                .or_else(|| ruling_cone_cylinder(cone, cyl, false)),
        ),
        (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::Cone(cone), AnalyticSurface::Sphere(sphere)) => {
            Ok(ruling_cone_sphere(cone, sphere, true))
        }
        (AnalyticSurface::Sphere(sphere), AnalyticSurface::Cone(cone)) => {
            Ok(ruling_cone_sphere(cone, sphere, false))
        }
        (AnalyticSurface::Torus(t), AnalyticSurface::Cylinder(c)) => {
            Ok(parallel_axis_torus_cylinder(t, c, true)
                .or_else(|| ruling_torus_cylinder(t, c, true)))
        }
        (AnalyticSurface::Cylinder(c), AnalyticSurface::Torus(t)) => {
            Ok(parallel_axis_torus_cylinder(t, c, false)
                .or_else(|| ruling_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)
    }))
}

/// Where two circles in a half-plane through an axis cross, as `(rho, z)`
/// pairs (distance from the axis, height along it): each sweeps a circle
/// about the axis. `None` (defer to the marcher) when the circles coincide
/// or touch, or a crossing lands on or past the axis; `Some` of none when
/// they miss.
fn meridian_crossings(
    first: (f64, f64, f64),
    second: (f64, f64, f64),
    scale: f64,
) -> Option<Vec<(f64, f64)>> {
    let ((x1, z1, r1), (x2, z2, r2)) = (first, second);
    let (dx, dz) = (x2 - x1, z2 - z1);
    let dist = dx.hypot(dz);
    let slack = 1e-9 * scale;
    if dist < slack || (dist - (r1 + r2)).abs() < slack || (dist - (r1 - r2).abs()).abs() < slack {
        return None;
    }
    if dist > r1 + r2 || dist < (r1 - r2).abs() {
        return Some(Vec::new());
    }
    let along = r2.mul_add(-r2, r1.mul_add(r1, dist * dist)) / (2.0 * dist);
    let across = r1.mul_add(r1, -(along * along)).max(0.0).sqrt();
    let (ux, uz) = (dx / dist, dz / dist);
    let mut crossings = Vec::with_capacity(2);
    for side in [1.0, -1.0] {
        let rho = x1 + along * ux - side * across * uz;
        if rho <= slack {
            return None;
        }
        crossings.push((rho, z1 + along * uz + side * across * ux));
    }
    Some(crossings)
}

/// Circles about an axis through `base`, at the given `(rho, z)` crossings.
fn circles_about_axis(
    base: Point3,
    axis: Vec3,
    crossings: &[(f64, f64)],
) -> Result<Vec<ExactIntersectionCurve>, MathError> {
    crossings
        .iter()
        .map(|&(rho, z)| {
            Circle3D::new(base + axis * z, axis, rho).map(ExactIntersectionCurve::Circle)
        })
        .collect()
}

/// Exact intersection of two tori sharing an axis: their tube cross-sections
/// in a half-plane through the axis cross in up to two points, and each sweeps
/// a circle about the axis.
///
/// `None` (defer to the marcher) unless the axes lie on one line, or when the
/// cross-sections coincide or touch.
///
/// # Errors
///
/// Returns an error if a section circle cannot be built.
pub fn exact_torus_torus(
    first: &ToroidalSurface,
    second: &ToroidalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let axis = first.z_axis();
    let scale = first.major_radius() + second.major_radius();
    let offset = second.center() - first.center();
    // A spindle torus's tube also crosses the far side of the axis.
    if first.minor_radius() >= first.major_radius()
        || second.minor_radius() >= second.major_radius()
        || axis.cross(second.z_axis()).length() > 1e-9
        || offset.cross(axis).length() > 1e-9 * scale
    {
        return Ok(None);
    }
    let Some(crossings) = meridian_crossings(
        (first.major_radius(), 0.0, first.minor_radius()),
        (
            second.major_radius(),
            offset.dot(axis),
            second.minor_radius(),
        ),
        scale,
    ) else {
        return Ok(None);
    };
    circles_about_axis(first.center(), axis, &crossings).map(Some)
}

/// Exact intersection of a torus with a cylinder sharing its axis.
///
/// The wall line and the tube's cross-section in a half-plane through the
/// axis cross in up to two points, each sweeping a circle about the axis.
///
/// `None` (defer to the marcher) unless the axes lie on one line, or when
/// the wall touches the tube.
///
/// # Errors
///
/// Returns an error if a section circle cannot be built.
pub fn exact_cylinder_torus(
    cylinder: &CylindricalSurface,
    torus: &ToroidalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let axis = torus.z_axis();
    let scale = torus.major_radius() + cylinder.radius();
    let offset = cylinder.origin() - torus.center();
    // A spindle torus's tube also crosses the far side of the axis.
    if torus.minor_radius() >= torus.major_radius()
        || axis.cross(cylinder.axis()).length() > 1e-9
        || offset.cross(axis).length() > 1e-9 * scale
    {
        return Ok(None);
    }
    let gap = cylinder.radius() - torus.major_radius();
    let small = torus.minor_radius();
    if (gap.abs() - small).abs() < 1e-9 * scale {
        return Ok(None);
    }
    if gap.abs() > small {
        return Ok(Some(Vec::new()));
    }
    let height = small.mul_add(small, -(gap * gap)).sqrt();
    circles_about_axis(
        torus.center(),
        axis,
        &[(cylinder.radius(), height), (cylinder.radius(), -height)],
    )
    .map(Some)
}

/// Exact intersection of a torus with a sphere centred on its axis.
///
/// The sphere's great circle and the tube's cross-section in a half-plane
/// through the axis cross in up to two points, and each sweeps a circle
/// about the axis.
///
/// `None` (defer to the marcher) unless the sphere's centre lies on the axis,
/// or when the two circles touch.
///
/// # Errors
///
/// Returns an error if a section circle cannot be built.
pub fn exact_sphere_torus(
    sphere: &SphericalSurface,
    torus: &ToroidalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let axis = torus.z_axis();
    let scale = torus.major_radius() + sphere.radius();
    let offset = sphere.center() - torus.center();
    // A spindle torus's tube also crosses the far side of the axis.
    if torus.minor_radius() >= torus.major_radius() || offset.cross(axis).length() > 1e-9 * scale {
        return Ok(None);
    }
    let Some(crossings) = meridian_crossings(
        (0.0, offset.dot(axis), sphere.radius()),
        (torus.major_radius(), 0.0, torus.minor_radius()),
        scale,
    ) else {
        return Ok(None);
    };
    circles_about_axis(torus.center(), axis, &crossings).map(Some)
}

/// 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, 0.0)? {
        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))
}

/// Exact coaxial cone-sphere intersection: the shared circles.
///
/// When the cone's axis runs through the sphere's centre `C`, every generator
/// `apex + v g` (`g` a unit direction, `v` the distance from the apex) has
/// the same `h = g·(apex − C)`, so all of them meet the sphere at the same
/// roots of `v² + 2hv + |apex − C|² − R² = 0`. Each root ahead of the apex is
/// a circle of radius `v cos a` at `v sin a` along the axis: two for a ball
/// the cone passes through, one for a ball that swallows the apex, passes
/// through it or touches the wall, none for a ball the cone misses or that
/// sits inside it clear of the wall.
///
/// Returns `Some(vec![..])` for a coaxial pair and `None` (the ruling trace
/// or the marcher) otherwise.
///
/// # Errors
///
/// Returns [`MathError`] if a shared `Circle3D` cannot be constructed.
pub fn exact_cone_sphere(
    cone: &ConicalSurface,
    sphere: &SphericalSurface,
) -> Result<Option<Vec<ExactIntersectionCurve>>, MathError> {
    let offset = cone.apex() - sphere.center();
    let along = offset.dot(cone.axis());
    if (offset - cone.axis() * along).length() > 1e-7 {
        return Ok(None);
    }
    let lin_tol = Tolerance::new().linear;
    let (sin_a, cos_a) = cone.half_angle().sin_cos();
    let (far_sq, radius_sq) = (offset.dot(offset), sphere.radius() * sphere.radius());
    let b = 2.0 * sin_a * along;
    let (disc, far, near) = ruling_quadratic(1.0, b, far_sq - radius_sq);
    // At a tangent ball `b² − 4c` cancels to a few ulps of its terms, which
    // can land below zero: that is a touch, not a miss.
    let noise = 16.0 * f64::EPSILON * 4.0f64.mul_add(far_sq + radius_sq, b * b);
    if disc < -noise {
        return Ok(Some(vec![]));
    }
    let roots: &[f64] = if far - near < lin_tol {
        &[far]
    } else {
        &[near, far]
    };
    let mut circles = Vec::new();
    for &v in roots {
        if v * cos_a > lin_tol {
            let centre = cone.apex() + cone.axis() * (v * sin_a);
            let circle = Circle3D::new(centre, cone.axis(), v * cos_a)?;
            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))
    })))
}

/// A cone and a cylinder whose axes are not parallel, traced along the
/// cylinder's rulings. A ruling `q + t w` meets the cone's double quadric
/// `|p - apex|^2 = h^2 / sin^2(half_angle)`, with `h` the offset along the
/// cone's axis, where a quadratic in `t` vanishes. `None` (the marcher's
/// case) when a ruling meets the far nappe, where no cone face lies, when
/// the rulings run along the cone's generators, or when no ruling meets it.
fn ruling_cone_cylinder(
    cone: &ConicalSurface,
    cyl: &CylindricalSurface,
    cone_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    let (sin_t, cos_t) = cone.half_angle().sin_cos();
    if sin_t < 1e-12 || cos_t < 1e-12 {
        return None;
    }
    let (apex, d, w) = (cone.apex(), cone.axis(), cyl.axis());
    let s = 1.0 / (sin_t * sin_t);
    let alpha = w.dot(d);
    let quad = 1.0 - s * alpha * alpha;
    if quad.abs() < 1e-9 {
        return None;
    }
    let roots = |u: f64| {
        let delta = cyl.evaluate(u, 0.0) - apex;
        let (dd, dw) = (delta.dot(d), delta.dot(w));
        let b = 2.0 * (dw - s * dd * alpha);
        let c = delta.dot(delta) - s * dd * dd;
        ruling_quadratic(quad, b, c)
    };
    let lin_tol = Tolerance::new().linear;
    let far_nappe = (0..WINDOW_SCAN * RULING_SAMPLES).any(|k| {
        #[allow(clippy::cast_precision_loss)]
        let u = TAU * (k as f64 + 0.5) / (WINDOW_SCAN * RULING_SAMPLES) as f64;
        let (disc, vp, vm) = roots(u);
        disc >= -lin_tol
            && [vp, vm]
                .iter()
                .any(|&t| (cyl.evaluate(u, t) - apex).dot(d) < -lin_tol)
    });
    if far_nappe {
        return None;
    }
    let samples = ruling_samples(cyl, &roots);
    // A window of meeting rulings narrower than the sampling would vanish
    // (a cone's tip just through the wall) while the others still made
    // loops; a finer scan finds every window, and any that the sampling
    // covers thinly goes to the marcher.
    let scan = WINDOW_SCAN * RULING_SAMPLES;
    // Ruling sample `i` lies midway between scan points
    // `WINDOW_SCAN i + 7` and `WINDOW_SCAN i + 8`.
    #[allow(clippy::cast_precision_loss)]
    let meets = |k: usize| roots(TAU * ((k % scan) as f64 + 0.5) / scan as f64).0 >= -lin_tol;
    if let Some(start) = (0..scan).find(|&k| !meets(k)) {
        let mut k = start;
        while k < start + scan {
            if !meets(k) {
                k += 1;
                continue;
            }
            let first = k;
            while k < start + scan && meets(k) {
                k += 1;
            }
            let covered = (first..k)
                .filter(|&j| j % WINDOW_SCAN == WINDOW_SCAN / 2 - 1 && meets(j + 1))
                .count();
            if covered < WINDOW_MIN_SAMPLES {
                return None;
            }
        }
    }
    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(cone.project_point(p), cyl.project_point(p), cone_first)
    }))
}

/// A torus and a cylinder whose axes are not parallel, where every ruling
/// of the cylinder meets the torus the same nonzero even number of times (a
/// rod through the ring's tube): a ruling meets the torus where a quartic in
/// its parameter vanishes, and each of its roots, taken in order, sweeps one
/// closed loop around the cylinder. `None` (the marcher's case) when the
/// count varies between rulings, where the curve turns back between them,
/// checked on a scan finer than the sampling so a narrow window of missing
/// rulings is not stepped over, and for a spindle torus, whose quartic also
/// holds its inner lemon.
fn ruling_torus_cylinder(
    torus: &ToroidalSurface,
    cyl: &CylindricalSurface,
    torus_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    if cyl.axis().dot(torus.z_axis()).abs() > 1.0 - 1e-9
        || torus.minor_radius() >= torus.major_radius()
    {
        return None;
    }
    let roots = |u: f64| intersect_line_torus(torus, cyl.evaluate(u, 0.0), cyl.axis());
    let rows: Vec<Vec<f64>> = (0..RULING_SAMPLES).map(|i| roots(ruling_u(i))).collect();
    let count = rows[0].len();
    let scan = WINDOW_SCAN * RULING_SAMPLES;
    #[allow(clippy::cast_precision_loss)]
    if count == 0
        || count % 2 == 1
        || (0..scan).any(|k| roots(TAU * (k as f64 + 0.5) / scan as f64).len() != count)
    {
        return None;
    }
    let loops: Vec<Vec<Point3>> = (0..count)
        .map(|j| {
            let mut pts: Vec<Point3> = rows
                .iter()
                .enumerate()
                .map(|(i, r)| cyl.evaluate(ruling_u(i), r[j]))
                .collect();
            pts.push(pts[0]);
            pts
        })
        .collect();
    Some(fit_ruling_loops(&loops, |p| {
        in_order(torus.project_point(p), cyl.project_point(p), torus_first)
    }))
}

/// A cone and a sphere whose centre is off the cone's axis, where every
/// generator of the cone crosses the sphere twice on the cone's nappe (a pin
/// through a ball): along a generator `apex + v g` the sphere is a quadratic
/// in `v`, and each of its two roots, taken in order, sweeps one closed loop
/// around the cone. A ball holding the apex is left once by every generator,
/// one loop. `None` when the centre lies on the axis (phase FF takes
/// [`exact_cone_sphere`]'s circles there), or the apex on the sphere. A
/// sphere only some generators cross goes to [`window_cone_sphere`].
fn ruling_cone_sphere(
    cone: &ConicalSurface,
    sphere: &SphericalSurface,
    cone_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    let (apex, centre, radius) = (cone.apex(), sphere.center(), sphere.radius());
    let offset = apex - centre;
    let lin_tol = Tolerance::new().linear;
    let along = offset.dot(cone.axis());
    let across = (offset - cone.axis() * along).length();
    if across < lin_tol {
        return None;
    }
    // With the apex inside the ball, `v² + 2hv + K` has roots of opposite
    // signs (`K = |offset|² − R² < 0`): every generator leaves the ball once
    // ahead of the apex, at `v = √(h² − K) − h`, one loop round the cone.
    // Where `h > 0` the root is taken as `−K / (h + √(h² − K))`, free of the
    // cancellation.
    let k = offset.dot(offset) - radius * radius;
    if radius - offset.length() > lin_tol {
        let exit = |u: f64| {
            let h = (cone.evaluate(u, 1.0) - apex).dot(offset);
            let root = h.mul_add(h, -k).sqrt();
            cone.evaluate(u, if h > 0.0 { -k / (h + root) } else { root - h })
        };
        // A ball that barely holds the apex turns the loop sharply where it
        // passes the apex, at the scale of `√−K`: a step whose midpoint
        // strays from its chord by a hundredth of the chord (a turn of about
        // 0.08 radians) is halved until the turn is resolved.
        let mut samples: Vec<(f64, Point3)> = (0..=RULING_SAMPLES)
            .map(|i| (ruling_u(i), exit(ruling_u(i))))
            .collect();
        for _ in 0..10 {
            let mut refined = Vec::with_capacity(2 * samples.len());
            for pair in samples.windows(2) {
                let ((u0, p0), (u1, p1)) = (pair[0], pair[1]);
                refined.push(pair[0]);
                let um = 0.5 * (u0 + u1);
                let pm = exit(um);
                let chord = (p1 - p0).length();
                if chord > lin_tol && (pm - (p0 + (p1 - p0) * 0.5)).length() > 0.01 * chord {
                    refined.push((um, pm));
                }
            }
            refined.extend(samples.last().copied());
            if refined.len() == samples.len() {
                break;
            }
            samples = refined;
        }
        let mut pts: Vec<Point3> = samples.iter().map(|&(_, p)| p).collect();
        if let Some(last) = pts.last_mut() {
            *last = samples[0].1;
        }
        return Some(fit_ruling_loops(&[pts], |p| {
            in_order(cone.project_point(p), sphere.project_point(p), cone_first)
        }));
    }
    // Both roots along a generator of unit direction `g` (`v` the distance
    // from the apex), nearer first, from `h = g·offset`, when it crosses the
    // sphere twice ahead of the apex.
    let crossing = |h: f64| {
        let (disc, vp, vm) = ruling_quadratic(1.0, 2.0 * h, k);
        (disc > lin_tol && vm >= lin_tol).then_some((vm, vp))
    };
    // Around the cone `h` runs between these bounds. Where both roots lie
    // ahead of the apex, raising `h` shrinks the discriminant and moves the
    // nearer root out, so every generator crosses when the two extreme ones
    // do.
    let (sin_a, cos_a) = cone.half_angle().sin_cos();
    if crossing(sin_a.mul_add(along, cos_a * across)).is_none()
        || crossing(sin_a.mul_add(along, -cos_a * across)).is_none()
    {
        return window_cone_sphere(cone, sphere, cone_first);
    }
    let rows: Vec<(f64, f64)> = (0..RULING_SAMPLES)
        .map(|i| crossing((cone.evaluate(ruling_u(i), 1.0) - apex).dot(offset)))
        .collect::<Option<_>>()?;
    let loops: Vec<Vec<Point3>> = [0, 1]
        .iter()
        .map(|&j| {
            let mut pts: Vec<Point3> = rows
                .iter()
                .enumerate()
                .map(|(i, &(near, far))| {
                    cone.evaluate(ruling_u(i), if j == 0 { near } else { far })
                })
                .collect();
            pts.push(pts[0]);
            pts
        })
        .collect();
    Some(fit_ruling_loops(&loops, |p| {
        in_order(cone.project_point(p), sphere.project_point(p), cone_first)
    }))
}

/// A cone and a sphere off its axis that only some generators cross (a ball
/// against the cone's side), with the apex outside the sphere: along the
/// generator at `u`, `h(u) = g·(apex − C)` is `c + A cos(u − φ)`, and the
/// generator crosses twice ahead of the apex where `h ≤ −√K`, with
/// `K = |apex − C|² − R²`. That is one arc of `u`, whose ends are where the
/// generator touches the sphere; the section is one loop, the nearer roots
/// out along the arc and the farther ones back. Sampled at
/// `u = mid − half·cos θ`, the roots' split `±√(h² − K)` changes sign with
/// `sin θ` and the loop stays smooth through both touching generators.
/// `Some(empty)` when no generator crosses ahead of the apex. `None`, the
/// marcher's case, when the apex lies within the sphere, when the centre is
/// too near the axis to place the arc, and when every generator reaches the
/// sphere (the extremes' test then failed on a touching generator or a root
/// at the apex).
fn window_cone_sphere(
    cone: &ConicalSurface,
    sphere: &SphericalSurface,
    cone_first: bool,
) -> Option<Vec<IntersectionCurve>> {
    let offset = cone.apex() - sphere.center();
    let lin_tol = Tolerance::new().linear;
    if offset.length() - sphere.radius() <= lin_tol {
        return None;
    }
    let k = offset.dot(offset) - sphere.radius() * sphere.radius();
    let (sin_a, cos_a) = cone.half_angle().sin_cos();
    let (ox, oy) = (offset.dot(cone.x_axis()), offset.dot(cone.y_axis()));
    let (c, a) = (sin_a * offset.dot(cone.axis()), cos_a * ox.hypot(oy));
    if a < lin_tol {
        return None;
    }
    let reach = (-k.sqrt() - c) / a;
    if reach <= -1.0 {
        return Some(Vec::new());
    }
    if reach >= 1.0 {
        return None;
    }
    let (mid, half) = (oy.atan2(ox) + std::f64::consts::PI, reach.acos());
    let half = std::f64::consts::PI - half;
    let n = RULING_SAMPLES;
    let mut pts: Vec<Point3> = (0..n)
        .map(|i| {
            #[allow(clippy::cast_precision_loss)]
            let theta = TAU * i as f64 / n as f64;
            let u = half.mul_add(-theta.cos(), mid);
            let h = a.mul_add((u - mid + std::f64::consts::PI).cos(), c);
            let split = h.mul_add(h, -k).max(0.0).sqrt();
            cone.evaluate(u, -h - split.copysign(theta.sin()))
        })
        .collect();
    pts.push(pts[0]);
    Some(fit_ruling_loops(&[pts], |p| {
        in_order(cone.project_point(p), sphere.project_point(p), cone_first)
    }))
}

/// Scan points per ruling sample when looking for windows of meeting
/// rulings, and the fewest samples a window needs to be fit.
const WINDOW_SCAN: usize = 16;
const WINDOW_MIN_SAMPLES: usize = 8;

/// 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 or oblique: `ruling_cone_cylinder` traces it.
    }

    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);
    }
    let (turn_lo, turn_hi) = (v_min, v_max);
    v_min = lo.max(v_min);
    v_max = hi.min(v_max);
    if v_max - v_min <= 1e-12 {
        return Ok(Some(vec![]));
    }
    // A cylinder beside the axis, with both turning points inside both faces,
    // meets the cone in a whole loop: traced along the cylinder's rulings,
    // each of which meets the nappe once, it comes back as one closed curve,
    // which the band split of the cylinder it winds needs. Around the axis
    // the loop winds the cone too and stays in two branches, and so does a
    // wall passing so near the axis that the loop's bend there, about
    // `(d - r) / sqrt(d r)` of a turn wide, spans fewer than a few rulings.
    let slack = Tolerance::new().linear;
    #[allow(clippy::cast_precision_loss)]
    let resolved = d - r > 3.0 * (d * r).sqrt() * TAU / RULING_SAMPLES as f64;
    if v_min <= turn_lo + slack && v_max >= turn_hi - slack && resolved {
        let mut pts: Vec<Point3> = (0..RULING_SAMPLES)
            .map(|i| {
                let (sin_u, cos_u) = ruling_u(i).sin_cos();
                let foot = cyl.origin() + (cyl.x_axis() * cos_u + cyl.y_axis() * sin_u) * r;
                let off = foot - apex;
                let across = off - axis * off.dot(axis);
                apex + across + axis * (across.length() * sin_t / cos_t)
            })
            .collect();
        pts.push(pts[0]);
        return Ok(Some(fit_ruling_loops(&[pts], |p| {
            (cone.project_point(p), cyl.project_point(p))
        })));
    }

    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,
    region: Option<Aabb3>,
) -> 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;
            }
            if region.is_some_and(|r| !r.contains_point(mid)) {
                points.push(mid);
                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;

    /// Arcs of a cone's hyperbola (a plane parallel to the axis) and parabola
    /// (a plane parallel to a ruling) between two of their sampled points
    /// stay on both the plane and the cone everywhere, not just at samples.
    #[test]
    fn plane_cone_conic_arcs_lie_on_both_surfaces() {
        let half_angle = 1.1_f64;
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            half_angle,
        )
        .unwrap();
        let ruling = Vec3::new(half_angle.sin(), 0.0, half_angle.cos());
        for (normal, d) in [(Vec3::new(1.0, 0.0, 0.0), 0.5), (ruling, 1.0)] {
            let chains =
                exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, d, 10.0)
                    .unwrap();
            let chain = chains
                .iter()
                .find_map(|c| match c {
                    ExactIntersectionCurve::Points(chain) => Some(chain),
                    _ => None,
                })
                .expect("a parabola or hyperbola section is sampled");
            let (from, to) = (chain[2], chain[chain.len() - 3]);
            let arc = plane_cone_conic_arc(&cone, normal, d, from, to)
                .unwrap()
                .expect("an exact arc");
            let (t0, t1) = arc.domain();
            assert!((arc.evaluate(t0) - from).length() < 1e-12);
            assert!((arc.evaluate(t1) - to).length() < 1e-12);
            for i in 0..=200 {
                let q = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0);
                let w = q - Point3::new(0.0, 0.0, 0.0);
                let off_plane = (normal.dot(w) - d).abs();
                let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
                assert!(off_plane < 1e-9, "off the plane by {off_plane}");
                assert!(off_cone < 1e-9, "off the cone by {off_cone}");
            }
        }
    }

    /// A plane a few 1e-10 short of parallel to a ruling cuts a vast ellipse
    /// that the parabola's closed form only approximates: the arc is either
    /// declined or on the cone, and an arc with coincident ends is declined.
    #[test]
    fn plane_cone_conic_arc_declines_a_near_parabolic_ellipse() {
        let half_angle = 1.1_f64;
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
            half_angle,
        )
        .unwrap();
        for shortfall in [1e-10, 3e-10, 8e-10] {
            let tilt = half_angle - shortfall / (2.0 * half_angle).sin();
            let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
            let chains =
                exact_plane_analytic_reaching(AnalyticSurface::Cone(&cone), normal, 1.0, 10.0)
                    .unwrap();
            let Some(chain) = chains.iter().find_map(|c| match c {
                ExactIntersectionCurve::Points(chain) => Some(chain),
                _ => None,
            }) else {
                continue;
            };
            let (from, to) = (chain[2], chain[chain.len() - 3]);
            assert!(
                plane_cone_conic_arc(&cone, normal, 1.0, from, from)
                    .unwrap()
                    .is_none(),
                "coincident ends"
            );
            let Some(arc) = plane_cone_conic_arc(&cone, normal, 1.0, from, to).unwrap() else {
                continue;
            };
            let (t0, t1) = arc.domain();
            for i in 0..=200 {
                let w = arc.evaluate(t0 + (t1 - t0) * f64::from(i) / 200.0)
                    - Point3::new(0.0, 0.0, 0.0);
                let off_cone = (w.z() - w.length() * half_angle.sin()).abs();
                assert!(off_cone < 1e-8, "{shortfall}: off the cone by {off_cone}");
            }
        }
    }

    #[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());
    }

    /// A 45 degree cone and a radius 0.1 tube tilted 40 degrees that pierces
    /// both nappes, so the ruling path declines and the general marcher runs:
    /// the tube's axis meets the near nappe at (0.1, -1.368, 1.370).
    fn cone_and_tilted_tube() -> (ConicalSurface, CylindricalSurface) {
        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 (s, c) = 40.0_f64.to_radians().sin_cos();
        let tube =
            CylindricalSurface::new(Point3::new(0.1, 0.0, 3.0), Vec3::new(0.0, s, c), 0.1).unwrap();
        (cone, tube)
    }

    #[test]
    fn marcher_keeps_to_its_region() {
        let (cone, tube) = cone_and_tilted_tube();
        let run = |region: Option<Aabb3>| {
            let (a, b) = (
                AnalyticSurface::Cone(&cone),
                AnalyticSurface::Cylinder(&tube),
            );
            let (va, vb) = (Some((0.5, 4.0)), Some((-5.0, 5.0)));
            match region {
                Some(r) => intersect_analytic_analytic_in_region(a, b, 32, va, vb, r),
                None => intersect_analytic_analytic_bounded(a, b, 32, va, vb),
            }
            .unwrap()
        };
        assert!(!run(None).is_empty());
        let near = Aabb3 {
            min: Point3::new(-0.5, -2.0, 0.8),
            max: Point3::new(0.7, -0.7, 2.0),
        };
        let curves = run(Some(near));
        assert!(!curves.is_empty(), "the loop through the region is kept");
        // The region's margin (two grid cells) and one marching step past it.
        let reach = near.expanded(0.6);
        for curve in &curves {
            assert!(curve.points.iter().all(|p| reach.contains_point(p.point)));
        }
        let away = Aabb3 {
            min: Point3::new(5.0, 5.0, 5.0),
            max: Point3::new(6.0, 6.0, 6.0),
        };
        assert!(run(Some(away)).is_empty(), "nothing is marched outside it");
    }

    #[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 oblique_cone_cylinder_traces_curves_on_both() {
        use crate::traits::ParametricCurve;
        // A pointed cone opening down from (0, 0, 3), radius half the depth,
        // and a rod along y through (x0, ., 1): one loop through the wall
        // when the rod pokes out, two when every ruling meets the cone.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 3.0),
            Vec3::new(0.0, 0.0, -1.0),
            2.0_f64.atan(),
        )
        .unwrap();
        for (x0, loops) in [(0.5, 1), (0.0, 2)] {
            let cyl =
                CylindricalSurface::new(Point3::new(x0, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
                    .unwrap();
            for cone_first in [true, false] {
                let (a, b) = if cone_first {
                    (
                        AnalyticSurface::Cone(&cone),
                        AnalyticSurface::Cylinder(&cyl),
                    )
                } else {
                    (
                        AnalyticSurface::Cylinder(&cyl),
                        AnalyticSurface::Cone(&cone),
                    )
                };
                let curves = intersect_analytic_analytic(a, b, 32).unwrap();
                assert_eq!(curves.len(), loops, "x0 {x0}: loops");
                for c in &curves {
                    let (t0, t1) = c.curve.domain();
                    for k in 0..=64 {
                        let t = (t1 - t0).mul_add(f64::from(k) / 64.0, t0);
                        let p = ParametricCurve::evaluate(&c.curve, t);
                        // A cubic through the ruling samples, bent most at the
                        // loop's branch points.
                        let rod = (p.x() - x0).hypot(p.z() - 1.0);
                        assert!(
                            (rod - 0.6).abs() < 1e-4,
                            "x0 {x0}: off the rod by {}",
                            rod - 0.6
                        );
                        let cone_r = p.x().hypot(p.y());
                        assert!(
                            (cone_r - 0.5 * (3.0 - p.z())).abs() < 1e-4,
                            "x0 {x0}: off the cone at {p:?}"
                        );
                    }
                }
            }
        }
    }

    #[test]
    fn a_rod_through_a_rings_tube_traces_four_loops() {
        use crate::traits::ParametricCurve;
        let ring = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
        // Along y through (0.5, ., 0.3): every ruling enters and leaves the
        // tube on either side of the hole, four roots, four loops.
        let rod =
            CylindricalSurface::new(Point3::new(0.5, 0.0, 0.3), Vec3::new(0.0, 1.0, 0.0), 0.6)
                .unwrap();
        let curves = ruling_torus_cylinder(&ring, &rod, true).unwrap();
        assert_eq!(curves.len(), 4);
        for c in &curves {
            let (t0, t1) = c.curve.domain();
            for k in 0..=64 {
                let p =
                    ParametricCurve::evaluate(&c.curve, (t1 - t0).mul_add(f64::from(k) / 64.0, t0));
                let on_rod = (p.x() - 0.5).hypot(p.z() - 0.3) - 0.6;
                let on_ring = (p.x().hypot(p.y()) - 4.0).hypot(p.z()) - 1.5;
                assert!(
                    on_rod.abs() < 1e-4 && on_ring.abs() < 1e-4,
                    "off by {on_rod}, {on_ring}"
                );
            }
        }
        // Higher, the rod's top rulings pass over the tube: the count varies.
        let high =
            CylindricalSurface::new(Point3::new(0.5, 0.0, 1.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
                .unwrap();
        assert!(ruling_torus_cylinder(&ring, &high, true).is_none());
        // Its top clears the tube only between two sampled rulings.
        let grazing =
            CylindricalSurface::new(Point3::new(0.5, 0.0, 0.9001), Vec3::new(0.0, 1.0, 0.0), 0.6)
                .unwrap();
        assert!(ruling_torus_cylinder(&ring, &grazing, true).is_none());
        // A spindle torus's quartic also holds its inner lemon.
        let spindle = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 1.0, 2.0).unwrap();
        let thin =
            CylindricalSurface::new(Point3::new(0.3, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), 0.2)
                .unwrap();
        assert!(ruling_torus_cylinder(&spindle, &thin, true).is_none());
    }

    #[test]
    fn a_pin_through_a_ball_traces_two_loops() {
        use crate::traits::ParametricCurve;
        let ball = SphericalSurface::new(Point3::new(0.0, 0.0, 0.0), 3.0).unwrap();
        // A pin tapering from radius 1.2 at (1, 0.5, -5) to 0.4 at 10 up:
        // radial 0.08 per unit of height.
        let half = 0.08_f64.atan();
        let apex = Point3::new(1.0, 0.5, -5.0 + 1.2 / 0.08);
        let pin = ConicalSurface::new(apex, Vec3::new(0.0, 0.0, -1.0), FRAC_PI_2 - half).unwrap();
        for cone_first in [true, false] {
            let (a, b) = if cone_first {
                (AnalyticSurface::Cone(&pin), AnalyticSurface::Sphere(&ball))
            } else {
                (AnalyticSurface::Sphere(&ball), AnalyticSurface::Cone(&pin))
            };
            let curves = intersect_analytic_analytic(a, b, 32).unwrap();
            assert_eq!(curves.len(), 2, "entry and exit loops");
            for c in &curves {
                let (t0, t1) = c.curve.domain();
                for k in 0..=64 {
                    let p = ParametricCurve::evaluate(
                        &c.curve,
                        (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
                    );
                    let on_ball = (p - Point3::new(0.0, 0.0, 0.0)).length() - 3.0;
                    let axial = apex.z() - p.z();
                    let on_pin = (p.x() - 1.0).hypot(p.y() - 0.5) - axial * half.tan();
                    assert!(
                        on_ball.abs() < 1e-4 && on_pin.abs() < 1e-4,
                        "off by {on_ball}, {on_pin}"
                    );
                }
            }
        }
        // On the ball's axis the loops are circles, which defer. A pin whose
        // apex is inside the ball leaves it once along every generator, and
        // one only partly through the ball meets it over the generators that
        // reach it: one loop each. One that opens away meets it not at all.
        let coaxial =
            ConicalSurface::new(Point3::new(0.0, 0.0, 10.0), Vec3::new(0.0, 0.0, -1.0), 1.4)
                .unwrap();
        assert!(ruling_cone_sphere(&coaxial, &ball, true).is_none());
        let aside = ConicalSurface::new(
            Point3::new(2.8, 0.0, 10.0),
            Vec3::new(0.0, 0.0, -1.0),
            FRAC_PI_2 - half,
        )
        .unwrap();
        assert_eq!(ruling_cone_sphere(&aside, &ball, true).unwrap().len(), 1);
        let holding = ConicalSurface::new(
            Point3::new(1.0, 0.5, 1.0),
            Vec3::new(0.0, 0.0, -1.0),
            FRAC_PI_2 - half,
        )
        .unwrap();
        assert_eq!(ruling_cone_sphere(&holding, &ball, true).unwrap().len(), 1);
        let away = ConicalSurface::new(
            Point3::new(1.0, 0.5, 10.0),
            Vec3::new(0.0, 0.0, 1.0),
            FRAC_PI_2 - half,
        )
        .unwrap();
        assert!(ruling_cone_sphere(&away, &ball, true).unwrap().is_empty());
        // Grazing: the generators that miss span about 0.002 radians,
        // narrower than a 2048-angle scan's step, and the near and far
        // roots join into one loop across them.
        let step = TAU / 2048.0;
        let grazed =
            SphericalSurface::new(Point3::new(step.cos(), step.sin(), 10.0), 9.255_250_971_8)
                .unwrap();
        let wide =
            ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.5).unwrap();
        assert_eq!(ruling_cone_sphere(&wide, &grazed, true).unwrap().len(), 1);
    }

    #[test]
    fn a_ball_beside_a_cone_meets_it_in_one_loop() {
        use crate::traits::ParametricCurve;
        // Apex 3 up, opening downward, the radius half the depth below it.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 3.0),
            Vec3::new(0.0, 0.0, -1.0),
            2.0_f64.atan(),
        )
        .unwrap();
        for (centre, radius) in [
            (Point3::new(1.0, 0.8, 1.2), 1.1),
            (Point3::new(1.5, 0.0, 0.0), 0.8),
        ] {
            let ball = SphericalSurface::new(centre, radius).unwrap();
            let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
            assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
            let (t0, t1) = curves[0].curve.domain();
            for k in 0..=64 {
                let p = ParametricCurve::evaluate(
                    &curves[0].curve,
                    (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
                );
                let on_ball = (p - centre).length() - radius;
                let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
                assert!(
                    on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5,
                    "ball at {centre:?}: off by {on_ball}, {on_cone}"
                );
            }
        }
        // A ball the nappe's generators miss, and one on the apex.
        let clear = SphericalSurface::new(Point3::new(4.0, 0.0, 0.0), 0.5).unwrap();
        assert!(ruling_cone_sphere(&cone, &clear, true).unwrap().is_empty());
        let on_apex = SphericalSurface::new(Point3::new(0.6, 0.0, 3.8), 1.0).unwrap();
        assert!(ruling_cone_sphere(&cone, &on_apex, true).is_none());
    }

    #[test]
    fn a_ball_holding_a_cones_apex_meets_it_in_one_loop() {
        use crate::traits::ParametricCurve;
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 3.0),
            Vec3::new(0.0, 0.0, -1.0),
            2.0_f64.atan(),
        )
        .unwrap();
        // Off the axis, and barely holding the apex: the generators pointing
        // away from the centre leave it just past the apex.
        for (centre, radius) in [
            (Point3::new(0.5, 0.0, 2.5), 2.0),
            (Point3::new(-0.4, 0.3, 2.0), 1.5),
            (Point3::new(0.0, 0.8, 3.0), 0.8001),
            (Point3::new(0.0, 0.8, 3.0), 0.800_001),
        ] {
            let ball = SphericalSurface::new(centre, radius).unwrap();
            let curves = ruling_cone_sphere(&cone, &ball, true).unwrap();
            assert_eq!(curves.len(), 1, "one loop for the ball at {centre:?}");
            let (t0, t1) = curves[0].curve.domain();
            for k in 0..=4096 {
                let p = ParametricCurve::evaluate(
                    &curves[0].curve,
                    (t1 - t0).mul_add(f64::from(k) / 4096.0, t0),
                );
                let on_ball = (p - centre).length() - radius;
                let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
                assert!(
                    on_ball.abs() < 1e-5 && on_cone.abs() < 1e-5 && p.z() < 3.0,
                    "ball at {centre:?}: off by {on_ball}, {on_cone} at {p:?}"
                );
            }
        }
    }

    #[test]
    fn oblique_cone_cylinder_defers_where_rulings_cannot_trace_it() {
        let t = 2.0_f64.atan();
        let cone =
            ConicalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 0.0, -1.0), t).unwrap();
        // A rod through the apex meets the far nappe.
        let through_apex =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 3.0), Vec3::new(0.0, 1.0, 0.0), 0.6)
                .unwrap();
        assert!(ruling_cone_cylinder(&cone, &through_apex, true).is_none());
        // A rod along a generator meets each ruling once.
        let generator = Vec3::new(t.cos(), 0.0, -t.sin());
        let along = CylindricalSurface::new(Point3::new(0.0, 0.3, 0.0), generator, 0.2).unwrap();
        assert!(ruling_cone_cylinder(&cone, &along, true).is_none());
        // A pin's tip just through a tube's wall: the tube's rulings that
        // meet it span a window narrower than the sampling.
        let pin =
            ConicalSurface::new(Point3::new(20.5, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), t).unwrap();
        let tube =
            CylindricalSurface::new(Point3::new(0.0, 0.0, -10.0), Vec3::new(0.0, 0.0, 1.0), 20.0)
                .unwrap();
        assert!(ruling_cone_cylinder(&pin, &tube, true).is_none());
    }

    #[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());
            }
        }
    }

    #[test]
    fn parallel_rod_through_a_cones_wall_closes_one_loop() {
        use crate::traits::ParametricCurve;
        // Apex 3 up, opening downward, the radius half the depth below it.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 3.0),
            Vec3::new(0.0, 0.0, -1.0),
            2.0_f64.atan(),
        )
        .unwrap();
        // Beside the axis: each ruling meets the nappe once.
        for (x, y) in [(0.0, 1.3), (1.2, 0.5)] {
            let rod =
                CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
                    .unwrap();
            let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
                .unwrap()
                .unwrap();
            assert_eq!(curves.len(), 1, "one closed loop at ({x}, {y})");
            let (t0, t1) = curves[0].curve.domain();
            let (first, last) = (
                ParametricCurve::evaluate(&curves[0].curve, t0),
                ParametricCurve::evaluate(&curves[0].curve, t1),
            );
            assert!((first - last).length() < 1e-9, "open at ({x}, {y})");
            for k in 0..=64 {
                let p = ParametricCurve::evaluate(
                    &curves[0].curve,
                    (t1 - t0).mul_add(f64::from(k) / 64.0, t0),
                );
                let on_rod = (p.x() - x).hypot(p.y() - y) - 0.6;
                let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
                assert!(
                    on_rod.abs() < 1e-5 && on_cone.abs() < 1e-5,
                    "({x}, {y}): off by {on_rod}, {on_cone}"
                );
            }
        }
        // Around the axis, and beside it by less than the bend the rulings
        // resolve: two branches.
        for (x, y) in [(0.3, 0.2), (0.65, 0.0)] {
            let rod =
                CylindricalSurface::new(Point3::new(x, y, -10.0), Vec3::new(0.0, 0.0, 1.0), 0.6)
                    .unwrap();
            let curves = algebraic_parallel_cone_cylinder(&cone, &rod, None, None)
                .unwrap()
                .unwrap();
            assert_eq!(curves.len(), 2, "two branches at ({x}, {y})");
        }
    }

    /// 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: its two branches touch at the node on
        // the inner equator. It must stay OPEN so a self-touching curve is
        // never sealed into a simple loop.
        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})"
        );
    }

    /// A wall parallel to the axis cuts one loop from the ring past the
    /// inner equator (its two turns where the branches meet), and two loops
    /// winding the tube inside it; each closes, at any scale, however near
    /// the wall runs to an equator.
    #[test]
    fn plane_torus_wall_sections_close_into_their_loops() {
        for (major, minor) in [(4.0, 1.5), (100.0, 30.0), (0.05, 0.01)] {
            let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
            for k in 1..200 {
                let (d, want) = match k.cmp(&100) {
                    // Within the inner equator: two loops winding the tube.
                    std::cmp::Ordering::Less => ((major - minor) * f64::from(k) / 100.0, 2),
                    // Past it, short of the outer equator: one loop.
                    std::cmp::Ordering::Greater => (
                        2.0f64.mul_add(minor * f64::from(k - 100) / 100.0, major - minor),
                        1,
                    ),
                    std::cmp::Ordering::Equal => continue,
                };
                let loops = plane_torus_loops(&torus, Vec3::new(1.0, 0.0, 0.0), d, 128);
                let closed = loops
                    .iter()
                    .filter(|l| (l[0].point - l[l.len() - 1].point).length() < 1e-12)
                    .count();
                assert_eq!(
                    (loops.len(), closed),
                    (want, want),
                    "R {major} r {minor}, wall at {d}"
                );
            }
        }
    }

    /// A plane through the centre tilted a little from the equator cuts two
    /// loops that each run all the way round the axis on a short run of `v`;
    /// their fitted curves stay on the torus.
    #[test]
    fn plane_torus_sections_round_the_axis_stay_on_the_torus() {
        let (major, minor) = (4.0, 1.5);
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), major, minor).unwrap();
        for tilt in [0.03_f64, 0.08, 0.2] {
            let normal = Vec3::new(tilt.sin(), 0.0, tilt.cos());
            let curves = intersect_plane_torus(&torus, normal, 0.0).unwrap();
            assert_eq!(curves.len(), 2, "tilt {tilt}");
            for c in &curves {
                let (t0, t1) = c.curve.domain();
                let off = (0..=400)
                    .map(|k| {
                        let p = c
                            .curve
                            .evaluate((t1 - t0).mul_add(f64::from(k) / 400.0, t0));
                        (p.x().hypot(p.y()) - major).hypot(p.z()) - minor
                    })
                    .fold(0.0_f64, |m, e| m.max(e.abs()));
                assert!(
                    off < 1e-6,
                    "tilt {tilt}: fitted section {off} off the torus"
                );
            }
        }
    }

    #[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"
        );
    }

    #[test]
    fn a_ball_on_a_cones_axis_meets_it_in_circles() {
        // Apex 3 up, opening downward, the radius half the depth below it.
        let cone = ConicalSurface::new(
            Point3::new(0.0, 0.0, 3.0),
            Vec3::new(0.0, 0.0, -1.0),
            2.0_f64.atan(),
        )
        .unwrap();
        for (height, radius, count) in [
            (0.0, 2.0, 2), // through the cone's middle
            (2.5, 1.3, 1), // swallowing the apex
            (2.5, 0.5, 1), // through the apex: the apex root is no circle
            (0.0, 1.0, 0), // inside, clear of the wall
            (5.0, 1.0, 0), // on the other nappe's side
        ] {
            let centre = Point3::new(0.0, 0.0, height);
            let ball = SphericalSurface::new(centre, radius).unwrap();
            let curves = exact_cone_sphere(&cone, &ball).unwrap().unwrap();
            let circles = circles_of(&curves);
            assert_eq!(circles.len(), count, "ball at {height}, radius {radius}");
            for circle in circles {
                for k in 0..16 {
                    let p = circle.evaluate(TAU * f64::from(k) / 16.0);
                    let on_ball = (p - centre).length() - radius;
                    let on_cone = p.x().hypot(p.y()) - (3.0 - p.z()) / 2.0;
                    assert!(
                        on_ball.abs() < 1e-9 && on_cone.abs() < 1e-9,
                        "ball at {height}: off by {on_ball}, {on_cone}"
                    );
                }
            }
        }
        let aside = SphericalSurface::new(Point3::new(0.5, 0.0, 0.0), 2.0).unwrap();
        assert!(exact_cone_sphere(&cone, &aside).unwrap().is_none());
        // A ball touching a wide cone a million units out, where the
        // discriminant cancels to a few ulps below zero.
        let wide =
            ConicalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), 0.3).unwrap();
        let far = SphericalSurface::new(Point3::new(0.0, 0.0, 1e6), 1e6 * 0.3_f64.cos()).unwrap();
        let curves = exact_cone_sphere(&wide, &far).unwrap().unwrap();
        let circles = circles_of(&curves);
        assert_eq!(circles.len(), 1, "the touch");
        let touch = 1e6 * 0.3_f64.sin() * 0.3_f64.cos();
        assert!(
            (circles[0].radius() - touch).abs() < 1e-3,
            "{}",
            circles[0].radius()
        );
    }

    /// The circles among exact section curves.
    fn circles_of(curves: &[ExactIntersectionCurve]) -> Vec<&Circle3D> {
        curves
            .iter()
            .filter_map(|c| match c {
                ExactIntersectionCurve::Circle(circle) => Some(circle),
                _ => None,
            })
            .collect()
    }

    /// Worst distance of a circle's points from a torus and from a second
    /// surface given by its own distance function.
    fn worst_off(
        circles: &[&Circle3D],
        torus: &ToroidalSurface,
        other: impl Fn(Point3) -> f64,
    ) -> f64 {
        let mut worst = 0.0_f64;
        for circle in circles {
            for k in 0..16 {
                let p = circle.evaluate(TAU * f64::from(k) / 16.0);
                let q = p - torus.center();
                let along = q.dot(torus.z_axis());
                let rho = (q - torus.z_axis() * along).length();
                let off = ((rho - torus.major_radius()).hypot(along) - torus.minor_radius()).abs();
                worst = worst.max(off).max(other(p).abs());
            }
        }
        worst
    }

    #[test]
    fn exact_sphere_torus_meets_a_ball_on_the_axis_in_circles() {
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
        for height in [0.0, 1.0] {
            let centre = Point3::new(0.0, 0.0, height);
            let sphere = SphericalSurface::new(centre, 3.0).unwrap();
            let curves = exact_sphere_torus(&sphere, &torus).unwrap().unwrap();
            let circles = circles_of(&curves);
            assert_eq!((curves.len(), circles.len()), (2, 2), "height {height}");
            let worst = worst_off(&circles, &torus, |p| (p - centre).length() - 3.0);
            assert!(worst < 1e-9, "height {height}: {worst}");
        }
    }

    #[test]
    fn exact_sphere_torus_misses_touches_and_defers() {
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
        let ball = |x: f64, r: f64| SphericalSurface::new(Point3::new(x, 0.0, 0.0), r).unwrap();
        assert!(
            exact_sphere_torus(&ball(0.0, 1.0), &torus)
                .unwrap()
                .unwrap()
                .is_empty(),
            "a small ball in the hole misses"
        );
        assert!(
            exact_sphere_torus(&ball(0.0, 2.5), &torus)
                .unwrap()
                .is_none(),
            "a ball touching the inner equator defers"
        );
        assert!(
            exact_sphere_torus(&ball(1.0, 3.0), &torus)
                .unwrap()
                .is_none(),
            "a ball off the axis defers"
        );
        let spindle = ToroidalSurface::with_axis_and_ref_dir(
            Point3::new(0.0, 0.0, 0.0),
            1.0,
            2.0,
            Vec3::new(0.0, 0.0, 1.0),
            Vec3::new(1.0, 0.0, 0.0),
        )
        .unwrap();
        assert!(
            exact_sphere_torus(&ball(0.0, 2.5), &spindle)
                .unwrap()
                .is_none()
        );
    }

    #[test]
    fn exact_cylinder_torus_meets_a_coaxial_rod_in_circles() {
        let torus = ToroidalSurface::new(Point3::new(0.0, 0.0, 0.0), 4.0, 1.5).unwrap();
        let z = Vec3::new(0.0, 0.0, 1.0);
        let rod = |r: f64| CylindricalSurface::new(Point3::new(0.0, 0.0, -5.0), z, r).unwrap();
        let curves = exact_cylinder_torus(&rod(4.2), &torus).unwrap().unwrap();
        let circles = circles_of(&curves);
        assert_eq!((curves.len(), circles.len()), (2, 2));
        let worst = worst_off(&circles, &torus, |p| p.x().hypot(p.y()) - 4.2);
        assert!(worst < 1e-9, "{worst}");
        assert!(
            exact_cylinder_torus(&rod(2.0), &torus)
                .unwrap()
                .unwrap()
                .is_empty(),
            "a rod clear in the hole misses"
        );
        assert!(
            exact_cylinder_torus(&rod(5.5), &torus).unwrap().is_none(),
            "a wall touching the outer equator defers"
        );
        let tilted =
            CylindricalSurface::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.1, 1.0), 4.2)
                .unwrap();
        let offset = CylindricalSurface::new(Point3::new(0.5, 0.0, 0.0), z, 4.2).unwrap();
        assert!(exact_cylinder_torus(&tilted, &torus).unwrap().is_none());
        assert!(exact_cylinder_torus(&offset, &torus).unwrap().is_none());
        let spindle = ToroidalSurface::with_axis_and_ref_dir(
            Point3::new(0.0, 0.0, 0.0),
            1.0,
            2.0,
            z,
            Vec3::new(1.0, 0.0, 0.0),
        )
        .unwrap();
        assert!(
            exact_cylinder_torus(&rod(0.5), &spindle).unwrap().is_none(),
            "a spindle torus's inner lemon also meets the rod"
        );
    }

    /// 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, 0.0).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, 0.0).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, 0.0).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, 0.0).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"
            );
        }
    }
}