use bone_types::{Angle, Length, Plane3, Point3, Tolerance, UnitVec3, Vec3};
use uom::si::angle::radian;
use uom::si::length::millimeter;

use crate::circle3::Circle3;

use super::{as_vec, basis, unit};

const PARALLEL_SIN_EPS: f64 = 1.0e-9;
const PERP_COS_EPS: f64 = 1.0e-9;
const SLOPE_EPS: f64 = 1.0e-9;

#[derive(Copy, Clone, Debug, PartialEq)]
pub struct InfiniteLine {
    point: Point3,
    direction: UnitVec3,
}

impl InfiniteLine {
    #[must_use]
    pub const fn new(point: Point3, direction: UnitVec3) -> Self {
        Self { point, direction }
    }

    #[must_use]
    pub const fn point(self) -> Point3 {
        self.point
    }

    #[must_use]
    pub const fn direction(self) -> UnitVec3 {
        self.direction
    }
}

#[derive(Copy, Clone, Debug, PartialEq)]
pub enum IntersectionCurve {
    Line(InfiniteLine),
    Circle(Circle3),
}

#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub enum SsiGap {
    Ellipse,
    GeneralQuadric,
}

#[derive(Clone, Debug, PartialEq)]
pub enum SurfaceIntersection {
    Empty,
    Coincident,
    Curves(Vec<IntersectionCurve>),
    Unsupported(SsiGap),
}

#[derive(Copy, Clone, Debug, PartialEq)]
pub enum AnalyticSurface {
    Plane {
        origin: Point3,
        normal: UnitVec3,
    },
    Cylinder {
        axis_origin: Point3,
        axis: UnitVec3,
        radius: Length,
    },
    Cone {
        apex: Point3,
        axis: UnitVec3,
        half_angle: Angle,
    },
    Sphere {
        center: Point3,
        radius: Length,
    },
    Torus {
        center: Point3,
        axis: UnitVec3,
        major_radius: Length,
        minor_radius: Length,
    },
}

#[must_use]
pub fn intersect_surfaces(
    a: AnalyticSurface,
    b: AnalyticSurface,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    use AnalyticSurface::Plane;
    if degenerate_cone(a) || degenerate_cone(b) {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let (a, b) = if rank(a) <= rank(b) { (a, b) } else { (b, a) };
    match a {
        Plane { origin, normal } => plane_against(origin, normal, b, tolerance),
        _ => quadric_against(a, b, tolerance),
    }
}

fn degenerate_cone(surface: AnalyticSurface) -> bool {
    matches!(
        surface,
        AnalyticSurface::Cone { half_angle, .. }
            if half_angle.get::<radian>().abs() >= core::f64::consts::FRAC_PI_2
    )
}

fn plane_against(
    origin: Point3,
    normal: UnitVec3,
    other: AnalyticSurface,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    use AnalyticSurface::{Cone, Cylinder, Plane, Sphere, Torus};
    match other {
        Plane {
            origin: o2,
            normal: n2,
        } => plane_plane(origin, normal, o2, n2, tolerance),
        Cylinder {
            axis_origin,
            axis,
            radius,
        } => plane_cylinder(origin, normal, axis_origin, axis, radius, tolerance),
        Cone {
            apex,
            axis,
            half_angle,
        } => plane_cone(origin, normal, apex, axis, half_angle, tolerance),
        Sphere { center, radius } => plane_sphere(origin, normal, center, radius, tolerance),
        Torus {
            center,
            axis,
            major_radius,
            minor_radius,
        } => plane_torus(
            origin,
            normal,
            center,
            axis,
            major_radius,
            minor_radius,
            tolerance,
        ),
    }
}

fn quadric_against(
    a: AnalyticSurface,
    b: AnalyticSurface,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    use AnalyticSurface::{Cone, Cylinder, Sphere};
    match (a, b) {
        (
            Cylinder {
                axis_origin: c1,
                axis: a1,
                radius: r1,
            },
            Cylinder {
                axis_origin: c2,
                axis: a2,
                radius: r2,
            },
        ) => cylinder_cylinder(c1, a1, r1, c2, a2, r2, tolerance),
        (
            Cylinder {
                axis_origin,
                axis,
                radius,
            },
            Cone {
                apex,
                axis: cone_axis,
                half_angle,
            },
        ) => cylinder_cone(
            axis_origin,
            axis,
            radius,
            apex,
            cone_axis,
            half_angle,
            tolerance,
        ),
        (
            Cylinder {
                axis_origin,
                axis,
                radius,
            },
            Sphere { center, radius: rs },
        ) => cylinder_sphere(axis_origin, axis, radius, center, rs, tolerance),
        (
            Cone {
                apex: p1,
                axis: a1,
                half_angle: h1,
            },
            Cone {
                apex: p2,
                axis: a2,
                half_angle: h2,
            },
        ) => cone_cone(p1, a1, h1, p2, a2, h2, tolerance),
        (
            Cone {
                apex,
                axis,
                half_angle,
            },
            Sphere { center, radius },
        ) => cone_sphere(apex, axis, half_angle, center, radius, tolerance),
        (
            Sphere {
                center: s1,
                radius: r1,
            },
            Sphere {
                center: s2,
                radius: r2,
            },
        ) => sphere_sphere(s1, r1, s2, r2, tolerance),
        _ => SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric),
    }
}

