geomcore 0.32.2

//! Hyperbolas in 3D: parametric evaluation, parameter inversion, and
//! center+two-points construction, a thin wrapper over
//! [`crate::curve_math::analytic`].

use crate::Intersect;
use crate::curve_math::analytic;
use crate::curves::{Curve2D, ParametrizeError};
use crate::frame::frame_error;
use crate::intersect::{ConicSurfaceIntersection, intersect_conic_quadric};
use crate::math::{polish_root, real_roots};
use crate::projection::{self, CurveProjection};
use crate::surfaces::Surface;
use crate::tol;
use crate::{Cone, Cylinder, Frame3D, Plane, Point3D, Sphere, Tolerance, Transform, Vector3D};
use std::fmt;

/// Error returned when a [`Hyperbola3D`] cannot be constructed from the
/// given inputs.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum HyperbolaConstructionError {
    /// The major or minor radius is zero or negative (a zero radius
    /// collapses the hyperbola onto lines).
    NonPositiveRadius,
    /// The normal (or the x direction) has zero length, or the x direction
    /// is parallel to the normal.
    NullNormal,
    /// Two (or more) of the three points given to build the hyperbola are
    /// coincident (or too close to distinguish).
    ConfusedPoints,
    /// The second point lies on the major-axis line (or is otherwise
    /// collinear with the center and the first point), so no minor axis can
    /// be derived.
    CollinearPoints,
    /// An input coordinate, radius, or angle is NaN or infinite.
    NonFinite,
}

impl fmt::Display for HyperbolaConstructionError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let message = match self {
            HyperbolaConstructionError::NonPositiveRadius => "radius is not positive",
            HyperbolaConstructionError::NullNormal => "normal has zero length",
            HyperbolaConstructionError::ConfusedPoints => "the points are confused",
            HyperbolaConstructionError::CollinearPoints => "the points are collinear",
            HyperbolaConstructionError::NonFinite => "an input is NaN or infinite",
        };
        f.write_str(message)
    }
}

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

/// A hyperbola in 3D: a plane [`Frame3D`] (origin at the center, plus local
/// x/y directions defining the plane, the major axis, and the transverse
/// direction), a semi-major radius, and a semi-minor radius, evaluated as
/// `center + major*cosh(u)*x_dir + minor*sinh(u)*y_dir`.
///
/// # Examples
///
/// ```
/// use geomcore::{Hyperbola3D, Point3D, Vector3D};
/// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
/// assert_eq!(hyperbola.eval_point(0.0), Point3D::new(2.0, 0.0, 0.0));
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Hyperbola3D {
    frame: Frame3D,
    major_radius: f64,
    minor_radius: f64,
}

impl Hyperbola3D {
    /// Creates a hyperbola from a center, a normal, an x-axis hint, a
    /// semi-major radius, and a semi-minor radius.
    ///
    /// The plane frame is derived from `normal` and `x_direction` via
    /// [`Frame3D::new`].
    ///
    /// # Errors
    ///
    /// Returns [`HyperbolaConstructionError::NullNormal`] if `normal` cannot
    /// be normalized, if `x_direction` cannot be normalized, or if
    /// `x_direction` is parallel to `normal`;
    /// [`HyperbolaConstructionError::NonPositiveRadius`] if
    /// `major_radius <= 0` or `minor_radius <= 0`; or
    /// [`HyperbolaConstructionError::NonFinite`] if any input is NaN or
    /// infinite.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.major_radius(), 2.0);
    /// assert_eq!(hyperbola.minor_radius(), 1.0);
    /// ```
    pub fn new(
        center: Point3D,
        normal: Vector3D,
        x_direction: Vector3D,
        major_radius: f64,
        minor_radius: f64,
    ) -> Result<Hyperbola3D, HyperbolaConstructionError> {
        let frame = Frame3D::new(center, normal, x_direction).map_err(frame_error(
            HyperbolaConstructionError::NonFinite,
            HyperbolaConstructionError::NullNormal,
        ))?;
        Hyperbola3D::from_frame(frame, major_radius, minor_radius)
    }

