geomcore 0.19.0

//! Cones in 3D: parametric evaluation, parameter inversion, and
//! construction from an apex angle or two circular sections, thin wrappers
//! over [`crate::surface_math::analytic`].

use crate::projection::{self, SurfaceProjection};
use crate::surface_math::analytic;
use crate::tol;
use crate::{
    Circle3D, ConeConeIntersection, ConeCylinderIntersection, Cylinder, Frame3D, Point3D,
    Tolerance, Vector3D,
};
use std::fmt;

/// Error returned when a [`Cone`] cannot be constructed from the given
/// inputs.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ConeConstructionError {
    /// The requested reference radius is negative.
    NegativeRadius,
    /// The semi-angle is out of the accepted range (see
    /// [`Cone::from_frame`]).
    BadAngle,
    /// The two points used to derive the axis are coincident (or too close
    /// to distinguish).
    NullAxis,
    /// The derived semi-angle is (numerically) zero or a right angle.
    NullAngle,
    /// The axis direction has zero length.
    NullNormal,
}

impl fmt::Display for ConeConstructionError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let message = match self {
            ConeConstructionError::NegativeRadius => "radius is negative",
            ConeConstructionError::BadAngle => "semi-angle is out of range",
            ConeConstructionError::NullAxis => "the two points are confused",
            ConeConstructionError::NullAngle => "semi-angle is null or a right angle",
            ConeConstructionError::NullNormal => "axis direction has zero length",
        };
        f.write_str(message)
    }
}

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

/// A cone in 3D: a [`Frame3D`] (origin plus local x/y/z directions, `z`
/// being the axis), a semi-angle, and a reference radius (the radius of the
/// circular section through the origin), evaluated as
/// `R = ref_radius + v*sin(semi_angle)`,
/// `origin + R*cos(u)*x_dir + R*sin(u)*y_dir + v*cos(semi_angle)*z_dir`.
///
/// # Examples
///
/// ```
/// use geomcore::{Cone, Point3D, Vector3D};
/// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
/// assert_eq!(cone.eval_point(0.0, 0.0), Point3D::new(2.0, 0.0, 0.0));
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Cone {
    frame: Frame3D,
    semi_angle: f64,
    ref_radius: f64,
}

impl Cone {
    /// Creates a cone from a frame, a semi-angle, and a reference radius.
    ///
    /// The semi-angle must satisfy `0 < |semi_angle| < PI/2` (both bounds
    /// excluded with `tol::ANGULAR` tolerance); a negative
    /// semi-angle flips the direction the cone opens in along `z`.
    ///
    /// # Errors
    ///
    /// Returns [`ConeConstructionError::BadAngle`] if `semi_angle` is
    /// outside the accepted range, or
    /// [`ConeConstructionError::NegativeRadius`] if `ref_radius < 0`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Frame3D};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// assert_eq!(cone.ref_radius(), 2.0);
    /// ```
    pub fn from_frame(
        frame: Frame3D,
        semi_angle: f64,
        ref_radius: f64,
    ) -> Result<Cone, ConeConstructionError> {
        if !is_valid_semi_angle(semi_angle) {
            return Err(ConeConstructionError::BadAngle);
        }
        if ref_radius < 0.0 {
            return Err(ConeConstructionError::NegativeRadius);
        }
        Ok(Cone {
            frame,
            semi_angle,
            ref_radius,
        })
    }

    /// Creates a cone from a center, an axis direction, a semi-angle, and a
    /// reference radius.
    ///
    /// The frame is derived from `axis_direction` via [`Frame3D::from_z`].
    ///
    /// # Errors
    ///
    /// Returns [`ConeConstructionError::NullNormal`] if `axis_direction`
    /// cannot be normalized (zero length), or any error from
    /// [`Cone::from_frame`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
    /// assert_eq!(cone.semi_angle(), 0.4);
    /// ```
    pub fn new(
        center: Point3D,
        axis_direction: Vector3D,
        semi_angle: f64,
        ref_radius: f64,
    ) -> Result<Cone, ConeConstructionError> {
        let frame = Frame3D::from_z(center, axis_direction)
            .map_err(|_| ConeConstructionError::NullNormal)?;
        Cone::from_frame(frame, semi_angle, ref_radius)
    }

