From 19028f3c02edba50e663df1809c676671877c463 Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 20:26:29 +0200 Subject: [PATCH 1/6] Implement cubic Hermite interpolation with error handling - Added `InterpolationError` enum to represent various errors in interpolation routines. - Introduced linear interpolation functionality in `linear.rs` with error checks for finite values, duplicate abscissae, and out-of-range queries. - Created `QuantityHermiteInterpolable` trait for handling cubic Hermite interpolation over typed quantities in `quantity.rs`. - Implemented cubic Hermite interpolation logic in `CubicHermiteQuantityTable` and related structures. - Added scalar interpolation primitives in `scalar.rs`, including cubic Hermite segment evaluation. - Defined `HermiteInterpolable` trait for combining complete typed values in `traits.rs`. - Updated `lib.rs` to expose new interpolation functionalities. - Created comprehensive tests for interpolation functionalities, covering various scenarios including cubic polynomial reproduction, exact node evaluation, and error handling for out-of-range and duplicate abscissae. --- src/interpolation/cubic_hermite.rs | 263 +++++++++++++++++++ src/interpolation/error.rs | 57 ++++ src/interpolation/linear.rs | 34 +++ src/interpolation/mod.rs | 23 ++ src/interpolation/quantity.rs | 409 +++++++++++++++++++++++++++++ src/interpolation/scalar.rs | 87 ++++++ src/interpolation/traits.rs | 247 +++++++++++++++++ src/lib.rs | 7 + tests/interpolation.rs | 271 +++++++++++++++++++ 9 files changed, 1398 insertions(+) create mode 100644 src/interpolation/cubic_hermite.rs create mode 100644 src/interpolation/error.rs create mode 100644 src/interpolation/linear.rs create mode 100644 src/interpolation/mod.rs create mode 100644 src/interpolation/quantity.rs create mode 100644 src/interpolation/scalar.rs create mode 100644 src/interpolation/traits.rs create mode 100644 tests/interpolation.rs diff --git a/src/interpolation/cubic_hermite.rs b/src/interpolation/cubic_hermite.rs new file mode 100644 index 0000000..71fd1bc --- /dev/null +++ b/src/interpolation/cubic_hermite.rs @@ -0,0 +1,263 @@ +//! Cubic Hermite interpolation tables. + +use super::error::InterpolationError; +use super::scalar::{cubic_hermite_segment, hermite_basis, HermiteEvaluation}; +use super::traits::HermiteInterpolable; + +/// A scalar Hermite table sample. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct HermiteSample { + /// Sample abscissa. + pub x: f64, + /// Sample value. + pub y: f64, + /// Sample derivative with respect to `x`. + pub dydx: f64, +} + +/// Piecewise scalar cubic Hermite spline. +/// +/// Derivative continuity is guaranteed only up to the first derivative supplied +/// at each sample. Acceleration or other second derivative quantities are not +/// generally continuous across sample boundaries. +#[derive(Debug, Clone, PartialEq)] +pub struct CubicHermiteSpline { + samples: Vec, +} + +impl CubicHermiteSpline { + /// Builds a spline from samples sorted by strictly increasing `x`. + pub fn new(samples: Vec) -> Result { + validate_len(samples.len())?; + for sample in &samples { + if !sample.x.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + if !sample.y.is_finite() || !sample.dydx.is_finite() { + return Err(InterpolationError::NonFiniteValue); + } + } + validate_sorted(samples.iter().map(|sample| sample.x))?; + Ok(Self { samples }) + } + + /// Returns the table samples. + pub fn samples(&self) -> &[HermiteSample] { + &self.samples + } + + /// Evaluates the spline without extrapolation. + pub fn evaluate(&self, x: f64) -> Result { + if !x.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + let (min, max) = self.range(); + if x < min || x > max { + return Err(InterpolationError::OutOfRange { x, min, max }); + } + + let segment = self.segment_index(x); + let s0 = self.samples[segment]; + let s1 = self.samples[segment + 1]; + if x == s0.x { + return Ok(HermiteEvaluation { + value: s0.y, + derivative: s0.dydx, + }); + } + if x == s1.x { + return Ok(HermiteEvaluation { + value: s1.y, + derivative: s1.dydx, + }); + } + cubic_hermite_segment(x, s0.x, s1.x, s0.y, s0.dydx, s1.y, s1.dydx) + } + + fn range(&self) -> (f64, f64) { + (self.samples[0].x, self.samples[self.samples.len() - 1].x) + } + + fn segment_index(&self, x: f64) -> usize { + match self + .samples + .binary_search_by(|sample| sample.x.total_cmp(&x)) + { + Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), + Err(index) => (index - 1).min(self.samples.len() - 2), + } + } +} + +/// A typed Hermite table node. +#[derive(Debug, Clone, PartialEq)] +pub struct HermiteNode +where + T: HermiteInterpolable, +{ + /// Sample abscissa. + pub x: f64, + /// Sample value. + pub value: T, + /// Sample derivative with respect to `x`. + pub derivative: T::Derivative, +} + +/// A typed Hermite table evaluation. +#[derive(Debug, Clone, PartialEq)] +pub struct HermiteTableEvaluation +where + T: HermiteInterpolable, +{ + /// Interpolated value. + pub value: T, + /// Interpolated derivative with respect to `x`. + pub derivative: T::Derivative, +} + +/// Piecewise cubic Hermite interpolation table for typed values. +/// +/// The abscissa is a monotonic scalar owned by the caller. For time-domain +/// ephemerides, downstream code should map epochs to a scalar such as elapsed +/// seconds before constructing this table. +/// +/// Derivative continuity is guaranteed only up to the first derivative supplied +/// at each sample. Acceleration or other second derivative quantities are not +/// generally continuous across sample boundaries. +pub struct CubicHermiteTable +where + T: HermiteInterpolable, +{ + samples: Vec>, +} + +impl CubicHermiteTable +where + T: HermiteInterpolable, +{ + /// Builds a typed table from nodes sorted by strictly increasing `x`. + pub fn new(samples: Vec>) -> Result { + validate_len(samples.len())?; + for sample in &samples { + if !sample.x.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + if !sample.value.hermite_value_is_finite() + || !T::hermite_derivative_is_finite(&sample.derivative) + { + return Err(InterpolationError::NonFiniteValue); + } + } + validate_sorted(samples.iter().map(|sample| sample.x))?; + Ok(Self { samples }) + } + + /// Returns the table samples. + pub fn samples(&self) -> &[HermiteNode] { + &self.samples + } +} + +impl CubicHermiteTable +where + T: HermiteInterpolable + Clone, + T::Derivative: Clone, +{ + /// Evaluates the table without extrapolation. + pub fn evaluate(&self, x: f64) -> Result, InterpolationError> { + if !x.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + let (min, max) = self.range(); + if x < min || x > max { + return Err(InterpolationError::OutOfRange { x, min, max }); + } + + let segment = self.segment_index(x); + let s0 = &self.samples[segment]; + let s1 = &self.samples[segment + 1]; + if x == s0.x { + return Ok(HermiteTableEvaluation { + value: s0.value.clone(), + derivative: s0.derivative.clone(), + }); + } + if x == s1.x { + return Ok(HermiteTableEvaluation { + value: s1.value.clone(), + derivative: s1.derivative.clone(), + }); + } + + let dt = s1.x - s0.x; + let t = (x - s0.x) / dt; + let basis = hermite_basis(t, dt); + + Ok(HermiteTableEvaluation { + value: T::hermite_combine( + basis.h00, + basis.h10_dt, + basis.h01, + basis.h11_dt, + s0.value.clone(), + s0.derivative.clone(), + s1.value.clone(), + s1.derivative.clone(), + ), + derivative: T::hermite_derivative_combine( + basis.dh00_over_dt, + basis.dh10_over_dt, + basis.dh01_over_dt, + basis.dh11_over_dt, + basis.inv_dt, + s0.value.clone(), + s0.derivative.clone(), + s1.value.clone(), + s1.derivative.clone(), + ), + }) + } + + fn range(&self) -> (f64, f64) { + (self.samples[0].x, self.samples[self.samples.len() - 1].x) + } + + fn segment_index(&self, x: f64) -> usize { + match self + .samples + .binary_search_by(|sample| sample.x.total_cmp(&x)) + { + Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), + Err(index) => (index - 1).min(self.samples.len() - 2), + } + } +} + +fn validate_len(len: usize) -> Result<(), InterpolationError> { + if len == 0 { + return Err(InterpolationError::EmptyTable); + } + if len < 2 { + return Err(InterpolationError::TooFewSamples { + required: 2, + actual: len, + }); + } + Ok(()) +} + +fn validate_sorted(xs: impl IntoIterator) -> Result<(), InterpolationError> { + let mut previous = None; + for x in xs { + if let Some(previous) = previous { + if x == previous { + return Err(InterpolationError::DuplicateAbscissa); + } + if x < previous { + return Err(InterpolationError::UnsortedAbscissa); + } + } + previous = Some(x); + } + Ok(()) +} diff --git a/src/interpolation/error.rs b/src/interpolation/error.rs new file mode 100644 index 0000000..40279de --- /dev/null +++ b/src/interpolation/error.rs @@ -0,0 +1,57 @@ +//! Error types for interpolation routines. + +use std::fmt; + +/// Errors returned by interpolation constructors and evaluators. +#[derive(Debug, Clone, PartialEq)] +pub enum InterpolationError { + /// The interpolation table contains no samples. + EmptyTable, + /// The interpolation table does not contain enough samples. + TooFewSamples { + /// Required sample count. + required: usize, + /// Actual sample count. + actual: usize, + }, + /// An abscissa value is not finite. + NonFiniteAbscissa, + /// A sample value or derivative is not finite. + NonFiniteValue, + /// Two adjacent samples have the same abscissa. + DuplicateAbscissa, + /// Input samples are not sorted by increasing abscissa. + UnsortedAbscissa, + /// The query abscissa is outside the interpolation table range. + OutOfRange { + /// Query abscissa. + x: f64, + /// Minimum supported abscissa. + min: f64, + /// Maximum supported abscissa. + max: f64, + }, +} + +impl fmt::Display for InterpolationError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::EmptyTable => write!(f, "interpolation table is empty"), + Self::TooFewSamples { required, actual } => write!( + f, + "interpolation table has too few samples: required {required}, got {actual}" + ), + Self::NonFiniteAbscissa => write!(f, "interpolation abscissa is not finite"), + Self::NonFiniteValue => write!(f, "interpolation value is not finite"), + Self::DuplicateAbscissa => { + write!(f, "interpolation samples contain duplicate abscissae") + } + Self::UnsortedAbscissa => write!(f, "interpolation samples are not sorted"), + Self::OutOfRange { x, min, max } => { + write!(f, "interpolation query {x} is outside range [{min}, {max}]") + } + } + } +} + +impl std::error::Error for InterpolationError {} diff --git a/src/interpolation/linear.rs b/src/interpolation/linear.rs new file mode 100644 index 0000000..f4a62c6 --- /dev/null +++ b/src/interpolation/linear.rs @@ -0,0 +1,34 @@ +//! Linear interpolation helpers. + +use super::error::InterpolationError; + +/// Evaluates scalar linear interpolation on one segment. +pub fn linear_segment( + x: f64, + x0: f64, + x1: f64, + y0: f64, + y1: f64, +) -> Result { + if !x.is_finite() || !x0.is_finite() || !x1.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + if !y0.is_finite() || !y1.is_finite() { + return Err(InterpolationError::NonFiniteValue); + } + if x1 <= x0 { + if x1 == x0 { + return Err(InterpolationError::DuplicateAbscissa); + } + return Err(InterpolationError::UnsortedAbscissa); + } + if x < x0 || x > x1 { + return Err(InterpolationError::OutOfRange { + x, + min: x0, + max: x1, + }); + } + let t = (x - x0) / (x1 - x0); + Ok((1.0 - t) * y0 + t * y1) +} diff --git a/src/interpolation/mod.rs b/src/interpolation/mod.rs new file mode 100644 index 0000000..5ef5209 --- /dev/null +++ b/src/interpolation/mod.rs @@ -0,0 +1,23 @@ +//! Interpolation primitives for typed affine geometry. +//! +//! This module is domain-agnostic: it works over a scalar abscissa and complete +//! typed values. Astronomy-specific epoch handling belongs in downstream crates. + +pub mod cubic_hermite; +pub mod error; +pub mod linear; +pub mod quantity; +pub mod scalar; +pub mod traits; + +pub use cubic_hermite::{ + CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, HermiteTableEvaluation, +}; +pub use error::InterpolationError; +pub use linear::linear_segment; +pub use quantity::{ + CubicHermiteQuantityTable, QuantityHermiteInterpolable, QuantityHermiteNode, + QuantityHermiteTableEvaluation, +}; +pub use scalar::{cubic_hermite_segment, HermiteEvaluation}; +pub use traits::HermiteInterpolable; diff --git a/src/interpolation/quantity.rs b/src/interpolation/quantity.rs new file mode 100644 index 0000000..0b7b490 --- /dev/null +++ b/src/interpolation/quantity.rs @@ -0,0 +1,409 @@ +//! Quantity-abscissa cubic Hermite interpolation tables. + +use super::error::InterpolationError; +use super::scalar::hermite_basis; +use crate::cartesian::{Position, Vector}; +use crate::centers::ReferenceCenter; +use crate::frames::ReferenceFrame; +use qtty::length::LengthUnit; +use qtty::{Quantity, Unit, UnitDiv, UnitMul}; + +/// A value that can be combined by cubic Hermite interpolation over a typed +/// quantity abscissa. +pub trait QuantityHermiteInterpolable: Sized { + /// Derivative type with respect to the abscissa unit `X`. + type Derivative; + + /// Combines endpoint values and derivatives into an interpolated value. + fn quantity_hermite_combine( + h00: f64, + h10_dx: Quantity, + h01: f64, + h11_dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self; + + /// Combines endpoint values and derivatives into an interpolated derivative. + fn quantity_hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative; + + /// Returns whether this value is finite. + fn quantity_hermite_value_is_finite(&self) -> bool { + true + } + + /// Returns whether this derivative is finite. + fn quantity_hermite_derivative_is_finite(_derivative: &Self::Derivative) -> bool { + true + } +} + +impl QuantityHermiteInterpolable for f64 { + type Derivative = f64; + + #[inline] + fn quantity_hermite_combine( + h00: f64, + h10_dx: Quantity, + h01: f64, + h11_dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + h00 * y0 + h10_dx.value() * dy0 + h01 * y1 + h11_dx.value() * dy1 + } + + #[inline] + fn quantity_hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + let inv_dx = 1.0 / dx.value(); + dh00_over_dt * inv_dx * y0 + + dh10_over_dt * dy0 + + dh01_over_dt * inv_dx * y1 + + dh11_over_dt * dy1 + } + + #[inline] + fn quantity_hermite_value_is_finite(&self) -> bool { + self.is_finite() + } + + #[inline] + fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.is_finite() + } +} + +impl QuantityHermiteInterpolable for Quantity +where + X: Unit, + U: Unit + UnitDiv, + >::Output: Unit + UnitMul, +{ + type Derivative = Quantity<>::Output>; + + #[inline] + fn quantity_hermite_combine( + h00: f64, + h10_dx: Quantity, + h01: f64, + h11_dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + y0 * h00 + (dy0 * h10_dx) + y1 * h01 + (dy1 * h11_dx) + } + + #[inline] + fn quantity_hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + (y0 / dx) * dh00_over_dt + + dy0 * dh10_over_dt + + (y1 / dx) * dh01_over_dt + + dy1 * dh11_over_dt + } + + #[inline] + fn quantity_hermite_value_is_finite(&self) -> bool { + self.value().is_finite() + } + + #[inline] + fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.value().is_finite() + } +} + +impl QuantityHermiteInterpolable for Position +where + X: Unit, + C: ReferenceCenter, + F: ReferenceFrame, + L: LengthUnit + UnitDiv, + >::Output: Unit + UnitMul, +{ + type Derivative = Vector>::Output>; + + #[inline] + fn quantity_hermite_combine( + _h00: f64, + h10_dx: Quantity, + h01: f64, + h11_dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + let chord = y1 - y0; + y0 + chord.scale(h01) + + displacement_from_derivative(dy0, h10_dx) + + displacement_from_derivative(dy1, h11_dx) + } + + #[inline] + fn quantity_hermite_derivative_combine( + _dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + dx: Quantity, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + let chord = y1 - y0; + Vector::>::Output>::new( + (chord.x() / dx) * dh01_over_dt + dy0.x() * dh10_over_dt + dy1.x() * dh11_over_dt, + (chord.y() / dx) * dh01_over_dt + dy0.y() * dh10_over_dt + dy1.y() * dh11_over_dt, + (chord.z() / dx) * dh01_over_dt + dy0.z() * dh10_over_dt + dy1.z() * dh11_over_dt, + ) + } + + #[inline] + fn quantity_hermite_value_is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + + #[inline] + fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.x().value().is_finite() + && derivative.y().value().is_finite() + && derivative.z().value().is_finite() + } +} + +#[inline] +fn displacement_from_derivative( + derivative: Vector>::Output>, + dx: Quantity, +) -> Vector +where + X: Unit, + F: ReferenceFrame, + L: LengthUnit + UnitDiv, + >::Output: Unit + UnitMul, +{ + Vector::::new( + derivative.x() * dx, + derivative.y() * dx, + derivative.z() * dx, + ) +} + +/// A typed Hermite table node with a quantity abscissa. +#[derive(Debug, Clone, PartialEq)] +pub struct QuantityHermiteNode +where + X: Unit, + T: QuantityHermiteInterpolable, +{ + /// Sample abscissa. + pub x: Quantity, + /// Sample value. + pub value: T, + /// Sample derivative with respect to `x`. + pub derivative: T::Derivative, +} + +/// A typed Hermite table evaluation with a quantity abscissa. +#[derive(Debug, Clone, PartialEq)] +pub struct QuantityHermiteTableEvaluation +where + X: Unit, + T: QuantityHermiteInterpolable, +{ + /// Interpolated value. + pub value: T, + /// Interpolated derivative with respect to `x`. + pub derivative: T::Derivative, + /// Evaluated abscissa. + pub x: Quantity, +} + +/// Piecewise cubic Hermite interpolation table over a typed quantity abscissa. +pub struct CubicHermiteQuantityTable +where + X: Unit, + T: QuantityHermiteInterpolable, +{ + samples: Vec>, +} + +impl CubicHermiteQuantityTable +where + X: Unit, + T: QuantityHermiteInterpolable, +{ + /// Builds a typed table from nodes sorted by strictly increasing `x`. + pub fn new(samples: Vec>) -> Result { + validate_len(samples.len())?; + for sample in &samples { + if !sample.x.value().is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + if !sample.value.quantity_hermite_value_is_finite() + || !T::quantity_hermite_derivative_is_finite(&sample.derivative) + { + return Err(InterpolationError::NonFiniteValue); + } + } + validate_sorted(samples.iter().map(|sample| sample.x.value()))?; + Ok(Self { samples }) + } + + /// Returns the table samples. + pub fn samples(&self) -> &[QuantityHermiteNode] { + &self.samples + } +} + +impl CubicHermiteQuantityTable +where + X: Unit, + T: QuantityHermiteInterpolable + Clone, + T::Derivative: Clone, +{ + /// Evaluates the table without extrapolation. + pub fn evaluate( + &self, + x: Quantity, + ) -> Result, InterpolationError> { + if !x.value().is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + let (min, max) = self.range_raw(); + let x_raw = x.value(); + if x_raw < min || x_raw > max { + return Err(InterpolationError::OutOfRange { x: x_raw, min, max }); + } + + let segment = self.segment_index(x_raw); + let s0 = &self.samples[segment]; + let s1 = &self.samples[segment + 1]; + if x_raw == s0.x.value() { + return Ok(QuantityHermiteTableEvaluation { + value: s0.value.clone(), + derivative: s0.derivative.clone(), + x, + }); + } + if x_raw == s1.x.value() { + return Ok(QuantityHermiteTableEvaluation { + value: s1.value.clone(), + derivative: s1.derivative.clone(), + x, + }); + } + + let dx_raw = s1.x.value() - s0.x.value(); + let normalized = (x_raw - s0.x.value()) / dx_raw; + let basis = hermite_basis(normalized, dx_raw); + + Ok(QuantityHermiteTableEvaluation { + value: T::quantity_hermite_combine( + basis.h00, + Quantity::::new(basis.h10_dt), + basis.h01, + Quantity::::new(basis.h11_dt), + s0.value.clone(), + s0.derivative.clone(), + s1.value.clone(), + s1.derivative.clone(), + ), + derivative: T::quantity_hermite_derivative_combine( + basis.dh00_over_dt, + basis.dh10_over_dt, + basis.dh01_over_dt, + basis.dh11_over_dt, + Quantity::::new(dx_raw), + s0.value.clone(), + s0.derivative.clone(), + s1.value.clone(), + s1.derivative.clone(), + ), + x, + }) + } + + fn range_raw(&self) -> (f64, f64) { + ( + self.samples[0].x.value(), + self.samples[self.samples.len() - 1].x.value(), + ) + } + + fn segment_index(&self, x: f64) -> usize { + match self + .samples + .binary_search_by(|sample| sample.x.value().total_cmp(&x)) + { + Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), + Err(index) => (index - 1).min(self.samples.len() - 2), + } + } +} + +fn validate_len(len: usize) -> Result<(), InterpolationError> { + if len == 0 { + return Err(InterpolationError::EmptyTable); + } + if len < 2 { + return Err(InterpolationError::TooFewSamples { + required: 2, + actual: len, + }); + } + Ok(()) +} + +fn validate_sorted(xs: impl IntoIterator) -> Result<(), InterpolationError> { + let mut previous = None; + for x in xs { + if let Some(previous) = previous { + if x == previous { + return Err(InterpolationError::DuplicateAbscissa); + } + if x < previous { + return Err(InterpolationError::UnsortedAbscissa); + } + } + previous = Some(x); + } + Ok(()) +} diff --git a/src/interpolation/scalar.rs b/src/interpolation/scalar.rs new file mode 100644 index 0000000..a76f100 --- /dev/null +++ b/src/interpolation/scalar.rs @@ -0,0 +1,87 @@ +//! Scalar interpolation primitives. + +use super::error::InterpolationError; + +/// Scalar Hermite interpolation result. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct HermiteEvaluation { + /// Interpolated value. + pub value: f64, + /// Derivative with respect to the original abscissa. + pub derivative: f64, +} + +/// Evaluates one scalar cubic Hermite segment. +/// +/// The returned derivative is with respect to `x`, not normalized `t`. +pub fn cubic_hermite_segment( + x: f64, + x0: f64, + x1: f64, + y0: f64, + dy0: f64, + y1: f64, + dy1: f64, +) -> Result { + if !x.is_finite() || !x0.is_finite() || !x1.is_finite() { + return Err(InterpolationError::NonFiniteAbscissa); + } + if !y0.is_finite() || !dy0.is_finite() || !y1.is_finite() || !dy1.is_finite() { + return Err(InterpolationError::NonFiniteValue); + } + if x1 <= x0 { + if x1 == x0 { + return Err(InterpolationError::DuplicateAbscissa); + } + return Err(InterpolationError::UnsortedAbscissa); + } + if x < x0 || x > x1 { + return Err(InterpolationError::OutOfRange { + x, + min: x0, + max: x1, + }); + } + + let dt = x1 - x0; + let t = (x - x0) / dt; + let basis = hermite_basis(t, dt); + + Ok(HermiteEvaluation { + value: basis.h00 * y0 + basis.h10_dt * dy0 + basis.h01 * y1 + basis.h11_dt * dy1, + derivative: basis.dh00_over_dt * basis.inv_dt * y0 + + basis.dh10_over_dt * dy0 + + basis.dh01_over_dt * basis.inv_dt * y1 + + basis.dh11_over_dt * dy1, + }) +} + +#[derive(Debug, Clone, Copy)] +pub(crate) struct HermiteBasis { + pub h00: f64, + pub h10_dt: f64, + pub h01: f64, + pub h11_dt: f64, + pub dh00_over_dt: f64, + pub dh10_over_dt: f64, + pub dh01_over_dt: f64, + pub dh11_over_dt: f64, + pub inv_dt: f64, +} + +pub(crate) fn hermite_basis(t: f64, dt: f64) -> HermiteBasis { + let t2 = t * t; + let t3 = t2 * t; + + HermiteBasis { + h00: 2.0 * t3 - 3.0 * t2 + 1.0, + h10_dt: (t3 - 2.0 * t2 + t) * dt, + h01: -2.0 * t3 + 3.0 * t2, + h11_dt: (t3 - t2) * dt, + dh00_over_dt: 6.0 * t2 - 6.0 * t, + dh10_over_dt: 3.0 * t2 - 4.0 * t + 1.0, + dh01_over_dt: -6.0 * t2 + 6.0 * t, + dh11_over_dt: 3.0 * t2 - 2.0 * t, + inv_dt: 1.0 / dt, + } +} diff --git a/src/interpolation/traits.rs b/src/interpolation/traits.rs new file mode 100644 index 0000000..166a2f6 --- /dev/null +++ b/src/interpolation/traits.rs @@ -0,0 +1,247 @@ +//! Traits for interpolating complete typed values. + +use crate::cartesian::{Position, Vector}; +use crate::centers::ReferenceCenter; +use crate::frames::ReferenceFrame; +use qtty::length::LengthUnit; +use qtty::{Quantity, Unit}; + +/// A value that can be combined by cubic Hermite interpolation. +/// +/// Implementations operate on complete typed values, preserving all frame, +/// center, and unit tags carried by the type. +pub trait HermiteInterpolable: Sized { + /// Derivative type with respect to the scalar interpolation abscissa. + type Derivative; + + /// Combines endpoint values and derivatives into an interpolated value. + fn hermite_combine( + h00: f64, + h10_dt: f64, + h01: f64, + h11_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self; + + /// Combines endpoint values and derivatives into an interpolated derivative. + fn hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + inv_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative; + + /// Returns whether this value is finite. + fn hermite_value_is_finite(&self) -> bool { + true + } + + /// Returns whether this derivative is finite. + fn hermite_derivative_is_finite(_derivative: &Self::Derivative) -> bool { + true + } +} + +impl HermiteInterpolable for f64 { + type Derivative = f64; + + #[inline] + fn hermite_combine( + h00: f64, + h10_dt: f64, + h01: f64, + h11_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + h00 * y0 + h10_dt * dy0 + h01 * y1 + h11_dt * dy1 + } + + #[inline] + fn hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + inv_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + dh00_over_dt * inv_dt * y0 + + dh10_over_dt * dy0 + + dh01_over_dt * inv_dt * y1 + + dh11_over_dt * dy1 + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.is_finite() + } +} + +impl HermiteInterpolable for Quantity { + type Derivative = Quantity; + + #[inline] + fn hermite_combine( + h00: f64, + h10_dt: f64, + h01: f64, + h11_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + y0 * h00 + dy0 * h10_dt + y1 * h01 + dy1 * h11_dt + } + + #[inline] + fn hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + inv_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + y0 * (dh00_over_dt * inv_dt) + + dy0 * dh10_over_dt + + y1 * (dh01_over_dt * inv_dt) + + dy1 * dh11_over_dt + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.value().is_finite() + } +} + +impl HermiteInterpolable for Vector +where + F: ReferenceFrame, + U: Unit, +{ + type Derivative = Vector; + + #[inline] + fn hermite_combine( + h00: f64, + h10_dt: f64, + h01: f64, + h11_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + y0.scale(h00) + dy0.scale(h10_dt) + y1.scale(h01) + dy1.scale(h11_dt) + } + + #[inline] + fn hermite_derivative_combine( + dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + inv_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + y0.scale(dh00_over_dt * inv_dt) + + dy0.scale(dh10_over_dt) + + y1.scale(dh01_over_dt * inv_dt) + + dy1.scale(dh11_over_dt) + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.x().value().is_finite() + && derivative.y().value().is_finite() + && derivative.z().value().is_finite() + } +} + +impl HermiteInterpolable for Position +where + C: ReferenceCenter, + F: ReferenceFrame, + U: LengthUnit, +{ + type Derivative = Vector; + + #[inline] + fn hermite_combine( + _h00: f64, + h10_dt: f64, + h01: f64, + h11_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + let chord = y1 - y0; + y0 + chord.scale(h01) + dy0.scale(h10_dt) + dy1.scale(h11_dt) + } + + #[inline] + fn hermite_derivative_combine( + _dh00_over_dt: f64, + dh10_over_dt: f64, + dh01_over_dt: f64, + dh11_over_dt: f64, + inv_dt: f64, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + let chord = y1 - y0; + chord.scale(dh01_over_dt * inv_dt) + dy0.scale(dh10_over_dt) + dy1.scale(dh11_over_dt) + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.x().value().is_finite() + && derivative.y().value().is_finite() + && derivative.z().value().is_finite() + } +} diff --git a/src/lib.rs b/src/lib.rs index f167a40..6f2b99d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -120,6 +120,7 @@ mod op_macros; // Coordinate type implementations pub mod cartesian; pub mod conic; +pub mod interpolation; pub mod spherical; // Core traits and marker types @@ -206,6 +207,12 @@ pub mod prelude { NonParabolicKindMarker, OrientedConic, Parabolic, ParabolicPeriapsis, PeriapsisParam, SemiMajorAxisParam, TypedPeriapsisParam, TypedSemiMajorAxisParam, }; + pub use crate::interpolation::{ + cubic_hermite_segment, CubicHermiteQuantityTable, CubicHermiteSpline, CubicHermiteTable, + HermiteEvaluation, HermiteInterpolable, HermiteNode, HermiteSample, HermiteTableEvaluation, + InterpolationError, QuantityHermiteInterpolable, QuantityHermiteNode, + QuantityHermiteTableEvaluation, + }; pub use crate::spherical::{Direction as SphericalDirection, Position as SphericalPosition}; // Operators diff --git a/tests/interpolation.rs b/tests/interpolation.rs new file mode 100644 index 0000000..f8c7a8b --- /dev/null +++ b/tests/interpolation.rs @@ -0,0 +1,271 @@ +use affn::cartesian::{Position, Velocity}; +use affn::interpolation::{ + CubicHermiteQuantityTable, CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, + InterpolationError, QuantityHermiteNode, +}; +use affn::{DeriveReferenceCenter as ReferenceCenter, DeriveReferenceFrame as ReferenceFrame}; +use qtty::unit::{Kilometer, Meter, Second}; +use qtty::{Per, Quantity}; + +#[derive(Debug, Copy, Clone, ReferenceFrame)] +struct TestFrame; + +#[derive(Debug, Copy, Clone, ReferenceCenter)] +struct TestCenter; + +type TestPosition = Position; +type TestVelocity = Velocity; +type TestKmPerSecond = Per; +type TestKilometerPosition = Position; +type TestKilometerVelocity = Velocity; + +fn cubic(x: f64) -> f64 { + x * x * x - 2.0 * x * x + x + 1.0 +} + +fn cubic_derivative(x: f64) -> f64 { + 3.0 * x * x - 4.0 * x + 1.0 +} + +fn cubic_position(t: f64) -> TestKilometerPosition { + TestKilometerPosition::new(cubic(t), -2.0 * cubic(t), 0.5 * cubic(t)) +} + +fn cubic_velocity(t: f64) -> TestKilometerVelocity { + TestKilometerVelocity::new( + Quantity::::new(cubic_derivative(t)), + Quantity::::new(-2.0 * cubic_derivative(t)), + Quantity::::new(0.5 * cubic_derivative(t)), + ) +} + +#[test] +fn scalar_cubic_polynomial_is_reproduced() { + let spline = CubicHermiteSpline::new(vec![ + HermiteSample { + x: -1.0, + y: cubic(-1.0), + dydx: cubic_derivative(-1.0), + }, + HermiteSample { + x: 0.5, + y: cubic(0.5), + dydx: cubic_derivative(0.5), + }, + HermiteSample { + x: 2.0, + y: cubic(2.0), + dydx: cubic_derivative(2.0), + }, + ]) + .unwrap(); + + for x in [-0.75, -0.1, 0.25, 1.25, 1.75] { + let evaluated = spline.evaluate(x).unwrap(); + assert!((evaluated.value - cubic(x)).abs() < 1e-12); + assert!((evaluated.derivative - cubic_derivative(x)).abs() < 1e-12); + } +} + +#[test] +fn exact_node_evaluation_returns_node_value_and_derivative() { + let spline = CubicHermiteSpline::new(vec![ + HermiteSample { + x: 0.0, + y: 10.0, + dydx: -3.0, + }, + HermiteSample { + x: 2.0, + y: 20.0, + dydx: 4.0, + }, + ]) + .unwrap(); + + let evaluated = spline.evaluate(2.0).unwrap(); + assert_eq!(evaluated.value, 20.0); + assert_eq!(evaluated.derivative, 4.0); +} + +#[test] +fn linear_motion_with_constant_velocity_is_exact() { + let table = CubicHermiteTable::new(vec![ + HermiteNode { + x: 0.0, + value: TestPosition::new(1.0, 2.0, 3.0), + derivative: TestVelocity::new(0.5, -1.0, 2.0), + }, + HermiteNode { + x: 4.0, + value: TestPosition::new(3.0, -2.0, 11.0), + derivative: TestVelocity::new(0.5, -1.0, 2.0), + }, + ]) + .unwrap(); + + let evaluated = table.evaluate(1.5).unwrap(); + assert!((evaluated.value.x().value() - 1.75).abs() < 1e-12); + assert!((evaluated.value.y().value() - 0.5).abs() < 1e-12); + assert!((evaluated.value.z().value() - 6.0).abs() < 1e-12); + assert!((evaluated.derivative.x().value() - 0.5).abs() < 1e-12); + assert!((evaluated.derivative.y().value() + 1.0).abs() < 1e-12); + assert!((evaluated.derivative.z().value() - 2.0).abs() < 1e-12); +} + +#[test] +fn quantity_table_accepts_position_over_seconds_with_velocity() { + let table = CubicHermiteQuantityTable::::new(vec![ + QuantityHermiteNode { + x: qtty::Second::new(0.0), + value: TestKilometerPosition::new(1.0, 2.0, 3.0), + derivative: TestKilometerVelocity::new( + Quantity::::new(0.5), + Quantity::::new(-1.0), + Quantity::::new(2.0), + ), + }, + QuantityHermiteNode { + x: qtty::Second::new(4.0), + value: TestKilometerPosition::new(3.0, -2.0, 11.0), + derivative: TestKilometerVelocity::new( + Quantity::::new(0.5), + Quantity::::new(-1.0), + Quantity::::new(2.0), + ), + }, + ]) + .unwrap(); + + let evaluated = table.evaluate(qtty::Second::new(1.5)).unwrap(); + assert!((evaluated.value.x().value() - 1.75).abs() < 1e-12); + assert!((evaluated.value.y().value() - 0.5).abs() < 1e-12); + assert!((evaluated.value.z().value() - 6.0).abs() < 1e-12); + assert!((evaluated.derivative.x().value() - 0.5).abs() < 1e-12); + assert!((evaluated.derivative.y().value() + 1.0).abs() < 1e-12); + assert!((evaluated.derivative.z().value() - 2.0).abs() < 1e-12); +} + +#[test] +fn quantity_table_reproduces_cubic_position_over_seconds() { + let table = CubicHermiteQuantityTable::::new(vec![ + QuantityHermiteNode { + x: qtty::Second::new(-1.0), + value: cubic_position(-1.0), + derivative: cubic_velocity(-1.0), + }, + QuantityHermiteNode { + x: qtty::Second::new(0.5), + value: cubic_position(0.5), + derivative: cubic_velocity(0.5), + }, + QuantityHermiteNode { + x: qtty::Second::new(2.0), + value: cubic_position(2.0), + derivative: cubic_velocity(2.0), + }, + ]) + .unwrap(); + + for t in [-0.75, -0.1, 0.25, 1.25, 1.75] { + let evaluated = table.evaluate(qtty::Second::new(t)).unwrap(); + let expected_position = cubic_position(t); + let expected_velocity = cubic_velocity(t); + assert!((evaluated.value.x().value() - expected_position.x().value()).abs() < 1e-12); + assert!((evaluated.value.y().value() - expected_position.y().value()).abs() < 1e-12); + assert!((evaluated.value.z().value() - expected_position.z().value()).abs() < 1e-12); + assert!((evaluated.derivative.x().value() - expected_velocity.x().value()).abs() < 1e-12); + assert!((evaluated.derivative.y().value() - expected_velocity.y().value()).abs() < 1e-12); + assert!((evaluated.derivative.z().value() - expected_velocity.z().value()).abs() < 1e-12); + } +} + +#[test] +fn non_uniform_sample_spacing_works() { + let spline = CubicHermiteSpline::new(vec![ + HermiteSample { + x: 0.0, + y: cubic(0.0), + dydx: cubic_derivative(0.0), + }, + HermiteSample { + x: 0.25, + y: cubic(0.25), + dydx: cubic_derivative(0.25), + }, + HermiteSample { + x: 2.5, + y: cubic(2.5), + dydx: cubic_derivative(2.5), + }, + ]) + .unwrap(); + + let evaluated = spline.evaluate(1.75).unwrap(); + assert!((evaluated.value - cubic(1.75)).abs() < 1e-12); + assert!((evaluated.derivative - cubic_derivative(1.75)).abs() < 1e-12); +} + +#[test] +fn out_of_range_queries_return_error() { + let spline = CubicHermiteSpline::new(vec![ + HermiteSample { + x: 0.0, + y: 0.0, + dydx: 1.0, + }, + HermiteSample { + x: 1.0, + y: 1.0, + dydx: 1.0, + }, + ]) + .unwrap(); + + assert_eq!( + spline.evaluate(2.0), + Err(InterpolationError::OutOfRange { + x: 2.0, + min: 0.0, + max: 1.0 + }) + ); +} + +#[test] +fn duplicate_abscissae_are_rejected() { + assert_eq!( + CubicHermiteSpline::new(vec![ + HermiteSample { + x: 0.0, + y: 0.0, + dydx: 0.0, + }, + HermiteSample { + x: 0.0, + y: 1.0, + dydx: 1.0, + }, + ]), + Err(InterpolationError::DuplicateAbscissa) + ); +} + +#[test] +fn unsorted_abscissae_are_rejected() { + assert_eq!( + CubicHermiteSpline::new(vec![ + HermiteSample { + x: 1.0, + y: 1.0, + dydx: 1.0, + }, + HermiteSample { + x: 0.0, + y: 0.0, + dydx: 0.0, + }, + ]), + Err(InterpolationError::UnsortedAbscissa) + ); +} From 996bdba550026bef990d7f1d13ab7c782ed341ed Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 20:56:57 +0200 Subject: [PATCH 2/6] Refactor interpolation module: remove quantity-based interpolation, enhance traits - Removed `quantity.rs` and its associated types and traits for quantity-based cubic Hermite interpolation. - Updated `interpolation/mod.rs` to reflect the removal of quantity-related modules and restructured public exports. - Enhanced `traits.rs` to introduce `HermiteBasis` for better handling of cubic Hermite interpolation across different abscissas. - Modified `HermiteInterpolable` trait to work with the new `HermiteBasis` structure. - Updated implementations for `f64`, `Quantity`, `Vector`, and `Position` to align with the new trait definitions. - Adjusted tests to remove references to quantity-based interpolation and ensure compatibility with the updated cubic Hermite table structure. --- src/cartesian/vector.rs | 42 ++- src/interpolation/abscissa.rs | 47 ++++ src/interpolation/cubic_hermite.rs | 110 ++++---- src/interpolation/mod.rs | 10 +- src/interpolation/quantity.rs | 409 ----------------------------- src/interpolation/traits.rs | 358 +++++++++++++++++++------ src/lib.rs | 7 +- tests/interpolation.rs | 112 +++++++- 8 files changed, 530 insertions(+), 565 deletions(-) create mode 100644 src/interpolation/abscissa.rs delete mode 100644 src/interpolation/quantity.rs diff --git a/src/cartesian/vector.rs b/src/cartesian/vector.rs index bd7915c..e6b48c2 100644 --- a/src/cartesian/vector.rs +++ b/src/cartesian/vector.rs @@ -51,10 +51,10 @@ use super::xyz::XYZ; use crate::frames::ReferenceFrame; use qtty::dimensionless::Ratio; use qtty::length::LengthUnit; -use qtty::{Quantity, Unit, UnitMul}; +use qtty::{Quantity, Unit, UnitDiv, UnitMul}; use std::marker::PhantomData; -use std::ops::{Add, Div, Neg, Sub}; +use std::ops::{Add, Div, Mul, Neg, Sub}; /// A free vector in 3D Cartesian coordinates. /// @@ -199,6 +199,23 @@ impl Vector { pub fn negate(&self) -> Self { Self::from_xyz(-self.xyz) } + + /// Divides this vector by a typed quantity, carrying the resulting unit in + /// the output vector. + /// + /// This is the dimensional counterpart to [`scale`](Self::scale). It is an + /// inherent method rather than a `/` operator so it can coexist with the + /// established same-unit division operator that returns a dimensionless + /// ratio vector. + #[inline] + pub fn div_quantity(&self, rhs: Quantity) -> Vector>::Output> + where + X: Unit, + U: UnitDiv, + >::Output: Unit, + { + Vector::>::Output>::new(self.x() / rhs, self.y() / rhs, self.z() / rhs) + } } // ============================================================================= @@ -321,14 +338,31 @@ impl Neg for Vector { forward_ref_unop! { impl[F: ReferenceFrame, U: Unit] Neg, neg for Vector } // ============================================================================= -// Scalar Division: Vector / Quantity → Vector +// Dimensional Operations: Vector * Quantity and / Quantity // ============================================================================= +impl Mul> for Vector +where + F: ReferenceFrame, + U: Unit + UnitMul, + X: Unit, +{ + type Output = Vector>::Output>; + + /// Multiplies every component by a typed quantity. + #[inline] + fn mul(self, rhs: Quantity) -> Self::Output { + Vector::>::Output>::new(self.x() * rhs, self.y() * rhs, self.z() * rhs) + } +} + +forward_ref_binop! { impl[F: ReferenceFrame, U: Unit + UnitMul, X: Unit] Mul, mul for Vector, Quantity } + impl Div> for Vector { type Output = Vector; /// Divides each component by a scalar with the same unit, producing a - /// dimensionless velocity vector (e.g. β = v/c for aberration). + /// dimensionless vector (e.g. beta = v/c for aberration). #[inline] fn div(self, rhs: Quantity) -> Vector { let c = rhs.value(); diff --git a/src/interpolation/abscissa.rs b/src/interpolation/abscissa.rs new file mode 100644 index 0000000..cc9afb6 --- /dev/null +++ b/src/interpolation/abscissa.rs @@ -0,0 +1,47 @@ +//! Interpolation abscissa types. + +use qtty::{Quantity, Unit}; + +/// An ordered interpolation abscissa. +/// +/// Implemented for raw scalar parameters and typed `qtty` quantities. New +/// interpolation APIs should use this trait instead of exposing separate raw +/// and quantity table types. +pub trait InterpolationAbscissa: Copy { + /// Difference between two abscissae. + type Delta: Copy; + + /// Returns the raw scalar used for ordering and normalized interpolation. + fn raw(self) -> f64; + + /// Builds a delta value from a raw difference in this abscissa' stored unit. + fn delta_from_raw(raw: f64) -> Self::Delta; +} + +impl InterpolationAbscissa for f64 { + type Delta = f64; + + #[inline] + fn raw(self) -> f64 { + self + } + + #[inline] + fn delta_from_raw(raw: f64) -> Self::Delta { + raw + } +} + +impl InterpolationAbscissa for Quantity { + type Delta = Quantity; + + #[inline] + fn raw(self) -> f64 { + self.value() + } + + #[inline] + fn delta_from_raw(raw: f64) -> Self::Delta { + Quantity::::new(raw) + } +} diff --git a/src/interpolation/cubic_hermite.rs b/src/interpolation/cubic_hermite.rs index 71fd1bc..b2dc17e 100644 --- a/src/interpolation/cubic_hermite.rs +++ b/src/interpolation/cubic_hermite.rs @@ -1,8 +1,9 @@ //! Cubic Hermite interpolation tables. use super::error::InterpolationError; -use super::scalar::{cubic_hermite_segment, hermite_basis, HermiteEvaluation}; -use super::traits::HermiteInterpolable; +use super::scalar::{cubic_hermite_segment, HermiteEvaluation}; +use super::traits::{HermiteBasis, HermiteInterpolable}; +use super::InterpolationAbscissa; /// A scalar Hermite table sample. #[derive(Debug, Clone, Copy, PartialEq)] @@ -89,57 +90,59 @@ impl CubicHermiteSpline { } } -/// A typed Hermite table node. +/// A Hermite table node. #[derive(Debug, Clone, PartialEq)] -pub struct HermiteNode +pub struct HermiteNode where - T: HermiteInterpolable, + X: InterpolationAbscissa, + T: HermiteInterpolable, { /// Sample abscissa. - pub x: f64, + pub x: X, /// Sample value. pub value: T, /// Sample derivative with respect to `x`. pub derivative: T::Derivative, } -/// A typed Hermite table evaluation. +/// A Hermite table evaluation. #[derive(Debug, Clone, PartialEq)] -pub struct HermiteTableEvaluation +pub struct HermiteTableEvaluation where - T: HermiteInterpolable, + X: InterpolationAbscissa, + T: HermiteInterpolable, { /// Interpolated value. pub value: T, /// Interpolated derivative with respect to `x`. pub derivative: T::Derivative, + /// Evaluated abscissa. + pub x: X, } /// Piecewise cubic Hermite interpolation table for typed values. /// -/// The abscissa is a monotonic scalar owned by the caller. For time-domain -/// ephemerides, downstream code should map epochs to a scalar such as elapsed -/// seconds before constructing this table. -/// -/// Derivative continuity is guaranteed only up to the first derivative supplied -/// at each sample. Acceleration or other second derivative quantities are not -/// generally continuous across sample boundaries. -pub struct CubicHermiteTable +/// `X` may be a raw scalar parameter (`f64`) or a typed `qtty::Quantity` such +/// as seconds or days. Use typed quantities for physical domains so derivatives +/// carry the expected units. +pub struct CubicHermiteTable where - T: HermiteInterpolable, + X: InterpolationAbscissa, + T: HermiteInterpolable, { - samples: Vec>, + samples: Vec>, } -impl CubicHermiteTable +impl CubicHermiteTable where - T: HermiteInterpolable, + X: InterpolationAbscissa, + T: HermiteInterpolable, { /// Builds a typed table from nodes sorted by strictly increasing `x`. - pub fn new(samples: Vec>) -> Result { + pub fn new(samples: Vec>) -> Result { validate_len(samples.len())?; for sample in &samples { - if !sample.x.is_finite() { + if !sample.x.raw().is_finite() { return Err(InterpolationError::NonFiniteAbscissa); } if !sample.value.hermite_value_is_finite() @@ -148,84 +151,85 @@ where return Err(InterpolationError::NonFiniteValue); } } - validate_sorted(samples.iter().map(|sample| sample.x))?; + validate_sorted(samples.iter().map(|sample| sample.x.raw()))?; Ok(Self { samples }) } /// Returns the table samples. - pub fn samples(&self) -> &[HermiteNode] { + pub fn samples(&self) -> &[HermiteNode] { &self.samples } } -impl CubicHermiteTable +impl CubicHermiteTable where - T: HermiteInterpolable + Clone, + X: InterpolationAbscissa, + T: HermiteInterpolable + Clone, T::Derivative: Clone, { /// Evaluates the table without extrapolation. - pub fn evaluate(&self, x: f64) -> Result, InterpolationError> { - if !x.is_finite() { + pub fn evaluate(&self, x: X) -> Result, InterpolationError> { + let x_raw = x.raw(); + if !x_raw.is_finite() { return Err(InterpolationError::NonFiniteAbscissa); } let (min, max) = self.range(); - if x < min || x > max { - return Err(InterpolationError::OutOfRange { x, min, max }); + if x_raw < min || x_raw > max { + return Err(InterpolationError::OutOfRange { x: x_raw, min, max }); } - let segment = self.segment_index(x); + let segment = self.segment_index(x_raw); let s0 = &self.samples[segment]; let s1 = &self.samples[segment + 1]; - if x == s0.x { + if x_raw == s0.x.raw() { return Ok(HermiteTableEvaluation { value: s0.value.clone(), derivative: s0.derivative.clone(), + x, }); } - if x == s1.x { + if x_raw == s1.x.raw() { return Ok(HermiteTableEvaluation { value: s1.value.clone(), derivative: s1.derivative.clone(), + x, }); } - let dt = s1.x - s0.x; - let t = (x - s0.x) / dt; - let basis = hermite_basis(t, dt); + let dx = s1.x.raw() - s0.x.raw(); + let t = (x_raw - s0.x.raw()) / dx; + let basis = HermiteBasis::::new(t, dx); Ok(HermiteTableEvaluation { - value: T::hermite_combine( - basis.h00, - basis.h10_dt, - basis.h01, - basis.h11_dt, + value: T::hermite_value( + basis, s0.value.clone(), s0.derivative.clone(), s1.value.clone(), s1.derivative.clone(), ), - derivative: T::hermite_derivative_combine( - basis.dh00_over_dt, - basis.dh10_over_dt, - basis.dh01_over_dt, - basis.dh11_over_dt, - basis.inv_dt, + derivative: T::hermite_derivative( + basis, s0.value.clone(), s0.derivative.clone(), s1.value.clone(), s1.derivative.clone(), ), + x, }) } fn range(&self) -> (f64, f64) { - (self.samples[0].x, self.samples[self.samples.len() - 1].x) + ( + self.samples[0].x.raw(), + self.samples[self.samples.len() - 1].x.raw(), + ) } fn segment_index(&self, x: f64) -> usize { match self .samples - .binary_search_by(|sample| sample.x.total_cmp(&x)) + .binary_search_by(|sample| sample.x.raw().total_cmp(&x)) { Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), Err(index) => (index - 1).min(self.samples.len() - 2), @@ -233,6 +237,12 @@ where } } +/// Explicit scalar table node alias. +pub type ScalarHermiteNode = HermiteNode; + +/// Explicit scalar table evaluation alias. +pub type ScalarHermiteTableEvaluation = HermiteTableEvaluation; + fn validate_len(len: usize) -> Result<(), InterpolationError> { if len == 0 { return Err(InterpolationError::EmptyTable); diff --git a/src/interpolation/mod.rs b/src/interpolation/mod.rs index 5ef5209..c93bc4f 100644 --- a/src/interpolation/mod.rs +++ b/src/interpolation/mod.rs @@ -3,21 +3,19 @@ //! This module is domain-agnostic: it works over a scalar abscissa and complete //! typed values. Astronomy-specific epoch handling belongs in downstream crates. +pub mod abscissa; pub mod cubic_hermite; pub mod error; pub mod linear; -pub mod quantity; pub mod scalar; pub mod traits; +pub use abscissa::InterpolationAbscissa; pub use cubic_hermite::{ CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, HermiteTableEvaluation, + ScalarHermiteNode, ScalarHermiteTableEvaluation, }; pub use error::InterpolationError; pub use linear::linear_segment; -pub use quantity::{ - CubicHermiteQuantityTable, QuantityHermiteInterpolable, QuantityHermiteNode, - QuantityHermiteTableEvaluation, -}; pub use scalar::{cubic_hermite_segment, HermiteEvaluation}; -pub use traits::HermiteInterpolable; +pub use traits::{HermiteBasis, HermiteInterpolable}; diff --git a/src/interpolation/quantity.rs b/src/interpolation/quantity.rs deleted file mode 100644 index 0b7b490..0000000 --- a/src/interpolation/quantity.rs +++ /dev/null @@ -1,409 +0,0 @@ -//! Quantity-abscissa cubic Hermite interpolation tables. - -use super::error::InterpolationError; -use super::scalar::hermite_basis; -use crate::cartesian::{Position, Vector}; -use crate::centers::ReferenceCenter; -use crate::frames::ReferenceFrame; -use qtty::length::LengthUnit; -use qtty::{Quantity, Unit, UnitDiv, UnitMul}; - -/// A value that can be combined by cubic Hermite interpolation over a typed -/// quantity abscissa. -pub trait QuantityHermiteInterpolable: Sized { - /// Derivative type with respect to the abscissa unit `X`. - type Derivative; - - /// Combines endpoint values and derivatives into an interpolated value. - fn quantity_hermite_combine( - h00: f64, - h10_dx: Quantity, - h01: f64, - h11_dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self; - - /// Combines endpoint values and derivatives into an interpolated derivative. - fn quantity_hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self::Derivative; - - /// Returns whether this value is finite. - fn quantity_hermite_value_is_finite(&self) -> bool { - true - } - - /// Returns whether this derivative is finite. - fn quantity_hermite_derivative_is_finite(_derivative: &Self::Derivative) -> bool { - true - } -} - -impl QuantityHermiteInterpolable for f64 { - type Derivative = f64; - - #[inline] - fn quantity_hermite_combine( - h00: f64, - h10_dx: Quantity, - h01: f64, - h11_dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self { - h00 * y0 + h10_dx.value() * dy0 + h01 * y1 + h11_dx.value() * dy1 - } - - #[inline] - fn quantity_hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self::Derivative { - let inv_dx = 1.0 / dx.value(); - dh00_over_dt * inv_dx * y0 - + dh10_over_dt * dy0 - + dh01_over_dt * inv_dx * y1 - + dh11_over_dt * dy1 - } - - #[inline] - fn quantity_hermite_value_is_finite(&self) -> bool { - self.is_finite() - } - - #[inline] - fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.is_finite() - } -} - -impl QuantityHermiteInterpolable for Quantity -where - X: Unit, - U: Unit + UnitDiv, - >::Output: Unit + UnitMul, -{ - type Derivative = Quantity<>::Output>; - - #[inline] - fn quantity_hermite_combine( - h00: f64, - h10_dx: Quantity, - h01: f64, - h11_dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self { - y0 * h00 + (dy0 * h10_dx) + y1 * h01 + (dy1 * h11_dx) - } - - #[inline] - fn quantity_hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self::Derivative { - (y0 / dx) * dh00_over_dt - + dy0 * dh10_over_dt - + (y1 / dx) * dh01_over_dt - + dy1 * dh11_over_dt - } - - #[inline] - fn quantity_hermite_value_is_finite(&self) -> bool { - self.value().is_finite() - } - - #[inline] - fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.value().is_finite() - } -} - -impl QuantityHermiteInterpolable for Position -where - X: Unit, - C: ReferenceCenter, - F: ReferenceFrame, - L: LengthUnit + UnitDiv, - >::Output: Unit + UnitMul, -{ - type Derivative = Vector>::Output>; - - #[inline] - fn quantity_hermite_combine( - _h00: f64, - h10_dx: Quantity, - h01: f64, - h11_dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self { - let chord = y1 - y0; - y0 + chord.scale(h01) - + displacement_from_derivative(dy0, h10_dx) - + displacement_from_derivative(dy1, h11_dx) - } - - #[inline] - fn quantity_hermite_derivative_combine( - _dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - dx: Quantity, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self::Derivative { - let chord = y1 - y0; - Vector::>::Output>::new( - (chord.x() / dx) * dh01_over_dt + dy0.x() * dh10_over_dt + dy1.x() * dh11_over_dt, - (chord.y() / dx) * dh01_over_dt + dy0.y() * dh10_over_dt + dy1.y() * dh11_over_dt, - (chord.z() / dx) * dh01_over_dt + dy0.z() * dh10_over_dt + dy1.z() * dh11_over_dt, - ) - } - - #[inline] - fn quantity_hermite_value_is_finite(&self) -> bool { - self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() - } - - #[inline] - fn quantity_hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.x().value().is_finite() - && derivative.y().value().is_finite() - && derivative.z().value().is_finite() - } -} - -#[inline] -fn displacement_from_derivative( - derivative: Vector>::Output>, - dx: Quantity, -) -> Vector -where - X: Unit, - F: ReferenceFrame, - L: LengthUnit + UnitDiv, - >::Output: Unit + UnitMul, -{ - Vector::::new( - derivative.x() * dx, - derivative.y() * dx, - derivative.z() * dx, - ) -} - -/// A typed Hermite table node with a quantity abscissa. -#[derive(Debug, Clone, PartialEq)] -pub struct QuantityHermiteNode -where - X: Unit, - T: QuantityHermiteInterpolable, -{ - /// Sample abscissa. - pub x: Quantity, - /// Sample value. - pub value: T, - /// Sample derivative with respect to `x`. - pub derivative: T::Derivative, -} - -/// A typed Hermite table evaluation with a quantity abscissa. -#[derive(Debug, Clone, PartialEq)] -pub struct QuantityHermiteTableEvaluation -where - X: Unit, - T: QuantityHermiteInterpolable, -{ - /// Interpolated value. - pub value: T, - /// Interpolated derivative with respect to `x`. - pub derivative: T::Derivative, - /// Evaluated abscissa. - pub x: Quantity, -} - -/// Piecewise cubic Hermite interpolation table over a typed quantity abscissa. -pub struct CubicHermiteQuantityTable -where - X: Unit, - T: QuantityHermiteInterpolable, -{ - samples: Vec>, -} - -impl CubicHermiteQuantityTable -where - X: Unit, - T: QuantityHermiteInterpolable, -{ - /// Builds a typed table from nodes sorted by strictly increasing `x`. - pub fn new(samples: Vec>) -> Result { - validate_len(samples.len())?; - for sample in &samples { - if !sample.x.value().is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - if !sample.value.quantity_hermite_value_is_finite() - || !T::quantity_hermite_derivative_is_finite(&sample.derivative) - { - return Err(InterpolationError::NonFiniteValue); - } - } - validate_sorted(samples.iter().map(|sample| sample.x.value()))?; - Ok(Self { samples }) - } - - /// Returns the table samples. - pub fn samples(&self) -> &[QuantityHermiteNode] { - &self.samples - } -} - -impl CubicHermiteQuantityTable -where - X: Unit, - T: QuantityHermiteInterpolable + Clone, - T::Derivative: Clone, -{ - /// Evaluates the table without extrapolation. - pub fn evaluate( - &self, - x: Quantity, - ) -> Result, InterpolationError> { - if !x.value().is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - let (min, max) = self.range_raw(); - let x_raw = x.value(); - if x_raw < min || x_raw > max { - return Err(InterpolationError::OutOfRange { x: x_raw, min, max }); - } - - let segment = self.segment_index(x_raw); - let s0 = &self.samples[segment]; - let s1 = &self.samples[segment + 1]; - if x_raw == s0.x.value() { - return Ok(QuantityHermiteTableEvaluation { - value: s0.value.clone(), - derivative: s0.derivative.clone(), - x, - }); - } - if x_raw == s1.x.value() { - return Ok(QuantityHermiteTableEvaluation { - value: s1.value.clone(), - derivative: s1.derivative.clone(), - x, - }); - } - - let dx_raw = s1.x.value() - s0.x.value(); - let normalized = (x_raw - s0.x.value()) / dx_raw; - let basis = hermite_basis(normalized, dx_raw); - - Ok(QuantityHermiteTableEvaluation { - value: T::quantity_hermite_combine( - basis.h00, - Quantity::::new(basis.h10_dt), - basis.h01, - Quantity::::new(basis.h11_dt), - s0.value.clone(), - s0.derivative.clone(), - s1.value.clone(), - s1.derivative.clone(), - ), - derivative: T::quantity_hermite_derivative_combine( - basis.dh00_over_dt, - basis.dh10_over_dt, - basis.dh01_over_dt, - basis.dh11_over_dt, - Quantity::::new(dx_raw), - s0.value.clone(), - s0.derivative.clone(), - s1.value.clone(), - s1.derivative.clone(), - ), - x, - }) - } - - fn range_raw(&self) -> (f64, f64) { - ( - self.samples[0].x.value(), - self.samples[self.samples.len() - 1].x.value(), - ) - } - - fn segment_index(&self, x: f64) -> usize { - match self - .samples - .binary_search_by(|sample| sample.x.value().total_cmp(&x)) - { - Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), - Err(index) => (index - 1).min(self.samples.len() - 2), - } - } -} - -fn validate_len(len: usize) -> Result<(), InterpolationError> { - if len == 0 { - return Err(InterpolationError::EmptyTable); - } - if len < 2 { - return Err(InterpolationError::TooFewSamples { - required: 2, - actual: len, - }); - } - Ok(()) -} - -fn validate_sorted(xs: impl IntoIterator) -> Result<(), InterpolationError> { - let mut previous = None; - for x in xs { - if let Some(previous) = previous { - if x == previous { - return Err(InterpolationError::DuplicateAbscissa); - } - if x < previous { - return Err(InterpolationError::UnsortedAbscissa); - } - } - previous = Some(x); - } - Ok(()) -} diff --git a/src/interpolation/traits.rs b/src/interpolation/traits.rs index 166a2f6..e519edb 100644 --- a/src/interpolation/traits.rs +++ b/src/interpolation/traits.rs @@ -1,25 +1,71 @@ //! Traits for interpolating complete typed values. +use super::abscissa::InterpolationAbscissa; use crate::cartesian::{Position, Vector}; use crate::centers::ReferenceCenter; use crate::frames::ReferenceFrame; use qtty::length::LengthUnit; -use qtty::{Quantity, Unit}; +use qtty::{Quantity, Unit, UnitDiv, UnitMul}; -/// A value that can be combined by cubic Hermite interpolation. +/// Basis coefficients for one cubic Hermite segment. +#[derive(Debug, Clone, Copy)] +pub struct HermiteBasis { + /// Value basis for the first endpoint. + pub h00: f64, + /// Derivative basis for the first endpoint, scaled by the abscissa delta. + pub h10_dx: X::Delta, + /// Value basis for the second endpoint. + pub h01: f64, + /// Derivative basis for the second endpoint, scaled by the abscissa delta. + pub h11_dx: X::Delta, + /// First derivative of `h00` with respect to normalized segment parameter. + pub dh00_over_dt: f64, + /// First derivative of `h10` with respect to normalized segment parameter. + pub dh10_over_dt: f64, + /// First derivative of `h01` with respect to normalized segment parameter. + pub dh01_over_dt: f64, + /// First derivative of `h11` with respect to normalized segment parameter. + pub dh11_over_dt: f64, + /// Full abscissa delta for the segment. + pub dx: X::Delta, + /// Raw inverse abscissa delta, for scalar-abscissa implementations. + pub inv_raw_dx: f64, +} + +impl HermiteBasis { + /// Constructs a basis from normalized `t` and raw segment width. + #[inline] + pub(crate) fn new(t: f64, raw_dx: f64) -> Self { + let t2 = t * t; + let t3 = t2 * t; + + Self { + h00: 2.0 * t3 - 3.0 * t2 + 1.0, + h10_dx: X::delta_from_raw((t3 - 2.0 * t2 + t) * raw_dx), + h01: -2.0 * t3 + 3.0 * t2, + h11_dx: X::delta_from_raw((t3 - t2) * raw_dx), + dh00_over_dt: 6.0 * t2 - 6.0 * t, + dh10_over_dt: 3.0 * t2 - 4.0 * t + 1.0, + dh01_over_dt: -6.0 * t2 + 6.0 * t, + dh11_over_dt: 3.0 * t2 - 2.0 * t, + dx: X::delta_from_raw(raw_dx), + inv_raw_dx: 1.0 / raw_dx, + } + } +} + +/// A value that can be combined by cubic Hermite interpolation over abscissa +/// type `X`. /// /// Implementations operate on complete typed values, preserving all frame, /// center, and unit tags carried by the type. -pub trait HermiteInterpolable: Sized { - /// Derivative type with respect to the scalar interpolation abscissa. +pub trait HermiteInterpolable: Sized { + /// Derivative type with respect to the interpolation abscissa. type Derivative; /// Combines endpoint values and derivatives into an interpolated value. - fn hermite_combine( - h00: f64, - h10_dt: f64, - h01: f64, - h11_dt: f64, + fn hermite_value( + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, @@ -27,12 +73,8 @@ pub trait HermiteInterpolable: Sized { ) -> Self; /// Combines endpoint values and derivatives into an interpolated derivative. - fn hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - inv_dt: f64, + fn hermite_derivative( + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, @@ -50,39 +92,71 @@ pub trait HermiteInterpolable: Sized { } } -impl HermiteInterpolable for f64 { +impl HermiteInterpolable for f64 { + type Derivative = f64; + + #[inline] + fn hermite_value( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + basis.h00 * y0 + basis.h10_dx * dy0 + basis.h01 * y1 + basis.h11_dx * dy1 + } + + #[inline] + fn hermite_derivative( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + basis.dh00_over_dt * basis.inv_raw_dx * y0 + + basis.dh10_over_dt * dy0 + + basis.dh01_over_dt * basis.inv_raw_dx * y1 + + basis.dh11_over_dt * dy1 + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.is_finite() + } +} + +impl HermiteInterpolable> for f64 { type Derivative = f64; #[inline] - fn hermite_combine( - h00: f64, - h10_dt: f64, - h01: f64, - h11_dt: f64, + fn hermite_value( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self { - h00 * y0 + h10_dt * dy0 + h01 * y1 + h11_dt * dy1 + basis.h00 * y0 + basis.h10_dx.value() * dy0 + basis.h01 * y1 + basis.h11_dx.value() * dy1 } #[inline] - fn hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - inv_dt: f64, + fn hermite_derivative( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - dh00_over_dt * inv_dt * y0 - + dh10_over_dt * dy0 - + dh01_over_dt * inv_dt * y1 - + dh11_over_dt * dy1 + basis.dh00_over_dt * basis.inv_raw_dx * y0 + + basis.dh10_over_dt * dy0 + + basis.dh01_over_dt * basis.inv_raw_dx * y1 + + basis.dh11_over_dt * dy1 } #[inline] @@ -96,39 +170,76 @@ impl HermiteInterpolable for f64 { } } -impl HermiteInterpolable for Quantity { +impl HermiteInterpolable for Quantity { type Derivative = Quantity; #[inline] - fn hermite_combine( - h00: f64, - h10_dt: f64, - h01: f64, - h11_dt: f64, + fn hermite_value( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + y0 * basis.h00 + dy0 * basis.h10_dx + y1 * basis.h01 + dy1 * basis.h11_dx + } + + #[inline] + fn hermite_derivative( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + y0 * (basis.dh00_over_dt * basis.inv_raw_dx) + + dy0 * basis.dh10_over_dt + + y1 * (basis.dh01_over_dt * basis.inv_raw_dx) + + dy1 * basis.dh11_over_dt + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.value().is_finite() + } +} + +impl HermiteInterpolable> for Quantity +where + X: Unit, + U: Unit + UnitDiv, + >::Output: Unit + UnitMul, +{ + type Derivative = Quantity<>::Output>; + + #[inline] + fn hermite_value( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self { - y0 * h00 + dy0 * h10_dt + y1 * h01 + dy1 * h11_dt + y0 * basis.h00 + (dy0 * basis.h10_dx) + y1 * basis.h01 + (dy1 * basis.h11_dx) } #[inline] - fn hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - inv_dt: f64, + fn hermite_derivative( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - y0 * (dh00_over_dt * inv_dt) - + dy0 * dh10_over_dt - + y1 * (dh01_over_dt * inv_dt) - + dy1 * dh11_over_dt + (y0 / basis.dx) * basis.dh00_over_dt + + dy0 * basis.dh10_over_dt + + (y1 / basis.dx) * basis.dh01_over_dt + + dy1 * basis.dh11_over_dt } #[inline] @@ -142,7 +253,7 @@ impl HermiteInterpolable for Quantity { } } -impl HermiteInterpolable for Vector +impl HermiteInterpolable for Vector where F: ReferenceFrame, U: Unit, @@ -150,35 +261,31 @@ where type Derivative = Vector; #[inline] - fn hermite_combine( - h00: f64, - h10_dt: f64, - h01: f64, - h11_dt: f64, + fn hermite_value( + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self { - y0.scale(h00) + dy0.scale(h10_dt) + y1.scale(h01) + dy1.scale(h11_dt) + y0.scale(basis.h00) + + dy0.scale(basis.h10_dx) + + y1.scale(basis.h01) + + dy1.scale(basis.h11_dx) } #[inline] - fn hermite_derivative_combine( - dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - inv_dt: f64, + fn hermite_derivative( + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - y0.scale(dh00_over_dt * inv_dt) - + dy0.scale(dh10_over_dt) - + y1.scale(dh01_over_dt * inv_dt) - + dy1.scale(dh11_over_dt) + y0.scale(basis.dh00_over_dt * basis.inv_raw_dx) + + dy0.scale(basis.dh10_over_dt) + + y1.scale(basis.dh01_over_dt * basis.inv_raw_dx) + + dy1.scale(basis.dh11_over_dt) } #[inline] @@ -194,7 +301,54 @@ where } } -impl HermiteInterpolable for Position +impl HermiteInterpolable> for Vector +where + X: Unit, + F: ReferenceFrame, + U: Unit + UnitDiv, + >::Output: Unit + UnitMul, +{ + type Derivative = Vector>::Output>; + + #[inline] + fn hermite_value( + basis: HermiteBasis>, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + y0.scale(basis.h00) + (dy0 * basis.h10_dx) + y1.scale(basis.h01) + (dy1 * basis.h11_dx) + } + + #[inline] + fn hermite_derivative( + basis: HermiteBasis>, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + y0.div_quantity(basis.dx).scale(basis.dh00_over_dt) + + dy0.scale(basis.dh10_over_dt) + + y1.div_quantity(basis.dx).scale(basis.dh01_over_dt) + + dy1.scale(basis.dh11_over_dt) + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.x().value().is_finite() + && derivative.y().value().is_finite() + && derivative.z().value().is_finite() + } +} + +impl HermiteInterpolable for Position where C: ReferenceCenter, F: ReferenceFrame, @@ -203,34 +357,78 @@ where type Derivative = Vector; #[inline] - fn hermite_combine( - _h00: f64, - h10_dt: f64, - h01: f64, - h11_dt: f64, + fn hermite_value( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self { + let chord = y1 - y0; + y0 + chord.scale(basis.h01) + dy0.scale(basis.h10_dx) + dy1.scale(basis.h11_dx) + } + + #[inline] + fn hermite_derivative( + basis: HermiteBasis, + y0: Self, + dy0: Self::Derivative, + y1: Self, + dy1: Self::Derivative, + ) -> Self::Derivative { + let chord = y1 - y0; + chord.scale(basis.dh01_over_dt * basis.inv_raw_dx) + + dy0.scale(basis.dh10_over_dt) + + dy1.scale(basis.dh11_over_dt) + } + + #[inline] + fn hermite_value_is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + + #[inline] + fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { + derivative.x().value().is_finite() + && derivative.y().value().is_finite() + && derivative.z().value().is_finite() + } +} + +impl HermiteInterpolable> for Position +where + X: Unit, + C: ReferenceCenter, + F: ReferenceFrame, + L: LengthUnit + UnitDiv, + >::Output: Unit + UnitMul, +{ + type Derivative = Vector>::Output>; + + #[inline] + fn hermite_value( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self { let chord = y1 - y0; - y0 + chord.scale(h01) + dy0.scale(h10_dt) + dy1.scale(h11_dt) + y0 + chord.scale(basis.h01) + (dy0 * basis.h10_dx) + (dy1 * basis.h11_dx) } #[inline] - fn hermite_derivative_combine( - _dh00_over_dt: f64, - dh10_over_dt: f64, - dh01_over_dt: f64, - dh11_over_dt: f64, - inv_dt: f64, + fn hermite_derivative( + basis: HermiteBasis>, y0: Self, dy0: Self::Derivative, y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { let chord = y1 - y0; - chord.scale(dh01_over_dt * inv_dt) + dy0.scale(dh10_over_dt) + dy1.scale(dh11_over_dt) + chord.div_quantity(basis.dx).scale(basis.dh01_over_dt) + + dy0.scale(basis.dh10_over_dt) + + dy1.scale(basis.dh11_over_dt) } #[inline] diff --git a/src/lib.rs b/src/lib.rs index 6f2b99d..7f63912 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -208,10 +208,9 @@ pub mod prelude { SemiMajorAxisParam, TypedPeriapsisParam, TypedSemiMajorAxisParam, }; pub use crate::interpolation::{ - cubic_hermite_segment, CubicHermiteQuantityTable, CubicHermiteSpline, CubicHermiteTable, - HermiteEvaluation, HermiteInterpolable, HermiteNode, HermiteSample, HermiteTableEvaluation, - InterpolationError, QuantityHermiteInterpolable, QuantityHermiteNode, - QuantityHermiteTableEvaluation, + cubic_hermite_segment, CubicHermiteSpline, CubicHermiteTable, HermiteEvaluation, + HermiteInterpolable, HermiteNode, HermiteSample, HermiteTableEvaluation, + InterpolationAbscissa, InterpolationError, ScalarHermiteNode, ScalarHermiteTableEvaluation, }; pub use crate::spherical::{Direction as SphericalDirection, Position as SphericalPosition}; diff --git a/tests/interpolation.rs b/tests/interpolation.rs index f8c7a8b..e9b4082 100644 --- a/tests/interpolation.rs +++ b/tests/interpolation.rs @@ -1,7 +1,6 @@ use affn::cartesian::{Position, Velocity}; use affn::interpolation::{ - CubicHermiteQuantityTable, CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, - InterpolationError, QuantityHermiteNode, + CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, InterpolationError, }; use affn::{DeriveReferenceCenter as ReferenceCenter, DeriveReferenceFrame as ReferenceFrame}; use qtty::unit::{Kilometer, Meter, Second}; @@ -39,6 +38,15 @@ fn cubic_velocity(t: f64) -> TestKilometerVelocity { ) } +fn expect_table_error( + result: Result, InterpolationError>, +) -> InterpolationError { + match result { + Ok(_) => panic!("table construction should fail"), + Err(err) => err, + } +} + #[test] fn scalar_cubic_polynomial_is_reproduced() { let spline = CubicHermiteSpline::new(vec![ @@ -90,7 +98,7 @@ fn exact_node_evaluation_returns_node_value_and_derivative() { #[test] fn linear_motion_with_constant_velocity_is_exact() { - let table = CubicHermiteTable::new(vec![ + let table = CubicHermiteTable::::new(vec![ HermiteNode { x: 0.0, value: TestPosition::new(1.0, 2.0, 3.0), @@ -114,9 +122,9 @@ fn linear_motion_with_constant_velocity_is_exact() { } #[test] -fn quantity_table_accepts_position_over_seconds_with_velocity() { - let table = CubicHermiteQuantityTable::::new(vec![ - QuantityHermiteNode { +fn typed_abscissa_table_accepts_position_over_seconds_with_velocity() { + let table = CubicHermiteTable::::new(vec![ + HermiteNode { x: qtty::Second::new(0.0), value: TestKilometerPosition::new(1.0, 2.0, 3.0), derivative: TestKilometerVelocity::new( @@ -125,7 +133,7 @@ fn quantity_table_accepts_position_over_seconds_with_velocity() { Quantity::::new(2.0), ), }, - QuantityHermiteNode { + HermiteNode { x: qtty::Second::new(4.0), value: TestKilometerPosition::new(3.0, -2.0, 11.0), derivative: TestKilometerVelocity::new( @@ -147,19 +155,19 @@ fn quantity_table_accepts_position_over_seconds_with_velocity() { } #[test] -fn quantity_table_reproduces_cubic_position_over_seconds() { - let table = CubicHermiteQuantityTable::::new(vec![ - QuantityHermiteNode { +fn typed_abscissa_table_reproduces_cubic_position_over_seconds() { + let table = CubicHermiteTable::::new(vec![ + HermiteNode { x: qtty::Second::new(-1.0), value: cubic_position(-1.0), derivative: cubic_velocity(-1.0), }, - QuantityHermiteNode { + HermiteNode { x: qtty::Second::new(0.5), value: cubic_position(0.5), derivative: cubic_velocity(0.5), }, - QuantityHermiteNode { + HermiteNode { x: qtty::Second::new(2.0), value: cubic_position(2.0), derivative: cubic_velocity(2.0), @@ -180,6 +188,86 @@ fn quantity_table_reproduces_cubic_position_over_seconds() { } } +#[test] +fn typed_abscissa_table_rejects_non_finite_abscissa() { + let err = expect_table_error( + CubicHermiteTable::::new(vec![ + HermiteNode { + x: qtty::Second::new(f64::NAN), + value: cubic_position(0.0), + derivative: cubic_velocity(0.0), + }, + HermiteNode { + x: qtty::Second::new(1.0), + value: cubic_position(1.0), + derivative: cubic_velocity(1.0), + }, + ]), + ); + + assert_eq!(err, InterpolationError::NonFiniteAbscissa); +} + +#[test] +fn typed_abscissa_table_rejects_non_finite_position_component() { + let err = expect_table_error( + CubicHermiteTable::::new(vec![ + HermiteNode { + x: qtty::Second::new(0.0), + value: TestKilometerPosition::new(f64::NAN, 0.0, 0.0), + derivative: cubic_velocity(0.0), + }, + HermiteNode { + x: qtty::Second::new(1.0), + value: cubic_position(1.0), + derivative: cubic_velocity(1.0), + }, + ]), + ); + + assert_eq!(err, InterpolationError::NonFiniteValue); +} + +#[test] +fn typed_abscissa_table_rejects_duplicate_seconds() { + let err = expect_table_error( + CubicHermiteTable::::new(vec![ + HermiteNode { + x: qtty::Second::new(0.0), + value: cubic_position(0.0), + derivative: cubic_velocity(0.0), + }, + HermiteNode { + x: qtty::Second::new(0.0), + value: cubic_position(1.0), + derivative: cubic_velocity(1.0), + }, + ]), + ); + + assert_eq!(err, InterpolationError::DuplicateAbscissa); +} + +#[test] +fn vector_dimensional_mul_and_div_quantity_work() { + let velocity = TestKilometerVelocity::new( + Quantity::::new(2.0), + Quantity::::new(-3.0), + Quantity::::new(4.0), + ); + let elapsed = qtty::Second::new(5.0); + + let displacement = velocity * elapsed; + assert!((displacement.x().value() - 10.0).abs() < 1e-12); + assert!((displacement.y().value() + 15.0).abs() < 1e-12); + assert!((displacement.z().value() - 20.0).abs() < 1e-12); + + let recovered_velocity = displacement.div_quantity(elapsed); + assert!((recovered_velocity.x().value() - 2.0).abs() < 1e-12); + assert!((recovered_velocity.y().value() + 3.0).abs() < 1e-12); + assert!((recovered_velocity.z().value() - 4.0).abs() < 1e-12); +} + #[test] fn non_uniform_sample_spacing_works() { let spline = CubicHermiteSpline::new(vec![ From d427e9d2edd772573630d5471a8d993f63a1fbc3 Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 21:15:30 +0200 Subject: [PATCH 3/6] Refactor interpolation module: remove linear interpolation and related structures, update cubic Hermite references --- src/interpolation/cubic_hermite.rs | 97 ++--------------------- src/interpolation/linear.rs | 34 -------- src/interpolation/mod.rs | 6 +- src/interpolation/scalar.rs | 87 -------------------- src/lib.rs | 6 +- tests/interpolation.rs | 123 ++++++++++++++++------------- 6 files changed, 77 insertions(+), 276 deletions(-) delete mode 100644 src/interpolation/linear.rs delete mode 100644 src/interpolation/scalar.rs diff --git a/src/interpolation/cubic_hermite.rs b/src/interpolation/cubic_hermite.rs index b2dc17e..9b246be 100644 --- a/src/interpolation/cubic_hermite.rs +++ b/src/interpolation/cubic_hermite.rs @@ -1,95 +1,9 @@ //! Cubic Hermite interpolation tables. use super::error::InterpolationError; -use super::scalar::{cubic_hermite_segment, HermiteEvaluation}; use super::traits::{HermiteBasis, HermiteInterpolable}; use super::InterpolationAbscissa; -/// A scalar Hermite table sample. -#[derive(Debug, Clone, Copy, PartialEq)] -pub struct HermiteSample { - /// Sample abscissa. - pub x: f64, - /// Sample value. - pub y: f64, - /// Sample derivative with respect to `x`. - pub dydx: f64, -} - -/// Piecewise scalar cubic Hermite spline. -/// -/// Derivative continuity is guaranteed only up to the first derivative supplied -/// at each sample. Acceleration or other second derivative quantities are not -/// generally continuous across sample boundaries. -#[derive(Debug, Clone, PartialEq)] -pub struct CubicHermiteSpline { - samples: Vec, -} - -impl CubicHermiteSpline { - /// Builds a spline from samples sorted by strictly increasing `x`. - pub fn new(samples: Vec) -> Result { - validate_len(samples.len())?; - for sample in &samples { - if !sample.x.is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - if !sample.y.is_finite() || !sample.dydx.is_finite() { - return Err(InterpolationError::NonFiniteValue); - } - } - validate_sorted(samples.iter().map(|sample| sample.x))?; - Ok(Self { samples }) - } - - /// Returns the table samples. - pub fn samples(&self) -> &[HermiteSample] { - &self.samples - } - - /// Evaluates the spline without extrapolation. - pub fn evaluate(&self, x: f64) -> Result { - if !x.is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - let (min, max) = self.range(); - if x < min || x > max { - return Err(InterpolationError::OutOfRange { x, min, max }); - } - - let segment = self.segment_index(x); - let s0 = self.samples[segment]; - let s1 = self.samples[segment + 1]; - if x == s0.x { - return Ok(HermiteEvaluation { - value: s0.y, - derivative: s0.dydx, - }); - } - if x == s1.x { - return Ok(HermiteEvaluation { - value: s1.y, - derivative: s1.dydx, - }); - } - cubic_hermite_segment(x, s0.x, s1.x, s0.y, s0.dydx, s1.y, s1.dydx) - } - - fn range(&self) -> (f64, f64) { - (self.samples[0].x, self.samples[self.samples.len() - 1].x) - } - - fn segment_index(&self, x: f64) -> usize { - match self - .samples - .binary_search_by(|sample| sample.x.total_cmp(&x)) - { - Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), - Err(index) => (index - 1).min(self.samples.len() - 2), - } - } -} - /// A Hermite table node. #[derive(Debug, Clone, PartialEq)] pub struct HermiteNode @@ -237,11 +151,14 @@ where } } -/// Explicit scalar table node alias. -pub type ScalarHermiteNode = HermiteNode; +/// Scalar cubic Hermite table alias. +pub type ScalarCubicHermiteTable = CubicHermiteTable; + +/// Scalar Hermite table node alias. +pub type ScalarHermiteNode = HermiteNode; -/// Explicit scalar table evaluation alias. -pub type ScalarHermiteTableEvaluation = HermiteTableEvaluation; +/// Scalar Hermite table evaluation alias. +pub type ScalarHermiteTableEvaluation = HermiteTableEvaluation; fn validate_len(len: usize) -> Result<(), InterpolationError> { if len == 0 { diff --git a/src/interpolation/linear.rs b/src/interpolation/linear.rs deleted file mode 100644 index f4a62c6..0000000 --- a/src/interpolation/linear.rs +++ /dev/null @@ -1,34 +0,0 @@ -//! Linear interpolation helpers. - -use super::error::InterpolationError; - -/// Evaluates scalar linear interpolation on one segment. -pub fn linear_segment( - x: f64, - x0: f64, - x1: f64, - y0: f64, - y1: f64, -) -> Result { - if !x.is_finite() || !x0.is_finite() || !x1.is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - if !y0.is_finite() || !y1.is_finite() { - return Err(InterpolationError::NonFiniteValue); - } - if x1 <= x0 { - if x1 == x0 { - return Err(InterpolationError::DuplicateAbscissa); - } - return Err(InterpolationError::UnsortedAbscissa); - } - if x < x0 || x > x1 { - return Err(InterpolationError::OutOfRange { - x, - min: x0, - max: x1, - }); - } - let t = (x - x0) / (x1 - x0); - Ok((1.0 - t) * y0 + t * y1) -} diff --git a/src/interpolation/mod.rs b/src/interpolation/mod.rs index c93bc4f..4a2ab04 100644 --- a/src/interpolation/mod.rs +++ b/src/interpolation/mod.rs @@ -6,16 +6,12 @@ pub mod abscissa; pub mod cubic_hermite; pub mod error; -pub mod linear; -pub mod scalar; pub mod traits; pub use abscissa::InterpolationAbscissa; pub use cubic_hermite::{ - CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, HermiteTableEvaluation, + CubicHermiteTable, HermiteNode, HermiteTableEvaluation, ScalarCubicHermiteTable, ScalarHermiteNode, ScalarHermiteTableEvaluation, }; pub use error::InterpolationError; -pub use linear::linear_segment; -pub use scalar::{cubic_hermite_segment, HermiteEvaluation}; pub use traits::{HermiteBasis, HermiteInterpolable}; diff --git a/src/interpolation/scalar.rs b/src/interpolation/scalar.rs deleted file mode 100644 index a76f100..0000000 --- a/src/interpolation/scalar.rs +++ /dev/null @@ -1,87 +0,0 @@ -//! Scalar interpolation primitives. - -use super::error::InterpolationError; - -/// Scalar Hermite interpolation result. -#[derive(Debug, Clone, Copy, PartialEq)] -pub struct HermiteEvaluation { - /// Interpolated value. - pub value: f64, - /// Derivative with respect to the original abscissa. - pub derivative: f64, -} - -/// Evaluates one scalar cubic Hermite segment. -/// -/// The returned derivative is with respect to `x`, not normalized `t`. -pub fn cubic_hermite_segment( - x: f64, - x0: f64, - x1: f64, - y0: f64, - dy0: f64, - y1: f64, - dy1: f64, -) -> Result { - if !x.is_finite() || !x0.is_finite() || !x1.is_finite() { - return Err(InterpolationError::NonFiniteAbscissa); - } - if !y0.is_finite() || !dy0.is_finite() || !y1.is_finite() || !dy1.is_finite() { - return Err(InterpolationError::NonFiniteValue); - } - if x1 <= x0 { - if x1 == x0 { - return Err(InterpolationError::DuplicateAbscissa); - } - return Err(InterpolationError::UnsortedAbscissa); - } - if x < x0 || x > x1 { - return Err(InterpolationError::OutOfRange { - x, - min: x0, - max: x1, - }); - } - - let dt = x1 - x0; - let t = (x - x0) / dt; - let basis = hermite_basis(t, dt); - - Ok(HermiteEvaluation { - value: basis.h00 * y0 + basis.h10_dt * dy0 + basis.h01 * y1 + basis.h11_dt * dy1, - derivative: basis.dh00_over_dt * basis.inv_dt * y0 - + basis.dh10_over_dt * dy0 - + basis.dh01_over_dt * basis.inv_dt * y1 - + basis.dh11_over_dt * dy1, - }) -} - -#[derive(Debug, Clone, Copy)] -pub(crate) struct HermiteBasis { - pub h00: f64, - pub h10_dt: f64, - pub h01: f64, - pub h11_dt: f64, - pub dh00_over_dt: f64, - pub dh10_over_dt: f64, - pub dh01_over_dt: f64, - pub dh11_over_dt: f64, - pub inv_dt: f64, -} - -pub(crate) fn hermite_basis(t: f64, dt: f64) -> HermiteBasis { - let t2 = t * t; - let t3 = t2 * t; - - HermiteBasis { - h00: 2.0 * t3 - 3.0 * t2 + 1.0, - h10_dt: (t3 - 2.0 * t2 + t) * dt, - h01: -2.0 * t3 + 3.0 * t2, - h11_dt: (t3 - t2) * dt, - dh00_over_dt: 6.0 * t2 - 6.0 * t, - dh10_over_dt: 3.0 * t2 - 4.0 * t + 1.0, - dh01_over_dt: -6.0 * t2 + 6.0 * t, - dh11_over_dt: 3.0 * t2 - 2.0 * t, - inv_dt: 1.0 / dt, - } -} diff --git a/src/lib.rs b/src/lib.rs index 7f63912..3c45d0f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -208,9 +208,9 @@ pub mod prelude { SemiMajorAxisParam, TypedPeriapsisParam, TypedSemiMajorAxisParam, }; pub use crate::interpolation::{ - cubic_hermite_segment, CubicHermiteSpline, CubicHermiteTable, HermiteEvaluation, - HermiteInterpolable, HermiteNode, HermiteSample, HermiteTableEvaluation, - InterpolationAbscissa, InterpolationError, ScalarHermiteNode, ScalarHermiteTableEvaluation, + CubicHermiteTable, HermiteInterpolable, HermiteNode, HermiteTableEvaluation, + InterpolationAbscissa, InterpolationError, ScalarCubicHermiteTable, ScalarHermiteNode, + ScalarHermiteTableEvaluation, }; pub use crate::spherical::{Direction as SphericalDirection, Position as SphericalPosition}; diff --git a/tests/interpolation.rs b/tests/interpolation.rs index e9b4082..a17d343 100644 --- a/tests/interpolation.rs +++ b/tests/interpolation.rs @@ -1,6 +1,6 @@ use affn::cartesian::{Position, Velocity}; use affn::interpolation::{ - CubicHermiteSpline, CubicHermiteTable, HermiteNode, HermiteSample, InterpolationError, + CubicHermiteTable, HermiteNode, InterpolationError, ScalarCubicHermiteTable, ScalarHermiteNode, }; use affn::{DeriveReferenceCenter as ReferenceCenter, DeriveReferenceFrame as ReferenceFrame}; use qtty::unit::{Kilometer, Meter, Second}; @@ -47,29 +47,38 @@ fn expect_table_error( } } +fn expect_scalar_table_error( + result: Result, +) -> InterpolationError { + match result { + Ok(_) => panic!("scalar table construction should fail"), + Err(err) => err, + } +} + #[test] fn scalar_cubic_polynomial_is_reproduced() { - let spline = CubicHermiteSpline::new(vec![ - HermiteSample { + let table = ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: -1.0, - y: cubic(-1.0), - dydx: cubic_derivative(-1.0), + value: cubic(-1.0), + derivative: cubic_derivative(-1.0), }, - HermiteSample { + ScalarHermiteNode { x: 0.5, - y: cubic(0.5), - dydx: cubic_derivative(0.5), + value: cubic(0.5), + derivative: cubic_derivative(0.5), }, - HermiteSample { + ScalarHermiteNode { x: 2.0, - y: cubic(2.0), - dydx: cubic_derivative(2.0), + value: cubic(2.0), + derivative: cubic_derivative(2.0), }, ]) .unwrap(); for x in [-0.75, -0.1, 0.25, 1.25, 1.75] { - let evaluated = spline.evaluate(x).unwrap(); + let evaluated = table.evaluate(x).unwrap(); assert!((evaluated.value - cubic(x)).abs() < 1e-12); assert!((evaluated.derivative - cubic_derivative(x)).abs() < 1e-12); } @@ -77,21 +86,21 @@ fn scalar_cubic_polynomial_is_reproduced() { #[test] fn exact_node_evaluation_returns_node_value_and_derivative() { - let spline = CubicHermiteSpline::new(vec![ - HermiteSample { + let table = ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: 0.0, - y: 10.0, - dydx: -3.0, + value: 10.0, + derivative: -3.0, }, - HermiteSample { + ScalarHermiteNode { x: 2.0, - y: 20.0, - dydx: 4.0, + value: 20.0, + derivative: 4.0, }, ]) .unwrap(); - let evaluated = spline.evaluate(2.0).unwrap(); + let evaluated = table.evaluate(2.0).unwrap(); assert_eq!(evaluated.value, 20.0); assert_eq!(evaluated.derivative, 4.0); } @@ -270,48 +279,48 @@ fn vector_dimensional_mul_and_div_quantity_work() { #[test] fn non_uniform_sample_spacing_works() { - let spline = CubicHermiteSpline::new(vec![ - HermiteSample { + let table = ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: 0.0, - y: cubic(0.0), - dydx: cubic_derivative(0.0), + value: cubic(0.0), + derivative: cubic_derivative(0.0), }, - HermiteSample { + ScalarHermiteNode { x: 0.25, - y: cubic(0.25), - dydx: cubic_derivative(0.25), + value: cubic(0.25), + derivative: cubic_derivative(0.25), }, - HermiteSample { + ScalarHermiteNode { x: 2.5, - y: cubic(2.5), - dydx: cubic_derivative(2.5), + value: cubic(2.5), + derivative: cubic_derivative(2.5), }, ]) .unwrap(); - let evaluated = spline.evaluate(1.75).unwrap(); + let evaluated = table.evaluate(1.75).unwrap(); assert!((evaluated.value - cubic(1.75)).abs() < 1e-12); assert!((evaluated.derivative - cubic_derivative(1.75)).abs() < 1e-12); } #[test] fn out_of_range_queries_return_error() { - let spline = CubicHermiteSpline::new(vec![ - HermiteSample { + let table = ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: 0.0, - y: 0.0, - dydx: 1.0, + value: 0.0, + derivative: 1.0, }, - HermiteSample { + ScalarHermiteNode { x: 1.0, - y: 1.0, - dydx: 1.0, + value: 1.0, + derivative: 1.0, }, ]) .unwrap(); assert_eq!( - spline.evaluate(2.0), + table.evaluate(2.0), Err(InterpolationError::OutOfRange { x: 2.0, min: 0.0, @@ -323,37 +332,37 @@ fn out_of_range_queries_return_error() { #[test] fn duplicate_abscissae_are_rejected() { assert_eq!( - CubicHermiteSpline::new(vec![ - HermiteSample { + expect_scalar_table_error(ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: 0.0, - y: 0.0, - dydx: 0.0, + value: 0.0, + derivative: 0.0, }, - HermiteSample { + ScalarHermiteNode { x: 0.0, - y: 1.0, - dydx: 1.0, + value: 1.0, + derivative: 1.0, }, - ]), - Err(InterpolationError::DuplicateAbscissa) + ])), + InterpolationError::DuplicateAbscissa ); } #[test] fn unsorted_abscissae_are_rejected() { assert_eq!( - CubicHermiteSpline::new(vec![ - HermiteSample { + expect_scalar_table_error(ScalarCubicHermiteTable::new(vec![ + ScalarHermiteNode { x: 1.0, - y: 1.0, - dydx: 1.0, + value: 1.0, + derivative: 1.0, }, - HermiteSample { + ScalarHermiteNode { x: 0.0, - y: 0.0, - dydx: 0.0, + value: 0.0, + derivative: 0.0, }, - ]), - Err(InterpolationError::UnsortedAbscissa) + ])), + InterpolationError::UnsortedAbscissa ); } From 85a93685f520f6ad313a9ff78cbe08fba6f8713b Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 21:53:08 +0200 Subject: [PATCH 4/6] Refactor interpolation traits and error handling - Introduced `AbscissaDelta` trait for scaling and diagnostics of abscissae. - Updated `InterpolationAbscissa` to use `AbscissaDelta` for delta calculations. - Enhanced error handling in `InterpolationError` to provide raw diagnostic values for out-of-range queries. - Refactored cubic Hermite interpolation to utilize typed `qtty::Quantity` for abscissae. - Improved documentation for clarity on usage and expectations. - Adjusted tests to reflect changes in abscissa handling and error reporting. --- src/cartesian/direction.rs | 2 +- src/cartesian/position.rs | 6 + src/cartesian/vector.rs | 6 + src/cartesian/xyz.rs | 4 +- src/ellipsoidal/position.rs | 4 +- src/frames/mod.rs | 2 +- src/interpolation/abscissa.rs | 116 +++++++++++-- src/interpolation/cubic_hermite.rs | 121 ++++++++------ src/interpolation/error.rs | 23 ++- src/interpolation/mod.rs | 2 +- src/interpolation/traits.rs | 260 ++++++++++++++--------------- src/lib.rs | 4 +- src/ops/rotation.rs | 2 +- tests/interpolation.rs | 68 ++++---- 14 files changed, 370 insertions(+), 250 deletions(-) diff --git a/src/cartesian/direction.rs b/src/cartesian/direction.rs index 2fc631c..2463089 100644 --- a/src/cartesian/direction.rs +++ b/src/cartesian/direction.rs @@ -75,7 +75,7 @@ use std::ops::Mul; /// # Invariants /// /// All public constructors ensure the direction is normalized. For unchecked -/// construction, use [`from_xyz_unchecked`](Self::from_xyz_unchecked). +/// construction, use the crate-internal unchecked constructor. /// /// # Zero-Cost Abstraction /// diff --git a/src/cartesian/position.rs b/src/cartesian/position.rs index e2e7974..c7a4cd8 100644 --- a/src/cartesian/position.rs +++ b/src/cartesian/position.rs @@ -276,6 +276,12 @@ impl Position { self.xyz.z() } + /// Returns whether all position components are finite. + #[inline] + pub fn is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + /// Returns a reference to the underlying `[Quantity; 3]` array. #[inline] pub fn as_array(&self) -> &[Quantity; 3] { diff --git a/src/cartesian/vector.rs b/src/cartesian/vector.rs index e6b48c2..778b905 100644 --- a/src/cartesian/vector.rs +++ b/src/cartesian/vector.rs @@ -200,6 +200,12 @@ impl Vector { Self::from_xyz(-self.xyz) } + /// Returns whether all vector components are finite. + #[inline] + pub fn is_finite(&self) -> bool { + self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + } + /// Divides this vector by a typed quantity, carrying the resulting unit in /// the output vector. /// diff --git a/src/cartesian/xyz.rs b/src/cartesian/xyz.rs index 1f25c8f..5f71b87 100644 --- a/src/cartesian/xyz.rs +++ b/src/cartesian/xyz.rs @@ -277,7 +277,7 @@ impl XYZ> { Quantity::new((x * x + y * y + z * z).sqrt()) } - /// Extracts raw f64 values as an XYZ. + /// Extracts raw f64 values as an `XYZ`. #[inline] pub fn to_raw(&self) -> XYZ { XYZ::new(self.0[0].value(), self.0[1].value(), self.0[2].value()) @@ -320,7 +320,7 @@ impl> XYZ> { /// This relies on `qtty`'s squared-unit support: for any `U: UnitMul` /// the output unit is `>::Output` (e.g. `Prod` /// for length). If you only need a raw `f64` for ordering or comparison, - /// use the deprecated [`magnitude_squared_raw`](Self::magnitude_squared_raw). + /// use [`to_raw`](Self::to_raw) and [`XYZ::magnitude_squared`]. #[inline] pub fn magnitude_squared(&self) -> Quantity<>::Output> { let x = self.0[0].value(); diff --git a/src/ellipsoidal/position.rs b/src/ellipsoidal/position.rs index f77dea4..8de1ceb 100644 --- a/src/ellipsoidal/position.rs +++ b/src/ellipsoidal/position.rs @@ -13,7 +13,7 @@ //! ## Type Parameters //! //! - `C`: Reference center (defines the origin) -//! - `F`: Reference frame (carries the ellipsoid via [`HasEllipsoid`](crate::ellipsoid::HasEllipsoid)) +//! - `F`: Reference frame (carries the ellipsoid via [`HasEllipsoid`]) //! - `U`: Length unit for the height (defaults to [`Meter`]) //! //! ## Conversion to Cartesian @@ -99,7 +99,7 @@ impl std::error::Error for GeodeticConvergenceError {} /// /// - `C`: The reference center (e.g., `Geocentric`) /// - `F`: The reference frame; determines the ellipsoid via -/// [`HasEllipsoid`](crate::ellipsoid::HasEllipsoid) +/// [`HasEllipsoid`] /// - `U`: The length unit for the height (defaults to [`Meter`]) /// /// # Field conventions diff --git a/src/frames/mod.rs b/src/frames/mod.rs index c94cdd8..d04f6f8 100644 --- a/src/frames/mod.rs +++ b/src/frames/mod.rs @@ -67,7 +67,7 @@ pub trait ReferenceFrame: Copy + Clone + std::fmt::Debug { /// /// When a frame implements [`SphericalNaming`], this returns /// `Some((polar_name, azimuth_name, distance_name))` so that Display - /// impls on [`spherical::Direction`] and [`spherical::Position`] can use + /// impls on [`crate::spherical::Direction`] and [`crate::spherical::Position`] can use /// domain-appropriate axis names (e.g. `"dec"`/`"ra"` for ICRS, /// `"alt"`/`"az"` for Horizontal) rather than generic Greek letters. /// diff --git a/src/interpolation/abscissa.rs b/src/interpolation/abscissa.rs index cc9afb6..b590c60 100644 --- a/src/interpolation/abscissa.rs +++ b/src/interpolation/abscissa.rs @@ -1,34 +1,99 @@ //! Interpolation abscissa types. +use super::InterpolationError; use qtty::{Quantity, Unit}; +use std::cmp::Ordering; + +mod sealed { + pub trait Sealed {} + impl Sealed for f64 {} + impl Sealed for qtty::Quantity {} +} + +/// Difference between two interpolation abscissae. +pub trait AbscissaDelta: Copy { + /// Scales this delta by a dimensionless Hermite basis coefficient. + fn scale(self, factor: f64) -> Self; + + /// Raw value used only for diagnostics and scalar-only formulas. + fn diagnostic_raw(self) -> f64; +} + +impl AbscissaDelta for f64 { + #[inline] + fn scale(self, factor: f64) -> Self { + self * factor + } + + #[inline] + fn diagnostic_raw(self) -> f64 { + self + } +} + +impl AbscissaDelta for Quantity { + #[inline] + fn scale(self, factor: f64) -> Self { + self * factor + } + + #[inline] + fn diagnostic_raw(self) -> f64 { + self.value() + } +} /// An ordered interpolation abscissa. /// -/// Implemented for raw scalar parameters and typed `qtty` quantities. New -/// interpolation APIs should use this trait instead of exposing separate raw -/// and quantity table types. -pub trait InterpolationAbscissa: Copy { +/// Physical domains should use typed [`qtty::Quantity`] abscissae. The raw +/// `f64` implementation is intended for explicitly scalar, dimensionless +/// interpolation such as `CubicHermiteTable`. +pub trait InterpolationAbscissa: sealed::Sealed + Copy { /// Difference between two abscissae. - type Delta: Copy; + type Delta: AbscissaDelta; + + /// Returns whether this abscissa is finite. + fn is_finite(self) -> bool; + + /// Compares two abscissae for table ordering. + fn cmp_abscissa(self, other: Self) -> Ordering; - /// Returns the raw scalar used for ordering and normalized interpolation. - fn raw(self) -> f64; + /// Returns `self - origin` as a typed delta. + fn delta_since(self, origin: Self) -> Self::Delta; - /// Builds a delta value from a raw difference in this abscissa' stored unit. - fn delta_from_raw(raw: f64) -> Self::Delta; + /// Returns the normalized segment coordinate in `[0, 1]`. + fn normalize_between(self, start: Self, end: Self) -> Result; + + /// Raw value used only for diagnostics. + fn diagnostic_raw(self) -> f64; } impl InterpolationAbscissa for f64 { type Delta = f64; #[inline] - fn raw(self) -> f64 { - self + fn is_finite(self) -> bool { + self.is_finite() + } + + #[inline] + fn cmp_abscissa(self, other: Self) -> Ordering { + self.total_cmp(&other) + } + + #[inline] + fn delta_since(self, origin: Self) -> Self::Delta { + self - origin + } + + #[inline] + fn normalize_between(self, start: Self, end: Self) -> Result { + Ok((self - start) / (end - start)) } #[inline] - fn delta_from_raw(raw: f64) -> Self::Delta { - raw + fn diagnostic_raw(self) -> f64 { + self } } @@ -36,12 +101,29 @@ impl InterpolationAbscissa for Quantity { type Delta = Quantity; #[inline] - fn raw(self) -> f64 { - self.value() + fn is_finite(self) -> bool { + self.value().is_finite() + } + + #[inline] + fn cmp_abscissa(self, other: Self) -> Ordering { + self.value().total_cmp(&other.value()) } #[inline] - fn delta_from_raw(raw: f64) -> Self::Delta { - Quantity::::new(raw) + fn delta_since(self, origin: Self) -> Self::Delta { + self - origin + } + + #[inline] + fn normalize_between(self, start: Self, end: Self) -> Result { + let offset = self - start; + let width = end - start; + Ok(offset.value() / width.value()) + } + + #[inline] + fn diagnostic_raw(self) -> f64 { + self.value() } } diff --git a/src/interpolation/cubic_hermite.rs b/src/interpolation/cubic_hermite.rs index 9b246be..add89cf 100644 --- a/src/interpolation/cubic_hermite.rs +++ b/src/interpolation/cubic_hermite.rs @@ -1,4 +1,17 @@ //! Cubic Hermite interpolation tables. +//! +//! Cubic Hermite interpolation is C1 continuous when neighboring samples carry +//! consistent first derivatives. Tables never extrapolate outside their +//! abscissa coverage. +//! +//! Use typed `qtty::Quantity` abscissae for physical domains such as time, +//! length, or angle. Raw `f64` abscissae are intended only for explicitly +//! scalar interpolation, for example [`ScalarCubicHermiteTable`]. +//! +//! `Position` interpolation is affine-safe: it interpolates from a segment-local +//! origin using a chord displacement and never adds two positions. Position +//! support currently requires centers with `Params = ()`; parameterized centers +//! need a future checked API. use super::error::InterpolationError; use super::traits::{HermiteBasis, HermiteInterpolable}; @@ -6,57 +19,57 @@ use super::InterpolationAbscissa; /// A Hermite table node. #[derive(Debug, Clone, PartialEq)] -pub struct HermiteNode +pub struct HermiteNode where - X: InterpolationAbscissa, - T: HermiteInterpolable, + A: InterpolationAbscissa, + T: HermiteInterpolable, { /// Sample abscissa. - pub x: X, + pub abscissa: A, /// Sample value. pub value: T, - /// Sample derivative with respect to `x`. + /// Sample derivative with respect to `abscissa`. pub derivative: T::Derivative, } /// A Hermite table evaluation. #[derive(Debug, Clone, PartialEq)] -pub struct HermiteTableEvaluation +pub struct HermiteTableEvaluation where - X: InterpolationAbscissa, - T: HermiteInterpolable, + A: InterpolationAbscissa, + T: HermiteInterpolable, { /// Interpolated value. pub value: T, - /// Interpolated derivative with respect to `x`. + /// Interpolated derivative with respect to `abscissa`. pub derivative: T::Derivative, /// Evaluated abscissa. - pub x: X, + pub abscissa: A, } /// Piecewise cubic Hermite interpolation table for typed values. /// -/// `X` may be a raw scalar parameter (`f64`) or a typed `qtty::Quantity` such +/// `A` may be a raw scalar parameter (`f64`) or a typed `qtty::Quantity` such /// as seconds or days. Use typed quantities for physical domains so derivatives /// carry the expected units. -pub struct CubicHermiteTable +pub struct CubicHermiteTable where - X: InterpolationAbscissa, - T: HermiteInterpolable, + A: InterpolationAbscissa, + T: HermiteInterpolable, { - samples: Vec>, + samples: Vec>, } -impl CubicHermiteTable +impl CubicHermiteTable where - X: InterpolationAbscissa, - T: HermiteInterpolable, + A: InterpolationAbscissa, + T: HermiteInterpolable, { - /// Builds a typed table from nodes sorted by strictly increasing `x`. - pub fn new(samples: Vec>) -> Result { + /// Builds a typed table from nodes sorted by strictly increasing abscissa. + pub fn new(samples: Vec>) -> Result { validate_len(samples.len())?; for sample in &samples { - if !sample.x.raw().is_finite() { + if !sample.abscissa.is_finite() { return Err(InterpolationError::NonFiniteAbscissa); } if !sample.value.hermite_value_is_finite() @@ -65,54 +78,60 @@ where return Err(InterpolationError::NonFiniteValue); } } - validate_sorted(samples.iter().map(|sample| sample.x.raw()))?; + validate_sorted(samples.iter().map(|sample| sample.abscissa))?; Ok(Self { samples }) } /// Returns the table samples. - pub fn samples(&self) -> &[HermiteNode] { + pub fn samples(&self) -> &[HermiteNode] { &self.samples } } -impl CubicHermiteTable +impl CubicHermiteTable where - X: InterpolationAbscissa, - T: HermiteInterpolable + Clone, + A: InterpolationAbscissa, + T: HermiteInterpolable + Clone, T::Derivative: Clone, { /// Evaluates the table without extrapolation. - pub fn evaluate(&self, x: X) -> Result, InterpolationError> { - let x_raw = x.raw(); - if !x_raw.is_finite() { + pub fn evaluate( + &self, + abscissa: A, + ) -> Result, InterpolationError> { + if !abscissa.is_finite() { return Err(InterpolationError::NonFiniteAbscissa); } let (min, max) = self.range(); - if x_raw < min || x_raw > max { - return Err(InterpolationError::OutOfRange { x: x_raw, min, max }); + if abscissa.cmp_abscissa(min).is_lt() || abscissa.cmp_abscissa(max).is_gt() { + return Err(InterpolationError::OutOfRange { + requested_raw: abscissa.diagnostic_raw(), + min_raw: min.diagnostic_raw(), + max_raw: max.diagnostic_raw(), + }); } - let segment = self.segment_index(x_raw); + let segment = self.segment_index(abscissa); let s0 = &self.samples[segment]; let s1 = &self.samples[segment + 1]; - if x_raw == s0.x.raw() { + if abscissa.cmp_abscissa(s0.abscissa).is_eq() { return Ok(HermiteTableEvaluation { value: s0.value.clone(), derivative: s0.derivative.clone(), - x, + abscissa, }); } - if x_raw == s1.x.raw() { + if abscissa.cmp_abscissa(s1.abscissa).is_eq() { return Ok(HermiteTableEvaluation { value: s1.value.clone(), derivative: s1.derivative.clone(), - x, + abscissa, }); } - let dx = s1.x.raw() - s0.x.raw(); - let t = (x_raw - s0.x.raw()) / dx; - let basis = HermiteBasis::::new(t, dx); + let dx = s1.abscissa.delta_since(s0.abscissa); + let tau = abscissa.normalize_between(s0.abscissa, s1.abscissa)?; + let basis = HermiteBasis::::new(tau, dx); Ok(HermiteTableEvaluation { value: T::hermite_value( @@ -129,21 +148,21 @@ where s1.value.clone(), s1.derivative.clone(), ), - x, + abscissa, }) } - fn range(&self) -> (f64, f64) { + fn range(&self) -> (A, A) { ( - self.samples[0].x.raw(), - self.samples[self.samples.len() - 1].x.raw(), + self.samples[0].abscissa, + self.samples[self.samples.len() - 1].abscissa, ) } - fn segment_index(&self, x: f64) -> usize { + fn segment_index(&self, abscissa: A) -> usize { match self .samples - .binary_search_by(|sample| sample.x.raw().total_cmp(&x)) + .binary_search_by(|sample| sample.abscissa.cmp_abscissa(abscissa)) { Ok(index) => index.saturating_sub(1).min(self.samples.len() - 2), Err(index) => (index - 1).min(self.samples.len() - 2), @@ -173,18 +192,20 @@ fn validate_len(len: usize) -> Result<(), InterpolationError> { Ok(()) } -fn validate_sorted(xs: impl IntoIterator) -> Result<(), InterpolationError> { +fn validate_sorted( + abscissae: impl IntoIterator, +) -> Result<(), InterpolationError> { let mut previous = None; - for x in xs { + for abscissa in abscissae { if let Some(previous) = previous { - if x == previous { + if abscissa.cmp_abscissa(previous).is_eq() { return Err(InterpolationError::DuplicateAbscissa); } - if x < previous { + if abscissa.cmp_abscissa(previous).is_lt() { return Err(InterpolationError::UnsortedAbscissa); } } - previous = Some(x); + previous = Some(abscissa); } Ok(()) } diff --git a/src/interpolation/error.rs b/src/interpolation/error.rs index 40279de..9a5b144 100644 --- a/src/interpolation/error.rs +++ b/src/interpolation/error.rs @@ -24,12 +24,12 @@ pub enum InterpolationError { UnsortedAbscissa, /// The query abscissa is outside the interpolation table range. OutOfRange { - /// Query abscissa. - x: f64, - /// Minimum supported abscissa. - min: f64, - /// Maximum supported abscissa. - max: f64, + /// Query abscissa, raw diagnostic value only. + requested_raw: f64, + /// Minimum supported abscissa, raw diagnostic value only. + min_raw: f64, + /// Maximum supported abscissa, raw diagnostic value only. + max_raw: f64, }, } @@ -47,8 +47,15 @@ impl fmt::Display for InterpolationError { write!(f, "interpolation samples contain duplicate abscissae") } Self::UnsortedAbscissa => write!(f, "interpolation samples are not sorted"), - Self::OutOfRange { x, min, max } => { - write!(f, "interpolation query {x} is outside range [{min}, {max}]") + Self::OutOfRange { + requested_raw, + min_raw, + max_raw, + } => { + write!( + f, + "interpolation query {requested_raw} is outside diagnostic range [{min_raw}, {max_raw}]" + ) } } } diff --git a/src/interpolation/mod.rs b/src/interpolation/mod.rs index 4a2ab04..9dae7c5 100644 --- a/src/interpolation/mod.rs +++ b/src/interpolation/mod.rs @@ -8,7 +8,7 @@ pub mod cubic_hermite; pub mod error; pub mod traits; -pub use abscissa::InterpolationAbscissa; +pub use abscissa::{AbscissaDelta, InterpolationAbscissa}; pub use cubic_hermite::{ CubicHermiteTable, HermiteNode, HermiteTableEvaluation, ScalarCubicHermiteTable, ScalarHermiteNode, ScalarHermiteTableEvaluation, diff --git a/src/interpolation/traits.rs b/src/interpolation/traits.rs index e519edb..5c8fca4 100644 --- a/src/interpolation/traits.rs +++ b/src/interpolation/traits.rs @@ -1,6 +1,6 @@ //! Traits for interpolating complete typed values. -use super::abscissa::InterpolationAbscissa; +use super::abscissa::{AbscissaDelta, InterpolationAbscissa}; use crate::cartesian::{Position, Vector}; use crate::centers::ReferenceCenter; use crate::frames::ReferenceFrame; @@ -9,49 +9,92 @@ use qtty::{Quantity, Unit, UnitDiv, UnitMul}; /// Basis coefficients for one cubic Hermite segment. #[derive(Debug, Clone, Copy)] -pub struct HermiteBasis { - /// Value basis for the first endpoint. - pub h00: f64, - /// Derivative basis for the first endpoint, scaled by the abscissa delta. - pub h10_dx: X::Delta, - /// Value basis for the second endpoint. - pub h01: f64, - /// Derivative basis for the second endpoint, scaled by the abscissa delta. - pub h11_dx: X::Delta, - /// First derivative of `h00` with respect to normalized segment parameter. - pub dh00_over_dt: f64, - /// First derivative of `h10` with respect to normalized segment parameter. - pub dh10_over_dt: f64, - /// First derivative of `h01` with respect to normalized segment parameter. - pub dh01_over_dt: f64, - /// First derivative of `h11` with respect to normalized segment parameter. - pub dh11_over_dt: f64, - /// Full abscissa delta for the segment. - pub dx: X::Delta, - /// Raw inverse abscissa delta, for scalar-abscissa implementations. - pub inv_raw_dx: f64, +#[non_exhaustive] +pub struct HermiteBasis { + h00: f64, + h10_dx: A::Delta, + h01: f64, + h11_dx: A::Delta, + dh00_dt: f64, + dh10_dt: f64, + dh01_dt: f64, + dh11_dt: f64, + dx: A::Delta, } -impl HermiteBasis { - /// Constructs a basis from normalized `t` and raw segment width. +impl HermiteBasis { + /// Constructs a basis from normalized `tau` and typed segment width. #[inline] - pub(crate) fn new(t: f64, raw_dx: f64) -> Self { - let t2 = t * t; - let t3 = t2 * t; + pub(crate) fn new(tau: f64, dx: A::Delta) -> Self { + let tau2 = tau * tau; + let tau3 = tau2 * tau; Self { - h00: 2.0 * t3 - 3.0 * t2 + 1.0, - h10_dx: X::delta_from_raw((t3 - 2.0 * t2 + t) * raw_dx), - h01: -2.0 * t3 + 3.0 * t2, - h11_dx: X::delta_from_raw((t3 - t2) * raw_dx), - dh00_over_dt: 6.0 * t2 - 6.0 * t, - dh10_over_dt: 3.0 * t2 - 4.0 * t + 1.0, - dh01_over_dt: -6.0 * t2 + 6.0 * t, - dh11_over_dt: 3.0 * t2 - 2.0 * t, - dx: X::delta_from_raw(raw_dx), - inv_raw_dx: 1.0 / raw_dx, + h00: 2.0 * tau3 - 3.0 * tau2 + 1.0, + h10_dx: dx.scale(tau3 - 2.0 * tau2 + tau), + h01: -2.0 * tau3 + 3.0 * tau2, + h11_dx: dx.scale(tau3 - tau2), + dh00_dt: 6.0 * tau2 - 6.0 * tau, + dh10_dt: 3.0 * tau2 - 4.0 * tau + 1.0, + dh01_dt: -6.0 * tau2 + 6.0 * tau, + dh11_dt: 3.0 * tau2 - 2.0 * tau, + dx, } } + + /// Value basis for the first endpoint. + #[inline] + pub fn h00(&self) -> f64 { + self.h00 + } + + /// Derivative basis for the first endpoint, scaled by segment width. + #[inline] + pub fn h10_dx(&self) -> A::Delta { + self.h10_dx + } + + /// Value basis for the second endpoint. + #[inline] + pub fn h01(&self) -> f64 { + self.h01 + } + + /// Derivative basis for the second endpoint, scaled by segment width. + #[inline] + pub fn h11_dx(&self) -> A::Delta { + self.h11_dx + } + + /// Normalized derivative of `h00`. + #[inline] + pub fn dh00_dt(&self) -> f64 { + self.dh00_dt + } + + /// Normalized derivative of `h10`. + #[inline] + pub fn dh10_dt(&self) -> f64 { + self.dh10_dt + } + + /// Normalized derivative of `h01`. + #[inline] + pub fn dh01_dt(&self) -> f64 { + self.dh01_dt + } + + /// Normalized derivative of `h11`. + #[inline] + pub fn dh11_dt(&self) -> f64 { + self.dh11_dt + } + + /// Typed segment width. + #[inline] + pub fn dx(&self) -> A::Delta { + self.dx + } } /// A value that can be combined by cubic Hermite interpolation over abscissa @@ -59,13 +102,13 @@ impl HermiteBasis { /// /// Implementations operate on complete typed values, preserving all frame, /// center, and unit tags carried by the type. -pub trait HermiteInterpolable: Sized { +pub trait HermiteInterpolable: Sized { /// Derivative type with respect to the interpolation abscissa. type Derivative; /// Combines endpoint values and derivatives into an interpolated value. fn hermite_value( - basis: HermiteBasis, + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, @@ -74,7 +117,7 @@ pub trait HermiteInterpolable: Sized { /// Combines endpoint values and derivatives into an interpolated derivative. fn hermite_derivative( - basis: HermiteBasis, + basis: HermiteBasis, y0: Self, dy0: Self::Derivative, y1: Self, @@ -103,7 +146,7 @@ impl HermiteInterpolable for f64 { y1: Self, dy1: Self::Derivative, ) -> Self { - basis.h00 * y0 + basis.h10_dx * dy0 + basis.h01 * y1 + basis.h11_dx * dy1 + basis.h00() * y0 + basis.h10_dx() * dy0 + basis.h01() * y1 + basis.h11_dx() * dy1 } #[inline] @@ -114,49 +157,9 @@ impl HermiteInterpolable for f64 { y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - basis.dh00_over_dt * basis.inv_raw_dx * y0 - + basis.dh10_over_dt * dy0 - + basis.dh01_over_dt * basis.inv_raw_dx * y1 - + basis.dh11_over_dt * dy1 - } - - #[inline] - fn hermite_value_is_finite(&self) -> bool { - self.is_finite() - } - - #[inline] - fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.is_finite() - } -} - -impl HermiteInterpolable> for f64 { - type Derivative = f64; - - #[inline] - fn hermite_value( - basis: HermiteBasis>, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self { - basis.h00 * y0 + basis.h10_dx.value() * dy0 + basis.h01 * y1 + basis.h11_dx.value() * dy1 - } - - #[inline] - fn hermite_derivative( - basis: HermiteBasis>, - y0: Self, - dy0: Self::Derivative, - y1: Self, - dy1: Self::Derivative, - ) -> Self::Derivative { - basis.dh00_over_dt * basis.inv_raw_dx * y0 - + basis.dh10_over_dt * dy0 - + basis.dh01_over_dt * basis.inv_raw_dx * y1 - + basis.dh11_over_dt * dy1 + (basis.dh00_dt() * y0 + basis.dh01_dt() * y1) / AbscissaDelta::diagnostic_raw(basis.dx()) + + basis.dh10_dt() * dy0 + + basis.dh11_dt() * dy1 } #[inline] @@ -181,7 +184,7 @@ impl HermiteInterpolable for Quantity { y1: Self, dy1: Self::Derivative, ) -> Self { - y0 * basis.h00 + dy0 * basis.h10_dx + y1 * basis.h01 + dy1 * basis.h11_dx + y0 * basis.h00() + dy0 * basis.h10_dx() + y1 * basis.h01() + dy1 * basis.h11_dx() } #[inline] @@ -192,10 +195,10 @@ impl HermiteInterpolable for Quantity { y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - y0 * (basis.dh00_over_dt * basis.inv_raw_dx) - + dy0 * basis.dh10_over_dt - + y1 * (basis.dh01_over_dt * basis.inv_raw_dx) - + dy1 * basis.dh11_over_dt + y0 * (basis.dh00_dt() / AbscissaDelta::diagnostic_raw(basis.dx())) + + dy0 * basis.dh10_dt() + + y1 * (basis.dh01_dt() / AbscissaDelta::diagnostic_raw(basis.dx())) + + dy1 * basis.dh11_dt() } #[inline] @@ -225,7 +228,7 @@ where y1: Self, dy1: Self::Derivative, ) -> Self { - y0 * basis.h00 + (dy0 * basis.h10_dx) + y1 * basis.h01 + (dy1 * basis.h11_dx) + y0 * basis.h00() + (dy0 * basis.h10_dx()) + y1 * basis.h01() + (dy1 * basis.h11_dx()) } #[inline] @@ -236,10 +239,10 @@ where y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - (y0 / basis.dx) * basis.dh00_over_dt - + dy0 * basis.dh10_over_dt - + (y1 / basis.dx) * basis.dh01_over_dt - + dy1 * basis.dh11_over_dt + (y0 / basis.dx()) * basis.dh00_dt() + + dy0 * basis.dh10_dt() + + (y1 / basis.dx()) * basis.dh01_dt() + + dy1 * basis.dh11_dt() } #[inline] @@ -268,10 +271,10 @@ where y1: Self, dy1: Self::Derivative, ) -> Self { - y0.scale(basis.h00) - + dy0.scale(basis.h10_dx) - + y1.scale(basis.h01) - + dy1.scale(basis.h11_dx) + y0.scale(basis.h00()) + + dy0.scale(basis.h10_dx()) + + y1.scale(basis.h01()) + + dy1.scale(basis.h11_dx()) } #[inline] @@ -282,22 +285,20 @@ where y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - y0.scale(basis.dh00_over_dt * basis.inv_raw_dx) - + dy0.scale(basis.dh10_over_dt) - + y1.scale(basis.dh01_over_dt * basis.inv_raw_dx) - + dy1.scale(basis.dh11_over_dt) + y0.scale(basis.dh00_dt() / AbscissaDelta::diagnostic_raw(basis.dx())) + + dy0.scale(basis.dh10_dt()) + + y1.scale(basis.dh01_dt() / AbscissaDelta::diagnostic_raw(basis.dx())) + + dy1.scale(basis.dh11_dt()) } #[inline] fn hermite_value_is_finite(&self) -> bool { - self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + self.is_finite() } #[inline] fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.x().value().is_finite() - && derivative.y().value().is_finite() - && derivative.z().value().is_finite() + derivative.is_finite() } } @@ -318,7 +319,10 @@ where y1: Self, dy1: Self::Derivative, ) -> Self { - y0.scale(basis.h00) + (dy0 * basis.h10_dx) + y1.scale(basis.h01) + (dy1 * basis.h11_dx) + y0.scale(basis.h00()) + + (dy0 * basis.h10_dx()) + + y1.scale(basis.h01()) + + (dy1 * basis.h11_dx()) } #[inline] @@ -329,22 +333,20 @@ where y1: Self, dy1: Self::Derivative, ) -> Self::Derivative { - y0.div_quantity(basis.dx).scale(basis.dh00_over_dt) - + dy0.scale(basis.dh10_over_dt) - + y1.div_quantity(basis.dx).scale(basis.dh01_over_dt) - + dy1.scale(basis.dh11_over_dt) + y0.div_quantity(basis.dx()).scale(basis.dh00_dt()) + + dy0.scale(basis.dh10_dt()) + + y1.div_quantity(basis.dx()).scale(basis.dh01_dt()) + + dy1.scale(basis.dh11_dt()) } #[inline] fn hermite_value_is_finite(&self) -> bool { - self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + self.is_finite() } #[inline] fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.x().value().is_finite() - && derivative.y().value().is_finite() - && derivative.z().value().is_finite() + derivative.is_finite() } } @@ -365,7 +367,7 @@ where dy1: Self::Derivative, ) -> Self { let chord = y1 - y0; - y0 + chord.scale(basis.h01) + dy0.scale(basis.h10_dx) + dy1.scale(basis.h11_dx) + y0 + chord.scale(basis.h01()) + dy0.scale(basis.h10_dx()) + dy1.scale(basis.h11_dx()) } #[inline] @@ -377,21 +379,19 @@ where dy1: Self::Derivative, ) -> Self::Derivative { let chord = y1 - y0; - chord.scale(basis.dh01_over_dt * basis.inv_raw_dx) - + dy0.scale(basis.dh10_over_dt) - + dy1.scale(basis.dh11_over_dt) + chord.scale(basis.dh01_dt() / AbscissaDelta::diagnostic_raw(basis.dx())) + + dy0.scale(basis.dh10_dt()) + + dy1.scale(basis.dh11_dt()) } #[inline] fn hermite_value_is_finite(&self) -> bool { - self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + self.is_finite() } #[inline] fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.x().value().is_finite() - && derivative.y().value().is_finite() - && derivative.z().value().is_finite() + derivative.is_finite() } } @@ -414,7 +414,7 @@ where dy1: Self::Derivative, ) -> Self { let chord = y1 - y0; - y0 + chord.scale(basis.h01) + (dy0 * basis.h10_dx) + (dy1 * basis.h11_dx) + y0 + chord.scale(basis.h01()) + (dy0 * basis.h10_dx()) + (dy1 * basis.h11_dx()) } #[inline] @@ -426,20 +426,18 @@ where dy1: Self::Derivative, ) -> Self::Derivative { let chord = y1 - y0; - chord.div_quantity(basis.dx).scale(basis.dh01_over_dt) - + dy0.scale(basis.dh10_over_dt) - + dy1.scale(basis.dh11_over_dt) + chord.div_quantity(basis.dx()).scale(basis.dh01_dt()) + + dy0.scale(basis.dh10_dt()) + + dy1.scale(basis.dh11_dt()) } #[inline] fn hermite_value_is_finite(&self) -> bool { - self.x().value().is_finite() && self.y().value().is_finite() && self.z().value().is_finite() + self.is_finite() } #[inline] fn hermite_derivative_is_finite(derivative: &Self::Derivative) -> bool { - derivative.x().value().is_finite() - && derivative.y().value().is_finite() - && derivative.z().value().is_finite() + derivative.is_finite() } } diff --git a/src/lib.rs b/src/lib.rs index 3c45d0f..a0a1573 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -15,12 +15,12 @@ //! //! ### Reference Centers //! -//! A [`ReferenceCenter`](centers::ReferenceCenter) defines the origin point of a coordinate system. +//! A [`ReferenceCenter`] defines the origin point of a coordinate system. //! Some centers require runtime parameters (stored in `ReferenceCenter::Params`). //! //! ### Reference Frames //! -//! A [`ReferenceFrame`](frames::ReferenceFrame) defines the orientation of coordinate axes. +//! A [`ReferenceFrame`] defines the orientation of coordinate axes. //! //! ### Coordinate Types //! diff --git a/src/ops/rotation.rs b/src/ops/rotation.rs index c11de89..71c6548 100644 --- a/src/ops/rotation.rs +++ b/src/ops/rotation.rs @@ -313,7 +313,7 @@ impl Rotation3 { /// /// Rotates the vector's components and re-tags the result with frame `F2`. /// The unit `U` is unchanged. This is a zero-overhead operation: the matrix - /// multiplication uses the same path as [`apply_array`], with a frame tag + /// multiplication uses the same path as [`Self::apply_array`], with a frame tag /// change at compile time. /// /// # Example diff --git a/tests/interpolation.rs b/tests/interpolation.rs index a17d343..07868a1 100644 --- a/tests/interpolation.rs +++ b/tests/interpolation.rs @@ -18,12 +18,12 @@ type TestKmPerSecond = Per; type TestKilometerPosition = Position; type TestKilometerVelocity = Velocity; -fn cubic(x: f64) -> f64 { - x * x * x - 2.0 * x * x + x + 1.0 +fn cubic(abscissa: f64) -> f64 { + abscissa * abscissa * abscissa - 2.0 * abscissa * abscissa + abscissa + 1.0 } -fn cubic_derivative(x: f64) -> f64 { - 3.0 * x * x - 4.0 * x + 1.0 +fn cubic_derivative(abscissa: f64) -> f64 { + 3.0 * abscissa * abscissa - 4.0 * abscissa + 1.0 } fn cubic_position(t: f64) -> TestKilometerPosition { @@ -60,17 +60,17 @@ fn expect_scalar_table_error( fn scalar_cubic_polynomial_is_reproduced() { let table = ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: -1.0, + abscissa: -1.0, value: cubic(-1.0), derivative: cubic_derivative(-1.0), }, ScalarHermiteNode { - x: 0.5, + abscissa: 0.5, value: cubic(0.5), derivative: cubic_derivative(0.5), }, ScalarHermiteNode { - x: 2.0, + abscissa: 2.0, value: cubic(2.0), derivative: cubic_derivative(2.0), }, @@ -88,12 +88,12 @@ fn scalar_cubic_polynomial_is_reproduced() { fn exact_node_evaluation_returns_node_value_and_derivative() { let table = ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: 10.0, derivative: -3.0, }, ScalarHermiteNode { - x: 2.0, + abscissa: 2.0, value: 20.0, derivative: 4.0, }, @@ -109,12 +109,12 @@ fn exact_node_evaluation_returns_node_value_and_derivative() { fn linear_motion_with_constant_velocity_is_exact() { let table = CubicHermiteTable::::new(vec![ HermiteNode { - x: 0.0, + abscissa: 0.0, value: TestPosition::new(1.0, 2.0, 3.0), derivative: TestVelocity::new(0.5, -1.0, 2.0), }, HermiteNode { - x: 4.0, + abscissa: 4.0, value: TestPosition::new(3.0, -2.0, 11.0), derivative: TestVelocity::new(0.5, -1.0, 2.0), }, @@ -134,7 +134,7 @@ fn linear_motion_with_constant_velocity_is_exact() { fn typed_abscissa_table_accepts_position_over_seconds_with_velocity() { let table = CubicHermiteTable::::new(vec![ HermiteNode { - x: qtty::Second::new(0.0), + abscissa: qtty::Second::new(0.0), value: TestKilometerPosition::new(1.0, 2.0, 3.0), derivative: TestKilometerVelocity::new( Quantity::::new(0.5), @@ -143,7 +143,7 @@ fn typed_abscissa_table_accepts_position_over_seconds_with_velocity() { ), }, HermiteNode { - x: qtty::Second::new(4.0), + abscissa: qtty::Second::new(4.0), value: TestKilometerPosition::new(3.0, -2.0, 11.0), derivative: TestKilometerVelocity::new( Quantity::::new(0.5), @@ -167,17 +167,17 @@ fn typed_abscissa_table_accepts_position_over_seconds_with_velocity() { fn typed_abscissa_table_reproduces_cubic_position_over_seconds() { let table = CubicHermiteTable::::new(vec![ HermiteNode { - x: qtty::Second::new(-1.0), + abscissa: qtty::Second::new(-1.0), value: cubic_position(-1.0), derivative: cubic_velocity(-1.0), }, HermiteNode { - x: qtty::Second::new(0.5), + abscissa: qtty::Second::new(0.5), value: cubic_position(0.5), derivative: cubic_velocity(0.5), }, HermiteNode { - x: qtty::Second::new(2.0), + abscissa: qtty::Second::new(2.0), value: cubic_position(2.0), derivative: cubic_velocity(2.0), }, @@ -202,12 +202,12 @@ fn typed_abscissa_table_rejects_non_finite_abscissa() { let err = expect_table_error( CubicHermiteTable::::new(vec![ HermiteNode { - x: qtty::Second::new(f64::NAN), + abscissa: qtty::Second::new(f64::NAN), value: cubic_position(0.0), derivative: cubic_velocity(0.0), }, HermiteNode { - x: qtty::Second::new(1.0), + abscissa: qtty::Second::new(1.0), value: cubic_position(1.0), derivative: cubic_velocity(1.0), }, @@ -222,12 +222,12 @@ fn typed_abscissa_table_rejects_non_finite_position_component() { let err = expect_table_error( CubicHermiteTable::::new(vec![ HermiteNode { - x: qtty::Second::new(0.0), + abscissa: qtty::Second::new(0.0), value: TestKilometerPosition::new(f64::NAN, 0.0, 0.0), derivative: cubic_velocity(0.0), }, HermiteNode { - x: qtty::Second::new(1.0), + abscissa: qtty::Second::new(1.0), value: cubic_position(1.0), derivative: cubic_velocity(1.0), }, @@ -242,12 +242,12 @@ fn typed_abscissa_table_rejects_duplicate_seconds() { let err = expect_table_error( CubicHermiteTable::::new(vec![ HermiteNode { - x: qtty::Second::new(0.0), + abscissa: qtty::Second::new(0.0), value: cubic_position(0.0), derivative: cubic_velocity(0.0), }, HermiteNode { - x: qtty::Second::new(0.0), + abscissa: qtty::Second::new(0.0), value: cubic_position(1.0), derivative: cubic_velocity(1.0), }, @@ -281,17 +281,17 @@ fn vector_dimensional_mul_and_div_quantity_work() { fn non_uniform_sample_spacing_works() { let table = ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: cubic(0.0), derivative: cubic_derivative(0.0), }, ScalarHermiteNode { - x: 0.25, + abscissa: 0.25, value: cubic(0.25), derivative: cubic_derivative(0.25), }, ScalarHermiteNode { - x: 2.5, + abscissa: 2.5, value: cubic(2.5), derivative: cubic_derivative(2.5), }, @@ -307,12 +307,12 @@ fn non_uniform_sample_spacing_works() { fn out_of_range_queries_return_error() { let table = ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: 0.0, derivative: 1.0, }, ScalarHermiteNode { - x: 1.0, + abscissa: 1.0, value: 1.0, derivative: 1.0, }, @@ -322,9 +322,9 @@ fn out_of_range_queries_return_error() { assert_eq!( table.evaluate(2.0), Err(InterpolationError::OutOfRange { - x: 2.0, - min: 0.0, - max: 1.0 + requested_raw: 2.0, + min_raw: 0.0, + max_raw: 1.0 }) ); } @@ -334,12 +334,12 @@ fn duplicate_abscissae_are_rejected() { assert_eq!( expect_scalar_table_error(ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: 0.0, derivative: 0.0, }, ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: 1.0, derivative: 1.0, }, @@ -353,12 +353,12 @@ fn unsorted_abscissae_are_rejected() { assert_eq!( expect_scalar_table_error(ScalarCubicHermiteTable::new(vec![ ScalarHermiteNode { - x: 1.0, + abscissa: 1.0, value: 1.0, derivative: 1.0, }, ScalarHermiteNode { - x: 0.0, + abscissa: 0.0, value: 0.0, derivative: 0.0, }, From 346ea999a182311c538c54903cdc164967ea8540 Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 21:56:26 +0200 Subject: [PATCH 5/6] chore(release): bump version to 0.8.0 for affn and affn-derive; update changelog with new interpolation module and features --- CHANGELOG.md | 15 +++++++++ Cargo.lock | 75 +++++++++++++++++++++--------------------- Cargo.toml | 2 +- affn-derive/Cargo.toml | 2 +- 4 files changed, 54 insertions(+), 40 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a545527..da83e59 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,21 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.8.0 - 2026-06-16] + +### Added +- New `affn::interpolation` module providing domain-agnostic interpolation primitives for typed affine geometry. +- `CubicHermiteTable` for cubic Hermite interpolation over any `HermiteInterpolable` type (including `Position`, `Vector`, `Direction`, and `XYZ`). +- `ScalarCubicHermiteTable` for efficient scalar-only interpolation. +- `HermiteNode` and `HermiteTableEvaluation` for representing nodes and evaluation results in Hermite interpolation. +- `InterpolationAbscissa` and `AbscissaDelta` for typed abscissa (e.g., time) with unit support via `qtty`. +- `InterpolationError` for comprehensive error handling during interpolation operations (out-of-bounds, insufficient nodes, etc.). +- Trait-based interpolation support via `HermiteInterpolable` (types that can be interpolated) and `HermiteBasis` (basis function implementations). +- `HermiteTableEvaluation::try_value()` for fallible evaluation with proper error propagation. + +### Changed +- Interpolation module is now domain-agnostic and focuses on scalar abscissae with strongly-typed values, enabling downstream crates to implement astronomy-specific epoch handling. + ## [0.7.3 - 2026-05-25] ### Changed diff --git a/Cargo.lock b/Cargo.lock index 9bcdc67..3ab694a 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4,7 +4,7 @@ version = 4 [[package]] name = "affn" -version = "0.7.3" +version = "0.8.0" dependencies = [ "affn-derive", "criterion", @@ -72,9 +72,9 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" [[package]] name = "cc" -version = "1.2.62" +version = "1.2.64" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a1dce859f0832a7d088c4f1119888ab94ef4b5d6795d1ce05afb7fe159d79f98" +checksum = "dad887fd958be91b5098c0248def011f4523ab786cd411be668777e55063501f" dependencies = [ "find-msvc-tools", "shlex", @@ -268,13 +268,12 @@ checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682" [[package]] name = "js-sys" -version = "0.3.99" +version = "0.3.102" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "142bc4740e452c1e57ade0cbc129f139c9093e354346f0872ef985f4f5cf5f11" +checksum = "03d04c30968dffe80775bd4d7fb676131cd04a1fb46d2686dbffbaec2d9dfd31" dependencies = [ "cfg-if", "futures-util", - "once_cell", "wasm-bindgen", ] @@ -292,9 +291,9 @@ checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981" [[package]] name = "memchr" -version = "2.8.0" +version = "2.8.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f8ca58f447f06ed17d5fc4043ce1b10dd205e060fb3ce5b979b8ed8e59ff3f79" +checksum = "88904434abc2901f197fe8cc55f0445e7ded921dba5911dad2e2b39b48e663c4" [[package]] name = "num-traits" @@ -372,9 +371,9 @@ dependencies = [ [[package]] name = "qtty" -version = "0.8.3" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "43953c1990c78d8a98a3018f9b51071f12b478c4c519f2eaf7fff954b315d08e" +checksum = "311988508597047d023edfb10cf4299b22a317be2b2c95745aaa4d489d9a9265" dependencies = [ "qtty-core", "qtty-derive", @@ -382,9 +381,9 @@ dependencies = [ [[package]] name = "qtty-core" -version = "0.8.3" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3abd20de3f1a4b6a4979f42f0032900b0ec030df95eaaa2b0294a2918107ffa6" +checksum = "dc07e03ed9d81d4b4494b661635dc87f0a1a38cdce3e51136c39af6887716b02" dependencies = [ "libm", "qtty-derive", @@ -394,9 +393,9 @@ dependencies = [ [[package]] name = "qtty-derive" -version = "0.8.3" +version = "0.8.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2ddc68ce74f35034435b61807f624f09001cbf68d3a0d4b1288d9cb68a40e1c2" +checksum = "7ee2574d42b017be334315c8cf4fdb241d728b093c5c7aba9bb1e7c5a3a708e8" dependencies = [ "proc-macro2", "quote", @@ -434,9 +433,9 @@ dependencies = [ [[package]] name = "regex" -version = "1.12.3" +version = "1.12.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e10754a14b9137dd7b1e3e5b0493cc9171fdd105e0ab477f51b72e7f3ac0e276" +checksum = "f1292b7759ae1cb9ec195452d1390a074f0cd8541ab7a5a8c31cd6db45d4a6ba" dependencies = [ "aho-corasick", "memchr", @@ -457,9 +456,9 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.8.10" +version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dc897dd8d9e8bd1ed8cdad82b5966c3e0ecae09fb1907d58efaa013543185d0a" +checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" [[package]] name = "rustversion" @@ -521,9 +520,9 @@ dependencies = [ [[package]] name = "shlex" -version = "1.3.0" +version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" +checksum = "f8fadd59c855ef2080decdef8ff161eb6661b86933c9d82e5ba29dc602a55aba" [[package]] name = "slab" @@ -533,9 +532,9 @@ checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" [[package]] name = "syn" -version = "2.0.117" +version = "2.0.118" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e665b8803e7b1d2a727f4023456bbbbe74da67099c585258af0ad9c5013b9b99" +checksum = "1b9ae57f904213ebb649ce6895b8a66c66f0203b9319718f69a5612a065b1422" dependencies = [ "proc-macro2", "quote", @@ -554,9 +553,9 @@ dependencies = [ [[package]] name = "typenum" -version = "1.20.0" +version = "1.20.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "40ce102ab67701b8526c123c1bab5cbe42d7040ccfd0f64af1a385808d2f43de" +checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20" [[package]] name = "unicode-ident" @@ -576,9 +575,9 @@ dependencies = [ [[package]] name = "wasm-bindgen" -version = "0.2.122" +version = "0.2.125" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3ed04576f974d2b2fba0f38c51dbc5518011e38c36bf1143164be765528fd409" +checksum = "8ddb3f79143bced6de84270411622a2699cee572fc0875aeaf1e7867cf9fca1a" dependencies = [ "cfg-if", "once_cell", @@ -589,9 +588,9 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro" -version = "0.2.122" +version = "0.2.125" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "916151b09da36bd82f6615cbf3a419e2f0ba23a03c6160e8e92eb6bd4aa1dec6" +checksum = "4e21a184b13fb19e157296e2c46056aec9092264fab83e4ba59e68c61b323c3d" dependencies = [ "quote", "wasm-bindgen-macro-support", @@ -599,9 +598,9 @@ dependencies = [ [[package]] name = "wasm-bindgen-macro-support" -version = "0.2.122" +version = "0.2.125" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "299047362ccbfce148b67ab7e73349f77748e00c8296f9542adfad2ad82c5c5e" +checksum = "fecefd9c35bd935a20fc3fc344b5f29138961e4f47fb03297d88f2587afb5ebd" dependencies = [ "bumpalo", "proc-macro2", @@ -612,18 +611,18 @@ dependencies = [ [[package]] name = "wasm-bindgen-shared" -version = "0.2.122" +version = "0.2.125" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9a929b2c61f11ba3e9bc35b50c1f25cb38e0e892c0c231ae2b8cf78d5dad4437" +checksum = "23939e44bb9a5d7576fa2b563dc2e136628f1224e88a8deed09e04858b77871f" dependencies = [ "unicode-ident", ] [[package]] name = "web-sys" -version = "0.3.99" +version = "0.3.102" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6d621441cfc37b84979402712047321980c178f299193a3589d05b99e8763436" +checksum = "a6430a72df5eb332242960fe84b3002a241163998241eb596d4f739b9757061d" dependencies = [ "js-sys", "wasm-bindgen", @@ -677,18 +676,18 @@ dependencies = [ [[package]] name = "zerocopy" -version = "0.8.48" +version = "0.8.52" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "eed437bf9d6692032087e337407a86f04cd8d6a16a37199ed57949d415bd68e9" +checksum = "ce1022995ff5ff5d841ad7d994facc23098cd40152f2c1d11cd607c6f530653f" dependencies = [ "zerocopy-derive", ] [[package]] name = "zerocopy-derive" -version = "0.8.48" +version = "0.8.52" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "70e3cd084b1788766f53af483dd21f93881ff30d7320490ec3ef7526d203bad4" +checksum = "1ae7f38b72ec2a254e2b87ef277cf2cd4fb97cbebf944faa6f33354da0867930" dependencies = [ "proc-macro2", "quote", diff --git a/Cargo.toml b/Cargo.toml index 476efa7..aaea099 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -2,7 +2,7 @@ [package] name = "affn" -version = "0.7.3" +version = "0.8.0" edition = "2021" authors = ["VPRamon "] description = "Affine geometry primitives: strongly-typed coordinate systems, reference frames, and centers for scientific computing." diff --git a/affn-derive/Cargo.toml b/affn-derive/Cargo.toml index 64dbfd4..2b381af 100644 --- a/affn-derive/Cargo.toml +++ b/affn-derive/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "affn-derive" -version = "0.7.3" +version = "0.8.0" edition = "2021" authors = ["VPRamon "] description = "Derive macros for affn: ReferenceFrame and ReferenceCenter" From 473bf662d40f189af01063d08c604a4572a7f6e4 Mon Sep 17 00:00:00 2001 From: VPRamon Date: Tue, 16 Jun 2026 22:00:42 +0200 Subject: [PATCH 6/6] update version --- Cargo.lock | 2 +- Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 3ab694a..5fc2dd0 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -15,7 +15,7 @@ dependencies = [ [[package]] name = "affn-derive" -version = "0.7.3" +version = "0.8.0" dependencies = [ "proc-macro2", "quote", diff --git a/Cargo.toml b/Cargo.toml index aaea099..1a5054f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,7 +21,7 @@ serde = ["dep:serde", "qtty/serde"] astro = ["qtty/astro"] [dependencies] -affn-derive = { version = "0.7", path = "affn-derive" } +affn-derive = { version = "0.8", path = "affn-derive" } qtty = "0.8" serde = { version = "1.0", default-features = false, features = ["derive"], optional = true }