    /// Creates a hyperbola from a plane frame, a semi-major radius, and a
    /// semi-minor radius directly.
    ///
    /// Unlike [`crate::Ellipse3D`], `minor_radius` may exceed `major_radius`
    /// (only positivity is checked).
    ///
    /// # Errors
    ///
    /// Returns [`HyperbolaConstructionError::NonPositiveRadius`] if
    /// `major_radius <= 0` or `minor_radius <= 0`, or
    /// [`HyperbolaConstructionError::NonFinite`] if either radius is NaN or
    /// infinite.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Frame3D};
    /// let hyperbola = Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.frame(), Frame3D::WORLD);
    /// ```
    pub fn from_frame(
        frame: Frame3D,
        major_radius: f64,
        minor_radius: f64,
    ) -> Result<Hyperbola3D, HyperbolaConstructionError> {
        if !(major_radius.is_finite() && minor_radius.is_finite()) {
            return Err(HyperbolaConstructionError::NonFinite);
        }
        if major_radius <= 0.0 || minor_radius <= 0.0 {
            return Err(HyperbolaConstructionError::NonPositiveRadius);
        }
        Ok(Hyperbola3D {
            frame,
            major_radius,
            minor_radius,
        })
    }

    /// Creates a hyperbola from a center, a major-axis end point `s1`, and a
    /// second point `s2` that fixes the minor axis and the plane.
    ///
    /// The major axis direction is `x_axis = normalize(s1 - center)`, with
    /// semi-major radius `major_radius = |s1 - center|`. The semi-minor
    /// radius `minor_radius` is the distance from `s2` to the line
    /// `(center, x_axis)`. The plane normal is
    /// `normalize(x_axis × (s2 - center))`, and the resulting frame is
    /// `Frame3D::new(center, normal, x_axis)`.
    ///
    /// # Errors
    ///
    /// Returns [`HyperbolaConstructionError::ConfusedPoints`] if any pair
    /// among `{s1, s2, center}` is within `tol::CONFUSION`; or
    /// [`HyperbolaConstructionError::CollinearPoints`] if `s2` lies on the
    /// line `(center, x_axis)`, or `x_axis` is parallel to `s2 - center`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D};
    /// let hyperbola = Hyperbola3D::from_center_and_points(
    ///     Point3D::ORIGIN,
    ///     Point3D::new(2.0, 0.0, 0.0),
    ///     Point3D::new(0.0, 1.0, 0.0),
    /// )
    /// .unwrap();
    /// assert_eq!(hyperbola.major_radius(), 2.0);
    /// assert_eq!(hyperbola.minor_radius(), 1.0);
    /// ```
    pub fn from_center_and_points(
        center: Point3D,
        s1: Point3D,
        s2: Point3D,
    ) -> Result<Hyperbola3D, HyperbolaConstructionError> {
        if center.distance(s1) < tol::CONFUSION
            || center.distance(s2) < tol::CONFUSION
            || s1.distance(s2) < tol::CONFUSION
        {
            return Err(HyperbolaConstructionError::ConfusedPoints);
        }

        let v1 = s1 - center;
        let d1 = v1.magnitude();
        let x_axis = v1 * (1.0 / d1);

        let v2 = s2 - center;
        let proj = v2.dot(x_axis);
        let perp = v2 - proj * x_axis;
        let d2 = perp.magnitude();

        if d2 < tol::CONFUSION {
            return Err(HyperbolaConstructionError::CollinearPoints);
        }

        let normal = x_axis.cross(v2);
        if normal.magnitude() < tol::CONFUSION {
            return Err(HyperbolaConstructionError::CollinearPoints);
        }

        let frame = Frame3D::new(center, normal, x_axis).map_err(frame_error(
            HyperbolaConstructionError::NonFinite,
            HyperbolaConstructionError::CollinearPoints,
        ))?;
        Hyperbola3D::from_frame(frame, d1, d2)
    }

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

    /// Returns the hyperbola's plane frame.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Frame3D};
    /// let hyperbola = Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.frame(), Frame3D::WORLD);
    /// ```
    pub fn frame(&self) -> Frame3D {
        self.frame
    }

    /// Returns the semi-major radius of the hyperbola.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Frame3D};
    /// let hyperbola = Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.major_radius(), 2.0);
    /// ```
    pub fn major_radius(&self) -> f64 {
        self.major_radius
    }

    /// Returns the semi-minor radius of the hyperbola.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Frame3D};
    /// let hyperbola = Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.minor_radius(), 1.0);
    /// ```
    pub fn minor_radius(&self) -> f64 {
        self.minor_radius
    }