    /// Creates a cone through two circular sections: `p1` with radius `r1`
    /// and `p2` with radius `r2`.
    ///
    /// The axis runs from `p1` to `p2`; the reference radius is `r1`
    /// (the section at `p1`) and the semi-angle is
    /// `atan((r2 - r1) / dist(p1, p2))` (negative when the cone narrows
    /// from `p1` to `p2`). The frame's origin is `p1`, its z direction is
    /// `normalize(p2 - p1)`, and its x direction is an arbitrary vector
    /// perpendicular to `z`, matching the reference implementation's
    /// convention (see `arbitrary_perpendicular_frame`).
    ///
    /// # Errors
    ///
    /// Returns [`ConeConstructionError::NullAxis`] if `p1` and `p2` are
    /// coincident (or too close to distinguish),
    /// [`ConeConstructionError::NegativeRadius`] if `r1 < 0 || r2 < 0`, or
    /// [`ConeConstructionError::NullAngle`] if the derived semi-angle is
    /// (numerically) zero or a right angle.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D};
    /// let cone = Cone::from_two_points_and_radii(
    ///     Point3D::ORIGIN,
    ///     Point3D::new(0.0, 0.0, 4.0),
    ///     2.0,
    ///     1.0,
    /// )
    /// .unwrap();
    /// assert_eq!(cone.ref_radius(), 2.0);
    /// ```
    pub fn from_two_points_and_radii(
        p1: Point3D,
        p2: Point3D,
        r1: f64,
        r2: f64,
    ) -> Result<Cone, ConeConstructionError> {
        let delta = p2 - p1;
        let dist = delta.magnitude();
        if dist < tol::CONFUSION {
            return Err(ConeConstructionError::NullAxis);
        }
        if r1 < 0.0 || r2 < 0.0 {
            return Err(ConeConstructionError::NegativeRadius);
        }

        let semi_angle = ((r2 - r1) / dist).atan();
        if !is_valid_semi_angle(semi_angle) {
            return Err(ConeConstructionError::NullAngle);
        }

        let z_dir = delta * (1.0 / dist);
        let frame = arbitrary_perpendicular_frame(p1, z_dir);
        Ok(Cone {
            frame,
            semi_angle,
            ref_radius: r1,
        })
    }

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

    /// Returns the cone's semi-angle.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Frame3D};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// assert_eq!(cone.semi_angle(), 0.4);
    /// ```
    pub fn semi_angle(&self) -> f64 {
        self.semi_angle
    }

