ogeom_geom/

curve.rs

//! Concrete space curves.
//!
//! Each analytic curve is a thin parameterization over the shape descriptions
//! in `ogeom-math`. The spline curve carries its own control points. All of them
//! are reachable through [`Curve`], an enum rather than a boxed trait object.
//!
//! # Why an enum
//!
//! Curves are stored in their millions in a real model, compared constantly,
//! and eventually serialized. An enum makes each of those cheap: no allocation,
//! no vtable, `Clone` and `PartialEq` derived, and (most usefully) exhaustive
//! matching, so adding a curve type produces a compile error at every site that
//! needs to know rather than a silent fallthrough.
//!
//! Deliberately *not* `#[non_exhaustive]`. Marking it so would force every
//! match outside this crate to carry a wildcard arm, which is exactly the
//! silent fallthrough the enum exists to prevent: a new curve type would then
//! compile everywhere and be mishandled everywhere. The cost is that adding a
//! variant is a breaking change, which for a kernel this size is the right
//! trade: a curve type nobody handles is worse than a version bump.
//!
//! [`CurveKind`] *is* non-exhaustive, because matching on it is for opting into
//! an analytic shortcut and a caller that does not recognise a kind should fall
//! back to the general path rather than fail to compile.
//!
//! The [`Curve3d`] trait is still the interface algorithms are written against;
//! [`Curve`] implements it and forwards.

use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
use ogeom_math::{
    Axis, Circle, Ellipse, Frame, Hyperbola, KnotVector, Parabola, Point, Transform, Vector,
    Weighted, bspline, elementary,
};

use crate::traits::{Continuity, Curve3d, CurveKind, Reversible, Transformable};

/// How far along a line the default domain reaches either side of its origin.
///
/// A line is unbounded, but every interface here works on a finite interval, so
/// an unbounded curve needs *some* domain. This is far beyond any real model
/// while staying well short of the range where `f64` spacing becomes coarse.
pub const LINE_EXTENT: f64 = 1.0e9;

/// A curve in space.
#[derive(Debug, Clone, PartialEq)]
pub enum Curve {
    /// A straight line.
    Line(LineCurve),
    /// A circle or arc.
    Circle(CircleCurve),
    /// An ellipse or arc.
    Ellipse(EllipseCurve),
    /// One branch of a hyperbola.
    Hyperbola(HyperbolaCurve),
    /// A parabola.
    Parabola(ParabolaCurve),
    /// A polynomial or rational B-spline.
    BSpline(BSplineCurve),
    /// A helix about an axis, the one transcendental the vocabulary keeps.
    Helix(HelixCurve),
    /// Another curve restricted to a sub-interval.
    Trimmed(Box<TrimmedCurve>),
    /// A curve at a constant distance from another, offset in the plane
    /// perpendicular to a reference direction.
    Offset(Box<OffsetCurve>),
    /// A surface curve: a pcurve composed with the surface it is drawn on.
    OnSurface(Box<CurveOnSurface>),
}

/// A straight line, parameterized by length from its origin.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct LineCurve {
    axis: Axis,
    domain: (f64, f64),
}

/// A circle, parameterized by angle.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CircleCurve {
    circle: Circle,
    reversed: bool,
}

/// An ellipse, parameterized by eccentric angle.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct EllipseCurve {
    ellipse: Ellipse,
    reversed: bool,
}

/// One branch of a hyperbola.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct HyperbolaCurve {
    hyperbola: Hyperbola,
    domain: (f64, f64),
    reversed: bool,
}

/// A parabola.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ParabolaCurve {
    parabola: Parabola,
    domain: (f64, f64),
    reversed: bool,
}

/// A helix about its frame's `z`, parameterized by turn angle.
///
/// The point at `t` sits at angle `t` around the axis, radius out along the
/// turned `x`, risen by `pitch·t/2π` along `z`, so one full turn advances
/// exactly one pitch, and a negative pitch winds the other hand. A non-zero
/// `taper` advances the radius the same way and winds a cone instead. A
/// helix is transcendental: no rational B-spline states it exactly, which
/// is why it is its own type rather than a conversion.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct HelixCurve {
    frame: Frame,
    radius: f64,
    pitch: f64,
    /// The radial advance per full turn: zero is the cylindrical helix,
    /// anything else winds a cone whose half-angle satisfies
    /// `tan(angle) = taper / pitch`.
    taper: f64,
    domain: (f64, f64),
    reversed: bool,
}

/// A B-spline curve, polynomial or rational.
#[derive(Debug, Clone, PartialEq)]
pub struct BSplineCurve {
    knots: KnotVector,
    control: Vec<Weighted<Point>>,
    rational: bool,
    periodic: bool,
}

/// Another curve restricted to a sub-interval of its domain.
#[derive(Debug, Clone, PartialEq)]
pub struct TrimmedCurve {
    basis: Curve,
    domain: (f64, f64),
    reversed: bool,
}

impl Curve {
    /// This curve over the wider parameter range `range`, tracing the same
    /// points at the same parameters over its own domain and running on as
    /// its own formula past it: a line, and a parabola, a hyperbola branch
    /// or a helix running forward, over any range; a forward trim over any
    /// range its basis holds. `None` for a curve whose formula does not
    /// reach past its domain (a B-spline, an offset or surface curve), one
    /// that repeats (a whole circle or ellipse has no end), one whose
    /// reversal mirrors its parameters across its domain, or a range that
    /// does not hold the domain.
    #[must_use]
    pub fn widened(&self, (lo, hi): (f64, f64)) -> Option<Self> {
        let (a, b) = self.domain();
        if !lo.is_finite() || !hi.is_finite() || lo > a || hi < b {
            return None;
        }
        match self {
            Self::Line(l) => Some(Self::Line(LineCurve {
                domain: (lo, hi),
                ..*l
            })),
            Self::Parabola(p) if !p.reversed => Some(Self::Parabola(ParabolaCurve {
                domain: (lo, hi),
                ..*p
            })),
            Self::Hyperbola(h) if !h.reversed => Some(Self::Hyperbola(HyperbolaCurve {
                domain: (lo, hi),
                ..*h
            })),
            Self::Helix(h) if !h.reversed => Some(Self::Helix(HelixCurve {
                domain: (lo, hi),
                ..*h
            })),
            Self::Trimmed(t) if !t.reversed => {
                let (lo, hi) = if t.basis.is_periodic() {
                    (lo, hi)
                } else {
                    let (ba, bb) = t.basis.domain();
                    (lo.max(ba), hi.min(bb))
                };
                (lo < a || hi > b).then(|| {
                    Self::Trimmed(Box::new(TrimmedCurve {
                        basis: t.basis.clone(),
                        domain: (lo, hi),
                        reversed: false,
                    }))
                })
            }
            _ => None,
        }
    }
}

impl LineCurve {
    /// A line along `axis`, spanning [`LINE_EXTENT`] either side of its origin.
    #[must_use]
    pub const fn new(axis: Axis) -> Self {
        Self {
            axis,
            domain: (-LINE_EXTENT, LINE_EXTENT),
        }
    }

    /// A line segment between two distinct points.
    ///
    /// The domain runs from zero to the distance between them, so the parameter
    /// is arc length, which makes every length query along the segment exact.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the points
    /// coincide.
    pub fn segment(from: Point, to: Point, tol: Tolerances) -> OgeomResult<Self> {
        let axis = Axis::through(from, to, tol)?;
        Ok(Self {
            axis,
            domain: (0.0, from.distance(to)),
        })
    }

    /// A line over an explicit parameter range.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the range
    /// is empty or non-finite.
    pub fn over(axis: Axis, start: f64, end: f64) -> OgeomResult<Self> {
        if !start.is_finite() || !end.is_finite() || end <= start {
            ogeom_bail!(Construction, "line range [{start}, {end}] is empty");
        }
        Ok(Self {
            axis,
            domain: (start, end),
        })
    }

    /// The underlying axis.
    #[must_use]
    pub const fn axis(&self) -> Axis {
        self.axis
    }
}

impl CircleCurve {
    /// A full circle, running counter-clockwise about its frame's `z`.
    #[must_use]
    pub const fn new(circle: Circle) -> Self {
        Self {
            circle,
            reversed: false,
        }
    }

    /// The underlying circle.
    #[must_use]
    pub const fn circle(&self) -> Circle {
        self.circle
    }

    /// Whether the curve runs backwards along its underlying circle.
    ///
    /// Part of the curve's state and not derivable from its circle, so
    /// anything that has to reproduce this curve exactly (the native format
    /// above all) needs to be able to read it.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }

    /// The same curve, the same point at every parameter, running forward:
    /// a reversed circle is its circle's frame turned half a turn about its
    /// `x` axis, the angle measured from the same `x` towards `-y`.
    ///
    /// # Errors
    ///
    /// None in practice: the turned frame is the given one's axes negated.
    pub fn forward(&self, tol: Tolerances) -> OgeomResult<Self> {
        if !self.reversed {
            return Ok(*self);
        }
        let frame = turned_about_x(self.circle.frame(), tol)?;
        Ok(Self::new(Circle::new(frame, self.circle.radius(), tol)?))
    }
}

/// A frame turned half a turn about its own `x` axis, of the same
/// handedness.
fn turned_about_x(frame: Frame, tol: Tolerances) -> OgeomResult<Frame> {
    Frame::from_axes(frame.origin(), frame.x(), -frame.y(), -frame.z(), tol)
}

impl EllipseCurve {
    /// A full ellipse.
    #[must_use]
    pub const fn new(ellipse: Ellipse) -> Self {
        Self {
            ellipse,
            reversed: false,
        }
    }

    /// The underlying ellipse.
    #[must_use]
    pub const fn ellipse(&self) -> Ellipse {
        self.ellipse
    }

    /// Whether the curve runs backwards along its underlying ellipse.
    ///
    /// Part of the curve's state and not derivable from its ellipse, so
    /// anything that has to reproduce this curve exactly (the native format
    /// above all) needs to be able to read it.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }

    /// The same curve, the same point at every parameter, running forward,
    /// as [`CircleCurve::forward`].
    ///
    /// # Errors
    ///
    /// None in practice: the turned frame is the given one's axes negated.
    pub fn forward(&self, tol: Tolerances) -> OgeomResult<Self> {
        if !self.reversed {
            return Ok(*self);
        }
        let e = self.ellipse;
        let frame = turned_about_x(e.frame(), tol)?;
        Ok(Self::new(Ellipse::new(
            frame,
            e.major_radius(),
            e.minor_radius(),
            tol,
        )?))
    }
}

impl HyperbolaCurve {
    /// A hyperbola branch over `[-extent, extent]` in its natural parameter.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `extent` is
    /// not finite and positive.
    pub fn new(hyperbola: Hyperbola, extent: f64) -> OgeomResult<Self> {
        if !extent.is_finite() || extent <= 0.0 {
            ogeom_bail!(
                Construction,
                "hyperbola extent {extent} must be finite and positive"
            );
        }
        Ok(Self {
            hyperbola,
            domain: (-extent, extent),
            reversed: false,
        })
    }

    /// A hyperbola branch over an arbitrary `[start, end]` in its natural
    /// parameter.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// interval is not finite and increasing.
    pub fn over(hyperbola: Hyperbola, start: f64, end: f64) -> OgeomResult<Self> {
        if !start.is_finite() || !end.is_finite() || start >= end {
            ogeom_bail!(
                Construction,
                "hyperbola domain [{start}, {end}] must be finite and increasing"
            );
        }
        Ok(Self {
            hyperbola,
            domain: (start, end),
            reversed: false,
        })
    }

    /// The underlying hyperbola.
    #[must_use]
    pub const fn hyperbola(&self) -> Hyperbola {
        self.hyperbola
    }

    /// Whether the curve runs backwards along its underlying hyperbola.
    ///
    /// Part of the curve's state and not derivable from its hyperbola, so
    /// anything that has to reproduce this curve exactly (the native format
    /// above all) needs to be able to read it.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }
}

impl ParabolaCurve {
    /// A parabola over `[-extent, extent]`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if `extent` is
    /// not finite and positive.
    pub fn new(parabola: Parabola, extent: f64) -> OgeomResult<Self> {
        if !extent.is_finite() || extent <= 0.0 {
            ogeom_bail!(
                Construction,
                "parabola extent {extent} must be finite and positive"
            );
        }
        Ok(Self {
            parabola,
            domain: (-extent, extent),
            reversed: false,
        })
    }

    /// A parabola over an arbitrary `[start, end]`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// interval is not finite and increasing.
    pub fn over(parabola: Parabola, start: f64, end: f64) -> OgeomResult<Self> {
        if !start.is_finite() || !end.is_finite() || start >= end {
            ogeom_bail!(
                Construction,
                "parabola domain [{start}, {end}] must be finite and increasing"
            );
        }
        Ok(Self {
            parabola,
            domain: (start, end),
            reversed: false,
        })
    }