    /// Returns the image of this hyperbola under `transform`: the center
    /// and frame directions map point-wise while both radii scale by the
    /// transform's uniform scale factor.
    ///
    /// A mirror flips handedness, so the image frame is re-derived
    /// right-handed: the point set is preserved but the transverse
    /// direction reverses. Fails with
    /// [`HyperbolaConstructionError::NullNormal`] only when the transform
    /// collapses the frame (a zero scale factor).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Transform, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let moved = hyperbola.transformed(Transform::translation(Vector3D::new(0.0, 3.0, 0.0))).unwrap();
    /// assert_eq!(moved.center(), Point3D::new(0.0, 3.0, 0.0));
    /// assert_eq!(moved.major_radius(), 2.0);
    /// ```
    pub fn transformed(
        &self,
        transform: Transform,
    ) -> Result<Hyperbola3D, HyperbolaConstructionError> {
        let (frame, scale) = transform.apply_frame(self.frame()).map_err(frame_error(
            HyperbolaConstructionError::NonFinite,
            HyperbolaConstructionError::NullNormal,
        ))?;
        Hyperbola3D::from_frame(
            frame,
            self.major_radius() * scale,
            self.minor_radius() * scale,
        )
    }

    /// Evaluates the point on the hyperbola at parameter `u`:
    /// `center + major*cosh(u)*x_dir + minor*sinh(u)*y_dir`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.eval_point(0.0), Point3D::new(2.0, 0.0, 0.0));
    /// ```
    pub fn eval_point(&self, u: f64) -> Point3D {
        analytic::hyperbola_d0(&self.frame, self.major_radius, self.minor_radius, u)
    }

    /// Evaluates the points on the hyperbola at each parameter in `us`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let points = hyperbola.eval_points(&[0.0, 1.0]);
    /// assert_eq!(points[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`.
    ///
    /// Unlike the trigonometric conics, hyperbolic derivatives do not cycle
    /// with period 4 (see `curve_math::analytic::hyperbola_dn`):
    /// odd orders follow the `sinh`/`cosh` pattern of the first derivative,
    /// even orders follow the `cosh`/`sinh` pattern of the second.
    ///
    /// # Panics
    ///
    /// Panics if `order == 0`; use [`Hyperbola3D::eval_point`] to evaluate
    /// the position itself.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.eval_derivative(0.0, 1), Vector3D::new(0.0, 1.0, 0.0));
    /// ```
    pub fn eval_derivative(&self, u: f64, order: u32) -> Vector3D {
        match order {
            0 => panic!("eval_derivative: order must be >= 1 (use eval_point for order 0)"),
            _ => {
                analytic::hyperbola_dn(&self.frame, self.major_radius, self.minor_radius, u, order)
            }
        }
    }

    /// Recovers the parameter of a point on (or near) the hyperbola:
    /// `asinh(((point - center) . y_dir) / minor_radius)`. Unbounded (no
    /// wrapping).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// assert!((hyperbola.parameter_of(Point3D::new(2.0, 0.0, 0.0)) - 0.0).abs() < 1e-9);
    /// ```
    pub fn parameter_of(&self, point: Point3D) -> f64 {
        analytic::hyperbola_parameter(&self.frame, self.minor_radius, point)
    }