fn rank(surface: AnalyticSurface) -> u8 {
    use AnalyticSurface::{Cone, Cylinder, Plane, Sphere, Torus};
    match surface {
        Plane { .. } => 0,
        Cylinder { .. } => 1,
        Cone { .. } => 2,
        Sphere { .. } => 3,
        Torus { .. } => 4,
    }
}

fn position(p: Point3) -> Vec3 {
    let (x, y, z) = p.coords_mm();
    Vec3::from_mm(x, y, z)
}

fn to_point(v: Vec3) -> Point3 {
    let (x, y, z) = v.coords_mm();
    Point3::from_mm(x, y, z)
}

fn ordered(a: Point3, b: Point3) -> (Point3, Point3) {
    if a.coords_mm() <= b.coords_mm() {
        (a, b)
    } else {
        (b, a)
    }
}

fn plane_plane(
    o1: Point3,
    n1: UnitVec3,
    o2: Point3,
    n2: UnitVec3,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let nv1 = as_vec(n1);
    let nv2 = as_vec(n2);
    let cross = nv1.cross(nv2);
    if cross.norm_mm() <= PARALLEL_SIN_EPS {
        let gap = (o2 - o1).dot_mm2(nv1);
        return if gap.abs() <= eps {
            SurfaceIntersection::Coincident
        } else {
            SurfaceIntersection::Empty
        };
    }
    let Some(direction) = unit(cross, tolerance) else {
        return SurfaceIntersection::Empty;
    };
    let cos = nv1.dot_mm2(nv2);
    let h1 = position(o1).dot_mm2(nv1);
    let h2 = position(o2).dot_mm2(nv2);
    let denom = 1.0 - cos * cos;
    let a = (h1 - h2 * cos) / denom;
    let b = (h2 - h1 * cos) / denom;
    let point = to_point(nv1 * a + nv2 * b);
    SurfaceIntersection::Curves(vec![IntersectionCurve::Line(InfiniteLine::new(
        point, direction,
    ))])
}

fn plane_cylinder(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    axis_origin: Point3,
    axis: UnitVec3,
    radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let np = as_vec(plane_normal);
    let av = as_vec(axis);
    let cos = np.dot_mm2(av);

    if cos.abs() >= 1.0 - PERP_COS_EPS {
        return plane_perpendicular_to_axis(
            plane_origin,
            plane_normal,
            axis_origin,
            av,
            cos,
            radius,
            tolerance,
        );
    }
    if cos.abs() <= PERP_COS_EPS {
        return plane_parallel_to_axis(
            plane_origin,
            plane_normal,
            axis_origin,
            axis,
            radius,
            tolerance,
        );
    }
    SurfaceIntersection::Unsupported(SsiGap::Ellipse)
}

fn plane_perpendicular_to_axis(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    axis_origin: Point3,
    av: Vec3,
    cos: f64,
    radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let t = (plane_origin - axis_origin).dot_mm2(as_vec(plane_normal)) / cos;
    let center = axis_origin + av * t;
    let Some((x, y)) = basis(plane_normal, tolerance) else {
        return SurfaceIntersection::Empty;
    };
    let Ok(plane) = Plane3::new(center, x, y, tolerance) else {
        return SurfaceIntersection::Empty;
    };
    match Circle3::new(plane, radius, tolerance) {
        Ok(circle) => SurfaceIntersection::Curves(vec![IntersectionCurve::Circle(circle)]),
        Err(_) => SurfaceIntersection::Empty,
    }
}

fn plane_parallel_to_axis(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    axis_origin: Point3,
    axis: UnitVec3,
    radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let np = as_vec(plane_normal);
    let av = as_vec(axis);
    let r = radius.get::<millimeter>();
    let dist = (axis_origin - plane_origin).dot_mm2(np);
    let d_abs = dist.abs();
    if d_abs > r + eps {
        return SurfaceIntersection::Empty;
    }
    let foot = axis_origin - np * dist;
    let Some(w) = unit(av.cross(np), tolerance) else {
        return SurfaceIntersection::Empty;
    };
    let half = (r * r - d_abs * d_abs).max(0.0).sqrt();
    if half <= eps {
        return SurfaceIntersection::Curves(vec![IntersectionCurve::Line(InfiniteLine::new(
            foot, axis,
        ))]);
    }
    let wv = as_vec(w);
    let (lo, hi) = ordered(foot + wv * half, foot - wv * half);
    SurfaceIntersection::Curves(vec![
        IntersectionCurve::Line(InfiniteLine::new(lo, axis)),
        IntersectionCurve::Line(InfiniteLine::new(hi, axis)),
    ])
}

