geomcore 0.32.2

//! Infinite lines in 2D and 3D: parametric evaluation and parameter
//! inversion, thin wrappers over [`crate::curve_math::analytic`].

use crate::Intersect;
use crate::curve_math::analytic;
use crate::curves::Curve2D;
use crate::curves::parametrize::{self, ParametrizeError};
use crate::frame::frame_error;
use crate::intersect::{
    ImplicitConic, LineCircle2DIntersection, LineCircle3DIntersection, LineConic3DIntersection,
    LineConicHit, LineLine2DIntersection, LineLine3DIntersection, LinePlaneIntersection,
    LineQuadricIntersection, LineTorusHit, LineTorusIntersection, PlacedConic,
};
use crate::math::{QuadraticSolution, quadratic_centre, roots_with_tangency, solve_quadratic};
use crate::projection::{self, CurveProjection};
use crate::surfaces::Surface;
use crate::{
    Axis2D, Axis3D, Circle2D, Circle3D, Cone, Cylinder, Ellipse3D, Hyperbola3D, Parabola3D, Plane,
    Point2D, Point3D, Sphere, Tolerance, Torus, Transform, Vector2D, Vector3D,
};
use std::fmt;

/// Error returned when a [`Line3D`] or [`Line2D`] cannot be constructed from
/// the given inputs.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LineConstructionError {
    /// A direction vector could not be normalized (zero length).
    NullDirection,
    /// The two points given to build the line are coincident (or too close
    /// to distinguish), so no direction can be derived from them.
    ConfusedPoints,
    /// An input coordinate, radius, or angle is NaN or infinite.
    NonFinite,
}

impl fmt::Display for LineConstructionError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let message = match self {
            LineConstructionError::NullDirection => "direction has zero length",
            LineConstructionError::ConfusedPoints => "the two points are confused",
            LineConstructionError::NonFinite => "an input is NaN or infinite",
        };
        f.write_str(message)
    }
}

impl std::error::Error for LineConstructionError {}

/// An infinite line in 3D: an origin point and a unit direction, evaluated
/// as `origin + u * direction`.
///
/// # Examples
///
/// ```
/// use geomcore::{Line3D, Point3D, Vector3D};
/// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
/// assert_eq!(line.eval_point(3.0), Point3D::new(3.0, 0.0, 0.0));
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Line3D {
    axis: Axis3D,
}

impl Line3D {
    /// Creates a new line from an origin and a direction.
    ///
    /// The direction is normalized. Returns
    /// [`LineConstructionError::NullDirection`] if `direction` cannot be
    /// normalized (zero length).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::new(2.0, 0.0, 0.0)).unwrap();
    /// assert_eq!(line.direction(), Vector3D::X);
    /// ```
    pub fn new(origin: Point3D, direction: Vector3D) -> Result<Line3D, LineConstructionError> {
        let axis = Axis3D::new(origin, direction).map_err(frame_error(
            LineConstructionError::NonFinite,
            LineConstructionError::NullDirection,
        ))?;
        Ok(Line3D { axis })
    }

    /// Creates a line from an axis directly.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Axis3D, Line3D, Point3D, Vector3D};
    /// let axis = Axis3D::new(Point3D::ORIGIN, Vector3D::Y).unwrap();
    /// let line = Line3D::from_axis(axis);
    /// assert_eq!(line.axis(), axis);
    /// ```
    pub fn from_axis(axis: Axis3D) -> Line3D {
        Line3D { axis }
    }

    /// Creates a line through two points; the direction points from `p1` to
    /// `p2`.
    ///
    /// The two points must be distinct. Returns
    /// [`LineConstructionError::ConfusedPoints`] if `p2 - p1` cannot be
    /// normalized (the points are coincident or too close to distinguish).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D};
    /// let line = Line3D::from_two_points(Point3D::ORIGIN, Point3D::new(2.0, 0.0, 0.0)).unwrap();
    /// assert_eq!(line.origin(), Point3D::ORIGIN);
    /// ```
    pub fn from_two_points(p1: Point3D, p2: Point3D) -> Result<Line3D, LineConstructionError> {
        let direction = (p2 - p1)
            .normalized()
            .ok_or(LineConstructionError::ConfusedPoints)?;
        Ok(Line3D {
            axis: Axis3D::new(p1, direction).map_err(frame_error(
                LineConstructionError::NonFinite,
                LineConstructionError::ConfusedPoints,
            ))?,
        })
    }

    /// Returns the origin point of the line.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::new(1.0, 2.0, 3.0), Vector3D::X).unwrap();
    /// assert_eq!(line.origin(), Point3D::new(1.0, 2.0, 3.0));
    /// ```
    pub fn origin(&self) -> Point3D {
        self.axis.origin()
    }

    /// Returns the unit direction of the line.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::new(0.0, 5.0, 0.0)).unwrap();
    /// assert_eq!(line.direction(), Vector3D::Y);
    /// ```
    pub fn direction(&self) -> Vector3D {
        self.axis.direction()
    }

    /// Returns the line's underlying axis.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Axis3D, Line3D, Point3D, Vector3D};
    /// let axis = Axis3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let line = Line3D::from_axis(axis);
    /// assert_eq!(line.axis(), axis);
    /// ```
    pub fn axis(&self) -> Axis3D {
        self.axis
    }

    /// Returns the image of this line under `transform`: the origin maps
    /// point-wise while the direction maps as a free vector (and is
    /// re-normalized).
    ///
    /// A uniform scale leaves the normalized direction unchanged. Fails
    /// with [`LineConstructionError::NullDirection`] only when the
    /// transform collapses the direction (a zero scale factor).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Transform, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let moved = line.transformed(Transform::translation(Vector3D::new(1.0, 2.0, 3.0))).unwrap();
    /// assert_eq!(moved.origin(), Point3D::new(1.0, 2.0, 3.0));
    /// assert_eq!(moved.direction(), Vector3D::X);
    /// ```
    pub fn transformed(&self, transform: Transform) -> Result<Line3D, LineConstructionError> {
        Line3D::new(
            transform.apply_point(self.origin()),
            transform.apply_vector(self.direction()),
        )
    }

    /// Evaluates the point on the line at parameter `u`: `origin + u * direction`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// assert_eq!(line.eval_point(3.0), Point3D::new(3.0, 0.0, 0.0));
    /// ```
    pub fn eval_point(&self, u: f64) -> Point3D {
        analytic::line_d0(&self.axis, u)
    }

    /// Evaluates the points on the line at each parameter in `us`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let points = line.eval_points(&[0.0, 1.0, 2.0]);
    /// assert_eq!(points, vec![Point3D::new(0.0, 0.0, 0.0), Point3D::new(1.0, 0.0, 0.0), Point3D::new(2.0, 0.0, 0.0)]);
    /// ```
    pub fn eval_points(&self, us: &[f64]) -> Vec<Point3D> {
        us.iter().map(|&u| self.eval_point(u)).collect()
    }

    /// Evaluates the derivative of the given `order` at parameter `u`.
    ///
    /// The line is affine in `u`, so the first derivative is the constant
    /// direction vector and every derivative of order 2 or higher is zero.
    /// `u` does not affect the result (included for API consistency with
    /// curved parametric types).
    ///
    /// # Panics
    ///
    /// Panics if `order == 0`; use [`Line3D::eval_point`] to evaluate the
    /// position itself.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::new(0.0, 3.0, 0.0)).unwrap();
    /// assert_eq!(line.eval_derivative(1.5, 1), Vector3D::Y);
    /// assert_eq!(line.eval_derivative(1.5, 2), Vector3D::ZERO);
    /// ```
    pub fn eval_derivative(&self, u: f64, order: u32) -> Vector3D {
        let _ = u;
        match order {
            0 => panic!("eval_derivative: order must be >= 1 (use eval_point for order 0)"),
            1 => analytic::line_d1(&self.axis),
            _ => Vector3D::ZERO,
        }
    }

    /// Recovers the parameter `u` of a point on (or near) the line: the
    /// signed projection `(point - origin) . direction`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// assert_eq!(line.parameter_of(Point3D::new(2.5, 0.0, 0.0)), 2.5);
    /// ```
    pub fn parameter_of(&self, point: Point3D) -> f64 {
        analytic::line_parameter(&self.axis, point)
    }

    /// Returns whether `point` lies on the line: the inverse parameter is
    /// recovered with [`Line3D::parameter_of`] and re-evaluated, and the
    /// point counts as contained when the re-evaluated point is within
    /// `tol.confusion` of it.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// assert!(line.contains(Point3D::new(2.5, 0.0, 0.0), tol));
    /// assert!(!line.contains(Point3D::new(2.5, 1.0, 0.0), tol));
    /// ```
    pub fn contains(&self, point: Point3D, tol: Tolerance) -> bool {
        let u = self.parameter_of(point);
        self.eval_point(u).distance(point) <= tol.confusion
    }

    /// Projects `point` onto the line, returning the parameter of the
    /// closest point and its distance. Distances within `tol.confusion`
    /// snap to `0.0`, matching [`Line3D::contains`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line3D, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let proj = line.project_point(Point3D::new(2.5, 1.0, 0.0), Tolerance::DEFAULT);
    /// assert_eq!(proj.parameter, 2.5);
    /// assert_eq!(proj.distance, 1.0);
    /// ```
    pub fn project_point(&self, point: Point3D, tol: Tolerance) -> CurveProjection {
        let u = self.parameter_of(point);
        let distance = self.eval_point(u).distance(point);
        CurveProjection {
            parameter: u,
            distance: projection::snap_distance(distance, tol.confusion),
        }
    }

    /// Projects each point in `points` onto the line.
    ///
    /// Default-style batch wrapper over [`Line3D::project_point`]: one
    /// native call per batch, mirroring [`Line3D::eval_points`].
    pub fn project_points(&self, points: &[Point3D], tol: Tolerance) -> Vec<CurveProjection> {
        points.iter().map(|&p| self.project_point(p, tol)).collect()
    }