    /// Returns whether `point` lies on the hyperbola: within `tol.confusion`
    /// of it, exactly when [`Hyperbola3D::project_point`] reports distance `0.0`.
    /// Its implicit equation decides most points at once: off the plane,
    /// clearly far from the curve or clearly on it, with rigorous bounds
    /// for a quadratic equation. Of the points left, within a few
    /// tolerances of the curve, the inverse parameter of
    /// [`Hyperbola3D::parameter_of`] (exact on the curve, cheap) accepts most; the
    /// projection settles the rest.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// assert!(hyperbola.contains(Point3D::new(2.0, 0.0, 0.0), tol));
    /// assert!(!hyperbola.contains(Point3D::new(2.0, 0.0, 5.0), tol));
    /// ```
    pub fn contains(&self, point: Point3D, tol: Tolerance) -> bool {
        let (a, b) = (self.major_radius, self.minor_radius);
        // `(x/a)² - (y/b)² = 1`: on it `|∇f| ≥ 2/a`. Whenever the bounds
        // apply, `a ≥ 2 tol`, so the branch (`x ≥ a`) is the only part of
        // the curve within `2 tol` of a point with `x ≥ a - tol`, and points
        // short of that are far from it.
        let near = analytic::conic_proximity(
            &self.frame,
            point,
            tol.confusion,
            2.0 / (a.min(b) * a.min(b)),
            2.0 / a,
            |x, y| {
                if x < a - tol.confusion {
                    return (f64::INFINITY, 1.0, 0.0);
                }
                let (sx, sy) = (x / (a * a), y / (b * b));
                (x * sx - y * sy - 1.0, 2.0 * sx, -2.0 * sy)
            },
        );
        if let Some(near) = near {
            return near;
        }
        let u = self.parameter_of(point);
        self.eval_point(u).distance(point) <= tol.confusion
            || self.project_point(point, tol).distance <= tol.confusion
    }

    /// All stationary points of the distance from `point` to the (single)
    /// hyperbola branch, ordered by ascending distance.
    ///
    /// The stationarity condition becomes a quartic in `u = e^t` (solved
    /// by `real_roots`; only positive roots give real parameters). Each
    /// root is then polished by Newton on the stationarity function in `t`,
    /// which the quartic's conditioning alone cannot guarantee. The
    /// first entry is the global closest point (see
    /// [`Hyperbola3D::project_point`]).
    ///
    /// A query point with a non-finite coordinate has no extrema.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let extrema = hyperbola.extrema(Point3D::new(2.0, 0.0, 5.0), Tolerance::DEFAULT);
    /// assert!(!extrema.is_empty());
    /// ```
    pub fn extrema(&self, point: Point3D, tol: Tolerance) -> Vec<CurveProjection> {
        if ![point.x, point.y, point.z].iter().all(|c| c.is_finite()) {
            return Vec::new();
        }
        let rel = point - self.center();
        let px = rel.dot(self.frame.x_direction());
        let py = rel.dot(self.frame.y_direction());
        let (a, b) = (self.major_radius, self.minor_radius);
        let s2 = a * a + b * b;
        if s2 == 0.0 {
            // Both radii zero: the branch is its centre, every parameter
            // equally close.
            let distance = self.center().distance(point);
            return vec![CurveProjection {
                parameter: 0.0,
                distance: projection::snap_distance(distance, tol.confusion),
            }];
        }
        let coeffs = [
            -s2,
            2.0 * (a * px - b * py),
            0.0,
            -2.0 * (a * px + b * py),
            s2,
        ];
        // The stationarity condition (H(t) - P).H'(t) in frame coordinates,
        // and its derivative: candidates from the quartic are polished on
        // it, since near the vertex a complex root pair crowds the real
        // root and the polynomial alone resolves it only to ~sqrt(eps).
        let g = |t: f64| s2 * t.sinh() * t.cosh() - a * px * t.sinh() - b * py * t.cosh();
        let dg = |t: f64| s2 * (2.0 * t).cosh() - a * px * t.cosh() - b * py * t.sinh();
        let mut out: Vec<CurveProjection> = real_roots(&coeffs, tol.confusion)
            .into_iter()
            .filter(|r| r.value > 0.0)
            .map(|r| {
                let t = polish_root(g, dg, r.value.ln(), 8);
                let distance = self.eval_point(t).distance(point);
                CurveProjection {
                    parameter: t,
                    distance: projection::snap_distance(distance, tol.confusion),
                }
            })
            .collect();
        out.sort_by(|x, y| x.distance.total_cmp(&y.distance));
        out
    }

    /// Projects `point` onto the hyperbola: the nearest of
    /// [`Hyperbola3D::extrema`]. Distances within `tol.confusion` snap to
    /// `0.0`, matching [`Hyperbola3D::contains`].
    ///
    /// A query point with a non-finite coordinate projects to a NaN
    /// parameter and distance.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Hyperbola3D, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let proj = hyperbola.project_point(Point3D::new(3.0, 0.0, 0.0), Tolerance::DEFAULT);
    /// assert!((proj.parameter.abs() - 0.62236).abs() < 1e-4);
    /// assert!((proj.distance - 0.8f64.sqrt()).abs() < 1e-9);
    /// ```
    pub fn project_point(&self, point: Point3D, tol: Tolerance) -> CurveProjection {
        self.extrema(point, tol)
            .into_iter()
            .next()
            .unwrap_or(CurveProjection::UNDEFINED)
    }