    /// Returns the cone's reference radius (the radius of the circular
    /// section through the origin, at `v = 0`).
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Frame3D};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// assert_eq!(cone.ref_radius(), 2.0);
    /// ```
    pub fn ref_radius(&self) -> f64 {
        self.ref_radius
    }

    /// Returns the apex of the cone: `origin + (-ref_radius /
    /// tan(semi_angle)) * z_dir`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, std::f64::consts::FRAC_PI_4, 1.0).unwrap();
    /// let apex = cone.apex();
    /// assert!((apex.z + 1.0).abs() < 1e-9);
    /// ```
    pub fn apex(&self) -> Point3D {
        let v = -self.ref_radius / self.semi_angle.tan();
        self.frame.point_at(0.0, 0.0, v)
    }

    /// Evaluates the point on the cone at `(u, v)`. See the type-level docs
    /// for the formula.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
    /// assert_eq!(cone.eval_point(0.0, 0.0), Point3D::new(2.0, 0.0, 0.0));
    /// ```
    pub fn eval_point(&self, u: f64, v: f64) -> Point3D {
        analytic::cone_d0(&self.frame, self.ref_radius, self.semi_angle, u, v)
    }

    /// Evaluates the points on the cone at each `(u, v)` in `uvs`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
    /// let points = cone.eval_points(&[(0.0, 0.0)]);
    /// assert_eq!(points[0], Point3D::new(2.0, 0.0, 0.0));
    /// ```
    pub fn eval_points(&self, uvs: &[(f64, f64)]) -> Vec<Point3D> {
        uvs.iter().map(|&(u, v)| self.eval_point(u, v)).collect()
    }

    /// Evaluates the derivative of order `(du, dv)` at `(u, v)`. See
    /// `surface_math::analytic::cone_derivative` for the formulas.
    ///
    /// # Panics
    ///
    /// Panics if `du + dv == 0` (use [`Cone::eval_point`] for the position
    /// itself) or if `du + dv > 2`.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
    /// let d1v = cone.eval_derivative(0.0, 0.0, 0, 1);
    /// assert!(d1v.x > 0.0);
    /// ```
    pub fn eval_derivative(&self, u: f64, v: f64, du: u32, dv: u32) -> Vector3D {
        match du + dv {
            0 => panic!(
                "eval_derivative: du + dv must be >= 1 (use eval_point for the (0, 0) order)"
            ),
            1..=2 => analytic::cone_derivative(
                &self.frame,
                self.ref_radius,
                self.semi_angle,
                u,
                v,
                du,
                dv,
            ),
            _ => panic!(
                "eval_derivative: order du={du}, dv={dv} is not supported (du + dv must be <= 2)"
            ),
        }
    }

    /// Recovers `(u, v)` of a point on (or near) the cone. See
    /// `surface_math::analytic::cone_parameters` for the formula,
    /// including the opposite-nappe handling below the apex.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Point3D, Vector3D};
    /// let cone = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
    /// let (u, v) = cone.parameters_of(Point3D::new(2.0, 0.0, 0.0));
    /// assert_eq!((u, v), (0.0, 0.0));
    /// ```
    pub fn parameters_of(&self, point: Point3D) -> (f64, f64) {
        analytic::cone_parameters(&self.frame, self.ref_radius, self.semi_angle, point)
    }

    /// Returns whether `point` lies on the cone: the inverse parameters are
    /// recovered with [`Cone::parameters_of`] and re-evaluated, and the
    /// point counts as contained when the re-evaluated point is within
    /// `tol.confusion` of it.
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Frame3D, Point3D, Tolerance};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// assert!(cone.contains(cone.eval_point(0.0, 0.0), tol));
    /// assert!(!cone.contains(Point3D::ORIGIN, tol));
    /// ```
    pub fn contains(&self, point: Point3D, tol: Tolerance) -> bool {
        let (u, v) = self.parameters_of(point);
        self.eval_point(u, v).distance(point) <= tol.confusion
    }

    /// Projects `point` onto the cone, returning the `(u, v)` parameters
    /// of the closest point and its distance. Distances within
    /// `tol.confusion` snap to `0.0`, matching [`Cone::contains`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, Frame3D, Point3D, Tolerance};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let p = cone.eval_point(1.0, 2.0);
    /// let proj = cone.project_point(p, Tolerance::DEFAULT);
    /// assert!((proj.u - 1.0).abs() < 1e-9);
    /// assert!((proj.v - 2.0).abs() < 1e-9);
    /// assert_eq!(proj.distance, 0.0);
    /// ```
    pub fn project_point(&self, point: Point3D, tol: Tolerance) -> SurfaceProjection {
        let (u, v) = self.parameters_of(point);
        let distance = self.eval_point(u, v).distance(point);
        SurfaceProjection {
            u,
            v,
            distance: projection::snap_distance(distance, tol.confusion),
        }
    }

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

    /// Intersects this cone with another cone.
    ///
    /// The analytic path needs coaxial axes: equal ring radii
    /// `alpha*tan(phi)` give one axial position
    /// `alpha* = -s*d*tan2 / (tan1 - s*tan2)` (`d` the signed apex
    /// separation, `s` the axis alignment), kept when both nappes cover
    /// it. Same apex, axis, and angle coincide; parallel distinct nappes
    /// and behind-apex roots miss. Anything else is a space quartic and
    /// reports [`ConeConeIntersection::NotAnalytic`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, ConeConeIntersection, Frame3D, Point3D, Tolerance, Vector3D};
    /// let c1 = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let c2 = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// assert_eq!(
    ///     c1.intersect_cone(&c2, tol),
    ///     ConeConeIntersection::Coincident
    /// );
    /// ```
    pub fn intersect_cone(&self, other: &Cone, tol: Tolerance) -> ConeConeIntersection {
        let a1 = self.frame().z_direction();
        let a2 = other.frame().z_direction();
        if a1.cross(a2).magnitude() > tol.angular {
            return ConeConeIntersection::NotAnalytic;
        }
        let w = other.apex() - self.apex();
        if (w - a1 * w.dot(a1)).magnitude() > tol.confusion {
            return ConeConeIntersection::NotAnalytic;
        }
        let s = a1.dot(a2);
        let d = w.dot(a1);
        let (t1, t2) = (self.semi_angle().tan(), other.semi_angle().tan());
        if d.abs() <= tol.confusion
            && (self.semi_angle() - other.semi_angle()).abs() <= tol.angular
            && s > 0.0
        {
            return ConeConeIntersection::Coincident;
        }
        let denom = t1 - s * t2;
        if denom.abs() <= tol.angular * (t1 + t2 + 1.0) {
            return ConeConeIntersection::Empty;
        }
        let alpha = -s * d * t2 / denom;
        let alpha2 = s * (alpha - d);
        if alpha < -tol.confusion || alpha2 < -tol.confusion {
            return ConeConeIntersection::Empty;
        }
        let center = self.apex() + a1 * alpha;
        let radius = (alpha * t1).max(0.0);
        if radius <= tol.confusion {
            ConeConeIntersection::ApexPoint(center)
        } else {
            let circle =
                Circle3D::new(center, a1, radius).expect("ring radius is positive by construction");
            ConeConeIntersection::Circle(circle)
        }
    }

    /// Intersects this cone with a cylinder.
    ///
    /// The analytic path needs coaxial axes: with the axis lines
    /// coincident (direction against `tol.angular`, separation against
    /// `tol.confusion`), the tube of radius `r` meets the nappe at axial
    /// `r/tan(semi_angle)`, always a single latitude circle — or just the
    /// apex for a degenerate (zero-radius) tube. Anything else is a space
    /// quartic and reports [`ConeCylinderIntersection::NotAnalytic`].
    ///
    /// # Examples
    ///
    /// ```
    /// use geomcore::{Cone, ConeCylinderIntersection, Cylinder, Frame3D, Point3D, Tolerance, Vector3D};
    /// let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
    /// let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
    /// let tol = Tolerance::DEFAULT;
    /// match cone.intersect_cylinder(&cylinder, tol) {
    ///     ConeCylinderIntersection::Circle(circle) => {
    ///         assert_eq!(circle.radius(), 2.0);
    ///         assert!(circle.center().distance(Point3D::ORIGIN) < 1e-12);
    ///     }
    ///     _ => panic!("expected a circle"),
    /// }
    /// ```
    pub fn intersect_cylinder(
        &self,
        cylinder: &Cylinder,
        tol: Tolerance,
    ) -> ConeCylinderIntersection {
        let a = self.frame().z_direction();
        let axis = cylinder.axis();
        if a.cross(axis.direction()).magnitude() > tol.angular {
            return ConeCylinderIntersection::NotAnalytic;
        }
        let w = axis.origin() - self.apex();
        if (w - a * w.dot(a)).magnitude() > tol.confusion {
            return ConeCylinderIntersection::NotAnalytic;
        }
        let r = cylinder.radius();
        if r <= tol.confusion {
            return ConeCylinderIntersection::ApexPoint(self.apex());
        }
        // Axial reach is always positive: exactly one latitude ring.
        let alpha = r / self.semi_angle().tan();
        let circle = Circle3D::new(self.apex() + a * alpha, a, r)
            .expect("tube radius is positive by construction");
        ConeCylinderIntersection::Circle(circle)
    }
}