    /// The underlying parabola.
    #[must_use]
    pub const fn parabola(&self) -> Parabola {
        self.parabola
    }

    /// Whether the curve runs backwards along its underlying parabola.
    ///
    /// Part of the curve's state and not derivable from its parabola, so
    /// anything that has to reproduce this curve exactly (the native format
    /// above all) needs to be able to read it.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }
}

impl HelixCurve {
    /// A helix over `turns` full revolutions from angle zero.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
    /// the radius is not finite and positive, the pitch is not finite and
    /// non-zero (a zero pitch is a circle, and there is a type for that),
    /// or `turns` is not finite and positive.
    pub fn new(frame: Frame, radius: f64, pitch: f64, turns: f64) -> OgeomResult<Self> {
        if !turns.is_finite() || turns <= 0.0 {
            ogeom_bail!(Construction, "a helix over {turns} turns is not a curve");
        }
        Self::over(frame, radius, pitch, 0.0, core::f64::consts::TAU * turns)
    }

    /// A helix over an arbitrary increasing angle interval.
    ///
    /// # Errors
    ///
    /// As [`HelixCurve::new`], with the interval checked instead of the turn
    /// count.
    pub fn over(frame: Frame, radius: f64, pitch: f64, start: f64, end: f64) -> OgeomResult<Self> {
        if !radius.is_finite() || radius <= 0.0 {
            ogeom_bail!(
                Construction,
                "helix radius {radius} must be finite and positive"
            );
        }
        if !pitch.is_finite() || pitch == 0.0 {
            ogeom_bail!(
                Construction,
                "helix pitch {pitch} must be finite and non-zero; a zero pitch is a circle"
            );
        }
        if !start.is_finite() || !end.is_finite() || start >= end {
            ogeom_bail!(
                Construction,
                "helix domain [{start}, {end}] must be finite and increasing"
            );
        }
        Ok(Self {
            frame,
            radius,
            pitch,
            taper: 0.0,
            domain: (start, end),
            reversed: false,
        })
    }

    /// A conical helix: the radius advances by `taper` per full turn while
    /// the point rises by `pitch`, winding the cone whose half-angle
    /// satisfies `tan(angle) = taper / pitch`. `radius` is the radius at
    /// angle zero, whether or not the interval holds it.
    ///
    /// # Errors
    ///
    /// As [`HelixCurve::over`], and additionally if `taper` is not finite
    /// or the radius runs non-positive anywhere on the interval. Past the
    /// apex there is no cone to wind.
    pub fn conical(
        frame: Frame,
        radius: f64,
        pitch: f64,
        taper: f64,
        start: f64,
        end: f64,
    ) -> OgeomResult<Self> {
        let mut helix = Self::over(frame, radius, pitch, start, end)?;
        if !taper.is_finite() {
            ogeom_bail!(Construction, "helix taper {taper} must be finite");
        }
        let slope = taper / core::f64::consts::TAU;
        let (ra, rb) = (slope.mul_add(start, radius), slope.mul_add(end, radius));
        if ra <= 0.0 || rb <= 0.0 {
            ogeom_bail!(
                Construction,
                "the helix radius runs non-positive on [{start}, {end}]; \
                 past the apex there is no cone to wind"
            );
        }
        helix.taper = taper;
        Ok(helix)
    }

    /// The frame the helix turns about.
    #[must_use]
    pub const fn frame(&self) -> &Frame {
        &self.frame
    }

    /// The radius.
    #[must_use]
    pub const fn radius(&self) -> f64 {
        self.radius
    }

    /// The advance along the axis per full turn. Negative winds left-handed.
    #[must_use]
    pub const fn pitch(&self) -> f64 {
        self.pitch
    }

    /// The radial advance per full turn. Zero is the cylindrical helix.
    #[must_use]
    pub const fn taper(&self) -> f64 {
        self.taper
    }

    /// Whether evaluation runs the domain backwards.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }

    /// The exact arc length between two parameters.
    ///
    /// The speed at angle `t` is `√(r(t)² + k²)`, where the radius
    /// `r(t) = radius + s·t` advances by `s = taper/2π` a radian and
    /// `k² = s² + (pitch/2π)²`. On a cylinder (`s` zero) that is constant, and
    /// the length is the speed times the swept angle. On a cone it is
    /// `∫√(r² + k²) dr / s`, in closed form.
    #[must_use]
    pub fn arc_length(&self, from: f64, to: f64) -> f64 {
        let tau = core::f64::consts::TAU;
        let (rise, slope) = (self.pitch / tau, self.taper / tau);
        if slope == 0.0 {
            return (to - from).abs() * self.radius.hypot(rise);
        }
        let k = slope.hypot(rise);
        let primitive = |r: f64| f64::midpoint(r * r.hypot(k), k * k * (r / k).asinh());
        let (r0, r1) = (
            slope.mul_add(from, self.radius),
            slope.mul_add(to, self.radius),
        );
        ((primitive(r1) - primitive(r0)) / slope).abs()
    }

    /// Point and first three derivatives at the raw (unreversed) angle.
    fn at(&self, t: f64) -> (Point, Vector, Vector, Vector) {
        let (sin, cos) = t.sin_cos();
        let x = self.frame.x().vector();
        let y = self.frame.y().vector();
        let z = self.frame.z().vector();
        let rise = self.pitch / core::f64::consts::TAU;
        let slope = self.taper / core::f64::consts::TAU;
        let r = slope.mul_add(t, self.radius);
        let point = self.frame.origin() + x * (r * cos) + y * (r * sin) + z * (rise * t);
        let d1 = x * slope.mul_add(cos, -(r * sin)) + y * slope.mul_add(sin, r * cos) + z * rise;
        let d2 = x * (2.0 * slope).mul_add(-sin, -(r * cos))
            + y * (2.0 * slope).mul_add(cos, -(r * sin));
        let d3 =
            x * (3.0 * slope).mul_add(-cos, r * sin) + y * (3.0 * slope).mul_add(-sin, -(r * cos));
        (point, d1, d2, d3)
    }
}

impl Curve3d for HelixCurve {
    fn domain(&self) -> (f64, f64) {
        self.domain
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        Ok(self.at(t).0)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        let (point, _, _, _) = self.at(t);
        // The radius runs linearly in the angle, so the `k`th derivative of
        // `r cos t` is `r cos⁽ᵏ⁾ t + k s cos⁽ᵏ⁻¹⁾ t`, and alike for the sine;
        // the rise survives only to first order. The chain rule for the
        // reversal flips odd orders.
        let (sin, cos) = t.sin_cos();
        let x = self.frame.x().vector();
        let y = self.frame.y().vector();
        let z = self.frame.z().vector();
        let rise = self.pitch / core::f64::consts::TAU;
        let slope = self.taper / core::f64::consts::TAU;
        let r = slope.mul_add(t, self.radius);
        let cos_k = |k: usize| [cos, -sin, -cos, sin][k % 4];
        let sin_k = |k: usize| [sin, cos, -sin, -cos][k % 4];
        let sign = if self.reversed { -1.0 } else { 1.0 };
        let mut out = vec![point.to_vector()];
        let mut factor = 1.0;
        for k in 1..=n {
            factor *= sign;
            #[allow(clippy::cast_precision_loss)]
            let ks = k as f64 * slope;
            let mut d = x * ks.mul_add(cos_k(k - 1), r * cos_k(k))
                + y * ks.mul_add(sin_k(k - 1), r * sin_k(k));
            if k == 1 {
                d += z * rise;
            }
            out.push(d * factor);
        }
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Helix
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        false
    }

    fn is_periodic(&self) -> bool {
        false
    }
}

/// A space curve displaced a constant signed distance perpendicular to a
/// reference direction: the offset direction at `t` is the unit vector of
/// `tangent x reference`, the classical spelling.
///
/// Point and first derivative are exact from the basis's first two
/// derivatives. The second would need the basis's third, which the
/// vocabulary does not carry, so `derivatives_at` beyond order one refuses
/// by name. Where the tangent runs along the reference the offset direction
/// is undefined and evaluation refuses.
#[derive(Debug, Clone, PartialEq)]
pub struct OffsetCurve {
    basis: Curve,
    distance: f64,
    reference: ogeom_math::Direction,
}

impl OffsetCurve {
    /// Offset `basis` by `distance` along `tangent x reference`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// distance is not finite and non-zero.
    pub fn new(basis: Curve, distance: f64, reference: ogeom_math::Direction) -> OgeomResult<Self> {
        if !distance.is_finite() || distance == 0.0 {
            ogeom_bail!(
                Construction,
                "an offset of {distance} is not a displacement"
            );
        }
        Ok(Self {
            basis,
            distance,
            reference,
        })
    }

    /// The curve being offset.
    #[must_use]
    pub const fn basis(&self) -> &Curve {
        &self.basis
    }

    /// The signed displacement.
    #[must_use]
    pub const fn distance(&self) -> f64 {
        self.distance
    }

    /// The direction the offset plane is perpendicular to.
    #[must_use]
    pub const fn reference(&self) -> ogeom_math::Direction {
        self.reference
    }

    /// The unit offset direction and its derivative at `t`.
    fn direction_and_slope(&self, t: f64, tol: Tolerances) -> OgeomResult<(Vector, Vector)> {
        let d = self.basis.derivatives_at(t, 2, tol)?;
        let v = self.reference.vector();
        let w = d[1].cross(v);
        let m = w.magnitude();
        if m <= tol.confusion() * d[1].magnitude().max(1.0) {
            ogeom_bail!(
                Construction,
                "the tangent at {t} runs along the reference; the offset direction is undefined"
            );
        }
        let n = w / m;
        let wp = d[2].cross(v);
        let np = (wp - n * n.dot(wp)) / m;
        Ok((n, np))
    }
}

impl Curve3d for OffsetCurve {
    fn domain(&self) -> (f64, f64) {
        self.basis.domain()
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let base = self.basis.point_at(u, tol)?;
        let (n, _) = self.direction_and_slope(u, tol)?;
        Ok(base + n * self.distance)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        let d1 = self.basis.d1_at(u, tol)?;
        let (_, np) = self.direction_and_slope(u, tol)?;
        Ok(d1 + np * self.distance)
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        if n >= 2 {
            ogeom_bail!(
                Construction,
                "an offset curve's second derivative needs its basis's third, which the \
                 vocabulary does not carry"
            );
        }
        let mut out = vec![self.point_at(u, tol)?.to_vector()];
        if n >= 1 {
            out.push(self.d1_at(u, tol)?);
        }
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Offset
    }

    fn continuity(&self) -> Continuity {
        match self.basis.continuity() {
            Continuity::CInfinity => Continuity::CInfinity,
            Continuity::C2 | Continuity::G2 => Continuity::C1,
            Continuity::C1 | Continuity::G1 | Continuity::C0 => Continuity::C0,
        }
    }

    fn is_closed(&self, tol: Tolerances) -> bool {
        self.basis.is_closed(tol)
    }

    fn is_periodic(&self) -> bool {
        self.basis.is_periodic()
    }
}

/// A pcurve composed with the surface it is drawn on: the space curve a
/// trimming boundary actually traces.
///
/// Everything is exact: the chain rule composes the pcurve's derivatives
/// with the surface's, both of which the vocabulary carries to second
/// order.
#[derive(Debug, Clone, PartialEq)]
pub struct CurveOnSurface {
    pcurve: crate::curve2d::PlanarCurve,
    surface: crate::surface::SurfaceGeometry,
}

impl CurveOnSurface {
    /// The composition of `pcurve` with `surface`.
    #[must_use]
    pub const fn new(
        pcurve: crate::curve2d::PlanarCurve,
        surface: crate::surface::SurfaceGeometry,
    ) -> Self {
        Self { pcurve, surface }
    }

    /// The curve in parameter space.
    #[must_use]
    pub const fn pcurve(&self) -> &crate::curve2d::PlanarCurve {
        &self.pcurve
    }

    /// The surface the pcurve is drawn on.
    #[must_use]
    pub const fn surface(&self) -> &crate::surface::SurfaceGeometry {
        &self.surface
    }
}