    /// Intersects with a [`Plane`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_plane(&self, plane: &Plane, tol: Tolerance) -> LinePlaneIntersection {
        self.intersect(plane, tol)
    }

    /// Intersects with a [`Sphere`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_sphere(&self, sphere: &Sphere, tol: Tolerance) -> LineQuadricIntersection {
        self.intersect(sphere, tol)
    }

    /// Intersects with a [`Cylinder`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_cylinder(
        &self,
        cylinder: &Cylinder,
        tol: Tolerance,
    ) -> LineQuadricIntersection {
        self.intersect(cylinder, tol)
    }

    /// Intersects with a [`Cone`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_cone(&self, cone: &Cone, tol: Tolerance) -> LineQuadricIntersection {
        self.intersect(cone, tol)
    }

    /// Intersects with a [`Torus`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_torus(&self, torus: &Torus, tol: Tolerance) -> LineTorusIntersection {
        self.intersect(torus, tol)
    }

    /// Computes the exact 2D representation of this line in a surface's
    /// parameter space: a [`Curve2D`] `q(t)` such that
    /// `surface.eval_point(q(t)) == self.eval_point(t)` for the same `t`.
    ///
    /// Only a few line/surface pairs admit a closed-form 2D image: a line on
    /// a plane (projects to a 2D line), a line parallel to a cylinder or cone
    /// axis / a cone generator (a vertical iso-`u` line in `(u, v)`). Every
    /// other pair — including a line on a sphere or torus — returns
    /// [`ParametrizeError::NotAnalytic`].
    ///
    /// The projection math assumes the line lies on the surface. As a
    /// geomcore safeguard, the candidate 2D image is verified against the
    /// surface at a few parameters; if `surface.eval_point(q(t))` is farther
    /// than `tol.confusion` from `self.eval_point(t)`,
    /// [`ParametrizeError::CurveNotOnSurface`] is returned. For periodic surfaces the result is normalized so `q(0)`
    /// lies in the canonical parameter window.
    ///
    /// # Examples
    ///
    /// A vertical line on a cylinder maps to a vertical line in `(u, v)`:
    ///
    /// ```
    /// use geomcore::curves::{Curve2D, ParametricCurve2D};
    /// use geomcore::{Cylinder, Line3D, Point3D, Tolerance, Vector3D};
    ///
    /// let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
    /// let line = Line3D::new(Point3D::new(2.0, 0.0, 0.0), Vector3D::Z).unwrap();
    /// let pcurve = line.parametrize_on(&cylinder, Tolerance::DEFAULT).unwrap();
    /// // q(0) sits at u = 0 (angle of x-axis), v = 0 (height of the origin).
    /// let q0 = pcurve.eval_point(0.0);
    /// assert!(q0.x.abs() < 1e-9 && q0.y.abs() < 1e-9);
    /// assert!(matches!(pcurve, Curve2D::Line(_)));
    /// ```
    pub fn parametrize_on(
        &self,
        surface: impl Into<Surface>,
        tol: Tolerance,
    ) -> Result<Curve2D, ParametrizeError> {
        parametrize::line_on_surface(self, &surface.into(), tol)
    }

    /// Intersects with a [`Line3D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_line(&self, other: &Line3D, tol: Tolerance) -> LineLine3DIntersection {
        self.intersect(other, tol)
    }

    /// Intersects with a [`Circle3D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_circle(&self, circle: &Circle3D, tol: Tolerance) -> LineCircle3DIntersection {
        self.intersect(circle, tol)
    }

    /// How this line meets the plane of a figure reaching `extent` from
    /// `centre`: in the plane when turning the line into it moves the line
    /// by at most `tol.confusion` across the figure, the offset measured
    /// at the line point closest to `centre` (ADR 0006).
    fn meet_plane_near(
        &self,
        plane: &Plane,
        centre: Point3D,
        extent: f64,
        tol: Tolerance,
    ) -> LinePlaneIntersection {
        let n = plane.normal();
        let b = n.dot(self.direction());
        let foot = self.parameter_of(centre);
        let c = n.dot(self.eval_point(foot) - plane.frame().origin());
        if b.abs() <= tol.angular_for_extent(extent) {
            if c.abs() <= tol.confusion {
                LinePlaneIntersection::Coincident
            } else {
                LinePlaneIntersection::Parallel
            }
        } else {
            let t = foot - c / b;
            LinePlaneIntersection::Point(t, self.eval_point(t))
        }
    }

    /// Shared solve for [`Line3D::intersect_ellipse`],
    /// [`Line3D::intersect_parabola`], and [`Line3D::intersect_hyperbola`]:
    /// a transversal line meets the conic's plane once (graze or miss); a
    /// coplanar line substitutes into the conic's implicit equation and the
    /// quadratic roots map back with this line's evaluation. Candidates
    /// off the conic (the unmodeled hyperbola branch) are dropped, so a
    /// lone hit reports [`LineConic3DIntersection::Tangent`] even when it
    /// is transversal.
    fn intersect_frame_conic(
        &self,
        conic: &PlacedConic,
        contains: impl Fn(Point3D) -> bool,
        conic_parameter_of: impl Fn(Point3D) -> f64,
        tol: Tolerance,
    ) -> LineConic3DIntersection {
        match self.meet_plane_near(&conic.plane, conic.origin, conic.extent, tol) {
            LinePlaneIntersection::Point(t, p) => {
                if contains(p) {
                    LineConic3DIntersection::Tangent(LineConicHit {
                        line_param: t,
                        point: p,
                        conic_param: conic_parameter_of(p),
                    })
                } else {
                    LineConic3DIntersection::Empty
                }
            }
            LinePlaneIntersection::Parallel => LineConic3DIntersection::Empty,
            LinePlaneIntersection::Coincident => {
                let d = self.direction();
                // Coefficients about the line point at parameter t0,
                // centred at the midpoint of the two roots.
                let coeffs = |t0: f64| {
                    let w = self.eval_point(t0) - conic.origin;
                    conic.implicit.line_substitute(
                        w.dot(conic.x),
                        w.dot(conic.y),
                        d.dot(conic.x),
                        d.dot(conic.y),
                    )
                };
                let (qa, qb, _) = coeffs(0.0);
                let t0 = quadratic_centre(qa, qb, tol);
                let (qa, qb, qc) = coeffs(t0);
                let mut hits = Vec::new();
                let mut push = |t: f64| {
                    let p = self.eval_point(t);
                    if contains(p) {
                        hits.push(LineConicHit {
                            line_param: t,
                            point: p,
                            conic_param: conic_parameter_of(p),
                        });
                    }
                };
                let sol = solve_quadratic(qa, qb, qc, tol)
                    .touching_at(-qb / (2.0 * qa), |s| contains(self.eval_point(t0 + s)));
                match sol.shifted(t0) {
                    QuadraticSolution::Two(t1, t2) => {
                        push(t1);
                        push(t2);
                    }
                    QuadraticSolution::One(t) | QuadraticSolution::Linear(t) => push(t),
                    QuadraticSolution::Degenerate | QuadraticSolution::Empty => {}
                }
                match hits.len() {
                    2 => LineConic3DIntersection::Points(hits[0], hits[1]),
                    1 => LineConic3DIntersection::Tangent(hits[0]),
                    _ => LineConic3DIntersection::Empty,
                }
            }
        }
    }

    /// Intersects with an [`Ellipse3D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_ellipse(
        &self,
        ellipse: &Ellipse3D,
        tol: Tolerance,
    ) -> LineConic3DIntersection {
        self.intersect(ellipse, tol)
    }

    /// Intersects with a [`Parabola3D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_parabola(
        &self,
        parabola: &Parabola3D,
        tol: Tolerance,
    ) -> LineConic3DIntersection {
        self.intersect(parabola, tol)
    }

    /// Intersects with a [`Hyperbola3D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_hyperbola(
        &self,
        hyperbola: &Hyperbola3D,
        tol: Tolerance,
    ) -> LineConic3DIntersection {
        self.intersect(hyperbola, tol)
    }
}

/// An infinite line in 2D: an origin point and a unit direction, evaluated
/// as `origin + u * direction`.
///
/// # Examples
///
/// ```
/// use geomcore::{Line2D, Point2D, Vector2D};
/// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
/// assert_eq!(line.eval_point(3.0), Point2D::new(3.0, 0.0));
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Line2D {
    axis: Axis2D,
}

impl Line2D {
    /// Creates a new line from an origin and a direction.
    ///
    /// The direction is normalized. Returns
    /// [`LineConstructionError::NullDirection`] if `direction` cannot be
    /// normalized (zero length).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::new(0.0, 5.0)).unwrap();
    /// assert_eq!(line.direction(), Vector2D::Y);
    /// ```
    pub fn new(origin: Point2D, direction: Vector2D) -> Result<Line2D, LineConstructionError> {
        let axis = Axis2D::new(origin, direction).map_err(frame_error(
            LineConstructionError::NonFinite,
            LineConstructionError::NullDirection,
        ))?;
        Ok(Line2D { axis })
    }

    /// Creates a line from an axis directly.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Axis2D, Line2D, Point2D, Vector2D};
    /// let axis = Axis2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let line = Line2D::from_axis(axis);
    /// assert_eq!(line.axis(), axis);
    /// ```
    pub fn from_axis(axis: Axis2D) -> Line2D {
        Line2D { axis }
    }

    /// Creates a line through two points; the direction points from `p1` to
    /// `p2`.
    ///
    /// The two points must be distinct. Returns
    /// [`LineConstructionError::ConfusedPoints`] if `p2 - p1` cannot be
    /// normalized (the points are coincident or too close to distinguish).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D};
    /// let line = Line2D::from_two_points(Point2D::ORIGIN, Point2D::new(0.0, 4.0)).unwrap();
    /// assert_eq!(line.origin(), Point2D::ORIGIN);
    /// ```
    pub fn from_two_points(p1: Point2D, p2: Point2D) -> Result<Line2D, LineConstructionError> {
        let direction = (p2 - p1)
            .normalized()
            .ok_or(LineConstructionError::ConfusedPoints)?;
        Ok(Line2D {
            axis: Axis2D::new(p1, direction).map_err(frame_error(
                LineConstructionError::NonFinite,
                LineConstructionError::ConfusedPoints,
            ))?,
        })
    }

    /// Returns the origin point of the line.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::new(1.0, 2.0), Vector2D::X).unwrap();
    /// assert_eq!(line.origin(), Point2D::new(1.0, 2.0));
    /// ```
    pub fn origin(&self) -> Point2D {
        self.axis.origin()
    }

    /// Returns the unit direction of the line.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::new(0.0, 5.0)).unwrap();
    /// assert_eq!(line.direction(), Vector2D::Y);
    /// ```
    pub fn direction(&self) -> Vector2D {
        self.axis.direction()
    }

    /// Returns the line's underlying axis.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Axis2D, Line2D, Point2D, Vector2D};
    /// let axis = Axis2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let line = Line2D::from_axis(axis);
    /// assert_eq!(line.axis(), axis);
    /// ```
    pub fn axis(&self) -> Axis2D {
        self.axis
    }

    /// Evaluates the point on the line at parameter `u`: `origin + u * direction`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// assert_eq!(line.eval_point(3.0), Point2D::new(3.0, 0.0));
    /// ```
    pub fn eval_point(&self, u: f64) -> Point2D {
        analytic::line2d_d0(&self.axis, u)
    }

    /// Evaluates the points on the line at each parameter in `us`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let points = line.eval_points(&[0.0, 1.0, 2.0]);
    /// assert_eq!(points, vec![Point2D::new(0.0, 0.0), Point2D::new(1.0, 0.0), Point2D::new(2.0, 0.0)]);
    /// ```
    pub fn eval_points(&self, us: &[f64]) -> Vec<Point2D> {
        us.iter().map(|&u| self.eval_point(u)).collect()
    }

    /// Evaluates the derivative of the given `order` at parameter `u`.
    ///
    /// The line is affine in `u`, so the first derivative is the constant
    /// direction vector and every derivative of order 2 or higher is zero.
    /// `u` does not affect the result (included for API consistency with
    /// curved parametric types).
    ///
    /// # Panics
    ///
    /// Panics if `order == 0`; use [`Line2D::eval_point`] to evaluate the
    /// position itself.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::new(0.0, 3.0)).unwrap();
    /// assert_eq!(line.eval_derivative(1.5, 1), Vector2D::Y);
    /// assert_eq!(line.eval_derivative(1.5, 2), Vector2D::ZERO);
    /// ```
    pub fn eval_derivative(&self, u: f64, order: u32) -> Vector2D {
        let _ = u;
        match order {
            0 => panic!("eval_derivative: order must be >= 1 (use eval_point for order 0)"),
            1 => self.axis.direction(),
            _ => Vector2D::ZERO,
        }
    }

    /// Recovers the parameter `u` of a point on (or near) the line: the
    /// signed projection `(point - origin) . direction`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// assert_eq!(line.parameter_of(Point2D::new(2.5, 0.0)), 2.5);
    /// ```
    pub fn parameter_of(&self, point: Point2D) -> f64 {
        analytic::line2d_parameter(&self.axis, point)
    }

    /// Returns whether `point` lies on the line: the inverse parameter is
    /// recovered with [`Line2D::parameter_of`] and re-evaluated, and the
    /// point counts as contained when the re-evaluated point is within
    /// `tol.confusion` of it.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Tolerance, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// assert!(line.contains(Point2D::new(2.5, 0.0), tol));
    /// assert!(!line.contains(Point2D::new(2.5, 1.0), tol));
    /// ```
    pub fn contains(&self, point: Point2D, tol: Tolerance) -> bool {
        let u = self.parameter_of(point);
        self.eval_point(u).distance(point) <= tol.confusion
    }