    /// Projects each point in `points` onto the hyperbola.
    ///
    /// Default-style batch wrapper over [`Hyperbola3D::project_point`]:
    /// one native call per batch, mirroring [`Hyperbola3D::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) -> ConicSurfaceIntersection {
        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) -> ConicSurfaceIntersection {
        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,
    ) -> ConicSurfaceIntersection {
        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) -> ConicSurfaceIntersection {
        self.intersect(cone, tol)
    }

    /// Computes the exact 2D representation of this hyperbola in a surface's
    /// parameter space.
    ///
    /// No hyperbola/surface pair has a closed-form 2D image in this release,
    /// so this always returns [`ParametrizeError::NotAnalytic`]. The `surface`
    /// and `tol` arguments are accepted for signature parity with the analytic
    /// [`crate::Line3D::parametrize_on`] and
    /// [`crate::Circle3D::parametrize_on`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::curves::ParametrizeError;
    /// use geomcore::{Hyperbola3D, Plane, Point3D, Tolerance, Vector3D};
    ///
    /// let plane = Plane::new(Point3D::ORIGIN, Vector3D::Z).unwrap();
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// assert_eq!(hyperbola.parametrize_on(&plane, Tolerance::DEFAULT), Err(ParametrizeError::NotAnalytic));
    /// ```
    pub fn parametrize_on(
        &self,
        surface: impl Into<Surface>,
        _tol: Tolerance,
    ) -> Result<Curve2D, ParametrizeError> {
        let _ = surface.into();
        Err(ParametrizeError::NotAnalytic)
    }
}

impl Intersect<Plane> for Hyperbola3D {
    type Output = ConicSurfaceIntersection;

    /// Intersects this hyperbola with a plane: quadratic in `u = e^t`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{ConicSurfaceIntersection, Hyperbola3D, Plane, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let plane = Plane::new(Point3D::new(2.0, 0.0, 0.0), Vector3D::X).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match hyperbola.intersect(&plane, tol) {
    ///     ConicSurfaceIntersection::Hits(hits) => assert_eq!(hits.len(), 1),
    ///     _ => panic!("expected hits"),
    /// }
    /// ```
    fn intersect(&self, plane: &Plane, tol: Tolerance) -> ConicSurfaceIntersection {
        intersect_conic_quadric(self, plane, tol)
    }
}

impl Intersect<Sphere> for Hyperbola3D {
    type Output = ConicSurfaceIntersection;

    /// Intersects this hyperbola with a sphere; hits verified on the
    /// sphere. See [`Hyperbola3D::intersect_plane`] for the calling shape.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{ConicSurfaceIntersection, Hyperbola3D, Point3D, Sphere, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let sphere = Sphere::new(Point3D::ORIGIN, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match hyperbola.intersect(&sphere, tol) {
    ///     ConicSurfaceIntersection::Hits(hits) => assert!(!hits.is_empty()),
    ///     _ => panic!("expected hits"),
    /// }
    /// ```
    fn intersect(&self, sphere: &Sphere, tol: Tolerance) -> ConicSurfaceIntersection {
        intersect_conic_quadric(self, sphere, tol)
    }
}

impl Intersect<Cylinder> for Hyperbola3D {
    type Output = ConicSurfaceIntersection;

    /// Intersects this hyperbola with a cylinder; hits verified on the
    /// cylinder. See
    /// [`Hyperbola3D::intersect_plane`] for the calling shape.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{ConicSurfaceIntersection, Cylinder, Hyperbola3D, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match hyperbola.intersect(&cylinder, tol) {
    ///     ConicSurfaceIntersection::Hits(hits) => assert!(!hits.is_empty()),
    ///     _ => panic!("expected hits"),
    /// }
    /// ```
    fn intersect(&self, cylinder: &Cylinder, tol: Tolerance) -> ConicSurfaceIntersection {
        intersect_conic_quadric(self, cylinder, tol)
    }
}

impl Intersect<Cone> for Hyperbola3D {
    type Output = ConicSurfaceIntersection;