impl Curve3d for CurveOnSurface {
    fn domain(&self) -> (f64, f64) {
        use crate::traits::Curve2d as _;
        self.pcurve.domain()
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        use crate::traits::{Curve2d as _, Surface as _};
        let p = self.pcurve.point_at(u, tol)?;
        self.surface.point_at(p.x, p.y, tol)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        use crate::traits::{Curve2d as _, Surface as _};
        if n >= 3 {
            ogeom_bail!(
                Construction,
                "a surface curve's third derivative needs the surface's, which the \
                 vocabulary does not carry"
            );
        }
        let d = self.pcurve.derivatives_at(u, n.max(2), tol)?;
        let at = d[0];
        let point = self.surface.point_at(at.x, at.y, tol)?;
        let mut out = vec![point.to_vector()];
        if n >= 1 {
            let (su, sv) = self.surface.d1_at(at.x, at.y, tol)?;
            out.push(su * d[1].x + sv * d[1].y);
            if n >= 2 {
                let (suu, suv, svv) = self.surface.d2_at(at.x, at.y, tol)?;
                // The chain rule's second order: quadratic in the pcurve's
                // slope, linear in its curvature.
                let second = suu * (d[1].x * d[1].x)
                    + suv * (2.0 * d[1].x * d[1].y)
                    + svv * (d[1].y * d[1].y)
                    + su * d[2].x
                    + sv * d[2].y;
                out.push(second);
            }
        }
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::OnSurface
    }

    fn continuity(&self) -> Continuity {
        use crate::traits::Surface as _;
        // The composition is as smooth as the rougher of the two.
        self.pcurve_continuity().min(self.surface.continuity())
    }

    fn is_closed(&self, tol: Tolerances) -> bool {
        use crate::traits::Curve2d as _;
        self.pcurve.is_closed(tol)
    }

    fn is_periodic(&self) -> bool {
        use crate::traits::Curve2d as _;
        self.pcurve.is_periodic()
    }
}

impl CurveOnSurface {
    fn pcurve_continuity(&self) -> Continuity {
        // Planar curves in the vocabulary are analytic or spline. The spline
        // reports through its own knots elsewhere, and C2 is the floor the
        // fitting machinery guarantees. Conservative either way.
        Continuity::C2
    }
}

impl BSplineCurve {
    /// A polynomial B-spline from a knot vector and control points.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) if the control point
    /// count disagrees with the knot vector.
    pub fn new(knots: KnotVector, control: Vec<Point>, tol: Tolerances) -> OgeomResult<Self> {
        let weighted = control
            .into_iter()
            .map(|p| Weighted::new(p, 1.0, tol))
            .collect::<OgeomResult<Vec<_>>>()?;
        Self::rational(knots, weighted)
    }

    /// A rational B-spline from a knot vector and weighted control points.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) if the control point
    /// count disagrees with the knot vector.
    pub fn rational(knots: KnotVector, control: Vec<Weighted<Point>>) -> OgeomResult<Self> {
        if control.len() != knots.control_point_count() {
            ogeom_bail!(
                Dimension,
                "knot vector describes {} control points, got {}",
                knots.control_point_count(),
                control.len()
            );
        }
        // A curve whose weights are all equal is polynomial regardless of what
        // that common value is, and saying so lets evaluation skip the divide.
        let first = control[0].weight;
        let rational = control
            .iter()
            .any(|w| (w.weight - first).abs() > 1e-12 * first.abs());
        Ok(Self {
            knots,
            control,
            rational,
            periodic: false,
        })
    }

    /// A smoothly periodic B-spline through a ring of control points.
    ///
    /// The ring is wrapped: the first `degree` controls repeat past the
    /// end over a uniform knot vector, and evaluation wraps its parameter,
    /// so the loop closes with `degree - 1` continuous derivatives and no
    /// clamped seam. The domain runs one knot step per ring point.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the ring
    /// has no more points than the degree, or the degree is zero.
    pub fn periodic(control: &[Point], degree: usize, tol: Tolerances) -> OgeomResult<Self> {
        if degree == 0 {
            ogeom_bail!(Construction, "a curve needs a degree of at least one");
        }
        if control.len() <= degree {
            ogeom_bail!(
                Construction,
                "a periodic ring of {} points cannot carry degree {degree}",
                control.len()
            );
        }
        let n = control.len();
        let mut wrapped: Vec<Point> = Vec::with_capacity(n + degree);
        wrapped.extend_from_slice(control);
        wrapped.extend_from_slice(&control[..degree]);
        #[allow(clippy::cast_precision_loss)]
        let knots: Vec<f64> = (0..wrapped.len() + degree + 1).map(|i| i as f64).collect();
        let mut built = Self::new(KnotVector::new(knots, degree)?, wrapped, tol)?;
        built.periodic = true;
        Ok(built)
    }

    /// The exact periodic representation, as a reader restores it: wrapped
    /// knots and control, with the wrap verified rather than assumed.
    ///
    /// # Errors
    ///
    /// As [`BSplineCurve::rational`], and additionally if the trailing
    /// `degree` controls do not repeat the leading ones. An unwrapped ring
    /// evaluated periodically would tear at the seam.
    pub fn periodic_from_parts(
        knots: KnotVector,
        control: Vec<Weighted<Point>>,
        tol: Tolerances,
    ) -> OgeomResult<Self> {
        let degree = knots.degree();
        if control.len() <= degree {
            ogeom_bail!(
                Construction,
                "a periodic curve of {} controls cannot carry degree {degree}",
                control.len()
            );
        }
        let n = control.len() - degree;
        for i in 0..degree {
            let (a, b) = (control[i], control[n + i]);
            if !a.point().is_equal(b.point(), tol) || (a.weight - b.weight).abs() > 1e-12 {
                ogeom_bail!(
                    Construction,
                    "a periodic curve's trailing controls must repeat its \
                     leading ones; control {} does not",
                    n + i
                );
            }
        }
        let mut built = Self::rational(knots, control)?;
        built.periodic = true;
        Ok(built)
    }

    /// The knot vector.
    #[must_use]
    pub const fn knots(&self) -> &KnotVector {
        &self.knots
    }

    /// The weighted control points.
    #[must_use]
    pub fn control_points(&self) -> &[Weighted<Point>] {
        &self.control
    }

    /// Whether the weights differ, so the curve is genuinely rational.
    #[must_use]
    pub const fn is_rational(&self) -> bool {
        self.rational
    }

    /// The degree.
    #[must_use]
    pub const fn degree(&self) -> usize {
        self.knots.degree()
    }

    /// Insert a knot without moving the curve.
    ///
    /// # Errors
    ///
    /// As [`bspline::insert_knot`].
    pub fn with_knot_inserted(&self, u: f64, count: usize, tol: Tolerances) -> OgeomResult<Self> {
        let (knots, control) = bspline::insert_knot(&self.knots, &self.control, u, count, tol)?;
        Ok(Self {
            knots,
            control,
            rational: self.rational,
            periodic: self.periodic,
        })
    }

    /// Raise the degree without moving the curve.
    ///
    /// # Errors
    ///
    /// As [`bspline::elevate_degree`].
    pub fn elevated(&self, tol: Tolerances) -> OgeomResult<Self> {
        let (knots, control) = bspline::elevate_degree(&self.knots, &self.control, tol)?;
        Ok(Self {
            knots,
            control,
            rational: self.rational,
            periodic: self.periodic,
        })
    }

    /// The same closed curve with its seam moved to `u`: what was the
    /// stretch from `u` to the end now comes first, and the stretch from
    /// the start to `u` follows it, joined where the old seam was. The
    /// domain keeps its length and begins at `u`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// curve is periodic (its seam is nowhere) or does not close, its
    /// two ends apart by more than a thousand confusions; as
    /// [`bspline::split`] if `u` is an end of the domain.
    pub fn reseamed_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Self> {
        if self.periodic {
            ogeom_bail!(Construction, "a periodic curve has no seam to move");
        }
        let (start, end) = self.domain();
        let (head, tail) = (self.point_at(start, tol)?, self.point_at(end, tol)?);
        // Closed to a thousand confusions: a marched section closes to
        // its own march's tolerance, and a loop whose ends sit a fraction
        // of a micron apart is closed for every purpose the seam serves.
        if head.distance(tail) > tol.confusion() * 1e3 {
            ogeom_bail!(
                Construction,
                "the curve does not close: its ends are {:.3e} apart",
                head.distance(tail)
            );
        }
        let (before, after) = self.split_at(u, tol)?;
        let (knots, control) = bspline::join(
            &(after.knots, after.control),
            &(before.knots, before.control),
        )?;
        Ok(Self {
            knots,
            control,
            rational: self.rational,
            periodic: self.periodic,
        })
    }

    /// This curve continued past one end by `length` in space: the
    /// polynomial continuation of its own end derivatives to the order
    /// `continuity`, joined on. A polynomial run of degree at most that
    /// order continues as itself, and so does a rational arc's homogeneous
    /// polynomial: a circle arc continued at order two stays on its circle.
    ///
    /// The continuation's parameter span is solved until its arc length is
    /// `length`. Extended at the start, the original run keeps its
    /// parameters and the domain grows downward.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// curve is periodic (it has no end to continue from) or stands still
    /// at that end; as [`bspline::extend`].
    pub fn extended(
        &self,
        at_end: bool,
        length: f64,
        continuity: usize,
        tol: Tolerances,
    ) -> OgeomResult<Self> {
        if !(length > 0.0 && length.is_finite()) {
            ogeom_bail!(
                Construction,
                "an extension needs a positive length; got {length}"
            );
        }
        let speed = self.end_speed(at_end, tol)?;
        let build = |span: f64| -> OgeomResult<Self> {
            let (knots, control) =
                bspline::extend(&self.knots, &self.control, at_end, span, continuity, tol)?;
            Ok(Self {
                knots,
                control,
                rational: self.rational,
                periodic: self.periodic,
            })
        };
        // The continuation over a span is the same polynomial whatever the
        // span: its length grows with the span, and the secant finds the
        // span that meets the length.
        let run = |curve: &Self, span: f64| -> OgeomResult<f64> {
            let (lo, hi) = curve.domain();
            let (a, b) = if at_end {
                (hi - span, hi)
            } else {
                (lo, lo + span)
            };
            curve.length_over((a, b), tol)
        };
        let mut s0 = length / speed;
        let mut c0 = build(s0)?;
        let mut l0 = run(&c0, s0)? - length;
        let mut s1 = s0 * 1.1;
        for _ in 0..40 {
            if l0.abs() <= tol.confusion() {
                break;
            }
            let c1 = build(s1)?;
            let l1 = run(&c1, s1)? - length;
            let next = if (l1 - l0).abs() > f64::EPSILON {
                s1 - l1 * (s1 - s0) / (l1 - l0)
            } else {
                s1
            };
            (s0, c0, l0) = (s1, c1, l1);
            s1 = if next > 0.0 { next } else { s0 * 0.5 };
        }
        Ok(c0)
    }

    /// This curve continued past one end to `target`: a piece carrying the
    /// curve's end derivatives to the order `continuity` and ending at the
    /// point, joined on. Extended at the start, the original run keeps its
    /// parameters and the domain grows downward.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// curve is periodic, stands still at that end, or already ends at
    /// `target`; as [`bspline::extend_to`].
    pub fn extended_to(
        &self,
        at_end: bool,
        target: Point,
        continuity: usize,
        tol: Tolerances,
    ) -> OgeomResult<Self> {
        let speed = self.end_speed(at_end, tol)?;
        let (lo, hi) = self.domain();
        let at = if at_end { hi } else { lo };
        let from = self.point_at(at, tol)?;
        let gap = from.distance(target);
        if gap <= tol.confusion() {
            ogeom_bail!(Construction, "the curve already ends at the point");
        }
        let end = if at_end {
            self.control[self.control.len() - 1]
        } else {
            self.control[0]
        };
        let weighted = Weighted::new(target, end.weight, tol)?;
        let (knots, control) = bspline::extend_to(
            &self.knots,
            &self.control,
            at_end,
            weighted,
            gap / speed,
            continuity,
            tol,
        )?;
        let rational = self.rational;
        Ok(Self {
            knots,
            control,
            rational,
            periodic: false,
        })
    }

    /// The speed at one end, refused where it is none.
    fn end_speed(&self, at_end: bool, tol: Tolerances) -> OgeomResult<f64> {
        if self.periodic {
            ogeom_bail!(Construction, "a periodic curve has no end to continue from");
        }
        let (lo, hi) = self.domain();
        let at = if at_end { hi } else { lo };
        let speed = self.d1_at(at, tol)?.magnitude();
        if speed <= tol.confusion() {
            ogeom_bail!(
                Construction,
                "the curve stands still at its end; there is no direction to continue in"
            );
        }
        Ok(speed)
    }