    /// Projects `point` onto the line, returning the parameter of the
    /// closest point and its distance. Distances within `tol.confusion`
    /// snap to `0.0`, matching [`Line2D::contains`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Line2D, Point2D, Tolerance, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let proj = line.project_point(Point2D::new(2.5, 1.0), Tolerance::DEFAULT);
    /// assert_eq!(proj.parameter, 2.5);
    /// assert_eq!(proj.distance, 1.0);
    /// ```
    pub fn project_point(&self, point: Point2D, tol: Tolerance) -> CurveProjection {
        let u = self.parameter_of(point);
        let distance = self.eval_point(u).distance(point);
        CurveProjection {
            parameter: u,
            distance: projection::snap_distance(distance, tol.confusion),
        }
    }

    /// Projects each point in `points` onto the line.
    ///
    /// Default-style batch wrapper over [`Line2D::project_point`]: one
    /// native call per batch, mirroring [`Line2D::eval_points`].
    pub fn project_points(&self, points: &[Point2D], tol: Tolerance) -> Vec<CurveProjection> {
        points.iter().map(|&p| self.project_point(p, tol)).collect()
    }

    /// Intersects with a [`Line2D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_line(&self, other: &Line2D, tol: Tolerance) -> LineLine2DIntersection {
        self.intersect(other, tol)
    }

    /// Intersects with a [`Circle2D`]; see its [`Intersect`] implementation.
    #[deprecated(since = "0.29.0", note = "use `Intersect::intersect`")]
    pub fn intersect_circle(&self, circle: &Circle2D, tol: Tolerance) -> LineCircle2DIntersection {
        self.intersect(circle, tol)
    }
}

impl Intersect<Plane> for Line3D {
    type Output = LinePlaneIntersection;

    /// Intersects this line with a plane.
    ///
    /// A transversal line meets the plane once; a line parallel to the
    /// plane (direction against the normal within `tol.angular`) is either
    /// offset ([`LinePlaneIntersection::Parallel`]) or contained
    /// ([`LinePlaneIntersection::Coincident`], decided by `tol.confusion`).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line3D, LinePlaneIntersection, Plane, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
    /// let plane = Plane::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&plane, tol) {
    ///     LinePlaneIntersection::Point(t, p) => {
    ///         assert_eq!(t, 0.0);
    ///         assert_eq!(p, Point3D::ORIGIN);
    ///     }
    ///     _ => panic!("expected a point"),
    /// }
    /// ```
    fn intersect(&self, plane: &Plane, tol: Tolerance) -> LinePlaneIntersection {
        let n = plane.normal();
        let b = n.dot(self.direction());
        let c = n.dot(self.origin() - plane.frame().origin());
        if b.abs() <= tol.angular {
            if c.abs() <= tol.confusion {
                LinePlaneIntersection::Coincident
            } else {
                LinePlaneIntersection::Parallel
            }
        } else {
            let t = -c / b;
            LinePlaneIntersection::Point(t, self.eval_point(t))
        }
    }
}

impl Intersect<Sphere> for Line3D {
    type Output = LineQuadricIntersection;

    /// Intersects this line with a sphere.
    ///
    /// Substituting the unit-speed parametrization into `|X - C|^2 = R^2`
    /// gives a quadratic (leading coefficient exactly 1), classified with
    /// `solve_quadratic`: two ordered hits, a grazing tangent, or empty.
    /// The quadratic is formed about the line's closest approach to the
    /// centre; the line is tangent when that closest point lies within
    /// `tol.confusion` of the sphere (`|dist(centre, line) - R| <= tol`),
    /// and the tangent hit reports there.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line3D, LineQuadricIntersection, Point3D, Sphere, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let sphere = Sphere::new(Point3D::ORIGIN, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&sphere, tol) {
    ///     LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
    ///         assert_eq!((t1, t2), (-2.0, 2.0));
    ///         assert_eq!(p1, Point3D::new(-2.0, 0.0, 0.0));
    ///         assert_eq!(p2, Point3D::new(2.0, 0.0, 0.0));
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, sphere: &Sphere, tol: Tolerance) -> LineQuadricIntersection {
        let d = self.direction();
        let t0 = -d.dot(self.origin() - sphere.center());
        let w = self.eval_point(t0) - sphere.center();
        // A is exactly 1: the direction is unit by construction.
        let qb = 2.0 * d.dot(w);
        let sol = solve_quadratic(1.0, qb, w.dot(w) - sphere.radius().powi(2), tol)
            .touching_at(-0.5 * qb, |s| sphere.contains(self.eval_point(t0 + s), tol));
        match sol.shifted(t0) {
            QuadraticSolution::Two(t1, t2) => LineQuadricIntersection::TwoPoints(
                (t1, self.eval_point(t1)),
                (t2, self.eval_point(t2)),
            ),
            QuadraticSolution::One(t) => LineQuadricIntersection::Tangent(t, self.eval_point(t)),
            QuadraticSolution::Empty => LineQuadricIntersection::Empty,
            QuadraticSolution::Linear(_) | QuadraticSolution::Degenerate => {
                unreachable!("unit direction gives a leading coefficient of 1")
            }
        }
    }
}

impl Intersect<Cylinder> for Line3D {
    type Output = LineQuadricIntersection;

    /// Intersects this line with a cylinder.
    ///
    /// The implicit equation `|X - C|^2 - ((X - C).a)^2 = r^2` along the
    /// line is quadratic; a direction parallel to the axis degrades to the
    /// linear case (one transversal hit) or, when fully degenerate, to a
    /// generator coincidence check.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Cylinder, Line3D, LineQuadricIntersection, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&cylinder, tol) {
    ///     LineQuadricIntersection::TwoPoints((t1, _), (t2, _)) => {
    ///         assert_eq!((t1, t2), (-2.0, 2.0));
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, cylinder: &Cylinder, tol: Tolerance) -> LineQuadricIntersection {
        let d = self.direction();
        let a = cylinder.axis().direction();
        let da = d.dot(a);
        let qa = 1.0 - da * da;
        // `qa` is the squared sine of the line's angle to the axis, so
        // the line is parallel to it below the squared angular tolerance.
        let tol = Tolerance {
            angular: tol.angular * tol.angular,
            ..tol
        };
        // Coefficients about the line point at parameter t0.
        let coeffs = |t0: f64| {
            let w = self.eval_point(t0) - cylinder.axis().origin();
            let wa = w.dot(a);
            (
                2.0 * (d.dot(w) - da * wa),
                w.dot(w) - wa * wa - cylinder.radius().powi(2),
            )
        };
        // Centred at the closest approach to the axis; tangent when that
        // point lies within tolerance of the cylinder.
        let t0 = quadratic_centre(qa, coeffs(0.0).0, tol);
        let (qb, qc) = coeffs(t0);
        let sol = solve_quadratic(qa, qb, qc, tol).touching_at(-qb / (2.0 * qa), |s| {
            cylinder.contains(self.eval_point(t0 + s), tol)
        });
        match sol.shifted(t0) {
            QuadraticSolution::Two(t1, t2) => LineQuadricIntersection::TwoPoints(
                (t1, self.eval_point(t1)),
                (t2, self.eval_point(t2)),
            ),
            QuadraticSolution::One(t) => LineQuadricIntersection::Tangent(t, self.eval_point(t)),
            QuadraticSolution::Linear(t) => {
                LineQuadricIntersection::OnePoint(t, self.eval_point(t))
            }
            QuadraticSolution::Empty => LineQuadricIntersection::Empty,
            QuadraticSolution::Degenerate => {
                if cylinder.contains(self.origin(), tol) {
                    LineQuadricIntersection::Coincident
                } else {
                    LineQuadricIntersection::Empty
                }
            }
        }
    }
}

impl Intersect<Cone> for Line3D {
    type Output = LineQuadricIntersection;

    /// Intersects this line with a cone.
    ///
    /// The half-angle equation `((X - A).a)^2 = |X - A|^2*cos^2(phi)`
    /// along the line is quadratic, and holds on both nappes of the
    /// (double) cone: two hits, a tangent, or empty. A direction parallel
    /// to a generator degrades to the linear case (one transversal hit);
    /// full degeneracy resolves to a generator coincidence check.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Cone, Frame3D, Line3D, LineQuadricIntersection, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::new(0.0, 0.0, 5.0), Vector3D::X).unwrap();
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&cone, tol) {
    ///     LineQuadricIntersection::TwoPoints((t1, _), (t2, _)) => {
    ///         assert!(t1 < t2);
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, cone: &Cone, tol: Tolerance) -> LineQuadricIntersection {
        let d = self.direction();
        let a = cone.frame().z_direction();
        let cos_phi = cone.semi_angle().cos();
        let da = d.dot(a);
        let qa = da * da - cos_phi * cos_phi;
        // Coefficients about the line point at parameter t0.
        let coeffs = |t0: f64| {
            let m = self.eval_point(t0) - cone.apex();
            let ma = m.dot(a);
            (
                2.0 * (da * ma - cos_phi * cos_phi * d.dot(m)),
                ma * ma - cos_phi * cos_phi * m.dot(m),
            )
        };
        // Centred at the midpoint of the two roots.
        let t0 = quadratic_centre(qa, coeffs(0.0).0, tol);
        let (qb, qc) = coeffs(t0);
        let hit = |t: f64| (t, self.eval_point(t));
        // Tangent when the vertex of the quadratic lies within tolerance
        // of the cone.
        let sol = solve_quadratic(qa, qb, qc, tol).touching_at(-qb / (2.0 * qa), |s| {
            cone.contains(self.eval_point(t0 + s), tol)
        });
        match sol.shifted(t0) {
            QuadraticSolution::Two(t1, t2) => LineQuadricIntersection::TwoPoints(hit(t1), hit(t2)),
            QuadraticSolution::One(t) => LineQuadricIntersection::Tangent(t, self.eval_point(t)),
            QuadraticSolution::Linear(t) => {
                LineQuadricIntersection::OnePoint(t, self.eval_point(t))
            }
            QuadraticSolution::Empty => LineQuadricIntersection::Empty,
            QuadraticSolution::Degenerate => {
                if cone.contains(self.origin(), tol) {
                    LineQuadricIntersection::Coincident
                } else {
                    LineQuadricIntersection::Empty
                }
            }
        }
    }
}

impl Intersect<Torus> for Line3D {
    type Output = LineTorusIntersection;