    /// Intersects this hyperbola with a cone; hits verified on
    /// either nappe.
    /// See [`Hyperbola3D::intersect_plane`] for the calling shape.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::Intersect;
    /// use geomcore::{Cone, ConicSurfaceIntersection, Frame3D, Hyperbola3D, Point3D, Tolerance, Vector3D};
    /// let hyperbola = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match hyperbola.intersect(&cone, tol) {
    ///     ConicSurfaceIntersection::Hits(hits) => assert!(!hits.is_empty()),
    ///     _ => panic!("expected hits"),
    /// }
    /// ```
    fn intersect(&self, cone: &Cone, tol: Tolerance) -> ConicSurfaceIntersection {
        intersect_conic_quadric(self, cone, tol)
    }
}

#[cfg(test)]
mod tests {
    use crate::Intersect;
    use crate::{
        Cone, ConicSurfaceIntersection, Frame3D, Hyperbola3D, HyperbolaConstructionError, Point3D,
        Tolerance, Vector3D, intersect_curve_surface,
    };

    // ---- construction ----

    #[test]
    fn test_hyperbola3d_new_ok() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        assert_eq!(h.center(), Point3D::ORIGIN);
        assert_eq!(h.major_radius(), 2.0);
        assert_eq!(h.minor_radius(), 1.0);
    }

    #[test]
    fn test_hyperbola3d_new_negative_major_radius_errors() {
        assert_eq!(
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, -2.0, 1.0),
            Err(HyperbolaConstructionError::NonPositiveRadius)
        );
    }

    #[test]
    fn test_hyperbola3d_new_negative_minor_radius_errors() {
        assert_eq!(
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, -1.0),
            Err(HyperbolaConstructionError::NonPositiveRadius)
        );
    }

    #[test]
    fn test_hyperbola3d_new_minor_greater_than_major_allowed() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 1.0, 2.0).unwrap();
        assert_eq!(h.major_radius(), 1.0);
        assert_eq!(h.minor_radius(), 2.0);
    }

    #[test]
    fn test_hyperbola3d_new_null_normal_errors() {
        assert_eq!(
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::ZERO, Vector3D::X, 2.0, 1.0),
            Err(HyperbolaConstructionError::NullNormal)
        );
    }

    #[test]
    fn test_hyperbola3d_new_parallel_x_direction_errors() {
        assert_eq!(
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::Z, 2.0, 1.0),
            Err(HyperbolaConstructionError::NullNormal)
        );
    }

    #[test]
    fn test_hyperbola3d_from_frame_ok() {
        let h = Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, 1.0).unwrap();
        assert_eq!(h.frame(), Frame3D::WORLD);
        assert_eq!(h.major_radius(), 2.0);
        assert_eq!(h.minor_radius(), 1.0);
    }

    #[test]
    fn test_hyperbola3d_from_frame_negative_major_radius_errors() {
        assert_eq!(
            Hyperbola3D::from_frame(Frame3D::WORLD, -2.0, 1.0),
            Err(HyperbolaConstructionError::NonPositiveRadius)
        );
    }

    #[test]
    fn test_hyperbola3d_from_frame_negative_minor_radius_errors() {
        assert_eq!(
            Hyperbola3D::from_frame(Frame3D::WORLD, 2.0, -1.0),
            Err(HyperbolaConstructionError::NonPositiveRadius)
        );
    }

    #[test]
    fn test_hyperbola3d_from_center_and_points_ok() {
        let h = Hyperbola3D::from_center_and_points(
            Point3D::ORIGIN,
            Point3D::new(2.0, 0.0, 0.0),
            Point3D::new(0.0, 1.0, 0.0),
        )
        .unwrap();
        assert_eq!(h.major_radius(), 2.0);
        assert_eq!(h.minor_radius(), 1.0);
    }

    #[test]
    fn test_hyperbola3d_from_center_and_points_confused_points_errors() {
        assert_eq!(
            Hyperbola3D::from_center_and_points(
                Point3D::ORIGIN,
                Point3D::ORIGIN,
                Point3D::new(0.0, 1.0, 0.0),
            ),
            Err(HyperbolaConstructionError::ConfusedPoints)
        );
    }

    #[test]
    fn test_hyperbola3d_from_center_and_points_collinear_points_errors() {
        assert_eq!(
            Hyperbola3D::from_center_and_points(
                Point3D::ORIGIN,
                Point3D::new(2.0, 0.0, 0.0),
                Point3D::new(1.0, 0.0, 0.0),
            ),
            Err(HyperbolaConstructionError::CollinearPoints)
        );
    }

    // ---- evaluation ----

    #[test]
    fn test_hyperbola3d_eval_point_zero() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let p = h.eval_point(0.0);
        assert!((p.x - 2.0).abs() < 1e-9);
        assert!(p.y.abs() < 1e-9);
        assert!(p.z.abs() < 1e-9);
    }

    #[test]
    fn test_hyperbola3d_eval_points_matches_loop() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let us = [0.0, 0.5, 1.5];
        let expected: Vec<Point3D> = us.iter().map(|&u| h.eval_point(u)).collect();
        assert_eq!(h.eval_points(&us), expected);
    }

    #[test]
    fn test_hyperbola3d_eval_derivative_order1() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let d1 = h.eval_derivative(0.0, 1);
        assert!(d1.x.abs() < 1e-9);
        assert!((d1.y - 1.0).abs() < 1e-9);
    }

    #[test]
    #[should_panic]
    fn test_hyperbola3d_eval_derivative_order0_panics() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        h.eval_derivative(0.0, 0);
    }

    #[test]
    fn test_hyperbola3d_parameter_of_round_trip() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        for u in [0.3, 2.0, -1.5] {
            let p = h.eval_point(u);
            assert!((h.parameter_of(p) - u).abs() < 1e-9);
        }
    }

    #[test]
    fn test_hyperbola3d_parameter_of_is_unbounded() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let p = h.eval_point(5.0);
        assert!((h.parameter_of(p) - 5.0).abs() < 1e-6);
    }

    // ---- HyperbolaConstructionError ----

    #[test]
    fn test_hyperbola_construction_error_display() {
        assert_eq!(
            HyperbolaConstructionError::NonPositiveRadius.to_string(),
            "radius is not positive"
        );
        assert_eq!(
            HyperbolaConstructionError::NullNormal.to_string(),
            "normal has zero length"
        );
        assert_eq!(
            HyperbolaConstructionError::ConfusedPoints.to_string(),
            "the points are confused"
        );
        assert_eq!(
            HyperbolaConstructionError::CollinearPoints.to_string(),
            "the points are collinear"
        );
    }

    #[test]
    fn test_hyperbola_construction_error_is_std_error() {
        fn takes_error(_e: &dyn std::error::Error) {}
        takes_error(&HyperbolaConstructionError::NonPositiveRadius);
    }

    #[test]
    fn test_hyperbola3d_contains() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let tol = Tolerance::DEFAULT;
        assert!(h.contains(Point3D::new(2.0, 0.0, 0.0), tol));
        assert!(h.contains(h.eval_point(1.0), tol));
        assert!(!h.contains(Point3D::ORIGIN, tol));
        assert!(!h.contains(Point3D::new(2.0, 0.0, 5.0), tol));
    }

    #[test]
    fn test_hyperbola3d_project_point() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let tol = Tolerance::DEFAULT;
        // On the transverse axis right of the vertex: the true minima are
        // a symmetric pair off the vertex, closer than the vertex itself.
        let proj = h.project_point(Point3D::new(3.0, 0.0, 0.0), tol);
        assert!((proj.parameter.abs() - 0.62236).abs() < 1e-4);
        assert!((proj.distance - 0.8f64.sqrt()).abs() < 1e-9);
        // Off-curve: verify against brute-force sampling (which can only
        // overestimate the true minimum).
        let query = Point3D::new(2.0, 0.0, 5.0);
        let off = h.project_point(query, tol);
        let mut best = f64::INFINITY;
        for i in -1000..=1000 {
            let d = h.eval_point(i as f64 * 0.01).distance(query);
            best = best.min(d);
        }
        assert!(off.distance <= best);
        assert!(best - off.distance < 1e-3);
    }

    #[test]
    fn test_hyperbola3d_project_point_polished() {
        // Thin hyperbola: the raw e^t quartic root left on-curve points
        // 3.4e-6 away before Newton polishing.
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 0.15, 5.7).unwrap();
        let proj = h.project_point(h.eval_point(0.0007), Tolerance::DEFAULT);
        assert!(proj.distance < 1e-12, "{}", proj.distance);
    }

    #[test]
    fn test_hyperbola3d_intersect_cone_crosscheck() {
        let tol = Tolerance::DEFAULT;
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        match h.intersect(&cone, tol) {
            ConicSurfaceIntersection::Hits(hits) => {
                assert!(!hits.is_empty());
                for h in &hits {
                    assert!(cone.contains(h.point, tol));
                }
                let generic = intersect_curve_surface(&h, &cone, tol);
                assert_eq!(generic.len(), hits.len());
            }
            _ => panic!("expected hits"),
        }
    }
}