    /// The arc length over a parameter range, by Gauss-Legendre on each of
    /// the knot spans it covers.
    fn length_over(&self, range: (f64, f64), tol: Tolerances) -> OgeomResult<f64> {
        const NODES: [(f64, f64); 5] = [
            (0.0, 0.568_888_888_888_888_9),
            (-0.538_469_310_105_683_1, 0.478_628_670_499_366_5),
            (0.538_469_310_105_683_1, 0.478_628_670_499_366_5),
            (-0.906_179_845_938_664, 0.236_926_885_056_189_1),
            (0.906_179_845_938_664, 0.236_926_885_056_189_1),
        ];
        let mut breaks: Vec<f64> = vec![range.0];
        breaks.extend(
            self.knots
                .distinct()
                .into_iter()
                .map(|(k, _)| k)
                .filter(|k| *k > range.0 && *k < range.1),
        );
        breaks.push(range.1);
        let mut total = 0.0;
        for w in breaks.windows(2) {
            // Each span in eight, for a rational run's uneven speed.
            for part in 0..8 {
                let a = w[0] + (w[1] - w[0]) * f64::from(part) / 8.0;
                let b = w[0] + (w[1] - w[0]) * f64::from(part + 1) / 8.0;
                let (mid, half) = (0.5 * (a + b), 0.5 * (b - a));
                for (x, weight) in NODES {
                    total += weight * half * self.d1_at(mid + half * x, tol)?.magnitude();
                }
            }
        }
        Ok(total)
    }

    /// The piece of this curve over `range`, exactly and keeping its
    /// parameters: the piece at `t` is this curve at `t`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if `range` is
    /// empty or leaves the domain.
    pub fn segment(&self, range: (f64, f64), tol: Tolerances) -> OgeomResult<Self> {
        let (a, b) = self.knots.domain();
        let eps = tol.parametric();
        if range.1 <= range.0 + eps || range.0 < a - eps || range.1 > b + eps {
            ogeom_bail!(
                Domain,
                "[{}, {}] is no piece of [{a}, {b}]",
                range.0,
                range.1
            );
        }
        let mut piece = Self {
            knots: self.knots.clone(),
            control: self.control.clone(),
            rational: self.rational,
            periodic: false,
        };
        if range.0 > a + eps {
            piece = piece.split_at(range.0, tol)?.1;
        }
        if range.1 < b - eps {
            piece = piece.split_at(range.1, tol)?.0;
        }
        Ok(piece)
    }

    /// Split into two curves meeting at `u`.
    ///
    /// # Errors
    ///
    /// As [`bspline::split`].
    pub fn split_at(&self, u: f64, tol: Tolerances) -> OgeomResult<(Self, Self)> {
        let ((lk, lc), (rk, rc)) = bspline::split(&self.knots, &self.control, u, tol)?;
        Ok((
            Self {
                knots: lk,
                control: lc,
                rational: self.rational,
                periodic: self.periodic,
            },
            Self {
                knots: rk,
                control: rc,
                rational: self.rational,
                periodic: self.periodic,
            },
        ))
    }
}

impl TrimmedCurve {
    /// Restrict `basis` to `[start, end]`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the range is empty or
    /// falls outside the basis curve's own domain.
    pub fn new(basis: Curve, start: f64, end: f64, tol: Tolerances) -> OgeomResult<Self> {
        let (a, b) = basis.domain();
        if !start.is_finite() || !end.is_finite() || end <= start + tol.parametric() {
            ogeom_bail!(Domain, "trim range [{start}, {end}] is empty");
        }
        if !basis.is_periodic() && (start < a - tol.parametric() || end > b + tol.parametric()) {
            ogeom_bail!(
                Domain,
                "trim range [{start}, {end}] leaves the basis domain [{a}, {b}]"
            );
        }
        Ok(Self {
            basis,
            domain: (start, end),
            reversed: false,
        })
    }

    /// The curve being trimmed.
    #[must_use]
    pub const fn basis(&self) -> &Curve {
        &self.basis
    }

    /// Whether the curve runs backwards along its underlying curve.
    ///
    /// Part of the curve's state and not derivable from its basis curve, so
    /// anything that has to reproduce this curve exactly (the native format
    /// above all) needs to be able to read it.
    #[must_use]
    pub const fn is_reversed(&self) -> bool {
        self.reversed
    }

    /// This curve's parameter mapped onto the basis curve's.
    fn basis_parameter(&self, u: f64, tol: Tolerances) -> OgeomResult<f64> {
        let u = self.normalize_parameter(u, tol)?;
        Ok(if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        })
    }
}

/// Reverse a parameter within `[a, b]`, so the curve runs the other way over
/// the same interval.
///
/// Preserving the domain matters: trimming ranges elsewhere refer to it, and a
/// reversal that also renumbered the parameters would invalidate them.
fn mirror(u: f64, a: f64, b: f64) -> f64 {
    a + b - u
}

/// A conic's derivatives to order `n` from its first two: a circle's and an
/// ellipse's cycle `d1, d2, -d1, -d2` (`period` four), a hyperbola's
/// `d1, d2` (`period` two). Each order is times `sign` to its power, the
/// reversal's chain rule.
fn cycled(
    point: Vector,
    d1: Vector,
    d2: Vector,
    n: usize,
    period: usize,
    sign: f64,
) -> Vec<Vector> {
    let mut out = vec![point];
    let mut factor = 1.0;
    for order in 1..=n {
        factor *= sign;
        let d = match (order - 1) % period {
            0 => d1,
            1 => d2,
            2 => -d1,
            _ => -d2,
        };
        out.push(d * factor);
    }
    out
}

impl Curve3d for LineCurve {
    fn domain(&self) -> (f64, f64) {
        self.domain
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        Ok(self.axis.point_at(u))
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        self.normalize_parameter(u, tol)?;
        Ok(self.axis.direction.vector())
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let p = self.point_at(u, tol)?;
        let mut out = vec![p.to_vector(), self.axis.direction.vector()];
        out.resize(n + 1, Vector::ZERO);
        out.truncate(n + 1);
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Line
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        false
    }

    fn is_periodic(&self) -> bool {
        false
    }
}

impl Curve3d for CircleCurve {
    fn domain(&self) -> (f64, f64) {
        (0.0, core::f64::consts::TAU)
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        let angle = if self.reversed { -u } else { u };
        Ok(elementary::circle_at(&self.circle, angle).point)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let angle = if self.reversed { -u } else { u };
        let c = elementary::circle_at(&self.circle, angle);
        let sign = if self.reversed { -1.0 } else { 1.0 };
        Ok(cycled(c.point.to_vector(), c.d1, c.d2, n, 4, sign))
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Circle
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        true
    }

    fn is_periodic(&self) -> bool {
        true
    }
}

impl Curve3d for EllipseCurve {
    fn domain(&self) -> (f64, f64) {
        (0.0, core::f64::consts::TAU)
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        let angle = if self.reversed { -u } else { u };
        Ok(elementary::ellipse_at(&self.ellipse, angle).point)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let angle = if self.reversed { -u } else { u };
        let c = elementary::ellipse_at(&self.ellipse, angle);
        let sign = if self.reversed { -1.0 } else { 1.0 };
        Ok(cycled(c.point.to_vector(), c.d1, c.d2, n, 4, sign))
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Ellipse
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        true
    }

    fn is_periodic(&self) -> bool {
        true
    }
}

impl Curve3d for HyperbolaCurve {
    fn domain(&self) -> (f64, f64) {
        self.domain
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        Ok(elementary::hyperbola_at(&self.hyperbola, t).point)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        let c = elementary::hyperbola_at(&self.hyperbola, t);
        let sign = if self.reversed { -1.0 } else { 1.0 };
        Ok(cycled(c.point.to_vector(), c.d1, c.d2, n, 2, sign))
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Hyperbola
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        false
    }

    fn is_periodic(&self) -> bool {
        false
    }
}

impl Curve3d for ParabolaCurve {
    fn domain(&self) -> (f64, f64) {
        self.domain
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        Ok(elementary::parabola_at(&self.parabola, t).point)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let t = if self.reversed {
            mirror(u, self.domain.0, self.domain.1)
        } else {
            u
        };
        let c = elementary::parabola_at(&self.parabola, t);
        let sign = if self.reversed { -1.0 } else { 1.0 };
        // Quadratic in its parameter: every order past the second is zero.
        let mut out = vec![c.point.to_vector(), c.d1 * sign, c.d2];
        out.resize(n.max(2) + 1, Vector::ZERO);
        out.truncate(n + 1);
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Parabola
    }

    fn continuity(&self) -> Continuity {
        Continuity::CInfinity
    }

    fn is_closed(&self, _tol: Tolerances) -> bool {
        false
    }

    fn is_periodic(&self) -> bool {
        false
    }
}

impl Curve3d for BSplineCurve {
    fn domain(&self) -> (f64, f64) {
        self.knots.domain()
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        let u = self.normalize_parameter(u, tol)?;
        if self.rational {
            bspline::evaluate_rational(&self.knots, &self.control, u, tol)
        } else {
            // All weights equal: the projection is a no-op up to that common
            // factor, so the cheaper polynomial path is exact here.
            Ok(bspline::evaluate(&self.knots, &self.control, u, tol)?.point())
        }
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        Ok(self.derivatives_at(u, 1, tol)?[1])
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let u = self.normalize_parameter(u, tol)?;
        let points = bspline::rational_derivatives(&self.knots, &self.control, u, n, tol)?;
        Ok(points.into_iter().map(Point::to_vector).collect())
    }

    fn kind(&self) -> CurveKind {
        CurveKind::BSpline
    }

    /// Continuity across the whole curve.
    ///
    /// A degree-`p` B-spline is `C^(p - m)` at an interior knot of multiplicity
    /// `m`, and the worst interior knot governs the curve. With no interior
    /// knots at all the curve is a single polynomial piece and so genuinely
    /// smooth to every order.
    ///
    /// Higher orders than `C2` report as `C2`, which is the highest
    /// [`Continuity`] names short of `CInfinity`. Reporting `CInfinity` for a
    /// merely-`C3` curve would be a claim that is false, and the distinction
    /// above `C2` is not one any algorithm here asks about.
    fn continuity(&self) -> Continuity {
        let degree = self.knots.degree();
        let (a, b) = self.knots.domain();
        let worst = self
            .knots
            .distinct()
            .into_iter()
            .filter(|(v, _)| *v > a && *v < b)
            .map(|(_, m)| m)
            .max();
        match worst {
            None => Continuity::CInfinity,
            Some(m) => match degree.saturating_sub(m) {
                0 => Continuity::C0,
                1 => Continuity::C1,
                _ => Continuity::C2,
            },
        }
    }

    /// A periodic spline closes by construction. Any other where its ends
    /// meet, which its end control points say only for clamped knots.
    fn is_closed(&self, tol: Tolerances) -> bool {
        self.periodic
            || match (self.start(tol), self.end(tol)) {
                (Ok(a), Ok(b)) => a.is_equal(b, tol),
                _ => false,
            }
    }

    fn is_periodic(&self) -> bool {
        self.periodic
    }
}

impl Curve3d for TrimmedCurve {
    fn domain(&self) -> (f64, f64) {
        self.domain
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        self.basis.point_at(self.basis_parameter(u, tol)?, tol)
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        let d = self.basis.d1_at(self.basis_parameter(u, tol)?, tol)?;
        Ok(if self.reversed { -d } else { d })
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        let t = self.basis_parameter(u, tol)?;
        let mut out = self.basis.derivatives_at(t, n, tol)?;
        if self.reversed {
            // Chain rule for u -> (s + e - u): each order picks up a factor of
            // -1, so odd orders flip.
            for (order, d) in out.iter_mut().enumerate() {
                if order % 2 == 1 {
                    *d = -*d;
                }
            }
        }
        Ok(out)
    }

    fn kind(&self) -> CurveKind {
        CurveKind::Trimmed
    }

    fn continuity(&self) -> Continuity {
        self.basis.continuity()
    }

    fn is_closed(&self, tol: Tolerances) -> bool {
        match (self.start(tol), self.end(tol)) {
            (Ok(a), Ok(b)) => a.is_equal(b, tol),
            _ => false,
        }
    }

    fn is_periodic(&self) -> bool {
        false
    }
}

/// Dispatch a method across every curve variant.
macro_rules! dispatch {
    ($self:ident, $c:ident => $body:expr) => {
        match $self {
            Self::Line($c) => $body,
            Self::Circle($c) => $body,
            Self::Ellipse($c) => $body,
            Self::Hyperbola($c) => $body,
            Self::Parabola($c) => $body,
            Self::BSpline($c) => $body,
            Self::Helix($c) => $body,
            Self::Trimmed($c) => $body,
            Self::Offset($c) => $body,
            Self::OnSurface($c) => $body,
        }
    };
}