    /// Intersects this line with a torus.
    ///
    /// The torus implicit equation
    /// `(|X - C|^2 + R^2 - r^2)^2 = 4*R^2*(|X - C|^2 - ((X - C).a)^2)`
    /// along the unit-speed line `X = O + t*d` is a quartic in `t` with
    /// leading coefficient exactly 1: with `w = O - C`, `b1 = 2*(w.d)`,
    /// `c1 = w.w + R^2 - r^2`, and the radial part
    /// `qa*t^2 + qb*t + qc`, the ascending coefficients are
    /// `[c1^2 - 4*R^2*qc, 2*b1*c1 - 4*R^2*qb,
    /// b1^2 + 2*c1 - 4*R^2*qa, 2*b1, 1]`, solved by `real_roots` with
    /// clustering inside `tol.confusion * (1 + |t|)` (the same band rule
    /// as the quadratic solver). The quartic is formed about the line's
    /// closest approach to the centre (`O` there, `t` measured from it),
    /// so neither the band nor the coefficients' accuracy depend on
    /// where the line's origin lies. A critical point of the quartic
    /// whose line point lies within `tol.confusion` of the torus is a
    /// grazing contact: the roots around it report as one tangent hit
    /// there (multiplicity 2 or more). Candidates verify with
    /// [`Torus::contains`], so spurious polynomial roots are dropped, and
    /// each hit also reports the torus `(u, v)` parameters.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line3D, LineTorusIntersection, Point3D, Tolerance, Torus, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let torus = Torus::new(Point3D::ORIGIN, Vector3D::Z, 4.0, 1.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&torus, tol) {
    ///     LineTorusIntersection::Hits(hits) => {
    ///         assert_eq!(hits.len(), 4);
    ///         assert_eq!(
    ///             [hits[0].line_param, hits[1].line_param, hits[2].line_param, hits[3].line_param],
    ///             [-5.0, -3.0, 3.0, 5.0]
    ///         );
    ///     }
    ///     _ => panic!("expected four hits"),
    /// }
    /// ```
    fn intersect(&self, torus: &Torus, tol: Tolerance) -> LineTorusIntersection {
        let a = torus.frame().z_direction();
        let (major, minor) = (torus.major_radius(), torus.minor_radius());
        let d = self.direction();
        let t0 = -d.dot(self.origin() - torus.center());
        let w = self.eval_point(t0) - torus.center();
        let wd = w.dot(d);
        let (da, wa) = (d.dot(a), w.dot(a));
        let b1 = 2.0 * wd;
        let c1 = w.dot(w) + major * major - minor * minor;
        let qa = 1.0 - da * da;
        let qb = 2.0 * (wd - wa * da);
        let qc = w.dot(w) - wa * wa;
        let rr = 4.0 * major * major;
        let coeffs = [
            c1 * c1 - rr * qc,
            2.0 * b1 * c1 - rr * qb,
            b1 * b1 + 2.0 * c1 - rr * qa,
            2.0 * b1,
            1.0,
        ];
        let touches = |s: f64| torus.contains(self.eval_point(t0 + s), tol);
        let mut hits: Vec<LineTorusHit> = roots_with_tangency(&coeffs, tol, touches)
            .iter()
            .filter_map(|r| {
                let t = t0 + r.value;
                let p = self.eval_point(t);
                if torus.contains(p, tol) {
                    Some(LineTorusHit {
                        line_param: t,
                        point: p,
                        surface_params: torus.parameters_of(p),
                        multiplicity: r.multiplicity,
                    })
                } else {
                    None
                }
            })
            .collect();
        hits.sort_by(|x, y| x.line_param.partial_cmp(&y.line_param).unwrap());
        if hits.is_empty() {
            LineTorusIntersection::Empty
        } else {
            LineTorusIntersection::Hits(hits)
        }
    }
}

impl Intersect<Line3D> for Line3D {
    type Output = LineLine3DIntersection;

    /// Intersects this line with another 3D line.
    ///
    /// Skew lines report both parameters and their distance (the classic
    /// `sc = (b*e - d)/(1 - b^2)` solution); parallel lines report their
    /// distance or coincidence. Transversal hits carry both parameters.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line3D, LineLine3DIntersection, Point3D, Tolerance, Vector3D};
    /// let l1 = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let l2 = Line3D::new(Point3D::ORIGIN, Vector3D::Y).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match l1.intersect(&l2, tol) {
    ///     LineLine3DIntersection::Point(s, p, t) => {
    ///         assert_eq!((s, t), (0.0, 0.0));
    ///         assert_eq!(p, Point3D::ORIGIN);
    ///     }
    ///     _ => panic!("expected a point"),
    /// }
    /// ```
    fn intersect(&self, other: &Line3D, tol: Tolerance) -> LineLine3DIntersection {
        let (d1, d2) = (self.direction(), other.direction());
        let w0 = self.origin() - other.origin();
        let b = d1.dot(d2);
        // The squared sine of the angle between the unit directions:
        // `1 - b²` cancels to zero below about 1e-8 rad.
        let denom = d1.cross(d2).square_magnitude();
        if denom <= tol.angular * tol.angular {
            // Parallel: distance from other's origin to self, or coincident.
            if other.contains(self.origin(), tol) {
                LineLine3DIntersection::Coincident
            } else {
                LineLine3DIntersection::Parallel(self.project_point(other.origin(), tol).distance)
            }
        } else {
            let d = d1.dot(w0);
            let e = d2.dot(w0);
            let s = (b * e - d) / denom;
            let t = (e - b * d) / denom;
            let (p1, p2) = (self.eval_point(s), other.eval_point(t));
            let distance = p1.distance(p2);
            if distance <= tol.confusion {
                LineLine3DIntersection::Point(s, p1, t)
            } else {
                LineLine3DIntersection::Skew { s, t, distance }
            }
        }
    }
}

impl Intersect<Circle3D> for Line3D {
    type Output = LineCircle3DIntersection;

    /// Intersects this line with a 3D circle.
    ///
    /// A transversal line meets the circle's plane once: on-circle hits
    /// graze, the rest miss. A line lying in the plane (turning it into
    /// the plane moves it by at most `tol.confusion` across the circle)
    /// reduces to the 2D problem in the circle's frame, mapping hits back
    /// with the plane evaluation and this line's inverse parameter.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Circle3D, Line3D, LineCircle3DIntersection, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let circle = Circle3D::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&circle, tol) {
    ///     LineCircle3DIntersection::Points((t1, _), (t2, _)) => {
    ///         assert_eq!((t1, t2), (-2.0, 2.0));
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, circle: &Circle3D, tol: Tolerance) -> LineCircle3DIntersection {
        let plane = Plane::from_frame(circle.frame());
        match self.meet_plane_near(&plane, circle.center(), circle.radius(), tol) {
            LinePlaneIntersection::Point(t, p) => {
                if circle.contains(p, tol) {
                    LineCircle3DIntersection::Tangent(t, p)
                } else {
                    LineCircle3DIntersection::Empty
                }
            }
            LinePlaneIntersection::Parallel => LineCircle3DIntersection::Empty,
            LinePlaneIntersection::Coincident => {
                // The line in the circle's frame coordinates, through its
                // point closest to the centre.
                let (x, y) = (plane.frame().x_direction(), plane.frame().y_direction());
                let (d, w) = (
                    self.direction(),
                    self.eval_point(self.parameter_of(circle.center())) - circle.center(),
                );
                let Ok(l2) = Line2D::new(
                    Point2D::new(w.dot(x), w.dot(y)),
                    Vector2D::new(d.dot(x), d.dot(y)),
                ) else {
                    return LineCircle3DIntersection::Empty;
                };
                let c2 = Circle2D::new(Point2D::ORIGIN, circle.radius())
                    .expect("radius is non-negative by construction");
                match l2.intersect(&c2, tol) {
                    LineCircle2DIntersection::Points((_, p1), (_, p2)) => {
                        let (q1, q2) = (plane.eval_point(p1.x, p1.y), plane.eval_point(p2.x, p2.y));
                        LineCircle3DIntersection::Points(
                            (self.parameter_of(q1), q1),
                            (self.parameter_of(q2), q2),
                        )
                    }
                    LineCircle2DIntersection::Tangent(_, p) => {
                        let q = plane.eval_point(p.x, p.y);
                        LineCircle3DIntersection::Tangent(self.parameter_of(q), q)
                    }
                    LineCircle2DIntersection::Empty => LineCircle3DIntersection::Empty,
                }
            }
        }
    }
}

impl Intersect<Ellipse3D> for Line3D {
    type Output = LineConic3DIntersection;

    /// Intersects this line with a 3D ellipse.
    ///
    /// A transversal line meets the ellipse's plane once: on-ellipse hits
    /// graze, the rest miss. A line lying in the plane substitutes into
    /// the ellipse equation (quadratic in the line parameter).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Ellipse3D, Line3D, LineConic3DIntersection, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let ellipse = Ellipse3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 3.0, 1.5).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&ellipse, tol) {
    ///     LineConic3DIntersection::Points(h1, h2) => {
    ///         assert_eq!((h1.line_param, h2.line_param), (-3.0, 3.0));
    ///         assert_eq!((h1.point, h2.point), (Point3D::new(-3.0, 0.0, 0.0), Point3D::new(3.0, 0.0, 0.0)));
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, ellipse: &Ellipse3D, tol: Tolerance) -> LineConic3DIntersection {
        let frame = ellipse.frame();
        let conic = PlacedConic {
            plane: Plane::from_frame(frame),
            x: frame.x_direction(),
            y: frame.y_direction(),
            origin: ellipse.center(),
            implicit: ImplicitConic::ellipse(ellipse.major_radius(), ellipse.minor_radius()),
            extent: ellipse.major_radius().max(ellipse.minor_radius()),
        };
        self.intersect_frame_conic(
            &conic,
            |p| ellipse.contains(p, tol),
            |p| ellipse.parameter_in_bounds(ellipse.parameter_of(p), tol),
            tol,
        )
    }
}

impl Intersect<Parabola3D> for Line3D {
    type Output = LineConic3DIntersection;

    /// Intersects this line with a 3D parabola.
    ///
    /// Same calling shape as [`Line3D::intersect_ellipse`]: a transversal
    /// line grazes or misses, a coplanar line solves the substituted
    /// quadratic (asymptote-style parallels report the linear root).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line3D, LineConic3DIntersection, Parabola3D, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
    /// let parabola = Parabola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 1.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&parabola, tol) {
    ///     LineConic3DIntersection::Tangent(hit) => {
    ///         assert_eq!(hit.line_param, 0.0);
    ///         assert_eq!(hit.point, Point3D::ORIGIN);
    ///     }
    ///     _ => panic!("expected a tangent"),
    /// }
    /// ```
    fn intersect(&self, parabola: &Parabola3D, tol: Tolerance) -> LineConic3DIntersection {
        let frame = parabola.frame();
        let conic = PlacedConic {
            plane: Plane::from_frame(frame),
            x: frame.x_direction(),
            y: frame.y_direction(),
            origin: parabola.apex(),
            implicit: ImplicitConic::parabola(parabola.focal()),
            extent: f64::INFINITY,
        };
        self.intersect_frame_conic(
            &conic,
            |p| parabola.contains(p, tol),
            |p| parabola.parameter_of(p),
            tol,
        )
    }
}

impl Intersect<Hyperbola3D> for Line3D {
    type Output = LineConic3DIntersection;

    /// Intersects this line with a 3D hyperbola.
    ///
    /// Same calling shape as [`Line3D::intersect_ellipse`]. Only the
    /// modeled branch participates ([`Hyperbola3D`] is single-branched):
    /// roots on the mirror branch are dropped, so a lone hit reports
    /// [`LineConic3DIntersection::Tangent`] even when it is transversal.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Hyperbola3D, Line3D, LineConic3DIntersection, Point3D, Tolerance, Vector3D};
    /// let line = Line3D::new(Point3D::new(2.0, -5.0, 0.0), Vector3D::Y).unwrap();
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&hyperbola, tol) {
    ///     LineConic3DIntersection::Tangent(hit) => {
    ///         assert_eq!(hit.point, Point3D::new(2.0, 0.0, 0.0));
    ///     }
    ///     _ => panic!("expected a tangent"),
    /// }
    /// ```
    fn intersect(&self, hyperbola: &Hyperbola3D, tol: Tolerance) -> LineConic3DIntersection {
        let frame = hyperbola.frame();
        let conic = PlacedConic {
            plane: Plane::from_frame(frame),
            x: frame.x_direction(),
            y: frame.y_direction(),
            origin: hyperbola.center(),
            implicit: ImplicitConic::hyperbola(hyperbola.major_radius(), hyperbola.minor_radius()),
            extent: f64::INFINITY,
        };
        self.intersect_frame_conic(
            &conic,
            |p| hyperbola.contains(p, tol),
            |p| hyperbola.parameter_of(p),
            tol,
        )
    }
}

impl Intersect<Line2D> for Line2D {
    type Output = LineLine2DIntersection;