/// Checks that `semi_angle` is within the accepted range for a cone:
/// `0 < |semi_angle| < PI/2`, both bounds excluded with
/// [`crate::tol::ANGULAR`] tolerance.
fn is_valid_semi_angle(semi_angle: f64) -> bool {
    let angle = semi_angle.abs();
    angle >= tol::ANGULAR && (std::f64::consts::FRAC_PI_2 - angle) >= tol::ANGULAR
}

/// Builds a frame at `origin` with `z_dir = normalize(axis_direction)` and
/// an arbitrary x direction perpendicular to it, matching the reference
/// implementation's convention for cones derived from two points and radii:
/// zero the smallest-magnitude component of `z_dir` and swap the other two.
/// This is the sign-mirror of the analogous helper in `plane.rs` (used for
/// `Plane::from_coefficients`); the two call paths were empirically found
/// (from `construction.json`'s golden cases) to use different underlying
/// "arbitrary perpendicular direction" conventions upstream.
///
/// `axis_direction` must already be a unit vector.
fn arbitrary_perpendicular_frame(origin: Point3D, axis_direction: Vector3D) -> Frame3D {
    let z = axis_direction;
    let (ax, ay, az) = (z.x.abs(), z.y.abs(), z.z.abs());
    let x_dir = if ax <= ay && ax <= az {
        Vector3D::new(0.0, -z.z, z.y)
    } else if ay <= ax && ay <= az {
        Vector3D::new(-z.z, 0.0, z.x)
    } else {
        Vector3D::new(-z.y, z.x, 0.0)
    }
    .normalized()
    .expect("z is a unit vector and nonzero on at least two axes, so the swap is nonzero");
    Frame3D::new(origin, z, x_dir).expect("x_dir constructed perpendicular to z by design")
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::Vector3D;
    use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};

    // ---- construction: from_frame ----

    #[test]
    fn test_from_frame_ok() {
        let c = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        assert_eq!(c.semi_angle(), 0.4);
        assert_eq!(c.ref_radius(), 2.0);
    }

    #[test]
    fn test_from_frame_negative_angle_ok() {
        let c = Cone::from_frame(Frame3D::WORLD, -0.4, 2.0).unwrap();
        assert_eq!(c.semi_angle(), -0.4);
    }

    #[test]
    fn test_from_frame_negative_radius_errors() {
        assert_eq!(
            Cone::from_frame(Frame3D::WORLD, 0.4, -1.0),
            Err(ConeConstructionError::NegativeRadius)
        );
    }

    #[test]
    fn test_from_frame_zero_angle_errors() {
        assert_eq!(
            Cone::from_frame(Frame3D::WORLD, 0.0, 1.0),
            Err(ConeConstructionError::BadAngle)
        );
    }

    #[test]
    fn test_from_frame_right_angle_errors() {
        assert_eq!(
            Cone::from_frame(Frame3D::WORLD, FRAC_PI_2, 1.0),
            Err(ConeConstructionError::BadAngle)
        );
        assert_eq!(
            Cone::from_frame(Frame3D::WORLD, -FRAC_PI_2, 1.0),
            Err(ConeConstructionError::BadAngle)
        );
    }

    #[test]
    fn test_from_frame_valid_boundary_ok() {
        assert!(Cone::from_frame(Frame3D::WORLD, FRAC_PI_4, 1.0).is_ok());
    }

    // ---- construction: new ----

    #[test]
    fn test_new_ok() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        assert_eq!(c.frame().z_direction(), Vector3D::Z);
    }

    #[test]
    fn test_new_null_normal_errors() {
        assert_eq!(
            Cone::new(Point3D::ORIGIN, Vector3D::ZERO, 0.4, 2.0),
            Err(ConeConstructionError::NullNormal)
        );
    }

    // ---- construction: from_two_points_and_radii ----

    #[test]
    fn test_from_two_points_and_radii_null_axis_errors() {
        let p = Point3D::new(1.0, 2.0, 3.0);
        assert_eq!(
            Cone::from_two_points_and_radii(p, p, 1.0, 2.0),
            Err(ConeConstructionError::NullAxis)
        );
    }

    #[test]
    fn test_from_two_points_and_radii_negative_radius_errors() {
        assert_eq!(
            Cone::from_two_points_and_radii(
                Point3D::ORIGIN,
                Point3D::new(0.0, 0.0, 1.0),
                -1.0,
                2.0
            ),
            Err(ConeConstructionError::NegativeRadius)
        );
        assert_eq!(
            Cone::from_two_points_and_radii(
                Point3D::ORIGIN,
                Point3D::new(0.0, 0.0, 1.0),
                1.0,
                -2.0
            ),
            Err(ConeConstructionError::NegativeRadius)
        );
    }

    #[test]
    fn test_from_two_points_and_radii_null_angle_errors_when_radii_equal() {
        // r1 == r2 gives a zero semi-angle: no valid cone (it's a cylinder).
        assert_eq!(
            Cone::from_two_points_and_radii(Point3D::ORIGIN, Point3D::new(0.0, 0.0, 1.0), 2.0, 2.0),
            Err(ConeConstructionError::NullAngle)
        );
    }

    #[test]
    fn test_from_two_points_and_radii_matches_golden_case_shrinking() {
        // From tests/fixtures/construction.json: cones_two_points_radii[0].
        let c =
            Cone::from_two_points_and_radii(Point3D::ORIGIN, Point3D::new(0.0, 0.0, 4.0), 2.0, 1.0)
                .unwrap();
        assert!((c.semi_angle() - (-0.24497866312686414)).abs() < 1e-9);
        assert_eq!(c.ref_radius(), 2.0);
        let frame = c.frame();
        assert_eq!(frame.origin(), Point3D::ORIGIN);
        assert_vector3_close(frame.x_direction(), Vector3D::new(0.0, -1.0, 0.0));
        assert_vector3_close(frame.y_direction(), Vector3D::new(1.0, 0.0, 0.0));
        assert_vector3_close(frame.z_direction(), Vector3D::Z);
    }

    #[test]
    fn test_from_two_points_and_radii_matches_golden_case_growing() {
        // From tests/fixtures/construction.json: cones_two_points_radii[1].
        let c = Cone::from_two_points_and_radii(
            Point3D::new(1.0, 1.0, 1.0),
            Point3D::new(2.0, 3.0, 1.0),
            1.0,
            2.5,
        )
        .unwrap();
        assert!((c.semi_angle() - 0.5908727501454191).abs() < 1e-9);
        assert_eq!(c.ref_radius(), 1.0);
        let frame = c.frame();
        assert_eq!(frame.origin(), Point3D::new(1.0, 1.0, 1.0));
        assert_vector3_close(
            frame.x_direction(),
            Vector3D::new(-0.894427190999916, 0.447213595499958, 0.0),
        );
        assert_vector3_close(frame.y_direction(), Vector3D::new(0.0, 0.0, 1.0));
        assert_vector3_close(
            frame.z_direction(),
            Vector3D::new(0.4472135954999579, 0.8944271909999159, 0.0),
        );
    }

    fn assert_vector3_close(actual: Vector3D, expected: Vector3D) {
        assert!(
            (actual.x - expected.x).abs() < 1e-9,
            "x: {actual:?} vs {expected:?}"
        );
        assert!(
            (actual.y - expected.y).abs() < 1e-9,
            "y: {actual:?} vs {expected:?}"
        );
        assert!(
            (actual.z - expected.z).abs() < 1e-9,
            "z: {actual:?} vs {expected:?}"
        );
    }

    // ---- evaluation ----

    #[test]
    fn test_eval_point() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        assert_eq!(c.eval_point(0.0, 0.0), Point3D::new(2.0, 0.0, 0.0));
    }

    #[test]
    fn test_eval_points_matches_loop() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        let uvs = [(0.0, 0.0), (0.5, 1.0)];
        let expected: Vec<Point3D> = uvs.iter().map(|&(u, v)| c.eval_point(u, v)).collect();
        assert_eq!(c.eval_points(&uvs), expected);
    }

    #[test]
    fn test_apex() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, FRAC_PI_4, 1.0).unwrap();
        let apex = c.apex();
        assert!(apex.x.abs() < 1e-9);
        assert!(apex.y.abs() < 1e-9);
        assert!((apex.z + 1.0).abs() < 1e-9);
    }

    #[test]
    #[should_panic(expected = "du + dv must be >= 1")]
    fn test_eval_derivative_zero_order_panics() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        c.eval_derivative(0.0, 0.0, 0, 0);
    }

    #[test]
    #[should_panic(expected = "du + dv must be <= 2")]
    fn test_eval_derivative_order_too_high_panics() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        c.eval_derivative(0.0, 0.0, 2, 1);
    }

    #[test]
    fn test_parameters_of_round_trip() {
        let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
        let (u, v) = c.parameters_of(Point3D::new(2.0, 0.0, 0.0));
        assert_eq!((u, v), (0.0, 0.0));
    }

    // ---- ConeConstructionError ----

    #[test]
    fn test_error_display() {
        assert_eq!(
            ConeConstructionError::NegativeRadius.to_string(),
            "radius is negative"
        );
        assert_eq!(
            ConeConstructionError::BadAngle.to_string(),
            "semi-angle is out of range"
        );
        assert_eq!(
            ConeConstructionError::NullAxis.to_string(),
            "the two points are confused"
        );
        assert_eq!(
            ConeConstructionError::NullAngle.to_string(),
            "semi-angle is null or a right angle"
        );
        assert_eq!(
            ConeConstructionError::NullNormal.to_string(),
            "axis direction has zero length"
        );
    }

    #[test]
    fn test_error_is_std_error() {
        fn takes_error(_e: &dyn std::error::Error) {}
        takes_error(&ConeConstructionError::NegativeRadius);
    }

    #[test]
    fn test_cone_contains() {
        let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        let tol = Tolerance::DEFAULT;
        assert!(cone.contains(cone.eval_point(0.0, 0.0), tol));
        assert!(cone.contains(cone.eval_point(1.0, 2.0), tol));
        assert!(!cone.contains(Point3D::ORIGIN, tol));
        assert!(!cone.contains(Point3D::new(0.0, 0.0, 10.0), tol));
    }

    #[test]
    fn test_cone_project_point() {
        let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        let tol = Tolerance::DEFAULT;
        let p = cone.eval_point(1.0, 2.0);
        let proj = cone.project_point(p, tol);
        assert!((proj.u - 1.0).abs() < 1e-9);
        assert!((proj.v - 2.0).abs() < 1e-9);
        assert_eq!(proj.distance, 0.0);
    }

    #[test]
    fn test_cone_cylinder_intersection() {
        let tol = Tolerance::DEFAULT;
        let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        let cylinder = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 2.0).unwrap();
        match cone.intersect_cylinder(&cylinder, tol) {
            ConeCylinderIntersection::Circle(circle) => {
                assert_eq!(circle.radius(), 2.0);
                assert!(circle.center().distance(Point3D::ORIGIN) < 1e-9);
                for p in circle.eval_points(&[0.0, 1.0, 3.0]) {
                    assert!(cone.contains(p, tol));
                    assert!(cylinder.contains(p, tol));
                }
            }
            _ => panic!("expected a circle"),
        }
        // Degenerate tube: just the apex.
        let axis_line = Cylinder::new(Point3D::ORIGIN, Vector3D::Z, 0.0).unwrap();
        assert_eq!(
            cone.intersect_cylinder(&axis_line, tol),
            ConeCylinderIntersection::ApexPoint(cone.apex())
        );
        // Offset axis: space quartic, no closed form.
        let off = Cylinder::new(Point3D::new(1.0, 0.0, 0.0), Vector3D::Z, 2.0).unwrap();
        assert_eq!(
            cone.intersect_cylinder(&off, tol),
            ConeCylinderIntersection::NotAnalytic
        );
        // Tilted axis: space quartic, no closed form.
        let tilted = Cylinder::new(
            Point3D::ORIGIN,
            Vector3D::new(0.0, 0.1, 1.0).normalized().unwrap(),
            2.0,
        )
        .unwrap();
        assert_eq!(
            cone.intersect_cylinder(&tilted, tol),
            ConeCylinderIntersection::NotAnalytic
        );
    }

    #[test]
    fn test_cone_cone_intersection() {
        let tol = Tolerance::DEFAULT;
        let c1 = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
        // Second coaxial cone: shifted frame, wider angle.
        let frame2 = Frame3D::new(Point3D::new(0.0, 0.0, 1.0), Vector3D::Z, Vector3D::X).unwrap();
        let c2 = Cone::from_frame(frame2, 0.6, 1.0).unwrap();
        match c1.intersect_cone(&c2, tol) {
            ConeConeIntersection::Circle(circle) => {
                for p in circle.eval_points(&[0.0, 1.0, 3.0]) {
                    assert!(c1.contains(p, tol));
                    assert!(c2.contains(p, tol));
                }
            }
            _ => panic!("expected a circle"),
        }
        // Parallel distinct nappes (same angle, offset apex): miss.
        let shifted = Cone::from_frame(frame2, 0.4, 1.0).unwrap();
        assert_eq!(
            c1.intersect_cone(&shifted, tol),
            ConeConeIntersection::Empty
        );
        // Tilted axis: space quartic, no closed form.
        let tilted = Cone::new(Point3D::ORIGIN, Vector3D::X, 0.4, 2.0).unwrap();
        assert_eq!(
            c1.intersect_cone(&tilted, tol),
            ConeConeIntersection::NotAnalytic
        );
    }
}