impl Curve3d for Curve {
    fn domain(&self) -> (f64, f64) {
        dispatch!(self, c => c.domain())
    }

    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point> {
        dispatch!(self, c => c.point_at(u, tol))
    }

    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector> {
        dispatch!(self, c => c.d1_at(u, tol))
    }

    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>> {
        dispatch!(self, c => c.derivatives_at(u, n, tol))
    }

    fn kind(&self) -> CurveKind {
        dispatch!(self, c => c.kind())
    }

    fn continuity(&self) -> Continuity {
        dispatch!(self, c => c.continuity())
    }

    fn is_closed(&self, tol: Tolerances) -> bool {
        dispatch!(self, c => c.is_closed(tol))
    }

    fn is_periodic(&self) -> bool {
        dispatch!(self, c => c.is_periodic())
    }
}

impl Transformable for Curve {
    fn transformed(&self, t: &Transform, tol: Tolerances) -> OgeomResult<Self> {
        if t.kind() == ogeom_math::TransformKind::Identity {
            return Ok(self.clone());
        }
        Ok(match self {
            Self::Line(c) => Self::Line(LineCurve {
                axis: Axis::new(
                    t.apply(c.axis.location),
                    t.apply_direction(c.axis.direction, tol)?,
                ),
                // The parameter is a length, so a scaling rescales the domain
                // with it. Otherwise the trimmed extent would silently change.
                domain: (
                    c.domain.0 * t.scale_factor().abs(),
                    c.domain.1 * t.scale_factor().abs(),
                ),
            }),
            // No parameter flip for a mirror: transforming the frame already
            // carries its x and y axes with it, so evaluating at the same angle
            // lands on the transformed point. Flipping as well would reverse
            // the arc twice.
            Self::Circle(c) => Self::Circle(CircleCurve {
                circle: c.circle.transformed(t, tol)?,
                ..*c
            }),
            Self::Ellipse(c) => Self::Ellipse(EllipseCurve {
                ellipse: c.ellipse.transformed(t, tol)?,
                ..*c
            }),
            Self::Hyperbola(c) => Self::Hyperbola(HyperbolaCurve {
                hyperbola: c.hyperbola.transformed(t, tol)?,
                ..*c
            }),
            // The parameter runs along the axis frame's `y`, a length, and a
            // scaling rescales it as it does a line's.
            Self::Parabola(c) => Self::Parabola(ParabolaCurve {
                parabola: c.parabola.transformed(t, tol)?,
                domain: (
                    c.domain.0 * t.scale_factor().abs(),
                    c.domain.1 * t.scale_factor().abs(),
                ),
                ..*c
            }),
            Self::BSpline(c) => {
                let control = c
                    .control
                    .iter()
                    .map(|w| Weighted::new(t.apply(w.point()), w.weight, tol))
                    .collect::<OgeomResult<Vec<_>>>()?;
                Self::BSpline(BSplineCurve {
                    control,
                    knots: c.knots.clone(),
                    rational: c.rational,
                    periodic: c.periodic,
                })
            }
            Self::Helix(c) => Self::Helix(HelixCurve {
                frame: t.apply_frame(&c.frame, tol)?,
                radius: c.radius * t.scale_factor().abs(),
                pitch: c.pitch * t.scale_factor().abs(),
                taper: c.taper * t.scale_factor().abs(),
                ..*c
            }),
            // The offset runs along `tangent x reference`, which a
            // reflection carries to the other side: the distance goes with it.
            Self::Offset(c) => Self::Offset(Box::new(OffsetCurve {
                basis: c.basis.transformed(t, tol)?,
                distance: c.distance
                    * t.scale_factor().abs()
                    * if t.preserves_handedness() { 1.0 } else { -1.0 },
                reference: t.apply_direction(c.reference, tol)?,
            })),
            // The surface's length directions rescale under a scaling (a
            // plane's both, a cylinder's height), and the chart curve
            // stretches with them so it names the same points' images.
            Self::OnSurface(c) => {
                use crate::traits::Surface as _;
                let surface = c.surface.transformed(t, tol)?;
                let (from, to) = (c.surface.domain(), surface.domain());
                let ratio = |a: (f64, f64), b: (f64, f64)| {
                    let span = a.1 - a.0;
                    if span.is_finite() && span != 0.0 && (b.1 - b.0).is_finite() {
                        (b.1 - b.0) / span
                    } else {
                        1.0
                    }
                };
                Self::OnSurface(Box::new(CurveOnSurface {
                    pcurve: c
                        .pcurve
                        .stretched(ratio(from.0, to.0), ratio(from.1, to.1), tol)?,
                    surface,
                }))
            }
            Self::Trimmed(c) => {
                let basis = c.basis.transformed(t, tol)?;
                // The trim lives in the basis's parameter and moves with it:
                // a line's is a length and rescales, and so does anything
                // built on one (a trim of a trim of a line). An angle or a
                // spline parameter does not, and the map is the identity.
                let (from, to) = (c.basis.domain(), basis.domain());
                Self::Trimmed(Box::new(TrimmedCurve {
                    domain: (
                        crate::surface::carried(c.domain.0, from, to),
                        crate::surface::carried(c.domain.1, from, to),
                    ),
                    basis,
                    reversed: c.reversed,
                }))
            }
        })
    }
}

impl Reversible for Curve {
    fn reversed(&self) -> Self {
        match self {
            Self::Line(c) => Self::Line(LineCurve {
                axis: Axis::new(
                    c.axis.point_at(c.domain.0 + c.domain.1),
                    c.axis.direction.reversed(),
                ),
                domain: c.domain,
            }),
            Self::Circle(c) => Self::Circle(CircleCurve {
                reversed: !c.reversed,
                ..*c
            }),
            Self::Ellipse(c) => Self::Ellipse(EllipseCurve {
                reversed: !c.reversed,
                ..*c
            }),
            Self::Hyperbola(c) => Self::Hyperbola(HyperbolaCurve {
                reversed: !c.reversed,
                ..*c
            }),
            Self::Parabola(c) => Self::Parabola(ParabolaCurve {
                reversed: !c.reversed,
                ..*c
            }),
            Self::BSpline(c) => {
                let (knots, control) = bspline::reverse(&c.knots, &c.control);
                Self::BSpline(BSplineCurve {
                    knots,
                    control,
                    rational: c.rational,
                    periodic: c.periodic,
                })
            }
            Self::Helix(c) => Self::Helix(HelixCurve {
                reversed: !c.reversed,
                ..*c
            }),
            // Reversing flips the tangent and with it the offset direction,
            // so the distance negates to keep the same point set.
            Self::Offset(c) => Self::Offset(Box::new(OffsetCurve {
                basis: c.basis.reversed(),
                distance: -c.distance,
                reference: c.reference,
            })),
            Self::OnSurface(c) => Self::OnSurface(Box::new(CurveOnSurface {
                pcurve: c.pcurve.reversed(),
                surface: c.surface.clone(),
            })),
            // A flag rather than reversing the basis: mirroring the trim range
            // within the basis domain would move this curve's own domain, and
            // trimming ranges held elsewhere refer to it.
            Self::Trimmed(c) => Self::Trimmed(Box::new(TrimmedCurve {
                reversed: !c.reversed,
                ..(**c).clone()
            })),
        }
    }
}

impl From<LineCurve> for Curve {
    fn from(c: LineCurve) -> Self {
        Self::Line(c)
    }
}
impl From<CircleCurve> for Curve {
    fn from(c: CircleCurve) -> Self {
        Self::Circle(c)
    }
}
impl From<EllipseCurve> for Curve {
    fn from(c: EllipseCurve) -> Self {
        Self::Ellipse(c)
    }
}
impl From<HelixCurve> for Curve {
    fn from(c: HelixCurve) -> Self {
        Self::Helix(c)
    }
}