    /// Intersects this line with another 2D line.
    ///
    /// The 2x2 system solves both parameters at once; directions parallel
    /// within `tol.angular` fall back to a coincidence check.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Line2D, LineLine2DIntersection, Point2D, Tolerance, Vector2D};
    /// let l1 = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let l2 = Line2D::new(Point2D::ORIGIN, Vector2D::Y).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match l1.intersect(&l2, tol) {
    ///     LineLine2DIntersection::Point(s, p, t) => {
    ///         assert_eq!((s, t), (0.0, 0.0));
    ///         assert_eq!(p, Point2D::ORIGIN);
    ///     }
    ///     _ => panic!("expected a point"),
    /// }
    /// ```
    fn intersect(&self, other: &Line2D, tol: Tolerance) -> LineLine2DIntersection {
        let (d1, d2) = (self.direction(), other.direction());
        let w = other.origin() - self.origin();
        // The sine of the angle between the unit directions.
        let cross = d1.x * d2.y - d1.y * d2.x;
        if cross.abs() <= tol.angular {
            return if other.contains(self.origin(), tol) {
                LineLine2DIntersection::Coincident
            } else {
                LineLine2DIntersection::Parallel
            };
        }
        let s = (w.x * d2.y - w.y * d2.x) / cross;
        let t = (w.x * d1.y - w.y * d1.x) / cross;
        LineLine2DIntersection::Point(s, self.eval_point(s), t)
    }
}

impl Intersect<Circle2D> for Line2D {
    type Output = LineCircle2DIntersection;

    /// Intersects this line with a 2D circle.
    ///
    /// Substituting the unit-speed parametrization into
    /// `|X - C|^2 = r^2` gives a quadratic with leading coefficient
    /// exactly 1, classified with `solve_quadratic` about the line's
    /// closest approach to the centre; the line is tangent when that
    /// closest point lies within `tol.confusion` of the circle.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Circle2D, Line2D, LineCircle2DIntersection, Point2D, Tolerance, Vector2D};
    /// let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
    /// let circle = Circle2D::new(Point2D::ORIGIN, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match line.intersect(&circle, tol) {
    ///     LineCircle2DIntersection::Points((t1, _), (t2, _)) => {
    ///         assert_eq!((t1, t2), (-2.0, 2.0));
    ///     }
    ///     _ => panic!("expected two points"),
    /// }
    /// ```
    fn intersect(&self, circle: &Circle2D, tol: Tolerance) -> LineCircle2DIntersection {
        let d = self.direction();
        let t0 = -d.dot(self.origin() - circle.center());
        let w = self.eval_point(t0) - circle.center();
        // A is exactly 1: the direction is unit by construction.
        let qb = 2.0 * d.dot(w);
        let sol = solve_quadratic(1.0, qb, w.dot(w) - circle.radius().powi(2), tol)
            .touching_at(-0.5 * qb, |s| circle.contains(self.eval_point(t0 + s), tol));
        match sol.shifted(t0) {
            QuadraticSolution::Two(t1, t2) => LineCircle2DIntersection::Points(
                (t1, self.eval_point(t1)),
                (t2, self.eval_point(t2)),
            ),
            QuadraticSolution::One(t) => LineCircle2DIntersection::Tangent(t, self.eval_point(t)),
            QuadraticSolution::Empty => LineCircle2DIntersection::Empty,
            QuadraticSolution::Linear(_) | QuadraticSolution::Degenerate => {
                unreachable!("unit direction gives a leading coefficient of 1")
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use crate::Intersect;
    use crate::{
        Axis2D, Axis3D, Circle2D, Circle3D, Cone, Cylinder, Ellipse3D, Frame3D, Hyperbola3D,
        Line2D, Line3D, LineCircle2DIntersection, LineCircle3DIntersection,
        LineConic3DIntersection, LineConstructionError, LineLine2DIntersection,
        LineLine3DIntersection, LinePlaneIntersection, LineQuadricIntersection,
        LineTorusIntersection, Parabola3D, Plane, Point2D, Point3D, Sphere, Tolerance, Torus,
        Transform, Vector2D, Vector3D, intersect_curve_surface,
    };

    // ---- Line3D construction ----

    #[test]
    fn test_line3d_new_normalizes_direction() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::new(2.0, 0.0, 0.0)).unwrap();
        assert_eq!(line.origin(), Point3D::ORIGIN);
        assert_eq!(line.direction(), Vector3D::X);
    }

    #[test]
    fn test_line3d_new_null_direction_errors() {
        assert_eq!(
            Line3D::new(Point3D::ORIGIN, Vector3D::ZERO),
            Err(LineConstructionError::NullDirection)
        );
    }

    #[test]
    fn test_line3d_from_axis() {
        let axis = Axis3D::new(Point3D::new(1.0, 2.0, 3.0), Vector3D::Y).unwrap();
        let line = Line3D::from_axis(axis);
        assert_eq!(line.origin(), Point3D::new(1.0, 2.0, 3.0));
        assert_eq!(line.direction(), Vector3D::Y);
        assert_eq!(line.axis(), axis);
    }

    #[test]
    fn test_line3d_from_two_points() {
        let p1 = Point3D::new(0.0, 0.0, 0.0);
        let p2 = Point3D::new(2.0, 0.0, 0.0);
        let line = Line3D::from_two_points(p1, p2).unwrap();
        assert_eq!(line.origin(), p1);
        assert_eq!(line.direction(), Vector3D::X);
    }

    #[test]
    fn test_line3d_from_two_points_confused_errors() {
        let p = Point3D::new(1.0, 2.0, 3.0);
        assert_eq!(
            Line3D::from_two_points(p, p),
            Err(LineConstructionError::ConfusedPoints)
        );
    }

    // ---- Line3D evaluation ----

    #[test]
    fn test_line3d_eval_point() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        assert_eq!(line.eval_point(3.0), Point3D::new(3.0, 0.0, 0.0));
    }

    #[test]
    fn test_line3d_eval_points_matches_loop() {
        let line = Line3D::new(Point3D::new(1.0, -2.0, 0.5), Vector3D::Y).unwrap();
        let us = [0.0, 1.5, -3.0];
        let expected: Vec<Point3D> = us.iter().map(|&u| line.eval_point(u)).collect();
        assert_eq!(line.eval_points(&us), expected);
    }