#[cfg(test)]
mod projection_input_tests {
    use super::*;

    #[test]
    fn test_hyperbola3d_project_non_finite_point_is_nan() {
        let h = Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let tol = Tolerance::DEFAULT;
        for p in [
            Point3D::new(f64::NAN, 0.0, 0.0),
            Point3D::new(f64::INFINITY, 0.0, 0.0),
        ] {
            assert!(h.extrema(p, tol).is_empty());
            let proj = h.project_point(p, tol);
            assert!(proj.parameter.is_nan() && proj.distance.is_nan());
        }
    }

    #[test]
    fn test_hyperbola3d_zero_radii_rejected() {
        let center = Point3D::new(1.0, 2.0, 3.0);
        for (a, b) in [(0.0, 0.0), (2.0, 0.0), (0.0, 1.0)] {
            assert_eq!(
                Hyperbola3D::new(center, Vector3D::Z, Vector3D::X, a, b),
                Err(HyperbolaConstructionError::NonPositiveRadius)
            );
        }
    }
}

#[cfg(test)]
mod transformed_tests {
    use super::*;
    use crate::{Axis3D, Point3D, Tolerance, Transform, Vector3D};

    #[test]
    fn test_hyperbola3d_transformed_translation_spot_check() {
        let hyperbola =
            Hyperbola3D::new(Point3D::ORIGIN, Vector3D::Z, Vector3D::X, 2.0, 1.0).unwrap();
        let t = Transform::translation(Vector3D::new(0.0, 0.0, 7.0));
        let moved = hyperbola.transformed(t).unwrap();
        assert_eq!(moved.center(), Point3D::new(0.0, 0.0, 7.0));
        assert_eq!(moved.major_radius(), 2.0);
        assert_eq!(moved.minor_radius(), 1.0);
        // Hand-computed: center + 2*cosh(0)*X + sinh(0)*Y, shifted.
        assert_eq!(moved.eval_point(0.0), Point3D::new(2.0, 0.0, 7.0));
    }