impl From<HyperbolaCurve> for Curve {
    fn from(c: HyperbolaCurve) -> Self {
        Self::Hyperbola(c)
    }
}
impl From<ParabolaCurve> for Curve {
    fn from(c: ParabolaCurve) -> Self {
        Self::Parabola(c)
    }
}
impl From<BSplineCurve> for Curve {
    fn from(c: BSplineCurve) -> Self {
        Self::BSpline(c)
    }
}
impl From<TrimmedCurve> for Curve {
    fn from(c: TrimmedCurve) -> Self {
        Self::Trimmed(Box::new(c))
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod reseam_tests {
    use super::*;
    use ogeom_core::Tolerances;
    use ogeom_math::{KnotVector, Point};

    #[test]
    fn a_reseamed_closed_curve_is_the_same_curve_from_a_new_start() {
        let tol = Tolerances::millimetres();
        // A closed cubic: a ring of control points ending where it began.
        let ring = [
            Point::new(1.0, 0.0, 0.0),
            Point::new(1.0, 1.0, 0.5),
            Point::new(-1.0, 1.0, 0.0),
            Point::new(-1.0, -1.0, -0.5),
            Point::new(1.0, -1.0, 0.0),
            Point::new(1.0, 0.0, 0.0),
        ];
        let knots = KnotVector::clamped_uniform(3, ring.len()).unwrap();
        let curve = BSplineCurve::new(knots, ring.to_vec(), tol).unwrap();
        let (start, end) = curve.domain();
        let seam = 0.35;
        let moved = curve.reseamed_at(seam, tol).unwrap();
        let (new_start, new_end) = moved.domain();
        assert!((new_start - seam).abs() < 1e-12, "begins at the new seam");
        assert!(
            ((new_end - new_start) - (end - start)).abs() < 1e-12,
            "keeps its length"
        );
        for i in 0..=40 {
            let s = (end - start) * f64::from(i) / 40.0;
            let old_u = if seam + s <= end {
                seam + s
            } else {
                seam + s - (end - start)
            };
            let a = curve.point_at(old_u, tol).unwrap();
            let b = moved.point_at(new_start + s, tol).unwrap();
            assert!(a.is_equal(b, tol), "at {s} along: {a:?} against {b:?}");
        }
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use ogeom_math::{Direction, Frame};

    const T: Tolerances = Tolerances::millimetres();

    /// A reversed circle or ellipse spelled forward is the same point at
    /// every parameter.
    #[test]
    fn a_reversed_conic_spelled_forward_is_the_same_curve() {
        let frame = Frame::new(Point::new(1.0, 2.0, 3.0), Direction::Y, Direction::Z, T).unwrap();
        let curves = [
            Curve::Circle(CircleCurve::new(Circle::new(frame, 2.0, T).unwrap())),
            Curve::Ellipse(EllipseCurve::new(Ellipse::new(frame, 5.0, 3.0, T).unwrap())),
        ];
        for reversed in curves.map(|c| c.reversed()) {
            let (forward, turned) = match &reversed {
                Curve::Circle(c) => {
                    let f = c.forward(T).unwrap();
                    (f.is_reversed(), Curve::Circle(f))
                }
                Curve::Ellipse(e) => {
                    let f = e.forward(T).unwrap();
                    (f.is_reversed(), Curve::Ellipse(f))
                }
                _ => unreachable!(),
            };
            assert!(!forward);
            for i in 0..8 {
                let u = 0.3 + f64::from(i) * 0.8;
                let (a, b) = (
                    reversed.point_at(u, T).unwrap(),
                    turned.point_at(u, T).unwrap(),
                );
                assert!(a.distance(b) < 1e-12, "{u}: {a:?} {b:?}");
            }
        }
    }

    /// A spline whose ends meet is closed, periodic or unclamped, and the
    /// surface it sweeps closes with it.
    #[test]
    fn a_periodic_spline_is_closed() {
        let ring = [
            Point::new(0.0, 0.0, 0.0),
            Point::new(10.0, 0.0, 1.0),
            Point::new(10.0, 10.0, 2.0),
            Point::new(0.0, 10.0, 0.0),
        ];
        let c: Curve = BSplineCurve::periodic(&ring, 3, T).unwrap().into();
        assert!(c.is_periodic() && c.is_closed(T));
        let sweep: crate::surface::SurfaceGeometry =
            crate::surface::ExtrusionSurface::new(c, Direction::Z, 10.0)
                .unwrap()
                .into();
        assert!(crate::traits::Surface::is_closed_u(&sweep, T));
    }

    /// A conical helix's arc length is the integral of its speed, which
    /// grows with its radius: the sum of fine chords agrees with it.
    #[test]
    fn a_conical_helix_measures_its_length() {
        for (taper, pitch) in [(0.0, 2.0), (3.0, 2.0), (-1.5, 0.5), (4.0, 0.1)] {
            let helix = HelixCurve::conical(Frame::WORLD, 10.0, pitch, taper, 0.0, 12.0).unwrap();
            let steps = 200_000;
            let mut chords = 0.0;
            let mut last = helix.point_at(0.0, T).unwrap();
            for k in 1..=steps {
                let next = helix
                    .point_at(12.0 * f64::from(k) / f64::from(steps), T)
                    .unwrap();
                chords += next.distance(last);
                last = next;
            }
            let exact = helix.arc_length(0.0, 12.0);
            assert!(
                (exact - chords).abs() < 1e-6 * exact,
                "{taper} {pitch}: {exact} against {chords}"
            );
        }
    }

    /// A scaled curve names, at every parameter, the image of the point it
    /// named: a parabola, whose parameter is a length; a trim of a trim of a
    /// line and a trim of an offset line in a chart; and a curve on a plane
    /// or on a cylinder, whose chart curve stretches with the surface's
    /// length directions.
    #[test]
    fn a_scaled_curve_is_the_image_of_the_original() {
        use crate::curve2d::{Line2d, Offset2d, PlanarCurve, Trimmed2d};
        use crate::surface::{CylinderSurface, PlaneSurface, SurfaceGeometry};
        use crate::traits::Curve2d as _;
        use ogeom_math::{Cylinder, Parabola, Plane, Point2, Transform, Transform2};
        let scaling = Transform::scaling(Point::ORIGIN, 2.0, T).unwrap();
        let parabola: Curve =
            ParabolaCurve::new(Parabola::new(Frame::WORLD, 5.0, T).unwrap(), 10.0)
                .unwrap()
                .into();
        let parabola: Curve = TrimmedCurve::new(parabola, -10.0, 10.0, T).unwrap().into();
        let segment: PlanarCurve =
            Line2d::segment(Point2::new(10.0, 0.0), Point2::new(10.0, 20.0), T)
                .unwrap()
                .into();
        let on_plane = Curve::OnSurface(Box::new(CurveOnSurface::new(
            segment.clone(),
            PlaneSurface::over(Plane::XY, (-100.0, 100.0), (-100.0, 100.0))
                .unwrap()
                .into(),
        )));
        let slant: PlanarCurve = Line2d::segment(Point2::new(0.0, 1.0), Point2::new(3.0, 4.0), T)
            .unwrap()
            .into();
        let drum: SurfaceGeometry =
            CylinderSurface::new(Cylinder::new(Frame::WORLD, 5.0, T).unwrap(), (-10.0, 10.0))
                .unwrap()
                .into();
        let on_drum = Curve::OnSurface(Box::new(CurveOnSurface::new(slant, drum)));
        for curve in [parabola, on_plane, on_drum] {
            let moved = curve.transformed(&scaling, T).unwrap();
            let (a, b) = curve.domain();
            let (c, d) = moved.domain();
            for k in 0..=8 {
                let f = f64::from(k) / 8.0;
                let want = scaling.apply(curve.point_at(a + (b - a) * f, T).unwrap());
                let got = moved.point_at(c + (d - c) * f, T).unwrap();
                assert!(got.distance(want) < 1e-9, "{got:?} against {want:?}");
            }
        }
        let flat = Transform2::scaling(Point2::ORIGIN, 2.0, T).unwrap();
        let line: PlanarCurve = Line2d::segment(Point2::new(0.0, 0.0), Point2::new(10.0, 0.0), T)
            .unwrap()
            .into();
        let nested: PlanarCurve = Trimmed2d::new(
            Trimmed2d::new(line.clone(), 1.0, 9.0, T).unwrap().into(),
            3.0,
            7.0,
            T,
        )
        .unwrap()
        .into();
        let beside: PlanarCurve = Trimmed2d::new(
            PlanarCurve::Offset(Box::new(Offset2d::new(line, 1.0).unwrap())),
            3.0,
            7.0,
            T,
        )
        .unwrap()
        .into();
        for curve in [nested, beside] {
            let moved = curve.transformed(&flat, T).unwrap();
            for (at, on) in [
                (curve.domain().0, moved.domain().0),
                (curve.domain().1, moved.domain().1),
            ] {
                let want = flat.apply(curve.point_at(at, T).unwrap());
                assert!(moved.point_at(on, T).unwrap().distance(want) < 1e-9);
            }
        }
    }

    /// An offset curve carried by a mirror or a point mirror names, at every
    /// parameter, the image of the point it named.
    #[test]
    fn a_mirrored_offset_curve_is_the_image_of_the_original() {
        use ogeom_math::{Circle, Transform};
        let circle: Curve = CircleCurve::new(Circle::new(Frame::WORLD, 20.0, T).unwrap()).into();
        let offset: Curve = Curve::Offset(Box::new(
            OffsetCurve::new(circle, 3.0, Direction::Z).unwrap(),
        ));
        for t in [
            Transform::plane_mirror(Point::ORIGIN, Direction::X),
            Transform::point_mirror(Point::new(1.0, 2.0, 3.0)),
            Transform::scaling(Point::ORIGIN, -2.0, T).unwrap(),
        ] {
            let moved = offset.transformed(&t, T).unwrap();
            for u in [0.0, 1.0, 2.5, 4.0] {
                let want = t.apply(offset.point_at(u, T).unwrap());
                let got = moved.point_at(u, T).unwrap();
                assert!(
                    got.distance(want) < 1e-9,
                    "{t:?} at {u}: {got:?} against {want:?}"
                );
            }
        }
    }

    fn tilted() -> Frame {
        Frame::new(
            Point::new(1.0, -2.0, 3.0),
            Direction::from_coords(1.0, 2.0, 3.0, T).unwrap(),
            Direction::X,
            T,
        )
        .unwrap()
    }

    fn every_curve() -> Vec<Curve> {
        let spline = {
            let control = vec![
                Point::new(0.0, 0.0, 0.0),
                Point::new(1.0, 2.0, 0.0),
                Point::new(3.0, 1.0, 1.0),
                Point::new(5.0, 0.0, 2.0),
                Point::new(6.0, -1.0, 0.0),
            ];
            let knots = KnotVector::clamped_uniform(3, control.len()).unwrap();
            BSplineCurve::new(knots, control, T).unwrap()
        };
        vec![
            LineCurve::segment(Point::ORIGIN, Point::new(3.0, 4.0, 0.0), T)
                .unwrap()
                .into(),
            CircleCurve::new(Circle::new(tilted(), 2.0, T).unwrap()).into(),
            EllipseCurve::new(Ellipse::new(tilted(), 5.0, 3.0, T).unwrap()).into(),
            HyperbolaCurve::new(Hyperbola::new(tilted(), 3.0, 4.0, T).unwrap(), 1.5)
                .unwrap()
                .into(),
            ParabolaCurve::new(Parabola::new(tilted(), 2.0, T).unwrap(), 4.0)
                .unwrap()
                .into(),
            spline.clone().into(),
            HelixCurve::new(tilted(), 2.5, 1.25, 2.0).unwrap().into(),
            TrimmedCurve::new(spline.into(), 0.2, 0.8, T)
                .unwrap()
                .into(),
        ]
    }

    #[test]
    fn a_helix_rises_one_pitch_per_turn_and_knows_its_length() {
        let helix = HelixCurve::new(Frame::WORLD, 3.0, 2.0, 2.0).unwrap();
        let tau = core::f64::consts::TAU;
        let start = helix.point_at(0.0, T).unwrap();
        let after_one_turn = helix.point_at(tau, T).unwrap();
        assert_relative_eq!(start.x, 3.0);
        assert_relative_eq!(after_one_turn.x, 3.0, epsilon = 1e-12);
        assert_relative_eq!(after_one_turn.y, 0.0, epsilon = 1e-12);
        assert_relative_eq!(after_one_turn.z - start.z, 2.0, epsilon = 1e-12);

        // Closed-form length: constant speed times swept angle, checked
        // against a fine chordal sum.
        let exact = helix.arc_length(0.0, 2.0 * tau);
        assert_relative_eq!(exact, 2.0 * tau * 3.0f64.hypot(2.0 / tau), epsilon = 1e-12);
        let mut chords = 0.0;
        let n = 20_000;
        for i in 0..n {
            let a = 2.0 * tau * f64::from(i) / f64::from(n);
            let b = 2.0 * tau * f64::from(i + 1) / f64::from(n);
            chords += helix
                .point_at(a, T)
                .unwrap()
                .distance(helix.point_at(b, T).unwrap());
        }
        assert!((exact - chords) / exact < 1e-6, "{exact} vs {chords}");

        // A negative pitch winds the other hand: same rise magnitude, the
        // quarter-turn point mirrored through the xz-plane... the y stays,
        // the z descends.
        let left = HelixCurve::new(Frame::WORLD, 3.0, -2.0, 2.0).unwrap();
        let q = left.point_at(tau / 4.0, T).unwrap();
        assert_relative_eq!(q.y, 3.0, epsilon = 1e-12);
        assert!(q.z < 0.0);
    }

    #[test]
    fn a_reversed_helix_swaps_its_ends_and_flips_its_tangent() {
        let helix: Curve = HelixCurve::new(tilted(), 2.0, 1.0, 1.5).unwrap().into();
        let (lo, hi) = helix.domain();
        let back = helix.reversed();
        assert_relative_eq!(
            helix
                .point_at(lo, T)
                .unwrap()
                .distance(back.point_at(hi, T).unwrap()),
            0.0,
            epsilon = 1e-12
        );
        let d_fwd = helix.d1_at(f64::midpoint(lo, hi), T).unwrap();
        let d_back = back.d1_at(f64::midpoint(lo, hi), T).unwrap();
        assert_relative_eq!((d_fwd + d_back).magnitude(), 0.0, epsilon = 1e-12);
    }

    #[test]
    fn an_offset_circle_is_the_larger_circle() {
        // The offset of a circle perpendicular to its own axis is the
        // concentric circle: for the counterclockwise traversal, tangent x z
        // points radially outward, so a positive distance grows the radius.
        let circle: Curve = CircleCurve::new(Circle::new(tilted(), 2.0, T).unwrap()).into();
        let bigger = Circle::new(tilted(), 3.0, T).unwrap();
        let offset = OffsetCurve::new(circle, 1.0, tilted().z()).unwrap();
        for i in 0..8 {
            let t = core::f64::consts::TAU * f64::from(i) / 8.0;
            let p = offset.point_at(t, T).unwrap();
            assert_relative_eq!(p.distance(bigger.centre()), 3.0, epsilon = 1e-12);
        }
        // The exact first derivative agrees with differencing the points.
        let h = 1e-6;
        let d = offset.d1_at(1.0, T).unwrap();
        let fd = (offset.point_at(1.0 + h, T).unwrap() - offset.point_at(1.0 - h, T).unwrap())
            / (2.0 * h);
        assert_relative_eq!((d - fd).magnitude(), 0.0, epsilon = 1e-5);
        // The second derivative is refused by name, not differenced quietly.
        assert!(offset.derivatives_at(1.0, 2, T).is_err());
    }

    #[test]
    fn a_sloped_line_on_a_cylinder_chart_is_a_helix() {
        use crate::curve2d::{Line2d, PlanarCurve};
        use crate::surface::{CylinderSurface, SurfaceGeometry};
        use ogeom_math::{Cylinder, Point2};

        // The pcurve u = t, v = pitch·t/2π on a cylinder chart lifts to
        // exactly the helix with that pitch: two constructions, no shared
        // code path, one curve.
        let radius = 3.0;
        let pitch = 2.0;
        let tau = core::f64::consts::TAU;
        let cylinder = SurfaceGeometry::Cylinder(
            CylinderSurface::new(Cylinder::new(Frame::WORLD, radius, T).unwrap(), (-1.0, 5.0))
                .unwrap(),
        );
        let rise = pitch / tau;
        let slope = (1.0 + rise * rise).sqrt();
        let line = Line2d::segment(Point2::new(0.0, 0.0), Point2::new(tau, pitch), T).unwrap();
        let on_surface = CurveOnSurface::new(PlanarCurve::Line(line), cylinder);
        let helix = HelixCurve::new(Frame::WORLD, radius, pitch, 1.0).unwrap();
        for i in 0..=8 {
            let angle = tau * f64::from(i) / 8.0;
            // The line is parameterized by 2D arc length; the helix by angle.
            let lifted = on_surface.point_at(angle * slope, T).unwrap();
            let wound = helix.point_at(angle, T).unwrap();
            assert_relative_eq!(lifted.distance(wound), 0.0, epsilon = 1e-9);
        }
        // Exact second derivative through the chain rule, checked against
        // differencing the exact first.
        let h = 1e-6;
        let d2 = on_surface.derivatives_at(2.0, 2, T).unwrap()[2];
        let fd = (on_surface.d1_at(2.0 + h, T).unwrap() - on_surface.d1_at(2.0 - h, T).unwrap())
            / (2.0 * h);
        assert_relative_eq!((d2 - fd).magnitude(), 0.0, epsilon = 1e-5);
    }

    #[test]
    fn a_helix_has_no_exact_spline_and_says_so() {
        let helix: Curve = HelixCurve::new(Frame::WORLD, 1.0, 1.0, 1.0).unwrap().into();
        assert!(helix.to_bspline(T).is_err());
    }

    /// Sample a curve evenly across its domain, avoiding the exact ends.
    fn interior(c: &Curve, n: usize) -> Vec<f64> {
        let (a, b) = c.domain();
        (1..n)
            .map(|i| {
                #[allow(clippy::cast_precision_loss)]
                let t = i as f64 / n as f64;
                a + (b - a) * t
            })
            .collect()
    }

    #[test]
    fn every_order_agrees_with_the_difference_of_the_one_below() {
        let h = 1e-5;
        let mut curves = every_curve();
        curves.push(
            HelixCurve::conical(tilted(), 4.0, 1.5, 2.0, 0.0, 6.0)
                .unwrap()
                .into(),
        );
        let reversed: Vec<Curve> = curves.iter().map(Reversible::reversed).collect();
        curves.extend(reversed);
        // The closed forms: a spline's orders jump at its knots.
        curves.retain(|c| !matches!(c.kind(), CurveKind::BSpline | CurveKind::Trimmed));
        for c in curves {
            for u in interior(&c, 6) {
                let at = c.derivatives_at(u, 5, T).unwrap();
                let below = c.derivatives_at(u + h, 4, T).unwrap();
                let above = c.derivatives_at(u - h, 4, T).unwrap();
                for order in 2..=5 {
                    let numeric = (below[order - 1] - above[order - 1]) * (1.0 / (2.0 * h));
                    let scale = numeric.magnitude().max(1.0);
                    assert!(
                        (at[order] - numeric).magnitude() <= 1e-5 * scale,
                        "{:?} order {order} at {u}: {:?} against {numeric:?}",
                        c.kind(),
                        at[order]
                    );
                }
            }
        }
    }

    #[test]
    fn every_curves_derivative_agrees_with_finite_differences() {
        let h = 1e-6;
        for c in every_curve() {
            for u in interior(&c, 8) {
                let d1 = c.d1_at(u, T).unwrap();
                let numeric = (c.point_at(u + h, T).unwrap() - c.point_at(u - h, T).unwrap())
                    * (1.0 / (2.0 * h));
                let scale = numeric.magnitude().max(1.0);
                assert!(
                    (d1 - numeric).magnitude() <= 1e-5 * scale,
                    "{:?} at {u}: {d1:?} vs {numeric:?}",
                    c.kind()
                );
            }
        }
    }

    #[test]
    fn derivatives_at_zero_returns_the_point_itself() {
        for c in every_curve() {
            for u in interior(&c, 4) {
                let d = c.derivatives_at(u, 0, T).unwrap();
                assert_eq!(d.len(), 1);
                assert!(Point::from_vector(d[0]).is_equal(c.point_at(u, T).unwrap(), T));
            }
        }
    }

    #[test]
    fn out_of_domain_parameters_are_refused_for_non_periodic_curves() {
        for c in every_curve() {
            let (a, b) = c.domain();
            if c.is_periodic() {
                // A periodic curve accepts anything and wraps it.
                assert!(c.point_at(b + 1.0, T).is_ok());
                assert!(c.point_at(a - 1.0, T).is_ok());
            } else {
                assert!(c.point_at(b + 1.0, T).is_err(), "{:?}", c.kind());
                assert!(c.point_at(a - 1.0, T).is_err(), "{:?}", c.kind());
            }
        }
    }

    #[test]
    fn a_periodic_curve_wraps_to_the_same_point() {
        let c: Curve = CircleCurve::new(Circle::new(tilted(), 2.0, T).unwrap()).into();
        assert!(c.is_periodic() && c.is_closed(T));
        let base = c.point_at(0.7, T).unwrap();
        for k in [-2.0_f64, -1.0, 1.0, 3.0] {
            let wrapped = c
                .point_at(k.mul_add(core::f64::consts::TAU, 0.7), T)
                .unwrap();
            assert!(base.is_equal(wrapped, T), "wrap by {k} moved the point");
        }
    }

    #[test]
    fn reversal_traverses_the_same_points_backwards() {
        for c in every_curve() {
            let r = c.reversed();
            let (a, b) = c.domain();
            assert_eq!(r.domain(), (a, b), "{:?} changed its domain", c.kind());
            for i in 0..=8 {
                let t = f64::from(i) / 8.0;
                let forward = c.point_at(a + (b - a) * t, T).unwrap();
                let backward = r.point_at(a + (b - a) * (1.0 - t), T).unwrap();
                assert!(
                    forward.is_equal(backward, T),
                    "{:?} at t = {t}: {forward:?} vs {backward:?}",
                    c.kind()
                );
            }
        }
    }

    #[test]
    fn reversing_twice_is_the_identity() {
        for c in every_curve() {
            let twice = c.reversed().reversed();
            for u in interior(&c, 8) {
                assert!(
                    c.point_at(u, T)
                        .unwrap()
                        .is_equal(twice.point_at(u, T).unwrap(), T),
                    "{:?}",
                    c.kind()
                );
            }
        }
    }

    #[test]
    fn a_reversed_curves_tangent_points_the_other_way() {
        for c in every_curve() {
            let r = c.reversed();
            let (a, b) = c.domain();
            let u = a + (b - a) * 0.4;
            let forward = c.tangent_at(u, T).unwrap();
            let backward = r.tangent_at(mirror(u, a, b), T).unwrap();
            assert!(forward.is_opposite(backward, T), "{:?}", c.kind());
        }
    }

    #[test]
    fn transforms_move_curves_and_preserve_their_shape() {
        let t =
            Transform::rotation(Axis::X, 0.7) * Transform::translation(Vector::new(1.0, 2.0, 3.0));
        for c in every_curve() {
            let moved = c.transformed(&t, T).unwrap();
            assert_eq!(moved.kind(), c.kind());
            for u in interior(&c, 8) {
                let expected = t.apply(c.point_at(u, T).unwrap());
                assert!(
                    moved.point_at(u, T).unwrap().is_equal(expected, T),
                    "{:?} at {u}",
                    c.kind()
                );
            }
        }
    }

    #[test]
    fn a_scaling_rescales_a_lines_arc_length_domain() {
        // A line's parameter is a length, so scaling must rescale the domain or
        // the segment silently changes extent.
        let line = LineCurve::segment(Point::ORIGIN, Point::new(3.0, 4.0, 0.0), T).unwrap();
        let c: Curve = line.into();
        assert_eq!(c.domain(), (0.0, 5.0));
        let scaled = c
            .transformed(&Transform::scaling(Point::ORIGIN, 2.0, T).unwrap(), T)
            .unwrap();
        assert_eq!(scaled.domain(), (0.0, 10.0));
        assert!(
            scaled
                .end(T)
                .unwrap()
                .is_equal(Point::new(6.0, 8.0, 0.0), T)
        );
    }

    #[test]
    fn mirroring_a_circle_moves_every_point_by_the_mirror() {
        // The frame carries its own axes through the transform, so evaluating
        // the mirrored circle at a parameter lands exactly where the mirror
        // sends the original point. No extra parameter flip is involved.
        let c: Curve = CircleCurve::new(Circle::new(Frame::WORLD, 2.0, T).unwrap()).into();
        let m = Transform::plane_mirror(Point::ORIGIN, Direction::X);
        let mirrored = c.transformed(&m, T).unwrap();
        for u in interior(&c, 8) {
            let expected = m.apply(c.point_at(u, T).unwrap());
            assert!(
                mirrored.point_at(u, T).unwrap().is_equal(expected, T),
                "at {u}"
            );
        }
    }

    #[test]
    fn a_line_segments_parameter_is_arc_length() {
        let line = LineCurve::segment(Point::ORIGIN, Point::new(3.0, 4.0, 0.0), T).unwrap();
        assert_eq!(line.domain(), (0.0, 5.0));
        assert!(line.point_at(0.0, T).unwrap().is_equal(Point::ORIGIN, T));
        assert!(
            line.point_at(5.0, T)
                .unwrap()
                .is_equal(Point::new(3.0, 4.0, 0.0), T)
        );
        assert!(
            line.point_at(2.5, T)
                .unwrap()
                .is_equal(Point::new(1.5, 2.0, 0.0), T)
        );
        assert_relative_eq!(
            line.d1_at(1.0, T).unwrap().magnitude(),
            1.0,
            epsilon = 1e-15
        );
    }

    #[test]
    fn degenerate_constructions_are_refused() {
        assert!(LineCurve::segment(Point::ORIGIN, Point::ORIGIN, T).is_err());
        assert!(LineCurve::over(Axis::X, 1.0, 1.0).is_err());
        assert!(LineCurve::over(Axis::X, 0.0, f64::NAN).is_err());
        assert!(HyperbolaCurve::new(Hyperbola::new(tilted(), 1.0, 1.0, T).unwrap(), 0.0).is_err());
        assert!(ParabolaCurve::new(Parabola::new(tilted(), 1.0, T).unwrap(), -1.0).is_err());
    }

    #[test]
    fn trimming_is_bounds_checked() {
        let base: Curve = LineCurve::over(Axis::X, 0.0, 10.0).unwrap().into();
        assert!(TrimmedCurve::new(base.clone(), 2.0, 8.0, T).is_ok());
        assert!(
            TrimmedCurve::new(base.clone(), 8.0, 2.0, T).is_err(),
            "empty"
        );
        assert!(
            TrimmedCurve::new(base.clone(), 5.0, 5.0, T).is_err(),
            "empty"
        );
        assert!(
            TrimmedCurve::new(base, -1.0, 5.0, T).is_err(),
            "outside the basis"
        );
    }

    #[test]
    fn a_trimmed_curve_agrees_with_its_basis() {
        let base: Curve = CircleCurve::new(Circle::new(tilted(), 2.0, T).unwrap()).into();
        let trimmed = TrimmedCurve::new(base.clone(), 0.5, 2.0, T).unwrap();
        assert_eq!(trimmed.domain(), (0.5, 2.0));
        for i in 0..=8 {
            let u = 0.5 + 1.5 * (f64::from(i) / 8.0);
            assert!(
                trimmed
                    .point_at(u, T)
                    .unwrap()
                    .is_equal(base.point_at(u, T).unwrap(), T)
            );
        }
        assert!(trimmed.point_at(0.4, T).is_err());
        assert!(trimmed.point_at(2.1, T).is_err());
    }

    #[test]
    fn a_rational_curve_is_recognized_and_a_uniformly_weighted_one_is_not() {
        let knots = KnotVector::clamped_uniform(2, 3).unwrap();
        let points = [
            Point::new(1.0, 0.0, 0.0),
            Point::new(1.0, 1.0, 0.0),
            Point::new(0.0, 1.0, 0.0),
        ];

        let uniform: Vec<_> = points
            .iter()
            .map(|p| Weighted::new(*p, 3.0, T).unwrap())
            .collect();
        assert!(
            !BSplineCurve::rational(knots.clone(), uniform)
                .unwrap()
                .is_rational(),
            "equal weights are polynomial whatever their value"
        );

        let w = core::f64::consts::FRAC_1_SQRT_2;
        let arc: Vec<_> = points
            .iter()
            .zip([1.0, w, 1.0])
            .map(|(p, w)| Weighted::new(*p, w, T).unwrap())
            .collect();
        let c = BSplineCurve::rational(knots, arc).unwrap();
        assert!(c.is_rational());
        // And it is an exact circular arc, which no polynomial curve is.
        for i in 0..=20 {
            let u = f64::from(i) / 20.0;
            assert_relative_eq!(
                c.point_at(u, T).unwrap().to_vector().magnitude(),
                1.0,
                epsilon = 1e-14
            );
        }
    }

    #[test]
    fn spline_continuity_follows_interior_knot_multiplicity() {
        let control = vec![
            Point::ORIGIN,
            Point::new(1.0, 1.0, 0.0),
            Point::new(2.0, 0.0, 0.0),
            Point::new(3.0, 1.0, 0.0),
            Point::new(4.0, 0.0, 0.0),
        ];
        // Degree 3 with simple interior knots is C2 there, not C-infinity,
        // which would be a false claim about a piecewise polynomial.
        let smooth = BSplineCurve::new(
            KnotVector::clamped_uniform(3, control.len()).unwrap(),
            control.clone(),
            T,
        )
        .unwrap();
        assert_eq!(smooth.continuity(), Continuity::C2);

        // A single Bezier piece has no interior knots and is smooth to every
        // order.
        let bezier = BSplineCurve::new(
            KnotVector::clamped_uniform(4, control.len()).unwrap(),
            control.clone(),
            T,
        )
        .unwrap();
        assert_eq!(bezier.continuity(), Continuity::CInfinity);

        // Degree 3, interior multiplicity 2: C1.
        let kinked = BSplineCurve::new(
            KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0], 3).unwrap(),
            {
                let mut c = control.clone();
                c.push(Point::new(5.0, 1.0, 0.0));
                c
            },
            T,
        )
        .unwrap();
        assert_eq!(kinked.continuity(), Continuity::C1);

        // An interior knot at full multiplicity is a corner.
        let corner = BSplineCurve::new(
            KnotVector::new(vec![0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0], 2).unwrap(),
            control,
            T,
        )
        .unwrap();
        assert_eq!(corner.continuity(), Continuity::C0);
    }

    #[test]
    fn spline_refinement_does_not_move_the_curve() {
        let control = vec![
            Point::ORIGIN,
            Point::new(1.0, 2.0, 0.0),
            Point::new(3.0, 1.0, 1.0),
            Point::new(5.0, 0.0, 2.0),
        ];
        let c = BSplineCurve::new(
            KnotVector::clamped_uniform(3, control.len()).unwrap(),
            control,
            T,
        )
        .unwrap();
        let refined = c.with_knot_inserted(0.5, 1, T).unwrap();
        let elevated = c.elevated(T).unwrap();
        assert_eq!(elevated.degree(), 4);
        for i in 0..=20 {
            let u = f64::from(i) / 20.0;
            let base = c.point_at(u, T).unwrap();
            assert!(refined.point_at(u, T).unwrap().is_equal(base, T));
            assert!(elevated.point_at(u, T).unwrap().is_equal(base, T));
        }
    }

    #[test]
    fn curvature_of_a_circle_is_the_reciprocal_of_its_radius() {
        for r in [0.5_f64, 2.0, 50.0] {
            let c: Curve = CircleCurve::new(Circle::new(tilted(), r, T).unwrap()).into();
            assert_relative_eq!(
                c.curvature_at(1.1, T).unwrap(),
                1.0 / r,
                max_relative = 1e-12
            );
        }
        // A line has no curvature.
        let line: Curve = LineCurve::segment(Point::ORIGIN, Point::new(1.0, 1.0, 1.0), T)
            .unwrap()
            .into();
        assert_relative_eq!(line.curvature_at(0.5, T).unwrap(), 0.0);
    }

    #[test]
    fn kinds_are_reported_for_dispatch() {
        let kinds: Vec<_> = every_curve().iter().map(Curve3d::kind).collect();
        assert_eq!(
            kinds,
            vec![
                CurveKind::Line,
                CurveKind::Circle,
                CurveKind::Ellipse,
                CurveKind::Hyperbola,
                CurveKind::Parabola,
                CurveKind::BSpline,
                CurveKind::Helix,
                CurveKind::Trimmed,
            ]
        );
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used, reason = "test code")]
mod conical_tests {
    use super::*;
    use ogeom_math::Vector;