    #[test]
    fn test_line3d_eval_derivative_order1_is_direction() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::new(0.0, 3.0, 0.0)).unwrap();
        assert_eq!(line.eval_derivative(1.5, 1), Vector3D::Y);
    }

    #[test]
    fn test_line3d_eval_derivative_order2_is_zero() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        assert_eq!(line.eval_derivative(1.5, 2), Vector3D::ZERO);
    }

    #[test]
    fn test_line3d_eval_derivative_order3_is_zero() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        assert_eq!(line.eval_derivative(1.5, 3), Vector3D::ZERO);
    }

    #[test]
    #[should_panic]
    fn test_line3d_eval_derivative_order0_panics() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        line.eval_derivative(1.5, 0);
    }

    #[test]
    fn test_line3d_parameter_of_round_trip() {
        let line = Line3D::new(Point3D::new(1.0, -2.0, 0.5), Vector3D::new(1.0, 2.0, 2.0)).unwrap();
        for u in [0.3, 2.0, -5.5] {
            let p = line.eval_point(u);
            assert!((line.parameter_of(p) - u).abs() < 1e-9);
        }
    }

    // ---- Line2D construction ----

    #[test]
    fn test_line2d_new_normalizes_direction() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::new(0.0, 5.0)).unwrap();
        assert_eq!(line.origin(), Point2D::ORIGIN);
        assert_eq!(line.direction(), Vector2D::Y);
    }

    #[test]
    fn test_line2d_new_null_direction_errors() {
        assert_eq!(
            Line2D::new(Point2D::ORIGIN, Vector2D::ZERO),
            Err(LineConstructionError::NullDirection)
        );
    }

    #[test]
    fn test_line2d_from_axis() {
        let axis = Axis2D::new(Point2D::new(1.0, 2.0), Vector2D::X).unwrap();
        let line = Line2D::from_axis(axis);
        assert_eq!(line.origin(), Point2D::new(1.0, 2.0));
        assert_eq!(line.direction(), Vector2D::X);
        assert_eq!(line.axis(), axis);
    }

    #[test]
    fn test_line2d_from_two_points() {
        let p1 = Point2D::new(0.0, 0.0);
        let p2 = Point2D::new(0.0, 4.0);
        let line = Line2D::from_two_points(p1, p2).unwrap();
        assert_eq!(line.origin(), p1);
        assert_eq!(line.direction(), Vector2D::Y);
    }

    #[test]
    fn test_line2d_from_two_points_confused_errors() {
        let p = Point2D::new(1.0, 2.0);
        assert_eq!(
            Line2D::from_two_points(p, p),
            Err(LineConstructionError::ConfusedPoints)
        );
    }

    // ---- Line2D evaluation ----

    #[test]
    fn test_line2d_eval_point() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        assert_eq!(line.eval_point(3.0), Point2D::new(3.0, 0.0));
    }

    #[test]
    fn test_line2d_eval_points_matches_loop() {
        let line = Line2D::new(Point2D::new(1.0, -2.0), Vector2D::Y).unwrap();
        let us = [0.0, 1.5, -3.0];
        let expected: Vec<Point2D> = us.iter().map(|&u| line.eval_point(u)).collect();
        assert_eq!(line.eval_points(&us), expected);
    }

    #[test]
    fn test_line2d_eval_derivative_order1_is_direction() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::new(0.0, 3.0)).unwrap();
        assert_eq!(line.eval_derivative(1.5, 1), Vector2D::Y);
    }

    #[test]
    fn test_line2d_eval_derivative_order2_is_zero() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        assert_eq!(line.eval_derivative(1.5, 2), Vector2D::ZERO);
    }

    #[test]
    #[should_panic]
    fn test_line2d_eval_derivative_order0_panics() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        line.eval_derivative(1.5, 0);
    }

    #[test]
    fn test_line2d_parameter_of_round_trip() {
        let line = Line2D::new(Point2D::new(1.0, -2.0), Vector2D::new(3.0, 4.0)).unwrap();
        for u in [0.3, 2.0, 5.5] {
            let p = line.eval_point(u);
            assert!((line.parameter_of(p) - u).abs() < 1e-9);
        }
    }

    // ---- LineConstructionError ----

    #[test]
    fn test_line_construction_error_display() {
        assert_eq!(
            LineConstructionError::NullDirection.to_string(),
            "direction has zero length"
        );
        assert_eq!(
            LineConstructionError::ConfusedPoints.to_string(),
            "the two points are confused"
        );
    }

    #[test]
    fn test_line_construction_error_is_std_error() {
        fn takes_error(_e: &dyn std::error::Error) {}
        takes_error(&LineConstructionError::NullDirection);
    }

    #[test]
    fn test_line3d_contains() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        let tol = Tolerance::DEFAULT;
        assert!(line.contains(Point3D::new(2.5, 0.0, 0.0), tol));
        assert!(line.contains(line.eval_point(-3.0), tol));
        assert!(!line.contains(Point3D::new(2.5, 1.0, 0.0), tol));
        // Tolerance boundary: 1e-6 off the line fails at default, passes loose.
        let near = Point3D::new(0.0, 1e-6, 0.0);
        assert!(!line.contains(near, tol));
        assert!(line.contains(
            near,
            Tolerance {
                confusion: 1e-5,
                ..tol
            }
        ));
    }

    #[test]
    fn test_line2d_contains() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        let tol = Tolerance::DEFAULT;
        assert!(line.contains(Point2D::new(2.5, 0.0), tol));
        assert!(line.contains(line.eval_point(-3.0), tol));
        assert!(!line.contains(Point2D::new(2.5, 1.0), tol));
    }

    #[test]
    fn test_line3d_project_point() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        let tol = Tolerance::DEFAULT;
        let on = line.project_point(Point3D::new(2.5, 0.0, 0.0), tol);
        assert_eq!(on.parameter, 2.5);
        assert_eq!(on.distance, 0.0);
        let off = line.project_point(Point3D::new(2.5, 1.0, 0.0), tol);
        assert_eq!(off.parameter, 2.5);
        assert_eq!(off.distance, 1.0);
        // contains agrees with zero distance.
        assert_eq!(
            line.contains(Point3D::new(2.5, 0.0, 0.0), tol),
            on.distance == 0.0
        );
    }

    #[test]
    fn test_line3d_project_points_batch() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        let tol = Tolerance::DEFAULT;
        let points = [
            Point3D::new(0.0, 0.0, 0.0),
            Point3D::new(1.0, 2.0, 0.0),
            Point3D::new(-1.0, 0.0, 3.0),
        ];
        let projs = line.project_points(&points, tol);
        assert_eq!(projs.len(), 3);
        assert_eq!(projs[0].distance, 0.0);
        assert_eq!(projs[1].parameter, 1.0);
        assert_eq!(projs[1].distance, 2.0);
        assert_eq!(projs[2].parameter, -1.0);
        assert_eq!(projs[2].distance, 3.0);
    }

    #[test]
    fn test_line2d_project_point() {
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        let tol = Tolerance::DEFAULT;
        let proj = line.project_point(Point2D::new(2.5, 1.0), tol);
        assert_eq!(proj.parameter, 2.5);
        assert_eq!(proj.distance, 1.0);
    }

    #[test]
    fn test_line3d_intersect_plane() {
        let tol = Tolerance::DEFAULT;
        let plane = Plane::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
        match line.intersect(&plane, tol) {
            LinePlaneIntersection::Point(t, p) => {
                assert_eq!(t, 0.0);
                assert_eq!(p, Point3D::ORIGIN);
            }
            _ => panic!("expected a point"),
        }
        // Parallel offset.
        let off = Line3D::new(Point3D::new(0.0, 0.0, 1.0), Vector3D::X).unwrap();
        assert_eq!(off.intersect(&plane, tol), LinePlaneIntersection::Parallel);
        // Contained.
        let flat = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        assert_eq!(
            flat.intersect(&plane, tol),
            LinePlaneIntersection::Coincident
        );
    }

    #[test]
    fn test_line3d_intersect_sphere() {
        let tol = Tolerance::DEFAULT;
        let sphere = Sphere::new(Point3D::ORIGIN, 2.0).unwrap();
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match line.intersect(&sphere, tol) {
            LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
                assert_eq!((t1, t2), (-2.0, 2.0));
                assert_eq!(p1, Point3D::new(-2.0, 0.0, 0.0));
                assert_eq!(p2, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected two points"),
        }
        // Tangent at the north pole.
        let tangent = Line3D::new(Point3D::new(0.0, 0.0, 2.0), Vector3D::X).unwrap();
        match tangent.intersect(&sphere, tol) {
            LineQuadricIntersection::Tangent(t, p) => {
                assert_eq!(t, 0.0);
                assert_eq!(p, Point3D::new(0.0, 0.0, 2.0));
            }
            _ => panic!("expected a tangent"),
        }
        // Clear miss.
        let miss = Line3D::new(Point3D::new(0.0, 0.0, 3.0), Vector3D::X).unwrap();
        assert_eq!(miss.intersect(&sphere, tol), LineQuadricIntersection::Empty);
    }

    #[test]
    fn test_line3d_intersect_cylinder() {
        let tol = Tolerance::DEFAULT;
        let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match line.intersect(&cylinder, tol) {
            LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
                assert_eq!((t1, t2), (-2.0, 2.0));
                assert_eq!(p1, Point3D::new(-2.0, 0.0, 0.0));
                assert_eq!(p2, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected two points"),
        }
        // Parallel offset: miss.
        let side = Line3D::new(Point3D::new(3.0, 0.0, 0.0), Vector3D::Z).unwrap();
        assert_eq!(
            side.intersect(&cylinder, tol),
            LineQuadricIntersection::Empty
        );
        // Generator: coincident.
        let generator = Line3D::new(Point3D::new(2.0, 0.0, 0.0), Vector3D::Z).unwrap();
        assert_eq!(
            generator.intersect(&cylinder, tol),
            LineQuadricIntersection::Coincident
        );
    }

    #[test]
    fn test_line3d_intersect_cone() {
        let tol = Tolerance::DEFAULT;
        let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        // Horizontal line at z = 5: radius there is (5 + 4.729...) * tan(0.4).
        let line = Line3D::new(Point3D::new(0.0, 0.0, 5.0), Vector3D::X).unwrap();
        match line.intersect(&cone, tol) {
            LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
                assert!(t1 < t2);
                assert!((t1.abs() - t2.abs()).abs() < 1e-9);
                assert!(cone.contains(p1, tol));
                assert!(cone.contains(p2, tol));
            }
            _ => panic!("expected two points"),
        }
        // Generator through the apex: coincident.
        let apex = cone.apex();
        let generator = Line3D::new(apex, Vector3D::new(0.4f64.sin(), 0.0, 0.4f64.cos())).unwrap();
        assert_eq!(
            generator.intersect(&cone, tol),
            LineQuadricIntersection::Coincident
        );
        // Vertical line at radial distance 20: one hit on each nappe.
        let side = Line3D::new(Point3D::new(20.0, 0.0, 0.0), Vector3D::Z).unwrap();
        match side.intersect(&cone, tol) {
            LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
                assert!(p1.z < apex.z && p2.z > apex.z);
                for (t, p) in [(t1, p1), (t2, p2)] {
                    assert!(cone.contains(p, tol));
                    assert_eq!(p, side.eval_point(t));
                }
            }
            _ => panic!("expected two points"),
        }
        // Horizontal line far from the axis: clear miss.
        let miss = Line3D::new(Point3D::new(0.0, 20.0, 0.0), Vector3D::X).unwrap();
        assert_eq!(miss.intersect(&cone, tol), LineQuadricIntersection::Empty);
        // Generator direction but offset: linear single hit.
        let slanted = Line3D::new(
            Point3D::new(0.0, 5.0, 0.0),
            Vector3D::new(0.4f64.sin(), 0.0, 0.4f64.cos()),
        )
        .unwrap();
        match slanted.intersect(&cone, tol) {
            LineQuadricIntersection::OnePoint(t, p) => {
                assert!(cone.contains(p, tol));
                assert_eq!(p, slanted.eval_point(t));
            }
            _ => panic!("expected one point"),
        }
    }

    #[test]
    fn test_line3d_intersect_line() {
        let tol = Tolerance::DEFAULT;
        let l1 = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        let l2 = Line3D::new(Point3D::ORIGIN, Vector3D::Y).unwrap();
        match l1.intersect(&l2, tol) {
            LineLine3DIntersection::Point(s, p, t) => {
                assert_eq!((s, t), (0.0, 0.0));
                assert_eq!(p, Point3D::ORIGIN);
            }
            _ => panic!("expected a point"),
        }
        // Skew: x-axis against a z-shifted y-line.
        let skew = Line3D::new(Point3D::new(0.0, 0.0, 1.0), Vector3D::Y).unwrap();
        match l1.intersect(&skew, tol) {
            LineLine3DIntersection::Skew { s, t, distance } => {
                assert_eq!((s, t), (0.0, 0.0));
                assert_eq!(distance, 1.0);
            }
            _ => panic!("expected skew"),
        }
        // Parallel distinct.
        let parallel = Line3D::new(Point3D::new(0.0, 1.0, 0.0), Vector3D::X).unwrap();
        match l1.intersect(&parallel, tol) {
            LineLine3DIntersection::Parallel(d) => assert_eq!(d, 1.0),
            _ => panic!("expected parallel"),
        }
        // Coincident.
        let same = Line3D::new(Point3D::new(2.0, 0.0, 0.0), Vector3D::X).unwrap();
        assert_eq!(l1.intersect(&same, tol), LineLine3DIntersection::Coincident);
    }

    #[test]
    fn test_line2d_intersect_line() {
        let tol = Tolerance::DEFAULT;
        let l1 = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        let l2 = Line2D::new(Point2D::ORIGIN, Vector2D::Y).unwrap();
        match l1.intersect(&l2, tol) {
            LineLine2DIntersection::Point(s, p, t) => {
                assert_eq!((s, t), (0.0, 0.0));
                assert_eq!(p, Point2D::ORIGIN);
            }
            _ => panic!("expected a point"),
        }
        let parallel = Line2D::new(Point2D::new(0.0, 1.0), Vector2D::X).unwrap();
        assert_eq!(
            l1.intersect(&parallel, tol),
            LineLine2DIntersection::Parallel
        );
        assert_eq!(l1.intersect(&l1, tol), LineLine2DIntersection::Coincident);
    }

    #[test]
    fn test_line2d_intersect_circle() {
        let tol = Tolerance::DEFAULT;
        let line = Line2D::new(Point2D::ORIGIN, Vector2D::X).unwrap();
        let circle = Circle2D::new(Point2D::ORIGIN, 2.0).unwrap();
        match line.intersect(&circle, tol) {
            LineCircle2DIntersection::Points((t1, p1), (t2, p2)) => {
                assert_eq!((t1, t2), (-2.0, 2.0));
                assert_eq!(p1, Point2D::new(-2.0, 0.0));
                assert_eq!(p2, Point2D::new(2.0, 0.0));
            }
            _ => panic!("expected two points"),
        }
        // Tangent at the top.
        let tangent = Line2D::new(Point2D::new(0.0, 2.0), Vector2D::X).unwrap();
        match tangent.intersect(&circle, tol) {
            LineCircle2DIntersection::Tangent(t, p) => {
                assert_eq!(t, 0.0);
                assert_eq!(p, Point2D::new(0.0, 2.0));
            }
            _ => panic!("expected a tangent"),
        }
        // Clear miss.
        let miss = Line2D::new(Point2D::new(0.0, 3.0), Vector2D::X).unwrap();
        assert_eq!(
            miss.intersect(&circle, tol),
            LineCircle2DIntersection::Empty
        );
    }

    #[test]
    fn test_line3d_intersect_circle() {
        let tol = Tolerance::DEFAULT;
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        let circle = Circle3D::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
        // In-plane line: two hits.
        match line.intersect(&circle, tol) {
            LineCircle3DIntersection::Points((t1, p1), (t2, p2)) => {
                assert_eq!((t1, t2), (-2.0, 2.0));
                assert_eq!(p1, Point3D::new(-2.0, 0.0, 0.0));
                assert_eq!(p2, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected two points"),
        }
        // Transversal graze through a circle point.
        let graze = Line3D::new(Point3D::new(2.0, 0.0, 0.0), Vector3D::Z).unwrap();
        match graze.intersect(&circle, tol) {
            LineCircle3DIntersection::Tangent(t, p) => {
                assert_eq!(t, 0.0);
                assert_eq!(p, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected a tangent"),
        }
        // Transversal miss.
        let miss = Line3D::new(Point3D::new(3.0, 0.0, 0.0), Vector3D::Z).unwrap();
        assert_eq!(
            miss.intersect(&circle, tol),
            LineCircle3DIntersection::Empty
        );
    }

    #[test]
    fn test_line3d_intersect_circle_in_plane_within_loose_tolerance() {
        // 1e-4 off the circle's plane is in the plane at confusion 1e-3, so
        // the line meets the circle twice (#26).
        let loose = Tolerance {
            confusion: 1e-3,
            ..Tolerance::DEFAULT
        };
        let circle = Circle3D::new(Point3D::ORIGIN, Vector3D::Z, 1.0).unwrap();
        for height in [1e-5, 1e-4] {
            let line = Line3D::new(Point3D::new(0.0, 0.0, height), Vector3D::X).unwrap();
            match line.intersect(&circle, loose) {
                LineCircle3DIntersection::Points((t1, p1), (t2, p2)) => {
                    assert!((t1 + 1.0).abs() < 1e-12 && (t2 - 1.0).abs() < 1e-12);
                    assert!(circle.contains(p1, loose) && circle.contains(p2, loose));
                    assert!(line.contains(p1, loose) && line.contains(p2, loose));
                }
                other => panic!("expected two points at height {height}, got {other:?}"),
            }
        }
    }

    #[test]
    fn test_line3d_intersect_ellipse() {
        let tol = Tolerance::DEFAULT;
        let ellipse = Ellipse3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 3.0, 1.5).unwrap();
        // In-plane line through both vertices.
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match line.intersect(&ellipse, tol) {
            LineConic3DIntersection::Points(h1, h2) => {
                assert_eq!((h1.line_param, h2.line_param), (-3.0, 3.0));
                assert_eq!(
                    (h1.point, h2.point),
                    (Point3D::new(-3.0, 0.0, 0.0), Point3D::new(3.0, 0.0, 0.0))
                );
                assert_eq!(h2.conic_param, 0.0);
            }
            _ => panic!("expected two points"),
        }
        // Transversal graze through a vertex.
        let graze = Line3D::new(Point3D::new(3.0, 0.0, 0.0), Vector3D::Z).unwrap();
        match graze.intersect(&ellipse, tol) {
            LineConic3DIntersection::Tangent(hit) => {
                assert_eq!(hit.line_param, 0.0);
                assert_eq!(hit.point, Point3D::new(3.0, 0.0, 0.0));
                assert_eq!(hit.conic_param, 0.0);
            }
            _ => panic!("expected a tangent"),
        }
        // Transversal miss.
        let miss = Line3D::new(Point3D::new(4.0, 0.0, 0.0), Vector3D::Z).unwrap();
        assert_eq!(
            miss.intersect(&ellipse, tol),
            LineConic3DIntersection::Empty
        );
    }

    #[test]
    fn test_line3d_intersect_parabola() {
        let tol = Tolerance::DEFAULT;
        let parabola = Parabola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 1.0).unwrap();
        // In-plane axis: the linear root at the apex.
        let axis = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match axis.intersect(&parabola, tol) {
            LineConic3DIntersection::Tangent(hit) => {
                assert_eq!(hit.line_param, 0.0);
                assert_eq!(hit.point, Point3D::ORIGIN);
            }
            _ => panic!("expected a tangent"),
        }
        // In-plane vertical line: two symmetric hits.
        let vertical = Line3D::new(Point3D::new(1.0, 0.0, 0.0), Vector3D::Y).unwrap();
        match vertical.intersect(&parabola, tol) {
            LineConic3DIntersection::Points(h1, h2) => {
                assert_eq!((h1.line_param, h2.line_param), (-2.0, 2.0));
                assert_eq!(
                    (h1.point, h2.point),
                    (Point3D::new(1.0, -2.0, 0.0), Point3D::new(1.0, 2.0, 0.0))
                );
            }
            _ => panic!("expected two points"),
        }
        // Parallel miss above the plane.
        let miss = Line3D::new(Point3D::new(0.0, 0.0, 1.0), Vector3D::X).unwrap();
        assert_eq!(
            miss.intersect(&parabola, tol),
            LineConic3DIntersection::Empty
        );
    }

    #[test]
    fn test_line3d_intersect_hyperbola() {
        let tol = Tolerance::DEFAULT;
        let hyperbola =
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        // Tangent at the right vertex.
        let tangent = Line3D::new(Point3D::new(2.0, -5.0, 0.0), Vector3D::Y).unwrap();
        match tangent.intersect(&hyperbola, tol) {
            LineConic3DIntersection::Tangent(hit) => {
                assert_eq!(hit.line_param, 5.0);
                assert_eq!(hit.point, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected a tangent"),
        }
        // Axis line: the mirror-branch root is dropped, leaving a lone hit.
        let axis = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match axis.intersect(&hyperbola, tol) {
            LineConic3DIntersection::Tangent(hit) => {
                assert_eq!(hit.line_param, 2.0);
                assert_eq!(hit.point, Point3D::new(2.0, 0.0, 0.0));
            }
            _ => panic!("expected a lone hit"),
        }
        // Parallel miss above the plane.
        let miss = Line3D::new(Point3D::new(0.0, 0.0, 1.0), Vector3D::X).unwrap();
        assert_eq!(
            miss.intersect(&hyperbola, tol),
            LineConic3DIntersection::Empty
        );
    }

    #[test]
    fn test_line3d_intersect_torus() {
        let tol = Tolerance::DEFAULT;
        let torus = Torus::new(Point3D::ORIGIN, Vector3D::Z, 4.0, 1.0).unwrap();
        // Equatorial line through the tube: four transversal hits.
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        match line.intersect(&torus, tol) {
            LineTorusIntersection::Hits(hits) => {
                assert_eq!(hits.len(), 4);
                let params: Vec<f64> = hits.iter().map(|h| h.line_param).collect();
                assert_eq!(params, vec![-5.0, -3.0, 3.0, 5.0]);
                for h in &hits {
                    assert_eq!(h.point, line.eval_point(h.line_param));
                    assert!(torus.contains(h.point, tol));
                    assert_eq!(h.multiplicity, 1);
                    let (u, v) = torus.parameters_of(h.point);
                    assert!((h.surface_params.0 - u).abs() < 1e-12);
                    assert!((h.surface_params.1 - v).abs() < 1e-12);
                }
            }
            _ => panic!("expected four hits"),
        }
        // Grazing the outer equator: a lone double hit.
        let graze = Line3D::new(Point3D::new(0.0, 5.0, 0.0), Vector3D::X).unwrap();
        match graze.intersect(&torus, tol) {
            LineTorusIntersection::Hits(hits) => {
                assert_eq!(hits.len(), 1);
                assert_eq!(hits[0].line_param, 0.0);
                assert_eq!(hits[0].point, Point3D::new(0.0, 5.0, 0.0));
                assert_eq!(hits[0].multiplicity, 2);
            }
            _ => panic!("expected a tangent hit"),
        }
        // Clear miss outside the tube.
        let miss = Line3D::new(Point3D::new(0.0, 6.0, 0.0), Vector3D::X).unwrap();
        assert_eq!(miss.intersect(&torus, tol), LineTorusIntersection::Empty);
        // Through the hole: no contact.
        let hole = Line3D::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
        assert_eq!(hole.intersect(&torus, tol), LineTorusIntersection::Empty);
    }

    #[test]
    fn test_line3d_intersect_torus_matches_generic() {
        // Deterministic pseudo-random lines (LCG) cross-checked against the
        // generic curve-surface solver: same count, same parameters.
        let tol = Tolerance::DEFAULT;
        let torus = Torus::new(Point3D::ORIGIN, Vector3D::Z, 4.0, 1.0).unwrap();
        let mut state: u64 = 0x1234_5678_9abc_def0;
        let mut next = || {
            state = state
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            (state >> 33) as f64 / 2147483648.0
        };
        for _ in 0..25 {
            let origin = Point3D::new(12.0 * next() - 6.0, 12.0 * next() - 6.0, 6.0 * next() - 3.0);
            let dir = Vector3D::new(2.0 * next() - 1.0, 2.0 * next() - 1.0, 2.0 * next() - 1.0);
            let Some(dir) = dir.normalized() else {
                continue;
            };
            let line = Line3D::new(origin, dir).unwrap();
            let analytic = line.intersect(&torus, tol);
            let generic = intersect_curve_surface(&line, &torus, tol);
            match analytic {
                LineTorusIntersection::Hits(hits) => {
                    assert_eq!(
                        hits.len(),
                        generic.len(),
                        "count mismatch for line {origin:?} {dir:?}"
                    );
                    for (h, g) in hits.iter().zip(generic.iter()) {
                        assert!(
                            (h.line_param - g.curve_param).abs() < 1e-6,
                            "param mismatch for line {origin:?} {dir:?}"
                        );
                        assert!(h.point.distance(g.point) < 1e-6);
                    }
                }
                LineTorusIntersection::Empty => assert!(
                    generic.is_empty(),
                    "generic found hits for line {origin:?} {dir:?}"
                ),
            }
        }
    }

    #[test]
    fn test_line3d_sphere_tangent_with_rounded_construction() {
        // A line built tangent to the sphere (its foot point placed at
        // distance r from the centre, then rounded): the closest approach is
        // within tolerance of the sphere, so the hit is a tangent even though
        // rounding pushed the double root apart in parameter space.
        let tol = Tolerance::DEFAULT;
        let sphere = Sphere::new(
            Point3D::new(-3.66014770679811, 2.1882136950508135, 0.8308364993480888),
            4.576151540371945,
        )
        .unwrap();
        let line = Line3D::new(
            Point3D::new(-7.265034215586097, 4.063943363277678, 3.160556308883224),
            Vector3D::new(
                -0.029375202317363944,
                -0.9552316968774335,
                0.2943968457192993,
            ),
        )
        .unwrap();
        match line.intersect(&sphere, tol) {
            LineQuadricIntersection::Tangent(t, p) => {
                assert_eq!(p, line.eval_point(t));
                assert!((p.distance(sphere.center()) - sphere.radius()).abs() < 1e-7);
            }
            other => panic!("expected a tangent, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_sphere_tangency_is_a_distance_criterion() {
        // Unit sphere, lines along X at height 1 + h: within tolerance of the
        // surface the contact is a tangent at the closest point, beyond it
        // the line misses or crosses.
        let tol = Tolerance::DEFAULT;
        let sphere = Sphere::new(Point3D::ORIGIN, 1.0).unwrap();
        for (h, expect) in [
            (5e-8, "tangent"),
            (-5e-8, "tangent"),
            (1e-6, "empty"),
            (-1e-6, "two"),
        ] {
            let line = Line3D::new(Point3D::new(-3.0, 1.0 + h, 0.0), Vector3D::X).unwrap();
            let got = line.intersect(&sphere, tol);
            match (expect, got) {
                ("tangent", LineQuadricIntersection::Tangent(t, _)) => {
                    assert!((t - 3.0).abs() < 1e-12)
                }
                ("empty", LineQuadricIntersection::Empty) => {}
                ("two", LineQuadricIntersection::TwoPoints(..)) => {}
                (_, other) => panic!("h = {h}: expected {expect}, got {other:?}"),
            }
        }
    }

    #[test]
    fn test_line3d_cylinder_tangent_with_rounded_construction() {
        // A line through a point at distance r from the axis, perpendicular to
        // both the axis and the radial direction (all rounded).
        let tol = Tolerance::DEFAULT;
        let frame = Frame3D::new(
            Point3D::new(-4.5491502150862, 3.488975760746989, 4.505226342025091),
            Vector3D::new(
                0.3378045276092872,
                -0.8749064394010623,
                -0.34702568121281996,
            ),
            Vector3D::new(
                -0.9411459179864347,
                -0.31849077016586064,
                -0.11317239228994093,
            ),
        )
        .unwrap();
        let cylinder = Cylinder::from_frame(frame, 4.790993298008813).unwrap();
        let line = Line3D::new(
            Point3D::new(-9.058174000607561, 1.9630886154046856, 3.96301816904436),
            Vector3D::new(
                -0.011509221699911018,
                0.36483194982582373,
                -0.9310022482261524,
            ),
        )
        .unwrap();
        match line.intersect(&cylinder, tol) {
            LineQuadricIntersection::Tangent(t, p) => {
                assert!(t.abs() < 1e-9);
                assert!(cylinder.contains(p, tol));
            }
            other => panic!("expected a tangent, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_torus_equator_tangents_with_rounded_construction() {
        // Lines tangent to the outer or inner equator of a tilted torus, three
        // units from their origin: one double hit at parameter 3 (plus the two
        // transversal crossings of the far side for the inner equator).
        let tol = Tolerance::DEFAULT;
        let cases = [
            (
                [-3.151766078083671, 1.7612911834316094, -3.8732580004320094],
                [0.591615246726186, -0.1992663849953292, 0.7812069557114858],
                [-0.11806191907195286, 0.9371156800331439, 0.3284441893245577],
                2.57900662578319,
                2.1693570817529384,
                [-1.3197794883641474, 7.070687415724088, -3.9063439891962752],
                [-0.7975291738276172, -0.2865433824866865, 0.5308861524347069],
                [-3.7123670098469996, 6.211057268264029, -2.313685531892154],
                1,
            ),
            (
                [
                    -0.42406925077418656,
                    -1.5070564593366398,
                    -1.7711640354794786,
                ],
                [0.8814613008971334, 0.29475868487218504, 0.3689760056062631],
                [-0.4162021275760656, 0.8540654733307351, 0.3120063401046585],
                5.377738053276436,
                1.313893854677209,
                [-1.4459605871845431, 3.249502898666323, -3.1297334451352716],
                [
                    -0.22316308839476426,
                    -0.4285901129946508,
                    0.8755048549387623,
                ],
                [-2.115449852368836, 1.9637325596823705, -0.5032188803189848],
                3,
            ),
        ];
        for (o, z, x, major, minor, lo, ld, contact, count) in cases {
            let frame = Frame3D::new(
                Point3D::new(o[0], o[1], o[2]),
                Vector3D::new(z[0], z[1], z[2]),
                Vector3D::new(x[0], x[1], x[2]),
            )
            .unwrap();
            let torus = Torus::from_frame(frame, major, minor).unwrap();
            let line = Line3D::new(
                Point3D::new(lo[0], lo[1], lo[2]),
                Vector3D::new(ld[0], ld[1], ld[2]),
            )
            .unwrap();
            let contact = Point3D::new(contact[0], contact[1], contact[2]);
            match line.intersect(&torus, tol) {
                LineTorusIntersection::Hits(hits) => {
                    assert_eq!(hits.len(), count, "{hits:?}");
                    let touch = hits
                        .iter()
                        .find(|h| h.multiplicity == 2)
                        .unwrap_or_else(|| panic!("no tangent hit in {hits:?}"));
                    assert!((touch.line_param - 3.0).abs() < 1e-6);
                    assert!(touch.point.distance(contact) < 1e-6);
                }
                other => panic!("expected hits, got {other:?}"),
            }
        }
    }

    #[test]
    fn test_line3d_far_origin_small_sphere_two_points() {
        // The line starts a million units away from a sphere of radius
        // 0.05: the chord is still two points, not a grazing contact.
        let tol = Tolerance::DEFAULT;
        let sphere = Sphere::new(Point3D::ORIGIN, 0.05).unwrap();
        let line = Line3D::new(Point3D::new(-1e6, 0.0, 0.0), Vector3D::X).unwrap();
        match line.intersect(&sphere, tol) {
            LineQuadricIntersection::TwoPoints((t1, p1), (t2, p2)) => {
                assert!((t1 - (1e6 - 0.05)).abs() < 1e-9);
                assert!((t2 - (1e6 + 0.05)).abs() < 1e-9);
                assert!((p1.x + 0.05).abs() < 1e-9 && (p2.x - 0.05).abs() < 1e-9);
            }
            other => panic!("expected two points, got {other:?}"),
        }
        // A genuine graze far from the origin stays a tangent.
        let graze = Line3D::new(Point3D::new(-1e6, 0.05, 0.0), Vector3D::X).unwrap();
        match graze.intersect(&sphere, tol) {
            LineQuadricIntersection::Tangent(t, p) => {
                assert!((t - 1e6).abs() < 1e-9);
                assert!(p.distance(Point3D::new(0.0, 0.05, 0.0)) < 1e-9);
            }
            other => panic!("expected a tangent, got {other:?}"),
        }
    }

    #[test]
    fn test_line2d_far_origin_small_circle_two_points() {
        let tol = Tolerance::DEFAULT;
        let circle = Circle2D::new(Point2D::ORIGIN, 0.05).unwrap();
        let line = Line2D::new(Point2D::new(-1e6, 0.0), Vector2D::X).unwrap();
        match line.intersect(&circle, tol) {
            LineCircle2DIntersection::Points((t1, _), (t2, _)) => {
                assert!((t1 - (1e6 - 0.05)).abs() < 1e-9);
                assert!((t2 - (1e6 + 0.05)).abs() < 1e-9);
            }
            other => panic!("expected two points, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_far_origin_thin_cylinder_two_points() {
        let tol = Tolerance::DEFAULT;
        let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 0.05).unwrap();
        let line = Line3D::new(Point3D::new(-1e6, 0.0, 3.0), Vector3D::X).unwrap();
        match line.intersect(&cylinder, tol) {
            LineQuadricIntersection::TwoPoints((t1, _), (t2, _)) => {
                assert!((t1 - (1e6 - 0.05)).abs() < 1e-9);
                assert!((t2 - (1e6 + 0.05)).abs() < 1e-9);
            }
            other => panic!("expected two points, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_torus_hits_from_offset_origin() {
        // The line's origin sits outside the torus: the parameters of the
        // four hits are offset by 300 but must stay exact.
        let tol = Tolerance::DEFAULT;
        let torus = Torus::new(Point3D::ORIGIN, Vector3D::Z, 100.0, 25.0).unwrap();
        let line = Line3D::new(Point3D::new(-300.0, 0.0, 0.0), Vector3D::X).unwrap();
        match line.intersect(&torus, tol) {
            LineTorusIntersection::Hits(hits) => {
                let params: Vec<f64> = hits.iter().map(|h| h.line_param).collect();
                assert_eq!(params, vec![175.0, 225.0, 375.0, 425.0]);
                assert!(hits.iter().all(|h| h.multiplicity == 1));
            }
            other => panic!("expected four hits, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_far_origin_small_torus_four_hits() {
        // A torus far smaller than the line's offset from it: four distinct
        // hits, none merged into tangencies by an offset-scaled band.
        let tol = Tolerance::DEFAULT;
        let torus = Torus::new(Point3D::ORIGIN, Vector3D::Z, 0.04, 0.01).unwrap();
        let line = Line3D::new(Point3D::new(-1e5, 0.0, 0.0), Vector3D::X).unwrap();
        match line.intersect(&torus, tol) {
            LineTorusIntersection::Hits(hits) => {
                assert_eq!(hits.len(), 4);
                let expected = [-0.05, -0.03, 0.03, 0.05];
                for (h, x) in hits.iter().zip(expected) {
                    assert!((h.line_param - (1e5 + x)).abs() < 1e-9);
                    assert_eq!(h.multiplicity, 1);
                }
            }
            other => panic!("expected four hits, got {other:?}"),
        }
    }

    #[test]
    fn test_line3d_in_plane_thin_ellipse_far_origin() {
        // A nearly vertical in-plane line crossing a thin, elongated
        // ellipse near the end of its major axis, the line's origin about
        // 30 units away, at a loose 1e-3 tolerance: the short chord is two
        // points, not a tangency.
        let tol = Tolerance {
            confusion: 1e-3,
            ..Tolerance::DEFAULT
        };
        let ellipse = Ellipse3D::new(
            Point3D::new(-610.348096534595, -710.720096056787, 0.0),
            Vector3D::Z,
            Vector3D::new(0.999999902285153, 0.000442074298181498, 0.0),
            15.0066332711999,
            0.291884102212871,
        )
        .unwrap();
        let line = Line3D::new(
            Point3D::new(-625.34430362036, -680.713463264921, 0.0),
            Vector3D::new(-0.000252749178714602, -0.999999968058926, 0.0),
        )
        .unwrap();
        match line.intersect(&ellipse, tol) {
            LineConic3DIntersection::Points(h1, h2) => {
                let expected = Point3D::new(-625.3518880129151, -710.7210476574684, 0.0);
                let near = h1.point.distance(expected).min(h2.point.distance(expected));
                assert!(near < 1e-6, "{h1:?} {h2:?}");
                assert!(h1.line_param < h2.line_param);
            }
            other => panic!("expected two points, got {other:?}"),
        }
    }

    // ---- Line3D transformed ----

    #[test]
    fn test_line3d_transformed_translation_spot_check() {
        let line = Line3D::new(Point3D::new(1.0, 0.0, 0.0), Vector3D::Y).unwrap();
        let t = Transform::translation(Vector3D::new(0.0, 5.0, -2.0));
        let moved = line.transformed(t).unwrap();
        assert_eq!(moved.origin(), Point3D::new(1.0, 5.0, -2.0));
        assert_eq!(moved.direction(), Vector3D::Y);
        // Hand-computed evaluation: origin + 2 * direction.
        assert_eq!(moved.eval_point(2.0), Point3D::new(1.0, 7.0, -2.0));
    }

    #[test]
    fn test_line3d_transformed_rotation_spot_check() {
        use std::f64::consts::PI;
        let line = Line3D::new(Point3D::new(1.0, 0.0, 0.0), Vector3D::X).unwrap();
        let axis = Axis3D::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
        let moved = line
            .transformed(Transform::rotation(axis, PI / 2.0))
            .unwrap();
        assert!(moved.origin().distance(Point3D::new(0.0, 1.0, 0.0)) < 1e-9);
        assert!((moved.direction() - Vector3D::Y).magnitude() < 1e-9);
    }

    #[test]
    fn test_line3d_transformed_round_trip() {
        let line = Line3D::new(Point3D::new(1.0, -2.0, 0.5), Vector3D::new(1.0, 2.0, 2.0)).unwrap();
        let axis = Axis3D::new(Point3D::new(0.5, 0.5, 0.0), Vector3D::Z).unwrap();
        let t = Transform::translation(Vector3D::new(3.0, -1.0, 2.0))
            .then(Transform::rotation(axis, 0.7))
            .then(Transform::scaling(Point3D::ORIGIN, 2.0));
        let moved = line.transformed(t).unwrap();
        // Transformed points stay on the transformed line.
        let tol = Tolerance::DEFAULT;
        for u in [0.0, 1.5, -2.5] {
            assert!(moved.contains(t.apply_point(line.eval_point(u)), tol));
        }
        // Round trip through the inverse restores the line.
        let inv = Transform::scaling(Point3D::ORIGIN, 0.5)
            .then(Transform::rotation(axis, -0.7))
            .then(Transform::translation(Vector3D::new(-3.0, 1.0, -2.0)));
        let back = moved.transformed(inv).unwrap();
        assert!(back.origin().distance(line.origin()) < 1e-9);
        assert!((back.direction().dot(line.direction()) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_line3d_transformed_mirror() {
        let line = Line3D::new(Point3D::new(1.0, 2.0, 3.0), Vector3D::X).unwrap();
        let m = Transform::mirror_plane(Frame3D::WORLD);
        let moved = line.transformed(m).unwrap();
        assert_eq!(moved.origin(), Point3D::new(1.0, 2.0, -3.0));
        assert_eq!(moved.direction(), Vector3D::X);
        // Mirroring twice is the identity.
        let back = moved.transformed(m).unwrap();
        assert_eq!(back, line);
    }

    #[test]
    fn test_line3d_transformed_zero_scale_errors() {
        let line = Line3D::new(Point3D::ORIGIN, Vector3D::X).unwrap();
        assert_eq!(
            line.transformed(Transform::scaling(Point3D::ORIGIN, 0.0)),
            Err(LineConstructionError::NullDirection)
        );
    }
}