fn cylinder_cylinder(
    c1: Point3,
    a1: UnitVec3,
    r1: Length,
    c2: Point3,
    a2: UnitVec3,
    r2: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let av1 = as_vec(a1);
    let av2 = as_vec(a2);
    if av1.cross(av2).norm_mm() > PARALLEL_SIN_EPS {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let r1m = r1.get::<millimeter>();
    let r2m = r2.get::<millimeter>();
    let delta = c2 - c1;
    let along = delta.dot_mm2(av1);
    let perp = delta - av1 * along;
    let dsep = perp.norm_mm();
    if dsep <= eps {
        return if (r1m - r2m).abs() <= eps {
            SurfaceIntersection::Coincident
        } else {
            SurfaceIntersection::Empty
        };
    }
    if dsep > r1m + r2m + eps || dsep < (r1m - r2m).abs() - eps {
        return SurfaceIntersection::Empty;
    }
    let Some(u) = unit(perp, tolerance) else {
        return SurfaceIntersection::Empty;
    };
    let uv = as_vec(u);
    let x = (dsep * dsep + r1m * r1m - r2m * r2m) / (2.0 * dsep);
    let half = (r1m * r1m - x * x).max(0.0).sqrt();
    let base = c1 + uv * x;
    if half <= eps {
        return SurfaceIntersection::Curves(vec![IntersectionCurve::Line(InfiniteLine::new(
            base, a1,
        ))]);
    }
    let wv = av1.cross(uv);
    let (lo, hi) = ordered(base + wv * half, base - wv * half);
    SurfaceIntersection::Curves(vec![
        IntersectionCurve::Line(InfiniteLine::new(lo, a1)),
        IntersectionCurve::Line(InfiniteLine::new(hi, a1)),
    ])
}

fn circle_curve(
    center: Point3,
    normal: UnitVec3,
    radius_mm: f64,
    tolerance: Tolerance,
) -> Option<IntersectionCurve> {
    let (x, y) = basis(normal, tolerance)?;
    let plane = Plane3::new(center, x, y, tolerance).ok()?;
    Circle3::new(plane, Length::new::<millimeter>(radius_mm), tolerance)
        .ok()
        .map(IntersectionCurve::Circle)
}

fn axial(origin: Point3, axis: UnitVec3, point: Point3) -> f64 {
    (point - origin).dot_mm2(as_vec(axis))
}

fn axial_center(origin: Point3, axis: UnitVec3, t: f64) -> Point3 {
    origin + as_vec(axis) * t
}

fn shared_axis(
    o1: Point3,
    a1: UnitVec3,
    o2: Point3,
    a2: UnitVec3,
    tolerance: Tolerance,
) -> Option<(Point3, UnitVec3)> {
    if as_vec(a1).cross(as_vec(a2)).norm_mm() > PARALLEL_SIN_EPS {
        return None;
    }
    let delta = o2 - o1;
    let along = delta.dot_mm2(as_vec(a1));
    if (delta - as_vec(a1) * along).norm_mm() > tolerance.value() {
        return None;
    }
    Some((o1, a1))
}

fn on_axis(origin: Point3, axis: UnitVec3, point: Point3, tolerance: Tolerance) -> bool {
    let delta = point - origin;
    let along = delta.dot_mm2(as_vec(axis));
    (delta - as_vec(axis) * along).norm_mm() <= tolerance.value()
}

fn cone_slope(axis: UnitVec3, half_angle: Angle, shared: UnitVec3) -> f64 {
    half_angle.get::<radian>().tan() * as_vec(axis).dot_mm2(as_vec(shared)).signum()
}

fn solve_quadratic(a: f64, b: f64, c: f64, eps: f64) -> Vec<f64> {
    if a.abs() <= eps {
        return if b.abs() <= eps { vec![] } else { vec![-c / b] };
    }
    let disc = b * b - 4.0 * a * c;
    if disc <= 0.0 {
        return vec![];
    }
    let root = disc.sqrt();
    let lo = (-b - root) / (2.0 * a);
    let hi = (-b + root) / (2.0 * a);
    if (hi - lo).abs() <= eps {
        vec![lo.midpoint(hi)]
    } else {
        vec![lo.min(hi), lo.max(hi)]
    }
}

fn axial_circles(
    origin: Point3,
    axis: UnitVec3,
    samples: impl IntoIterator<Item = (f64, f64)>,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let curves: Vec<IntersectionCurve> = samples
        .into_iter()
        .filter_map(|(t, radius)| {
            circle_curve(axial_center(origin, axis, t), axis, radius, tolerance)
        })
        .collect();
    if curves.is_empty() {
        SurfaceIntersection::Empty
    } else {
        SurfaceIntersection::Curves(curves)
    }
}

fn plane_sphere(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    center: Point3,
    radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let n = as_vec(plane_normal);
    let signed = (center - plane_origin).dot_mm2(n);
    let r = radius.get::<millimeter>();
    let half = (r * r - signed * signed).max(0.0).sqrt();
    if half <= eps {
        return SurfaceIntersection::Empty;
    }
    let foot = center - n * signed;
    circle_curve(foot, plane_normal, half, tolerance).map_or(SurfaceIntersection::Empty, |curve| {
        SurfaceIntersection::Curves(vec![curve])
    })
}

fn sphere_sphere(
    c1: Point3,
    r1: Length,
    c2: Point3,
    r2: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let r1m = r1.get::<millimeter>();
    let r2m = r2.get::<millimeter>();
    let delta = c2 - c1;
    let d = delta.norm_mm();
    if d <= eps {
        return if (r1m - r2m).abs() <= eps {
            SurfaceIntersection::Coincident
        } else {
            SurfaceIntersection::Empty
        };
    }
    if d > r1m + r2m + eps || d < (r1m - r2m).abs() - eps {
        return SurfaceIntersection::Empty;
    }
    let Some(axis) = unit(delta, tolerance) else {
        return SurfaceIntersection::Empty;
    };
    let x = (d * d + r1m * r1m - r2m * r2m) / (2.0 * d);
    let half = (r1m * r1m - x * x).max(0.0).sqrt();
    if half <= eps {
        return SurfaceIntersection::Empty;
    }
    let center = c1 + as_vec(axis) * x;
    circle_curve(center, axis, half, tolerance).map_or(SurfaceIntersection::Empty, |curve| {
        SurfaceIntersection::Curves(vec![curve])
    })
}

fn plane_cone(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    apex: Point3,
    axis: UnitVec3,
    half_angle: Angle,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let cos = as_vec(plane_normal).dot_mm2(as_vec(axis));
    if cos.abs() < 1.0 - PERP_COS_EPS {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let tq = (plane_origin - apex).dot_mm2(as_vec(plane_normal)) / cos;
    let radius = tq * half_angle.get::<radian>().tan();
    if radius <= tolerance.value() {
        return SurfaceIntersection::Empty;
    }
    let center = axial_center(apex, axis, tq);
    circle_curve(center, axis, radius, tolerance).map_or(SurfaceIntersection::Empty, |curve| {
        SurfaceIntersection::Curves(vec![curve])
    })
}

fn plane_torus(
    plane_origin: Point3,
    plane_normal: UnitVec3,
    center: Point3,
    axis: UnitVec3,
    major_radius: Length,
    minor_radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let eps = tolerance.value();
    let cos = as_vec(plane_normal).dot_mm2(as_vec(axis));
    if cos.abs() < 1.0 - PERP_COS_EPS {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let major = major_radius.get::<millimeter>();
    let minor = minor_radius.get::<millimeter>();
    let z = axial(center, axis, plane_origin);
    if z.abs() > minor + eps {
        return SurfaceIntersection::Empty;
    }
    let half = (minor * minor - z * z).max(0.0).sqrt();
    let ring_center = axial_center(center, axis, z);
    let radii = if half <= eps {
        vec![major]
    } else {
        vec![major + half, major - half]
    };
    let curves: Vec<IntersectionCurve> = radii
        .into_iter()
        .filter(|radius| *radius > eps)
        .filter_map(|radius| circle_curve(ring_center, axis, radius, tolerance))
        .collect();
    if curves.is_empty() {
        SurfaceIntersection::Empty
    } else {
        SurfaceIntersection::Curves(curves)
    }
}

fn cylinder_sphere(
    axis_origin: Point3,
    axis: UnitVec3,
    radius: Length,
    center: Point3,
    sphere_radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    if !on_axis(axis_origin, axis, center, tolerance) {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let eps = tolerance.value();
    let rc = radius.get::<millimeter>();
    let rs = sphere_radius.get::<millimeter>();
    if rs < rc - eps {
        return SurfaceIntersection::Empty;
    }
    let tc = axial(axis_origin, axis, center);
    let half = (rs * rs - rc * rc).max(0.0).sqrt();
    let samples = if half <= eps {
        vec![(tc, rc)]
    } else {
        vec![(tc - half, rc), (tc + half, rc)]
    };
    axial_circles(axis_origin, axis, samples, tolerance)
}

fn cylinder_cone(
    axis_origin: Point3,
    axis: UnitVec3,
    radius: Length,
    apex: Point3,
    cone_axis: UnitVec3,
    half_angle: Angle,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let Some((origin, shared)) = shared_axis(axis_origin, axis, apex, cone_axis, tolerance) else {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    };
    let slope = cone_slope(cone_axis, half_angle, shared);
    if slope.abs() <= SLOPE_EPS {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let rc = radius.get::<millimeter>();
    let t_apex = axial(origin, shared, apex);
    let t = t_apex + rc / slope;
    axial_circles(origin, shared, [(t, rc)], tolerance)
}

fn cone_cone(
    apex1: Point3,
    axis1: UnitVec3,
    half1: Angle,
    apex2: Point3,
    axis2: UnitVec3,
    half2: Angle,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    let Some((origin, shared)) = shared_axis(apex1, axis1, apex2, axis2, tolerance) else {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    };
    let eps = tolerance.value();
    let slope1 = cone_slope(axis1, half1, shared);
    let slope2 = cone_slope(axis2, half2, shared);
    let t1 = axial(origin, shared, apex1);
    let t2 = axial(origin, shared, apex2);
    if (slope1 - slope2).abs() <= SLOPE_EPS {
        return if (t1 - t2).abs() <= tolerance.value() {
            SurfaceIntersection::Coincident
        } else {
            SurfaceIntersection::Empty
        };
    }
    let t = (t1 * slope1 - t2 * slope2) / (slope1 - slope2);
    let radius = (t - t1) * slope1;
    if radius < eps {
        return SurfaceIntersection::Empty;
    }
    axial_circles(origin, shared, [(t, radius)], tolerance)
}

fn cone_sphere(
    apex: Point3,
    axis: UnitVec3,
    half_angle: Angle,
    center: Point3,
    radius: Length,
    tolerance: Tolerance,
) -> SurfaceIntersection {
    if !on_axis(apex, axis, center, tolerance) {
        return SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric);
    }
    let eps = tolerance.value();
    let tan = half_angle.get::<radian>().tan();
    let tc = axial(apex, axis, center);
    let rs = radius.get::<millimeter>();
    let a = tan * tan + 1.0;
    let b = -2.0 * tc;
    let c = tc * tc - rs * rs;
    let samples: Vec<(f64, f64)> = solve_quadratic(a, b, c, eps)
        .into_iter()
        .filter_map(|t| {
            let ring = t * tan;
            (ring > eps).then_some((t, ring))
        })
        .collect();
    axial_circles(apex, axis, samples, tolerance)
}

#[cfg(test)]
mod tests {
    use super::{
        AnalyticSurface, IntersectionCurve, SsiGap, SurfaceIntersection, intersect_surfaces,
    };
    use crate::circle3::Circle3;
    use bone_types::{Angle, Length, Point3, Tolerance, UnitVec3};
    use uom::si::angle::radian;
    use uom::si::length::millimeter;

    const TOL: Tolerance = Tolerance::new(1.0e-7);

    fn mm(value: f64) -> Length {
        Length::new::<millimeter>(value)
    }

    fn rad(value: f64) -> Angle {
        Angle::new::<radian>(value)
    }

    fn point(x: f64, y: f64, z: f64) -> Point3 {
        Point3::from_mm(x, y, z)
    }

    fn dir(x: f64, y: f64, z: f64) -> UnitVec3 {
        let Ok(unit) = UnitVec3::try_from_components(x, y, z, TOL) else {
            panic!("{x},{y},{z} is a direction");
        };
        unit
    }

    fn plane(origin: Point3, normal: UnitVec3) -> AnalyticSurface {
        AnalyticSurface::Plane { origin, normal }
    }

    fn cylinder(axis_origin: Point3, axis: UnitVec3, radius: f64) -> AnalyticSurface {
        AnalyticSurface::Cylinder {
            axis_origin,
            axis,
            radius: mm(radius),
        }
    }

    fn sphere(center: Point3, radius: f64) -> AnalyticSurface {
        AnalyticSurface::Sphere {
            center,
            radius: mm(radius),
        }
    }

    fn cone(apex: Point3, axis: UnitVec3, half_angle: f64) -> AnalyticSurface {
        AnalyticSurface::Cone {
            apex,
            axis,
            half_angle: rad(half_angle),
        }
    }

    fn torus(center: Point3, axis: UnitVec3, major: f64, minor: f64) -> AnalyticSurface {
        AnalyticSurface::Torus {
            center,
            axis,
            major_radius: mm(major),
            minor_radius: mm(minor),
        }
    }

    fn single_circle(result: &SurfaceIntersection) -> Circle3 {
        let SurfaceIntersection::Curves(curves) = result else {
            panic!("expected intersection curves, got {result:?}");
        };
        let [IntersectionCurve::Circle(circle)] = curves.as_slice() else {
            panic!("expected a single circle, got {curves:?}");
        };
        *circle
    }

    fn on_plane(p: Point3, origin: Point3, normal: UnitVec3) -> bool {
        (p - origin).dot_mm2(vec(normal)).abs() <= 1.0e-5
    }

    fn vec(u: UnitVec3) -> bone_types::Vec3 {
        let (x, y, z) = u.components();
        bone_types::Vec3::from_mm(x, y, z)
    }

    fn radial_distance(p: Point3, c: Point3, a: UnitVec3) -> f64 {
        let d = p - c;
        let along = d.dot_mm2(vec(a));
        (d - vec(a) * along).norm_mm()
    }

    #[test]
    fn perpendicular_planes_meet_in_a_line() {
        let xy = plane(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0));
        let xz = plane(point(0.0, 0.0, 0.0), dir(0.0, 1.0, 0.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(xy, xz, TOL) else {
            panic!("two transversal planes meet");
        };
        let [IntersectionCurve::Line(line)] = curves.as_slice() else {
            panic!("plane-plane is a single line");
        };
        assert!(on_plane(
            line.point(),
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0)
        ));
        assert!(on_plane(
            line.point(),
            point(0.0, 0.0, 0.0),
            dir(0.0, 1.0, 0.0)
        ));
        let (dx, dy, dz) = line.direction().components();
        assert!(dx.abs() > 0.99 && dy.abs() < 1.0e-6 && dz.abs() < 1.0e-6);
    }

    #[test]
    fn offset_planes_meet_off_origin() {
        let xy = plane(point(0.0, 0.0, 2.0), dir(0.0, 0.0, 1.0));
        let xz = plane(point(0.0, 3.0, 0.0), dir(0.0, 1.0, 0.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(xy, xz, TOL) else {
            panic!("two transversal planes meet");
        };
        let [IntersectionCurve::Line(line)] = curves.as_slice() else {
            panic!("single line");
        };
        let (_, py, pz) = line.point().coords_mm();
        assert!((py - 3.0).abs() <= 1.0e-5);
        assert!((pz - 2.0).abs() <= 1.0e-5);
    }

    #[test]
    fn parallel_coincident_planes_are_coincident() {
        let a = plane(point(0.0, 0.0, 1.0), dir(0.0, 0.0, 1.0));
        let b = plane(point(5.0, -2.0, 1.0), dir(0.0, 0.0, 1.0));
        assert_eq!(
            intersect_surfaces(a, b, TOL),
            SurfaceIntersection::Coincident
        );
    }

    #[test]
    fn parallel_distinct_planes_are_empty() {
        let a = plane(point(0.0, 0.0, 1.0), dir(0.0, 0.0, 1.0));
        let b = plane(point(0.0, 0.0, 4.0), dir(0.0, 0.0, 1.0));
        assert_eq!(intersect_surfaces(a, b, TOL), SurfaceIntersection::Empty);
    }

    #[test]
    fn plane_perpendicular_to_axis_is_a_circle() {
        let cyl = cylinder(point(1.0, 1.0, 0.0), dir(0.0, 0.0, 1.0), 4.0);
        let cut = plane(point(0.0, 0.0, 5.0), dir(0.0, 0.0, 1.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(cyl, cut, TOL) else {
            panic!("cutting a cylinder square gives a circle");
        };
        let [IntersectionCurve::Circle(circle)] = curves.as_slice() else {
            panic!("expected a circle");
        };
        assert!((circle.radius_mm() - 4.0).abs() <= 1.0e-6);
        let (cx, cy, cz) = circle.center().coords_mm();
        assert!((cx - 1.0).abs() <= 1.0e-6);
        assert!((cy - 1.0).abs() <= 1.0e-6);
        assert!((cz - 5.0).abs() <= 1.0e-6);
    }

    #[test]
    fn plane_parallel_to_axis_cuts_two_lines() {
        let cyl = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        let cut = plane(point(0.0, 3.0, 0.0), dir(0.0, 1.0, 0.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(cyl, cut, TOL) else {
            panic!("a chord plane cuts two generator lines");
        };
        assert_eq!(curves.len(), 2);
        curves.iter().for_each(|curve| {
            let IntersectionCurve::Line(line) = curve else {
                panic!("expected a line");
            };
            assert!(on_plane(
                line.point(),
                point(0.0, 3.0, 0.0),
                dir(0.0, 1.0, 0.0)
            ));
            let r = radial_distance(line.point(), point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0));
            assert!((r - 5.0).abs() <= 1.0e-5, "generator lies on the cylinder");
            let (dx, dy, dz) = line.direction().components();
            assert!(dx.abs() < 1.0e-6 && dy.abs() < 1.0e-6 && dz.abs() > 0.99);
        });
    }

    #[test]
    fn plane_tangent_to_axis_is_one_line() {
        let cyl = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        let cut = plane(point(0.0, 5.0, 0.0), dir(0.0, 1.0, 0.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(cyl, cut, TOL) else {
            panic!("a tangent plane touches one generator");
        };
        assert_eq!(curves.len(), 1);
    }

    #[test]
    fn plane_clear_of_axis_is_empty() {
        let cyl = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        let cut = plane(point(0.0, 9.0, 0.0), dir(0.0, 1.0, 0.0));
        assert_eq!(
            intersect_surfaces(cyl, cut, TOL),
            SurfaceIntersection::Empty
        );
    }

    #[test]
    fn oblique_plane_cylinder_is_unsupported_ellipse() {
        let cyl = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        let cut = plane(point(0.0, 0.0, 0.0), dir(0.0, 1.0, 1.0));
        assert_eq!(
            intersect_surfaces(cyl, cut, TOL),
            SurfaceIntersection::Unsupported(SsiGap::Ellipse)
        );
    }

    #[test]
    fn parallel_cylinders_cut_two_lines() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(4.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(a, b, TOL) else {
            panic!("overlapping parallel cylinders cross in two lines");
        };
        assert_eq!(curves.len(), 2);
        curves.iter().for_each(|curve| {
            let IntersectionCurve::Line(line) = curve else {
                panic!("expected a line");
            };
            let ra = radial_distance(line.point(), point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0));
            let rb = radial_distance(line.point(), point(4.0, 0.0, 0.0), dir(0.0, 0.0, 1.0));
            assert!((ra - 3.0).abs() <= 1.0e-5);
            assert!((rb - 3.0).abs() <= 1.0e-5);
        });
    }

    #[test]
    fn coaxial_equal_cylinders_are_coincident() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(0.0, 0.0, 7.0), dir(0.0, 0.0, 1.0), 3.0);
        assert_eq!(
            intersect_surfaces(a, b, TOL),
            SurfaceIntersection::Coincident
        );
    }

    #[test]
    fn coaxial_unequal_cylinders_are_empty() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        assert_eq!(intersect_surfaces(a, b, TOL), SurfaceIntersection::Empty);
    }

    #[test]
    fn separated_parallel_cylinders_are_empty() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(20.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        assert_eq!(intersect_surfaces(a, b, TOL), SurfaceIntersection::Empty);
    }

    #[test]
    fn crossed_cylinders_are_unsupported_quadric() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(0.0, 0.0, 0.0), dir(1.0, 0.0, 0.0), 3.0);
        assert_eq!(
            intersect_surfaces(a, b, TOL),
            SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric)
        );
    }

    #[test]
    fn intersection_is_deterministic() {
        let a = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let b = cylinder(point(4.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        assert_eq!(intersect_surfaces(a, b, TOL), intersect_surfaces(a, b, TOL));
    }

    #[test]
    fn plane_cuts_sphere_in_a_circle() {
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        let cut = plane(point(0.0, 0.0, 3.0), dir(0.0, 0.0, 1.0));
        let circle = single_circle(&intersect_surfaces(ball, cut, TOL));
        assert!((circle.radius_mm() - 4.0).abs() <= 1.0e-6);
        let (cx, cy, cz) = circle.center().coords_mm();
        assert!(cx.abs() <= 1.0e-6 && cy.abs() <= 1.0e-6 && (cz - 3.0).abs() <= 1.0e-6);
    }

    #[test]
    fn plane_tangent_to_sphere_is_empty() {
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        let cut = plane(point(0.0, 0.0, 5.0), dir(0.0, 0.0, 1.0));
        assert_eq!(
            intersect_surfaces(ball, cut, TOL),
            SurfaceIntersection::Empty
        );
    }

    #[test]
    fn plane_clear_of_sphere_is_empty() {
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        let cut = plane(point(0.0, 0.0, 9.0), dir(0.0, 0.0, 1.0));
        assert_eq!(
            intersect_surfaces(ball, cut, TOL),
            SurfaceIntersection::Empty
        );
    }

    #[test]
    fn two_spheres_meet_in_a_circle() {
        let a = sphere(point(0.0, 0.0, 0.0), 5.0);
        let b = sphere(point(8.0, 0.0, 0.0), 5.0);
        let circle = single_circle(&intersect_surfaces(a, b, TOL));
        assert!((circle.radius_mm() - 3.0).abs() <= 1.0e-6);
        let (cx, _, _) = circle.center().coords_mm();
        assert!((cx - 4.0).abs() <= 1.0e-6);
    }

    #[test]
    fn concentric_equal_spheres_are_coincident() {
        let a = sphere(point(1.0, 2.0, 3.0), 5.0);
        let b = sphere(point(1.0, 2.0, 3.0), 5.0);
        assert_eq!(
            intersect_surfaces(a, b, TOL),
            SurfaceIntersection::Coincident
        );
    }

    #[test]
    fn separated_spheres_are_empty() {
        let a = sphere(point(0.0, 0.0, 0.0), 3.0);
        let b = sphere(point(20.0, 0.0, 0.0), 3.0);
        assert_eq!(intersect_surfaces(a, b, TOL), SurfaceIntersection::Empty);
    }

    #[test]
    fn plane_perpendicular_to_cone_axis_is_a_circle() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let cut = plane(point(0.0, 0.0, 4.0), dir(0.0, 0.0, 1.0));
        let circle = single_circle(&intersect_surfaces(tip, cut, TOL));
        assert!((circle.radius_mm() - 4.0).abs() <= 1.0e-6);
        let (_, _, cz) = circle.center().coords_mm();
        assert!((cz - 4.0).abs() <= 1.0e-6);
    }

    #[test]
    fn plane_below_cone_apex_is_empty() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let cut = plane(point(0.0, 0.0, -4.0), dir(0.0, 0.0, 1.0));
        assert_eq!(
            intersect_surfaces(tip, cut, TOL),
            SurfaceIntersection::Empty
        );
    }

    #[test]
    fn oblique_plane_cone_is_unsupported() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let cut = plane(point(0.0, 0.0, 4.0), dir(0.0, 1.0, 1.0));
        assert_eq!(
            intersect_surfaces(tip, cut, TOL),
            SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric)
        );
    }

    #[test]
    fn plane_perpendicular_to_torus_axis_cuts_two_circles() {
        let ring = torus(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0, 2.0);
        let cut = plane(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0));
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(ring, cut, TOL) else {
            panic!("the equatorial plane cuts a torus in two circles");
        };
        assert_eq!(curves.len(), 2);
        let radii: Vec<f64> = curves
            .iter()
            .map(|curve| {
                let IntersectionCurve::Circle(circle) = curve else {
                    panic!("expected circles");
                };
                circle.radius_mm()
            })
            .collect();
        assert!((radii[0] - 7.0).abs() <= 1.0e-6);
        assert!((radii[1] - 3.0).abs() <= 1.0e-6);
    }

    #[test]
    fn oblique_plane_torus_is_unsupported() {
        let ring = torus(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0, 2.0);
        let cut = plane(point(0.0, 0.0, 0.0), dir(0.0, 1.0, 1.0));
        assert_eq!(
            intersect_surfaces(ring, cut, TOL),
            SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric)
        );
    }

    #[test]
    fn coaxial_sphere_through_cylinder_cuts_two_circles() {
        let tube = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(tube, ball, TOL) else {
            panic!("a sphere wider than the bore meets the tube in two circles");
        };
        assert_eq!(curves.len(), 2);
        curves.iter().for_each(|curve| {
            let IntersectionCurve::Circle(circle) = curve else {
                panic!("expected circles");
            };
            assert!((circle.radius_mm() - 3.0).abs() <= 1.0e-6);
            let (_, _, cz) = circle.center().coords_mm();
            assert!((cz.abs() - 4.0).abs() <= 1.0e-6);
        });
    }

    #[test]
    fn sphere_inside_cylinder_bore_is_empty() {
        let tube = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0);
        let ball = sphere(point(0.0, 0.0, 0.0), 3.0);
        assert_eq!(
            intersect_surfaces(tube, ball, TOL),
            SurfaceIntersection::Empty
        );
    }

    #[test]
    fn offset_axis_cylinder_sphere_is_unsupported() {
        let tube = cylinder(point(10.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        assert_eq!(
            intersect_surfaces(tube, ball, TOL),
            SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric)
        );
    }

    #[test]
    fn coaxial_cone_meets_cylinder_in_one_circle() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let tube = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 4.0);
        let circle = single_circle(&intersect_surfaces(tip, tube, TOL));
        assert!((circle.radius_mm() - 4.0).abs() <= 1.0e-6);
        let (_, _, cz) = circle.center().coords_mm();
        assert!((cz - 4.0).abs() <= 1.0e-6);
    }

    #[test]
    fn coaxial_cones_meet_in_one_circle() {
        let up = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let down = cone(
            point(0.0, 0.0, 8.0),
            dir(0.0, 0.0, 1.0),
            -core::f64::consts::FRAC_PI_4,
        );
        let circle = single_circle(&intersect_surfaces(up, down, TOL));
        let (_, _, cz) = circle.center().coords_mm();
        assert!((cz - 4.0).abs() <= 1.0e-6);
        assert!((circle.radius_mm() - 4.0).abs() <= 1.0e-6);
    }

    #[test]
    fn identical_coaxial_cones_are_coincident() {
        let a = cone(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 0.5);
        let b = cone(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 0.5);
        assert_eq!(
            intersect_surfaces(a, b, TOL),
            SurfaceIntersection::Coincident
        );
    }

    #[test]
    fn coaxial_cone_meets_sphere_in_a_circle() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        let SurfaceIntersection::Curves(curves) = intersect_surfaces(tip, ball, TOL) else {
            panic!("an apex-centered sphere meets the cone in one upward circle");
        };
        assert_eq!(curves.len(), 1);
        let IntersectionCurve::Circle(circle) = &curves[0] else {
            panic!("expected a circle");
        };
        let expected = 5.0 / 2.0_f64.sqrt();
        assert!((circle.radius_mm() - expected).abs() <= 1.0e-6);
    }

    #[test]
    fn torus_against_quadric_is_unsupported() {
        let ring = torus(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 5.0, 2.0);
        let tube = cylinder(point(0.0, 0.0, 0.0), dir(0.0, 0.0, 1.0), 3.0);
        assert_eq!(
            intersect_surfaces(ring, tube, TOL),
            SurfaceIntersection::Unsupported(SsiGap::GeneralQuadric)
        );
    }

    #[test]
    fn quadric_intersection_is_deterministic() {
        let tip = cone(
            point(0.0, 0.0, 0.0),
            dir(0.0, 0.0, 1.0),
            core::f64::consts::FRAC_PI_4,
        );
        let ball = sphere(point(0.0, 0.0, 0.0), 5.0);
        assert_eq!(
            intersect_surfaces(tip, ball, TOL),
            intersect_surfaces(tip, ball, TOL)
        );
    }
}