    /// A curve continued to a point ends there, keeps its own run, and
    /// meets the continuation with its tangent and curvature, at either
    /// end, polynomial and rational alike.
    #[test]
    fn a_curve_extended_to_a_point_ends_there_smoothly() {
        use crate::Curve3d as _;
        let cubic = BSplineCurve::rational(
            KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0], 3).unwrap(),
            [
                Point::new(0.0, 0.0, 0.0),
                Point::new(1.0, 2.0, 0.0),
                Point::new(3.0, 2.5, 1.0),
                Point::new(4.0, 0.5, 1.0),
                Point::new(6.0, 1.0, 0.0),
            ]
            .iter()
            .map(|p| Weighted::new(*p, 1.0, T).unwrap())
            .collect(),
        )
        .unwrap();
        let circle = Circle::new(Frame::WORLD, 5.0, T).unwrap();
        let arc = Curve::Trimmed(Box::new(
            TrimmedCurve::new(
                Curve::Circle(CircleCurve::new(circle)),
                0.0,
                core::f64::consts::FRAC_PI_2,
                T,
            )
            .unwrap(),
        ))
        .to_bspline(T)
        .unwrap();
        for curve in [cubic, arc] {
            let (lo, hi) = curve.domain();
            for at_end in [true, false] {
                let target = if at_end {
                    Point::new(8.0, 3.0, -1.0)
                } else {
                    Point::new(-2.0, -1.0, 0.5)
                };
                let longer = curve.extended_to(at_end, target, 2, T).unwrap();
                let (elo, ehi) = longer.domain();
                let end = longer.point_at(if at_end { ehi } else { elo }, T).unwrap();
                assert!(end.distance(target) < 1e-9, "ends at the point: {end:?}");
                for i in 0..=8 {
                    let u = lo + (hi - lo) * f64::from(i) / 8.0;
                    let (was, now) = (
                        curve.point_at(u, T).unwrap(),
                        longer.point_at(u, T).unwrap(),
                    );
                    assert!(was.distance(now) < 1e-9, "the run itself at {u}");
                }
                let join = if at_end { hi } else { lo };
                let step = if at_end { 1e-6 } else { -1e-6 };
                let inside = curve.derivatives_at(join - step, 2, T).unwrap();
                let outside = longer.derivatives_at(join + step, 2, T).unwrap();
                for order in 1..=2 {
                    let gap = (inside[order] - outside[order]).magnitude();
                    let scale = inside[order].magnitude().max(1.0);
                    assert!(gap < scale * 1e-3, "order {order} across the join: {gap}");
                }
            }
        }
    }

    /// A rational quarter circle continued at order two stays on its
    /// circle: the homogeneous polynomial continues as itself. The original
    /// run keeps its points, at either end, and the continuation reaches
    /// the length asked.
    #[test]
    fn a_rational_arc_extended_stays_on_its_circle() {
        use crate::Curve3d as _;
        let circle = Circle::new(Frame::WORLD, 5.0, T).unwrap();
        let arc = TrimmedCurve::new(
            Curve::Circle(CircleCurve::new(circle)),
            0.0,
            core::f64::consts::FRAC_PI_2,
            T,
        )
        .unwrap();
        let spline = Curve::Trimmed(Box::new(arc)).to_bspline(T).unwrap();
        let (lo, hi) = spline.domain();
        for at_end in [true, false] {
            let longer = spline.extended(at_end, 4.0, 2, T).unwrap();
            let (elo, ehi) = longer.domain();
            for i in 0..=8 {
                let u = lo + (hi - lo) * f64::from(i) / 8.0;
                let (was, now) = (
                    spline.point_at(u, T).unwrap(),
                    longer.point_at(u, T).unwrap(),
                );
                assert!(was.distance(now) < 1e-9, "the arc itself at {u}");
            }
            let (from, to) = if at_end { (hi, ehi) } else { (elo, lo) };
            let mut swept = 0.0;
            let mut last = longer.point_at(from, T).unwrap();
            for i in 1..=64 {
                let u = from + (to - from) * f64::from(i) / 64.0;
                let p = longer.point_at(u, T).unwrap();
                assert!(
                    (p.distance(Point::ORIGIN) - 5.0).abs() < 1e-9,
                    "off the circle at {u}: {p:?}"
                );
                swept += last.distance(p);
                last = p;
            }
            assert!((swept - 4.0).abs() < 1e-3, "the length asked: {swept}");
        }
        let ring = BSplineCurve::periodic(
            &[
                Point::new(0.0, 0.0, 0.0),
                Point::new(1.0, 0.0, 0.0),
                Point::new(1.0, 1.0, 0.0),
                Point::new(0.0, 1.0, 0.0),
            ],
            3,
            T,
        )
        .unwrap();
        assert!(ring.extended(true, 1.0, 2, T).is_err(), "a ring has no end");
    }

    const T: Tolerances = Tolerances::millimetres();

    #[test]
    fn a_conical_helix_winds_its_cone_and_speaks_its_derivatives() {
        let tau = core::f64::consts::TAU;
        // Radius 5 at angle zero, growing 2 per turn, rising 3 per turn.
        let helix = HelixCurve::conical(Frame::WORLD, 5.0, 3.0, 2.0, 0.0, 2.0 * tau).unwrap();
        for i in 0..=8 {
            let t = 2.0 * tau * f64::from(i) / 8.0;
            let p = helix.point_at(t, T).unwrap();
            let r = 2.0f64.mul_add(t / tau, 5.0);
            assert!((p.to_vector().dot(Vector::Z) - 3.0 * t / tau).abs() < 1e-9);
            assert!((p.x.hypot(p.y) - r).abs() < 1e-9, "radius at {t}");
        }
        // The derivative against finite differences.
        let t = 1.234;
        let h = 1e-6;
        let d1 = helix.d1_at(t, T).unwrap();
        let fwd = helix.point_at(t + h, T).unwrap();
        let bwd = helix.point_at(t - h, T).unwrap();
        let fd = (fwd - bwd) / (2.0 * h);
        assert!((d1 - fd).magnitude() < 1e-6, "d1 {d1:?} against {fd:?}");
    }

    #[test]
    fn a_helix_past_its_apex_is_refused_by_name() {
        let tau = core::f64::consts::TAU;
        let err = HelixCurve::conical(Frame::WORLD, 1.0, 3.0, -2.0, 0.0, 2.0 * tau).unwrap_err();
        assert!(err.to_string().contains("apex"), "{err}");
    }

    #[test]
    fn a_periodic_bspline_wraps_smoothly_and_says_so() {
        let ring: Vec<Point> = (0..8)
            .map(|i| {
                let a = core::f64::consts::TAU * f64::from(i) / 8.0;
                Point::new(a.cos() * 4.0, a.sin() * 4.0, 0.0)
            })
            .collect();
        let curve = BSplineCurve::periodic(&ring, 3, T).unwrap();
        assert!(Curve3d::is_periodic(&curve));
        let (lo, hi) = Curve3d::domain(&curve);
        // The loop closes with matching tangents across the seam.
        let p_lo = curve.point_at(lo, T).unwrap();
        let p_hi = curve.point_at(hi, T).unwrap();
        assert!(p_lo.distance(p_hi) < 1e-9);
        let d_lo = curve.d1_at(lo, T).unwrap();
        let d_hi = curve.d1_at(hi, T).unwrap();
        assert!((d_lo - d_hi).magnitude() < 1e-9, "C1 across the seam");
        // And a wrapped parameter lands where its image does.
        let inside = curve.point_at(lo + 0.4, T).unwrap();
        let wrapped = curve.point_at(hi + 0.4, T).unwrap();
        assert!(inside.distance(wrapped) < 1e-9);
    }
}