    #[test]
    fn test_hyperbola3d_transformed_round_trip_and_mirror() {
        let hyperbola = Hyperbola3D::new(
            Point3D::new(1.0, 0.0, 0.0),
            Vector3D::Z,
            Vector3D::X,
            2.0,
            1.0,
        )
        .unwrap();
        let axis = Axis3D::new(Point3D::ORIGIN, Vector3D::Y).unwrap();
        let t = Transform::translation(Vector3D::new(0.0, 4.0, 0.0))
            .then(Transform::rotation(axis, std::f64::consts::FRAC_PI_2));
        let moved = hyperbola.transformed(t).unwrap();
        let tol = Tolerance::DEFAULT;
        for u in [0.0, 0.5, -1.0] {
            assert!(moved.contains(t.apply_point(hyperbola.eval_point(u)), tol));
        }
        let inv = Transform::rotation(axis, -std::f64::consts::FRAC_PI_2)
            .then(Transform::translation(Vector3D::new(0.0, -4.0, 0.0)));
        let back = moved.transformed(inv).unwrap();
        for u in [0.0, 0.5, -1.0] {
            assert!(back.eval_point(u).distance(hyperbola.eval_point(u)) < 1e-9);
        }
        // Mirror across the xy-plane preserves the point set.
        let m = Transform::mirror_plane(Frame3D::WORLD);
        let mirrored = hyperbola.transformed(m).unwrap();
        for u in [0.0, 0.5, -1.0] {
            assert!(mirrored.contains(m.apply_point(hyperbola.eval_point(u)), tol));
        }
        assert_eq!(mirrored.transformed(m).unwrap(), hyperbola);
    }
}