#[cfg(test)]
#[allow(clippy::unwrap_used, clippy::expect_used, reason = "test code")]
mod widened_tests {
    use super::*;
    use ogeom_math::{Direction, Frame, Parabola};

    const T: Tolerances = Tolerances::millimetres();

    /// A line and a forward parabola widened trace the same points at the
    /// same parameters over their old domain, and their own formula past
    /// it; a reversed one, whose parameters mirror across its domain, and
    /// a whole circle are not widened.
    #[test]
    fn a_curve_widens_as_its_own_formula_or_not_at_all() {
        let line: Curve = LineCurve::segment(Point::ORIGIN, Point::new(3.0, 4.0, 0.0), T)
            .unwrap()
            .into();
        let frame = Frame::new(Point::new(1.0, 2.0, 3.0), Direction::Z, Direction::X, T).unwrap();
        let parabola: Curve = ParabolaCurve::new(Parabola::new(frame, 0.7, T).unwrap(), 2.0)
            .unwrap()
            .into();
        for curve in [line, parabola] {
            let (a, b) = curve.domain();
            let wide = curve.widened((a - 1.5, b + 2.5)).expect("widens");
            assert_eq!(wide.domain(), (a - 1.5, b + 2.5));
            for k in 0..=8 {
                let t = a + (b - a) * f64::from(k) / 8.0;
                let gap = wide
                    .point_at(t, T)
                    .unwrap()
                    .distance(curve.point_at(t, T).unwrap());
                assert!(gap <= 1e-12, "{gap} at {t}");
            }
            assert!(curve.widened((a + 0.1, b)).is_none(), "a narrower range");
            // A reversed line is a line along the other way, and widens.
            let reversed = curve.reversed().widened((a - 1.0, b + 1.0));
            assert_eq!(reversed.is_some(), matches!(curve, Curve::Line(_)));
        }
        let circle: Curve =
            CircleCurve::new(ogeom_math::Circle::new(frame, 2.0, T).unwrap()).into();
        let (a, b) = circle.domain();
        assert!(circle.widened((a - 1.0, b + 1.0)).is_none());
    }
}