diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b2a5ace..2485fa0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -24,6 +24,8 @@ jobs: - run: cargo clippy --all-targets --all-features -- -D warnings - run: cargo check --no-default-features - run: cargo check --no-default-features --features alloc + - run: cargo check --no-default-features --features serde + - run: cargo check --no-default-features --features alloc,serde test: name: Test diff --git a/CHANGELOG.md b/CHANGELOG.md index 39bbf0a..a65a39c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,31 +7,73 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] +### Added + +- `MeanAnomaly`, `TrueAnomaly`, `EccentricAnomaly`, `HyperbolicAnomaly`, + `ParabolicAnomaly` typed anomaly newtypes; all anomaly solver functions now + accept and return these typed values instead of raw `Radians` or `f64`. +- `AnomalyOptions::try_new(max_iter, tol)` fallible constructor; fields are now + private; `max_iter()` and `tol()` accessors added. +- `AnomalyError::InvalidOptions` variant for invalid solver configuration. +- `ConicRegime` and `EccentricityError` moved to `eccentricity` module; + `elements::ConicRegime` is a re-export for backwards compatibility. +- `Eccentricity::try_new`, `new_elliptic`, `new_hyperbolic`, `parabolic`, + `classify` constructors. +- `parabolic_from_mean`, `true_from_parabolic`, `parabolic_from_true` parabolic + anomaly helpers (replaces the old `kepler_parabolic` raw-f64 function). +- `KeplerianElements::try_to_cartesian` fallible Cartesian conversion (replaces + infallible `to_cartesian`). +- `ConversionError::IncoherentRegime` — returned when the semi-major axis sign + is inconsistent with the eccentricity regime. +- `a`/`e` coherence check in `KeplerianElements::new`: elliptic requires `a > 0`; + hyperbolic requires `a < 0`. +- `TransferError` enum for fallible transfer helpers. +- `try_orbital_period`, `try_hohmann_delta_v`, `try_vis_viva_speed`, + `try_escape_speed` fallible variants of transfer functions. +- `prelude` module re-exporting the most commonly used public items. +- `#![cfg_attr(not(feature = "std"), no_std)]` — crate is now no-std capable + with `alloc` feature for heap-using APIs. +- `Clone` derive on `LambertError`; `PartialEq` and serde derives on + `LambertDiagnostics`, `LambertBranch`, `NRevBranch`, `TypedLambertSolution`. +- CI: additional checks for `--features serde` and `--features alloc,serde` + no-default-features combinations. + +### Changed + +- `KeplerProblem.mu` field is now private; use `mu()` accessor. +- `TransferCandidate.total_dv` renamed to `endpoint_speed_sum`. +- `kepler_parabolic` renamed to `parabolic_from_mean`; accepts and returns + `ParabolicAnomaly` instead of raw `f64`. + +### Removed + +- Infallible `KeplerianElements::to_cartesian`; replaced by + `try_to_cartesian`. +- `AnomalyOptions` struct literal construction (fields privatised); use + `AnomalyOptions::try_new`. + ## [0.1.0] - 2026-05-22 ### Added -- `Eccentricity` newtype with `new_unchecked`, `new_elliptic`, `new_hyperbolic`, - and `classify` constructors; compile-time validation enforced. -- `AnomalyOptions` control struct (max iterations, tolerance) and - `AnomalyError` (non-convergence) for all iterative solvers. +- `Eccentricity` newtype with `new_unchecked` constructor; compile-time + validation enforced. +- `AnomalyOptions` control struct and `AnomalyError` for all iterative solvers. - Elliptic anomaly solvers: `eccentric_from_mean`, `true_from_eccentric`, `eccentric_from_true`, `mean_from_eccentric`. -- Parabolic anomaly solver: `kepler_parabolic` (Barker's equation). - Hyperbolic anomaly solvers: `hyperbolic_from_mean`, `true_from_hyperbolic`. - `KeplerianElements` — six classical elements typed over an `affn` frame; - conversions `to_cartesian` / `from_cartesian`. + conversions `from_cartesian`. - `KeplerProblem` — two-body IVP solver; `new` + `propagate`. - Lambert boundary-value solver (`lambert`, `lambert_n_rev`) with multi-rev support. - Transfer invariants: `specific_orbital_energy`, `specific_angular_momentum`, `vis_viva_speed`, `orbital_period`, `escape_speed`, `hohmann_delta_v`. -- `OrbitalState` and `StateDerivative` Cartesian state types. +- `CartesianState` Cartesian state type. - `KeplerError` crate-level error family. - `alloc`-gated `search` module: Lambert transfer-search grids. -- Optional `serde` feature for all public data types. -- CI workflow with fmt, Clippy (default + all features), no-std/alloc checks, - tests (default + all features), and doc-test jobs. +- Optional `serde` feature for public data types. +- CI workflow with fmt, Clippy, no-std/alloc checks, tests, and doc-test jobs. - Audit, deny, coverage (llvm-cov ≥ 90 % line rate), and Miri jobs. - Publish workflow triggered on `v*.*.*` tags. diff --git a/Cargo.lock b/Cargo.lock index 75b9437..22bc1b9 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4,9 +4,9 @@ version = 4 [[package]] name = "affn" -version = "0.7.2" +version = "0.7.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3e64f45b201739dea0c81b5fd4314f9cbc35f3d6ed5163e9bfec17be47178d89" +checksum = "c87b9ab9380ca5e99f103892edf0d608fd693aee15445e76d9669500ae6560a4" dependencies = [ "affn-derive", "qtty", @@ -15,9 +15,9 @@ dependencies = [ [[package]] name = "affn-derive" -version = "0.7.2" +version = "0.7.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8e59162518fab751e2af6b0e7ec7394f04cef094bd667cc5b09b3f2ecd6f0b5a" +checksum = "0b0de78baee7e6ec28ef4f556d50a710a2fde9dcb89da1bd80e4cc4a36643984" dependencies = [ "proc-macro2", "quote", @@ -214,9 +214,9 @@ checksum = "32a66949e030da00e8c7d4434b251670a91556f4144941d37452769c25d58a53" [[package]] name = "log" -version = "0.4.29" +version = "0.4.30" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5e5032e24019045c762d3c0f28f5b6b8bbf38563a65908389bf7978758920897" +checksum = "616ec5685824bcc94416c6d4a7a446eea774a31efd7062c8480ba6fd06d7a6e5" [[package]] name = "memchr" @@ -288,9 +288,9 @@ dependencies = [ [[package]] name = "qtty" -version = "0.8.2" +version = "0.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9e37724c555d5c157a10522b8a09f4ab1f0fa4eb10d3107b1270b92ad56ba5aa" +checksum = "43953c1990c78d8a98a3018f9b51071f12b478c4c519f2eaf7fff954b315d08e" dependencies = [ "qtty-core", "qtty-derive", @@ -298,9 +298,9 @@ dependencies = [ [[package]] name = "qtty-core" -version = "0.8.2" +version = "0.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "40d5237d4b90129f128f7d115ab7be10fc2dd6957a52900d75400253bd7c724b" +checksum = "3abd20de3f1a4b6a4979f42f0032900b0ec030df95eaaa2b0294a2918107ffa6" dependencies = [ "libm", "qtty-derive", @@ -310,9 +310,9 @@ dependencies = [ [[package]] name = "qtty-derive" -version = "0.8.2" +version = "0.8.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6b53a35f7c86504e9f8db9273ec58bf541e35a34ef424cde04b4fbc8bb73bf63" +checksum = "2ddc68ce74f35034435b61807f624f09001cbf68d3a0d4b1288d9cb68a40e1c2" dependencies = [ "proc-macro2", "quote", diff --git a/examples/anomaly_solve.rs b/examples/anomaly_solve.rs index 1ab7701..b22ca2e 100644 --- a/examples/anomaly_solve.rs +++ b/examples/anomaly_solve.rs @@ -1,14 +1,14 @@ //! Example binary for elliptic Kepler solves. #![allow(clippy::print_stdout)] -use keplerian::anomaly::{kepler_elliptic, AnomalyOptions}; +use keplerian::anomaly::{kepler_elliptic, AnomalyOptions, MeanAnomaly}; use keplerian::Eccentricity; -use qtty::angular::Radians; fn main() { let ecc = Eccentricity::new(0.2).unwrap(); for m in [0.0, 0.5, 1.0, 1.5] { - let e = kepler_elliptic(Radians::new(m), ecc, AnomalyOptions::default()).unwrap(); + let e = + kepler_elliptic(MeanAnomaly::from_value(m), ecc, AnomalyOptions::default()).unwrap(); println!("M = {m:.3} rad -> E = {:.6} rad", e.value()); } } diff --git a/examples/state_elements_roundtrip.rs b/examples/state_elements_roundtrip.rs index e4800b3..4a79216 100644 --- a/examples/state_elements_roundtrip.rs +++ b/examples/state_elements_roundtrip.rs @@ -32,7 +32,7 @@ fn main() { Velocity::::new(0.0, 7.2, 1.0), ); let elements = KeplerianElements::from_cartesian(&state, mu).unwrap(); - let back = elements.to_cartesian::
(mu); + let back = elements.try_to_cartesian::
(mu).unwrap(); assert!((back.position().x().value() - state.position().x().value()).abs() < 1e-8); println!( "a = {:.3} km, e = {:.6}", diff --git a/src/anomaly.rs b/src/anomaly.rs index 1b0b9f3..9de7cf0 100644 --- a/src/anomaly.rs +++ b/src/anomaly.rs @@ -1,7 +1,7 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon -//! Typed anomaly solvers and anomaly conversions. +//! Typed anomaly newtypes, solvers, and anomaly conversions. //! //! ## Scientific scope //! This module implements the standard elliptic, parabolic, and hyperbolic @@ -9,7 +9,7 @@ //! target ordinary astrodynamics workloads and do not model perturbations. //! //! ## Technical scope -//! Public APIs use [`qtty::angular::Radians`] for angular quantities and +//! Public APIs use the five anomaly newtypes defined here and //! [`crate::eccentricity::Eccentricity`] for the conic parameter. Private //! numeric kernels remain raw `f64`. //! @@ -24,13 +24,244 @@ use qtty::angular::Radians; use crate::eccentricity::Eccentricity; +// ── Anomaly newtypes ───────────────────────────────────────────────────────── + +/// Mean anomaly for elliptic motion (angular, wraps a [`Radians`]). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::MeanAnomaly; +/// use qtty::angular::Radians; +/// let m = MeanAnomaly::new(Radians::new(1.0)); +/// assert_eq!(m.value(), 1.0); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct MeanAnomaly(Radians); + +impl MeanAnomaly { + /// Creates a mean anomaly from a typed [`Radians`]. + #[must_use] + pub const fn new(r: Radians) -> Self { + Self(r) + } + + /// Creates a mean anomaly from a raw radian value. + #[must_use] + pub fn from_value(v: f64) -> Self { + Self(Radians::new(v)) + } + + /// Returns the typed radians value. + #[must_use] + pub const fn radians(self) -> Radians { + self.0 + } + + /// Returns the raw radian value. + #[must_use] + pub fn value(self) -> f64 { + self.0.value() + } +} + +/// True anomaly (angular, wraps a [`Radians`]). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::TrueAnomaly; +/// use qtty::angular::Radians; +/// let nu = TrueAnomaly::new(Radians::new(0.5)); +/// assert_eq!(nu.value(), 0.5); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct TrueAnomaly(Radians); + +impl TrueAnomaly { + /// Creates a true anomaly from a typed [`Radians`]. + #[must_use] + pub const fn new(r: Radians) -> Self { + Self(r) + } + + /// Creates a true anomaly from a raw radian value. + #[must_use] + pub fn from_value(v: f64) -> Self { + Self(Radians::new(v)) + } + + /// Returns the typed radians value. + #[must_use] + pub const fn radians(self) -> Radians { + self.0 + } + + /// Returns the raw radian value. + #[must_use] + pub fn value(self) -> f64 { + self.0.value() + } +} + +/// Eccentric anomaly for elliptic motion (angular, wraps a [`Radians`]). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::EccentricAnomaly; +/// use qtty::angular::Radians; +/// let ea = EccentricAnomaly::new(Radians::new(1.2)); +/// assert_eq!(ea.value(), 1.2); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct EccentricAnomaly(Radians); + +impl EccentricAnomaly { + /// Creates an eccentric anomaly from a typed [`Radians`]. + #[must_use] + pub const fn new(r: Radians) -> Self { + Self(r) + } + + /// Creates an eccentric anomaly from a raw radian value. + #[must_use] + pub fn from_value(v: f64) -> Self { + Self(Radians::new(v)) + } + + /// Returns the typed radians value. + #[must_use] + pub const fn radians(self) -> Radians { + self.0 + } + + /// Returns the raw radian value. + #[must_use] + pub fn value(self) -> f64 { + self.0.value() + } +} + +/// Hyperbolic anomaly (dimensionless Barker-like parameter; **not** an angle). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::HyperbolicAnomaly; +/// let f = HyperbolicAnomaly::new(0.5); +/// assert_eq!(f.value(), 0.5); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct HyperbolicAnomaly(f64); + +impl HyperbolicAnomaly { + /// Creates a hyperbolic anomaly from a raw dimensionless value. + #[must_use] + pub const fn new(v: f64) -> Self { + Self(v) + } + + /// Returns the raw dimensionless value. + #[must_use] + pub const fn value(self) -> f64 { + self.0 + } +} + +/// Parabolic anomaly — Barker's parameter `D = tan(ν/2)` (dimensionless). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::ParabolicAnomaly; +/// let d = ParabolicAnomaly::new(0.5); +/// assert_eq!(d.value(), 0.5); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq, PartialOrd)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub struct ParabolicAnomaly(f64); + +impl ParabolicAnomaly { + /// Creates a parabolic anomaly from a raw dimensionless value. + #[must_use] + pub const fn new(v: f64) -> Self { + Self(v) + } + + /// Returns the raw dimensionless value. + #[must_use] + pub const fn value(self) -> f64 { + self.0 + } +} + +// ── AnomalyOptions ─────────────────────────────────────────────────────────── + /// Options controlling iterative Kepler-equation solves. +/// +/// Construct via [`AnomalyOptions::try_new`] or use [`Default`] (64 iterations, +/// tolerance `1e-12`). +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::AnomalyOptions; +/// let opts = AnomalyOptions::try_new(50, 1e-14).unwrap(); +/// assert_eq!(opts.max_iter(), 50); +/// ``` #[derive(Debug, Clone, Copy, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct AnomalyOptions { - /// Maximum Newton or bisection iterations. - pub max_iter: u32, - /// Absolute residual tolerance in raw radians. - pub tol: f64, + max_iter: u32, + tol: f64, +} + +impl AnomalyOptions { + /// Creates solver options. + /// + /// # Errors + /// + /// Returns [`AnomalyError::InvalidOptions`] if `max_iter == 0` or `tol` is + /// not a strictly positive finite value. + /// + /// # Examples + /// + /// ``` + /// use keplerian::anomaly::{AnomalyOptions, AnomalyError}; + /// assert!(AnomalyOptions::try_new(64, 1e-12).is_ok()); + /// assert!(matches!(AnomalyOptions::try_new(0, 1e-12), Err(AnomalyError::InvalidOptions(_)))); + /// assert!(matches!(AnomalyOptions::try_new(10, 0.0), Err(AnomalyError::InvalidOptions(_)))); + /// ``` + pub fn try_new(max_iter: u32, tol: f64) -> Result { + if max_iter == 0 { + return Err(AnomalyError::InvalidOptions( + "max_iter must be greater than zero", + )); + } + if !tol.is_finite() || tol <= 0.0 { + return Err(AnomalyError::InvalidOptions( + "tol must be a strictly positive finite value", + )); + } + Ok(Self { max_iter, tol }) + } + + /// Returns the maximum iteration count. + #[must_use] + pub const fn max_iter(self) -> u32 { + self.max_iter + } + + /// Returns the convergence tolerance. + #[must_use] + pub const fn tol(self) -> f64 { + self.tol + } } impl Default for AnomalyOptions { @@ -42,6 +273,8 @@ impl Default for AnomalyOptions { } } +// ── AnomalyError ───────────────────────────────────────────────────────────── + /// Errors returned by anomaly solvers. #[derive(Debug, Clone, Copy, PartialEq, thiserror::Error)] pub enum AnomalyError { @@ -61,9 +294,14 @@ pub enum AnomalyError { /// Mean anomaly is not finite. #[error("invalid mean anomaly {0}")] InvalidMeanAnomaly(f64), + /// Invalid solver options. + #[error("invalid AnomalyOptions: {0}")] + InvalidOptions(&'static str), } -/// Solves the elliptic Kepler equation `M = E - e sin(E)`. +// ── Elliptic solvers ───────────────────────────────────────────────────────── + +/// Solves the elliptic Kepler equation `M = E − e sin(E)`. /// /// # Errors /// @@ -73,17 +311,20 @@ pub enum AnomalyError { /// # Examples /// /// ``` -/// use keplerian::anomaly::{kepler_elliptic, AnomalyOptions}; +/// use keplerian::anomaly::{kepler_elliptic, AnomalyOptions, MeanAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; -/// let e = kepler_elliptic(Radians::new(0.5), Eccentricity::new_unchecked(0.1), AnomalyOptions::default()).unwrap(); +/// let e = kepler_elliptic( +/// MeanAnomaly::from_value(0.5), +/// Eccentricity::new_unchecked(0.1), +/// AnomalyOptions::default(), +/// ).unwrap(); /// assert!((e.value() - 0.5524799869).abs() < 1e-9); /// ``` pub fn kepler_elliptic( - mean_anomaly: Radians, + mean_anomaly: MeanAnomaly, ecc: Eccentricity, opts: AnomalyOptions, -) -> Result { +) -> Result { let ecc_value = ecc.value(); let mean_value = mean_anomaly.value(); if !(0.0..1.0).contains(&ecc_value) || !ecc_value.is_finite() { @@ -93,13 +334,13 @@ pub fn kepler_elliptic( return Err(AnomalyError::InvalidMeanAnomaly(mean_value)); } if ecc_value == 0.0 { - return Ok(mean_anomaly); + return Ok(EccentricAnomaly::new(Radians::new(mean_value))); } let mut e_anom = mean_value + ecc_value * mean_value.sin(); for i in 0..opts.max_iter { let mut residual = elliptic_residual(e_anom, ecc_value, mean_value); if residual.abs() <= opts.tol { - return Ok(Radians::new(e_anom)); + return Ok(EccentricAnomaly::new(Radians::new(e_anom))); } let fp = 1.0 - ecc_value * e_anom.cos(); e_anom -= residual / fp; @@ -109,14 +350,14 @@ pub fn kepler_elliptic( if i + 1 == opts.max_iter { residual = elliptic_residual(e_anom, ecc_value, mean_value); if residual.abs() <= opts.tol { - return Ok(Radians::new(e_anom)); + return Ok(EccentricAnomaly::new(Radians::new(e_anom))); } } } - elliptic_bisection(mean_value, ecc_value, opts).map(Radians::new) + elliptic_bisection(mean_value, ecc_value, opts).map(|v| EccentricAnomaly::new(Radians::new(v))) } -/// Solves the hyperbolic Kepler equation `M = e sinh(F) - F`. +/// Solves the hyperbolic Kepler equation `M = e sinh(F) − F`. /// /// # Errors /// @@ -126,18 +367,17 @@ pub fn kepler_elliptic( /// # Examples /// /// ``` -/// use keplerian::anomaly::{kepler_hyperbolic, mean_from_hyperbolic, AnomalyOptions}; +/// use keplerian::anomaly::{kepler_hyperbolic, mean_from_hyperbolic, AnomalyOptions, MeanAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(1.5); -/// let f = kepler_hyperbolic(Radians::new(1.0), ecc, AnomalyOptions::default()).unwrap(); -/// assert!((mean_from_hyperbolic(f.value(), ecc).value() - 1.0).abs() < 1e-12); +/// let f = kepler_hyperbolic(MeanAnomaly::from_value(1.0), ecc, AnomalyOptions::default()).unwrap(); +/// assert!((mean_from_hyperbolic(f, ecc).value() - 1.0).abs() < 1e-12); /// ``` pub fn kepler_hyperbolic( - mean_anomaly: Radians, + mean_anomaly: MeanAnomaly, ecc: Eccentricity, opts: AnomalyOptions, -) -> Result { +) -> Result { let ecc_value = ecc.value(); let mean_value = mean_anomaly.value(); if ecc_value <= 1.0 || !ecc_value.is_finite() { @@ -147,7 +387,7 @@ pub fn kepler_hyperbolic( return Err(AnomalyError::InvalidMeanAnomaly(mean_value)); } if mean_value == 0.0 { - return Ok(Radians::new(0.0)); + return Ok(HyperbolicAnomaly::new(0.0)); } let abs_mean = mean_value.abs(); let mut f_anom = if abs_mean > 50.0 * ecc_value { @@ -158,7 +398,7 @@ pub fn kepler_hyperbolic( for _ in 0..opts.max_iter { let residual = hyperbolic_residual(f_anom, ecc_value, mean_value); if residual.abs() <= opts.tol { - return Ok(Radians::new(f_anom)); + return Ok(HyperbolicAnomaly::new(f_anom)); } let fp = ecc_value * f_anom.cosh() - 1.0; f_anom -= residual / fp; @@ -166,62 +406,66 @@ pub fn kepler_hyperbolic( break; } } - hyperbolic_bisection(mean_value, ecc_value, opts).map(Radians::new) + hyperbolic_bisection(mean_value, ecc_value, opts).map(HyperbolicAnomaly::new) } -/// Solves Barker's parabolic equation `M = D + D³/3` analytically. +/// Solves Barker's parabolic equation `M_D = D + D³/3` analytically for D. +/// +/// Both input and output are [`ParabolicAnomaly`] (Barker's dimensionless +/// parameter). The input represents the mean-like quantity +/// `M_D = √(μ/2p³)·(t − t_p)`, and the output is `D = tan(ν/2)`. /// /// # Examples /// /// ``` -/// let d = keplerian::anomaly::kepler_parabolic(0.25); -/// assert!((d + d.powi(3) / 3.0 - 0.25).abs() < 1e-12); +/// use keplerian::anomaly::{parabolic_from_mean, ParabolicAnomaly}; +/// let m_d = ParabolicAnomaly::new(0.25); +/// let d = parabolic_from_mean(m_d); +/// assert!((d.value() + d.value().powi(3) / 3.0 - 0.25).abs() < 1e-12); /// ``` #[must_use] -pub fn kepler_parabolic(mean_anomaly_d: f64) -> f64 { - let a = 1.5 * mean_anomaly_d; - (a + (a * a + 1.0).sqrt()).cbrt() - ((a * a + 1.0).sqrt() - a).cbrt() +pub fn parabolic_from_mean(m_d: ParabolicAnomaly) -> ParabolicAnomaly { + let a = 1.5 * m_d.value(); + ParabolicAnomaly::new((a + (a * a + 1.0).sqrt()).cbrt() - ((a * a + 1.0).sqrt() - a).cbrt()) } -/// Converts eccentric anomaly to true anomaly for elliptic motion. +/// Converts true anomaly to eccentric anomaly for elliptic motion. /// /// # Examples /// /// ``` -/// use keplerian::anomaly::true_from_eccentric; +/// use keplerian::anomaly::{eccentric_from_true, TrueAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; -/// let nu = true_from_eccentric(Radians::new(0.5), Eccentricity::new_unchecked(0.1)); -/// assert!(nu.value().is_finite()); +/// let ea = eccentric_from_true(TrueAnomaly::from_value(0.5), Eccentricity::new_unchecked(0.1)); +/// assert!(ea.value().is_finite()); /// ``` #[must_use] -pub fn true_from_eccentric(ea: Radians, ecc: Eccentricity) -> Radians { - let ea_value = ea.value(); +pub fn eccentric_from_true(nu: TrueAnomaly, ecc: Eccentricity) -> EccentricAnomaly { + let nu_value = nu.value(); let ecc_value = ecc.value(); - let s = ((1.0 + ecc_value).sqrt() * (0.5 * ea_value).sin()) - .atan2((1.0 - ecc_value).sqrt() * (0.5 * ea_value).cos()); - Radians::new(wrap_two_pi_raw(2.0 * s)) + let e = 2.0 + * (((1.0 - ecc_value).sqrt() * (0.5 * nu_value).sin()) + .atan2((1.0 + ecc_value).sqrt() * (0.5 * nu_value).cos())); + EccentricAnomaly::new(Radians::new(wrap_two_pi_raw(e))) } -/// Converts true anomaly to eccentric anomaly for elliptic motion. +/// Converts eccentric anomaly to true anomaly for elliptic motion. /// /// # Examples /// /// ``` -/// use keplerian::anomaly::eccentric_from_true; +/// use keplerian::anomaly::{true_from_eccentric, EccentricAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; -/// let ea = eccentric_from_true(Radians::new(0.5), Eccentricity::new_unchecked(0.1)); -/// assert!(ea.value().is_finite()); +/// let nu = true_from_eccentric(EccentricAnomaly::from_value(0.5), Eccentricity::new_unchecked(0.1)); +/// assert!(nu.value().is_finite()); /// ``` #[must_use] -pub fn eccentric_from_true(nu: Radians, ecc: Eccentricity) -> Radians { - let nu_value = nu.value(); +pub fn true_from_eccentric(ea: EccentricAnomaly, ecc: Eccentricity) -> TrueAnomaly { + let ea_value = ea.value(); let ecc_value = ecc.value(); - let e = 2.0 - * (((1.0 - ecc_value).sqrt() * (0.5 * nu_value).sin()) - .atan2((1.0 + ecc_value).sqrt() * (0.5 * nu_value).cos())); - Radians::new(wrap_two_pi_raw(e)) + let s = ((1.0 + ecc_value).sqrt() * (0.5 * ea_value).sin()) + .atan2((1.0 - ecc_value).sqrt() * (0.5 * ea_value).cos()); + TrueAnomaly::new(Radians::new(wrap_two_pi_raw(2.0 * s))) } /// Converts eccentric anomaly to mean anomaly for elliptic motion. @@ -229,16 +473,15 @@ pub fn eccentric_from_true(nu: Radians, ecc: Eccentricity) -> Radians { /// # Examples /// /// ``` -/// use keplerian::anomaly::mean_from_eccentric; +/// use keplerian::anomaly::{mean_from_eccentric, EccentricAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; -/// let m = mean_from_eccentric(Radians::new(0.5), Eccentricity::new_unchecked(0.1)); +/// let m = mean_from_eccentric(EccentricAnomaly::from_value(0.5), Eccentricity::new_unchecked(0.1)); /// assert!((m.value() - (0.5 - 0.1_f64 * 0.5_f64.sin())).abs() < 1e-12); /// ``` #[must_use] -pub fn mean_from_eccentric(ea: Radians, ecc: Eccentricity) -> Radians { +pub fn mean_from_eccentric(ea: EccentricAnomaly, ecc: Eccentricity) -> MeanAnomaly { let ea_value = ea.value(); - Radians::new(ea_value - ecc.value() * ea_value.sin()) + MeanAnomaly::new(Radians::new(ea_value - ecc.value() * ea_value.sin())) } /// Converts mean anomaly to eccentric anomaly for elliptic motion. @@ -250,20 +493,19 @@ pub fn mean_from_eccentric(ea: Radians, ecc: Eccentricity) -> Radians { /// # Examples /// /// ``` -/// use keplerian::anomaly::{eccentric_from_mean, mean_from_eccentric, AnomalyOptions}; +/// use keplerian::anomaly::{eccentric_from_mean, mean_from_eccentric, AnomalyOptions, MeanAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(0.2); -/// let m = Radians::new(1.0); +/// let m = MeanAnomaly::from_value(1.0); /// let ea = eccentric_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); /// let m2 = mean_from_eccentric(ea, ecc); /// assert!((m2.value() - m.value()).abs() < 1e-12); /// ``` pub fn eccentric_from_mean( - m: Radians, + m: MeanAnomaly, ecc: Eccentricity, opts: AnomalyOptions, -) -> Result { +) -> Result { kepler_elliptic(m, ecc, opts) } @@ -276,20 +518,19 @@ pub fn eccentric_from_mean( /// # Examples /// /// ``` -/// use keplerian::anomaly::{true_from_mean, mean_from_true, AnomalyOptions}; +/// use keplerian::anomaly::{true_from_mean, mean_from_true, AnomalyOptions, MeanAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(0.2); -/// let m = Radians::new(1.0); +/// let m = MeanAnomaly::from_value(1.0); /// let nu = true_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); /// let m2 = mean_from_true(nu, ecc); /// assert!((m2.value() - m.value()).abs() < 1e-12); /// ``` pub fn true_from_mean( - m: Radians, + m: MeanAnomaly, ecc: Eccentricity, opts: AnomalyOptions, -) -> Result { +) -> Result { Ok(true_from_eccentric(eccentric_from_mean(m, ecc, opts)?, ecc)) } @@ -298,61 +539,61 @@ pub fn true_from_mean( /// # Examples /// /// ``` -/// use keplerian::anomaly::{mean_from_true, true_from_mean, AnomalyOptions}; +/// use keplerian::anomaly::{mean_from_true, true_from_mean, AnomalyOptions, TrueAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(0.15); -/// let nu = Radians::new(0.8); +/// let nu = TrueAnomaly::from_value(0.8); /// let m = mean_from_true(nu, ecc); /// let nu2 = true_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); /// assert!((nu2.value() - nu.value()).abs() < 1e-12); /// ``` #[must_use] -pub fn mean_from_true(nu: Radians, ecc: Eccentricity) -> Radians { +pub fn mean_from_true(nu: TrueAnomaly, ecc: Eccentricity) -> MeanAnomaly { mean_from_eccentric(eccentric_from_true(nu, ecc), ecc) } -/// Converts hyperbolic anomaly to true anomaly. +// ── Hyperbolic conversions ──────────────────────────────────────────────────── + +/// Converts true anomaly to hyperbolic anomaly. /// /// # Examples /// /// ``` -/// use keplerian::anomaly::{true_from_hyperbolic, hyperbolic_from_true}; +/// use keplerian::anomaly::{hyperbolic_from_true, true_from_hyperbolic, TrueAnomaly}; /// use keplerian::Eccentricity; -/// let ecc = keplerian::Eccentricity::new_unchecked(2.0); -/// let f = 0.5_f64; -/// let nu = true_from_hyperbolic(f, ecc); -/// let f2 = hyperbolic_from_true(nu, ecc); -/// assert!((f2 - f).abs() < 1e-12); +/// let ecc = Eccentricity::new_unchecked(2.0); +/// let nu = TrueAnomaly::from_value(0.4); +/// let f = hyperbolic_from_true(nu, ecc); +/// let nu2 = true_from_hyperbolic(f, ecc); +/// assert!((nu2.value() - nu.value()).abs() < 1e-12); /// ``` #[must_use] -pub fn true_from_hyperbolic(fa: f64, ecc: Eccentricity) -> Radians { +pub fn hyperbolic_from_true(nu: TrueAnomaly, ecc: Eccentricity) -> HyperbolicAnomaly { let ecc_value = ecc.value(); - Radians::new( - 2.0 * (((ecc_value + 1.0).sqrt() * (0.5 * fa).sinh()) - .atan2((ecc_value - 1.0).sqrt() * (0.5 * fa).cosh())), - ) + let t = (0.5 * nu.value()).tan() * ((ecc_value - 1.0) / (ecc_value + 1.0)).sqrt(); + HyperbolicAnomaly::new(2.0 * t.atanh()) } -/// Converts true anomaly to hyperbolic anomaly. +/// Converts hyperbolic anomaly to true anomaly. /// /// # Examples /// /// ``` -/// use keplerian::anomaly::{hyperbolic_from_true, true_from_hyperbolic}; +/// use keplerian::anomaly::{true_from_hyperbolic, hyperbolic_from_true, TrueAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(2.0); -/// let nu = Radians::new(0.4); -/// let f = hyperbolic_from_true(nu, ecc); -/// let nu2 = true_from_hyperbolic(f, ecc); -/// assert!((nu2.value() - nu.value()).abs() < 1e-12); +/// let f = keplerian::anomaly::HyperbolicAnomaly::new(0.5); +/// let nu = true_from_hyperbolic(f, ecc); +/// let f2 = hyperbolic_from_true(nu, ecc); +/// assert!((f2.value() - f.value()).abs() < 1e-12); /// ``` #[must_use] -pub fn hyperbolic_from_true(nu: Radians, ecc: Eccentricity) -> f64 { +pub fn true_from_hyperbolic(fa: HyperbolicAnomaly, ecc: Eccentricity) -> TrueAnomaly { let ecc_value = ecc.value(); - let t = (0.5 * nu.value()).tan() * ((ecc_value - 1.0) / (ecc_value + 1.0)).sqrt(); - 2.0 * t.atanh() + TrueAnomaly::new(Radians::new( + 2.0 * (((ecc_value + 1.0).sqrt() * (0.5 * fa.value()).sinh()) + .atan2((ecc_value - 1.0).sqrt() * (0.5 * fa.value()).cosh())), + )) } /// Converts hyperbolic anomaly to mean anomaly. @@ -360,16 +601,16 @@ pub fn hyperbolic_from_true(nu: Radians, ecc: Eccentricity) -> f64 { /// # Examples /// /// ``` -/// use keplerian::anomaly::mean_from_hyperbolic; +/// use keplerian::anomaly::{mean_from_hyperbolic, HyperbolicAnomaly}; /// use keplerian::Eccentricity; -/// let ecc = keplerian::Eccentricity::new_unchecked(1.5); -/// let f = 0.3_f64; +/// let ecc = Eccentricity::new_unchecked(1.5); +/// let f = HyperbolicAnomaly::new(0.3); /// let m = mean_from_hyperbolic(f, ecc); /// assert!((m.value() - (1.5 * 0.3_f64.sinh() - 0.3)).abs() < 1e-12); /// ``` #[must_use] -pub fn mean_from_hyperbolic(fa: f64, ecc: Eccentricity) -> Radians { - Radians::new(ecc.value() * fa.sinh() - fa) +pub fn mean_from_hyperbolic(fa: HyperbolicAnomaly, ecc: Eccentricity) -> MeanAnomaly { + MeanAnomaly::new(Radians::new(ecc.value() * fa.value().sinh() - fa.value())) } /// Converts mean anomaly to hyperbolic anomaly. @@ -381,22 +622,56 @@ pub fn mean_from_hyperbolic(fa: f64, ecc: Eccentricity) -> Radians { /// # Examples /// /// ``` -/// use keplerian::anomaly::{hyperbolic_from_mean, mean_from_hyperbolic, AnomalyOptions}; +/// use keplerian::anomaly::{hyperbolic_from_mean, mean_from_hyperbolic, AnomalyOptions, MeanAnomaly}; /// use keplerian::Eccentricity; -/// use qtty::angular::Radians; /// let ecc = Eccentricity::new_unchecked(1.5); -/// let m = Radians::new(1.0); +/// let m = MeanAnomaly::from_value(1.0); /// let f = hyperbolic_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); /// assert!((mean_from_hyperbolic(f, ecc).value() - m.value()).abs() < 1e-12); /// ``` pub fn hyperbolic_from_mean( - m: Radians, + m: MeanAnomaly, ecc: Eccentricity, opts: AnomalyOptions, -) -> Result { - kepler_hyperbolic(m, ecc, opts).map(|value| value.value()) +) -> Result { + kepler_hyperbolic(m, ecc, opts) +} + +// ── Parabolic conversions ───────────────────────────────────────────────────── + +/// Converts parabolic anomaly `D = tan(ν/2)` to true anomaly. +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::{true_from_parabolic, parabolic_from_true, TrueAnomaly}; +/// let nu = TrueAnomaly::from_value(0.6); +/// let d = parabolic_from_true(nu); +/// let nu2 = true_from_parabolic(d); +/// assert!((nu2.value() - nu.value()).abs() < 1e-12); +/// ``` +#[must_use] +pub fn true_from_parabolic(dp: ParabolicAnomaly) -> TrueAnomaly { + TrueAnomaly::new(Radians::new(2.0 * dp.value().atan())) } +/// Converts true anomaly to parabolic anomaly `D = tan(ν/2)`. +/// +/// # Examples +/// +/// ``` +/// use keplerian::anomaly::{parabolic_from_true, true_from_parabolic, TrueAnomaly}; +/// let nu = TrueAnomaly::from_value(0.6); +/// let d = parabolic_from_true(nu); +/// assert!((d.value() - (0.6_f64 / 2.0).tan()).abs() < 1e-12); +/// ``` +#[must_use] +pub fn parabolic_from_true(nu: TrueAnomaly) -> ParabolicAnomaly { + ParabolicAnomaly::new((nu.value() / 2.0).tan()) +} + +// ── Internal helpers ────────────────────────────────────────────────────────── + /// Wraps a raw angle to `[0, 2π)` for internal numeric use. #[must_use] pub(crate) fn wrap_two_pi_raw(x: f64) -> f64 { @@ -496,22 +771,36 @@ mod tests { fn solvers_converge() { let zero = Eccentricity::new(0.0).unwrap(); assert_eq!( - kepler_elliptic(Radians::new(0.4), zero, AnomalyOptions::default()).unwrap(), - Radians::new(0.4) + kepler_elliptic( + MeanAnomaly::from_value(0.4), + zero, + AnomalyOptions::default() + ) + .unwrap() + .value(), + 0.4 ); let elliptic_ecc = Eccentricity::new(0.4).unwrap(); - let e = - kepler_elliptic(Radians::new(1.0), elliptic_ecc, AnomalyOptions::default()).unwrap(); + let e = kepler_elliptic( + MeanAnomaly::from_value(1.0), + elliptic_ecc, + AnomalyOptions::default(), + ) + .unwrap(); assert!((mean_from_eccentric(e, elliptic_ecc).value() - 1.0).abs() < 1e-12); let hyperbolic_ecc = Eccentricity::new(1.4).unwrap(); - let f = kepler_hyperbolic(Radians::new(1.0), hyperbolic_ecc, AnomalyOptions::default()) - .unwrap(); - assert!((mean_from_hyperbolic(f.value(), hyperbolic_ecc).value() - 1.0).abs() < 1e-12); + let f = kepler_hyperbolic( + MeanAnomaly::from_value(1.0), + hyperbolic_ecc, + AnomalyOptions::default(), + ) + .unwrap(); + assert!((mean_from_hyperbolic(f, hyperbolic_ecc).value() - 1.0).abs() < 1e-12); } #[test] fn round_trip_elliptic_anomalies() { - let nu = Radians::new(1.2); + let nu = TrueAnomaly::from_value(1.2); let ecc = Eccentricity::new(0.3).unwrap(); let ea = eccentric_from_true(nu, ecc); let m = mean_from_eccentric(ea, ecc); @@ -523,19 +812,19 @@ mod tests { #[test] fn round_trip_hyperbolic_anomalies() { - let nu = Radians::new(0.8); + let nu = TrueAnomaly::from_value(0.8); let ecc = Eccentricity::new(1.7).unwrap(); let f = hyperbolic_from_true(nu, ecc); let m = mean_from_hyperbolic(f, ecc); let f2 = hyperbolic_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); - assert!((f2 - f).abs() < 1e-12); + assert!((f2.value() - f.value()).abs() < 1e-12); } #[test] fn invalid_eccentricity() { assert!(matches!( kepler_elliptic( - Radians::new(0.0), + MeanAnomaly::from_value(0.0), Eccentricity::new_unchecked(1.0), AnomalyOptions::default() ), @@ -543,7 +832,7 @@ mod tests { )); assert!(matches!( kepler_hyperbolic( - Radians::new(0.0), + MeanAnomaly::from_value(0.0), Eccentricity::new_unchecked(0.9), AnomalyOptions::default() ), @@ -551,15 +840,22 @@ mod tests { )); } + #[test] + fn try_new_rejects_invalid_options() { + assert!(AnomalyOptions::try_new(0, 1e-12).is_err()); + assert!(AnomalyOptions::try_new(10, 0.0).is_err()); + assert!(AnomalyOptions::try_new(10, -1.0).is_err()); + assert!(AnomalyOptions::try_new(10, f64::NAN).is_err()); + } + #[test] fn detects_non_convergence() { + // 1 iteration with extreme tolerance forces non-convergence. + let opts = AnomalyOptions::try_new(1, 1e-20).unwrap(); let err = kepler_elliptic( - Radians::new(1.0), + MeanAnomaly::from_value(1.0), Eccentricity::new(0.9).unwrap(), - AnomalyOptions { - max_iter: 0, - tol: 1e-16, - }, + opts, ) .unwrap_err(); assert!(matches!(err, AnomalyError::NotConverged { .. })); @@ -568,18 +864,112 @@ mod tests { #[test] fn near_parabolic_hyperbolic_solver_converges() { let ecc = Eccentricity::new(1.0000001).unwrap(); - let mean = Radians::new(0.001_f64.to_radians() + 0.01720209895 * 0.5); - let f = kepler_hyperbolic( - mean, - ecc, - AnomalyOptions { - max_iter: 100, - tol: 1e-14, - }, - ) - .unwrap(); + let mean = MeanAnomaly::from_value(0.001_f64.to_radians() + 0.01720209895 * 0.5); + let f = kepler_hyperbolic(mean, ecc, AnomalyOptions::try_new(100, 1e-14).unwrap()).unwrap(); assert!(f.value().is_finite()); assert!(hyperbolic_residual(f.value(), ecc.value(), mean.value()).abs() < 1e-13); } + + #[test] + fn elliptic_circular_orbit_returns_mean_anomaly() { + let zero = Eccentricity::new(0.0).unwrap(); + let m = MeanAnomaly::from_value(1.23); + let e = eccentric_from_mean(m, zero, AnomalyOptions::default()).unwrap(); + assert!((e.value() - m.value()).abs() < 1e-14); + } + + #[test] + fn hyperbolic_mean_anomaly_zero_returns_zero() { + let ecc = Eccentricity::new_unchecked(1.5); + let f = hyperbolic_from_mean(MeanAnomaly::from_value(0.0), ecc, AnomalyOptions::default()) + .unwrap(); + assert!(f.value().abs() < 1e-14); + } + + #[test] + fn hyperbolic_mean_anomaly_large_branch() { + let ecc = Eccentricity::new_unchecked(1.2); + let m = MeanAnomaly::from_value(100.0); + assert!(hyperbolic_from_mean(m, ecc, AnomalyOptions::default()).is_ok()); + } + + #[test] + fn parabolic_round_trips() { + let m = 0.5_f64; + let d = parabolic_from_mean(ParabolicAnomaly::new(m)); + let m_back = d.value() + d.value().powi(3) / 3.0; + assert!((m_back - m).abs() < 1e-12); + } + + #[test] + fn hyperbolic_anomaly_round_trips_all_conversions() { + let ecc = Eccentricity::new_unchecked(2.0); + let nu = TrueAnomaly::from_value(0.8); + let f = hyperbolic_from_true(nu, ecc); + let nu2 = true_from_hyperbolic(f, ecc); + assert!((nu2.value() - nu.value()).abs() < 1e-12); + let m = mean_from_hyperbolic(f, ecc); + let f2 = hyperbolic_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); + assert!((f2.value() - f.value()).abs() < 1e-12); + } + + #[test] + fn hyperbolic_kepler_rejects_nan_mean_anomaly() { + let ecc = Eccentricity::new_unchecked(1.5); + let err = hyperbolic_from_mean( + MeanAnomaly::from_value(f64::NAN), + ecc, + AnomalyOptions::default(), + ); + assert!(err.is_err()); + } + + #[test] + fn elliptic_kepler_rejects_nan_mean_anomaly() { + let ecc = Eccentricity::new_unchecked(0.5); + let err = eccentric_from_mean( + MeanAnomaly::from_value(f64::NAN), + ecc, + AnomalyOptions::default(), + ); + assert!(err.is_err()); + } + + #[test] + fn elliptic_bisection_fallback_path_is_exercised() { + let ecc = Eccentricity::new_unchecked(0.5); + let opts = AnomalyOptions::try_new(2, 1e-300).unwrap(); + let _ = kepler_elliptic(MeanAnomaly::from_value(1.0), ecc, opts); + } + + #[test] + fn hyperbolic_bisection_fallback_path_is_exercised() { + let ecc = Eccentricity::new_unchecked(1.5); + let opts = AnomalyOptions::try_new(2, 1e-300).unwrap(); + let _ = kepler_hyperbolic(MeanAnomaly::from_value(1.0), ecc, opts); + } + + #[test] + fn hyperbolic_bisection_large_mean_anomaly_path_is_exercised() { + let ecc = Eccentricity::new_unchecked(1.1); + let opts = AnomalyOptions::try_new(2, 1e-300).unwrap(); + let _ = kepler_hyperbolic(MeanAnomaly::from_value(50.0), ecc, opts); + } + + #[test] + fn anomaly_newtypes_round_trip_values() { + let ea = EccentricAnomaly::from_value(1.2); + assert_eq!(ea.value(), 1.2); + assert_eq!(ea.radians().value(), 1.2); + + let ta = TrueAnomaly::from_value(0.5); + assert_eq!(ta.value(), 0.5); + + let ma = MeanAnomaly::from_value(0.8); + assert_eq!(ma.value(), 0.8); + + let pa = ParabolicAnomaly::new(0.3); + assert_eq!(pa.value(), 0.3); + } } diff --git a/src/eccentricity.rs b/src/eccentricity.rs index 6769628..94f0737 100644 --- a/src/eccentricity.rs +++ b/src/eccentricity.rs @@ -1,23 +1,53 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon -//! Domain-semantic eccentricity scalar. +//! Domain-semantic eccentricity scalar and conic regime classifier. //! //! ## Scientific scope //! Eccentricity classifies conic sections in central-force motion. This module -//! models only the scalar `e` itself; it does not attach any frame, epoch, or -//! body semantics. +//! models the scalar `e` itself and the resulting [`ConicRegime`]; it does not +//! attach any frame, epoch, or body semantics. //! //! ## Technical scope //! [`Eccentricity`] wraps a raw `f64` so public APIs can distinguish orbital //! eccentricity from unrelated dimensionless diagnostics such as tolerances or -//! residuals. +//! residuals. [`ConicRegime`] is derived via [`Eccentricity::classify`]. //! //! ## References //! - Battin, R. H. (1999). *An Introduction to the Mathematics and Methods of //! Astrodynamics*. //! - Vallado, D. A. (2013). *Fundamentals of Astrodynamics and Applications*. +/// Conic regime classified by eccentricity. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum ConicRegime { + /// Bounded ellipse, `0 ≤ e < 1 − ε`. + Elliptic, + /// Parabola, `|e − 1| ≤ ε`. + Parabolic, + /// Unbounded hyperbola, `e > 1 + ε`. + Hyperbolic, +} + +/// Errors returned by fallible eccentricity constructors. +#[derive(Debug, Clone, Copy, PartialEq, thiserror::Error)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] +pub enum EccentricityError { + /// The value is negative. + #[error("eccentricity must be non-negative, got {0}")] + Negative(f64), + /// The value is not finite. + #[error("eccentricity must be finite, got {0}")] + NotFinite(f64), + /// An elliptic constructor received `e ≥ 1`. + #[error("elliptic eccentricity requires e < 1, got {0}")] + NotElliptic(f64), + /// A hyperbolic constructor received `e ≤ 1`. + #[error("hyperbolic eccentricity requires e > 1, got {0}")] + NotHyperbolic(f64), +} + /// Dimensionless eccentricity scalar (domain-semantic newtype). /// /// Prevents accidentally mixing eccentricity with tolerances, residuals, @@ -36,9 +66,7 @@ pub struct Eccentricity(f64); impl Eccentricity { - /// Creates a new eccentricity. - /// - /// Returns `None` if `e` is negative or non-finite. + /// Creates a new eccentricity, returning `None` if `e` is negative or non-finite. /// /// # Examples /// @@ -52,6 +80,85 @@ impl Eccentricity { (e.is_finite() && e >= 0.0).then_some(Self(e)) } + /// Creates a new eccentricity, returning [`EccentricityError`] if invalid. + /// + /// # Errors + /// + /// Returns [`EccentricityError::NotFinite`] for NaN/±∞, or + /// [`EccentricityError::Negative`] for `e < 0`. + /// + /// # Examples + /// + /// ``` + /// use keplerian::eccentricity::{Eccentricity, EccentricityError}; + /// assert!(Eccentricity::try_new(0.5).is_ok()); + /// assert!(matches!(Eccentricity::try_new(-0.1), Err(EccentricityError::Negative(_)))); + /// ``` + pub fn try_new(e: f64) -> Result { + if !e.is_finite() { + return Err(EccentricityError::NotFinite(e)); + } + if e < 0.0 { + return Err(EccentricityError::Negative(e)); + } + Ok(Self(e)) + } + + /// Creates an elliptic eccentricity (`0 ≤ e < 1`). + /// + /// # Errors + /// + /// Returns [`EccentricityError`] if `e < 0`, non-finite, or `e ≥ 1`. + /// + /// # Examples + /// + /// ``` + /// use keplerian::eccentricity::{Eccentricity, EccentricityError}; + /// assert!(Eccentricity::new_elliptic(0.5).is_ok()); + /// assert!(matches!(Eccentricity::new_elliptic(1.5), Err(EccentricityError::NotElliptic(_)))); + /// ``` + pub fn new_elliptic(e: f64) -> Result { + let ec = Self::try_new(e)?; + if e >= 1.0 { + return Err(EccentricityError::NotElliptic(e)); + } + Ok(ec) + } + + /// Creates a hyperbolic eccentricity (`e > 1`). + /// + /// # Errors + /// + /// Returns [`EccentricityError`] if `e < 0`, non-finite, or `e ≤ 1`. + /// + /// # Examples + /// + /// ``` + /// use keplerian::eccentricity::{Eccentricity, EccentricityError}; + /// assert!(Eccentricity::new_hyperbolic(1.5).is_ok()); + /// assert!(matches!(Eccentricity::new_hyperbolic(0.5), Err(EccentricityError::NotHyperbolic(_)))); + /// ``` + pub fn new_hyperbolic(e: f64) -> Result { + let ec = Self::try_new(e)?; + if e <= 1.0 { + return Err(EccentricityError::NotHyperbolic(e)); + } + Ok(ec) + } + + /// Returns the parabolic eccentricity `e = 1.0` exactly. + /// + /// # Examples + /// + /// ``` + /// use keplerian::eccentricity::Eccentricity; + /// assert_eq!(Eccentricity::parabolic().value(), 1.0); + /// ``` + #[must_use] + pub const fn parabolic() -> Self { + Self(1.0) + } + /// Creates a new eccentricity without checking. /// /// # Examples @@ -79,6 +186,30 @@ impl Eccentricity { self.0 } + /// Classifies this eccentricity into a [`ConicRegime`] using tolerance `eps`. + /// + /// `|e − 1| ≤ eps` → [`ConicRegime::Parabolic`]; `e < 1 − eps` → Elliptic; + /// `e > 1 + eps` → Hyperbolic. + /// + /// # Examples + /// + /// ``` + /// use keplerian::eccentricity::{ConicRegime, Eccentricity}; + /// assert_eq!(Eccentricity::new_unchecked(0.5).classify(1e-10), ConicRegime::Elliptic); + /// assert_eq!(Eccentricity::new_unchecked(1.0).classify(1e-10), ConicRegime::Parabolic); + /// assert_eq!(Eccentricity::new_unchecked(1.5).classify(1e-10), ConicRegime::Hyperbolic); + /// ``` + #[must_use] + pub fn classify(self, eps: f64) -> ConicRegime { + if (self.0 - 1.0).abs() <= eps { + ConicRegime::Parabolic + } else if self.0 < 1.0 { + ConicRegime::Elliptic + } else { + ConicRegime::Hyperbolic + } + } + /// Returns `true` for an elliptic orbit (`0 ≤ e < 1`). /// /// # Examples @@ -119,3 +250,98 @@ impl Eccentricity { self.0 > 1.0 } } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn is_hyperbolic_classifies_strictly_above_one() { + assert!(Eccentricity::new_unchecked(1.5).is_hyperbolic()); + assert!(!Eccentricity::new_unchecked(0.5).is_hyperbolic()); + } + + #[test] + fn try_new_accepts_valid_values() { + let e = Eccentricity::try_new(0.5).unwrap(); + assert!((e.value() - 0.5).abs() < 1e-15); + assert_eq!(Eccentricity::try_new(0.0).unwrap().value(), 0.0); + } + + #[test] + fn try_new_rejects_invalid_values() { + assert!(matches!( + Eccentricity::try_new(-0.1), + Err(EccentricityError::Negative(_)) + )); + assert!(matches!( + Eccentricity::try_new(f64::NAN), + Err(EccentricityError::NotFinite(_)) + )); + assert!(matches!( + Eccentricity::try_new(f64::INFINITY), + Err(EccentricityError::NotFinite(_)) + )); + } + + #[test] + fn new_elliptic_validates_regime() { + let e = Eccentricity::new_elliptic(0.3).unwrap(); + assert!((e.value() - 0.3).abs() < 1e-15); + assert!(matches!( + Eccentricity::new_elliptic(1.0), + Err(EccentricityError::NotElliptic(_)) + )); + assert!(matches!( + Eccentricity::new_elliptic(1.5), + Err(EccentricityError::NotElliptic(_)) + )); + assert!(matches!( + Eccentricity::new_elliptic(-0.1), + Err(EccentricityError::Negative(_)) + )); + } + + #[test] + fn new_hyperbolic_validates_regime() { + let e = Eccentricity::new_hyperbolic(1.5).unwrap(); + assert!((e.value() - 1.5).abs() < 1e-15); + assert!(matches!( + Eccentricity::new_hyperbolic(1.0), + Err(EccentricityError::NotHyperbolic(_)) + )); + assert!(matches!( + Eccentricity::new_hyperbolic(0.5), + Err(EccentricityError::NotHyperbolic(_)) + )); + assert!(matches!( + Eccentricity::new_hyperbolic(f64::NAN), + Err(EccentricityError::NotFinite(_)) + )); + } + + #[test] + fn parabolic_constructor_is_exactly_one() { + assert_eq!(Eccentricity::parabolic().value(), 1.0); + } + + #[test] + fn classify_uses_tolerance_band() { + assert_eq!( + Eccentricity::new_unchecked(0.5).classify(1e-10), + ConicRegime::Elliptic + ); + assert_eq!( + Eccentricity::new_unchecked(1.0).classify(1e-10), + ConicRegime::Parabolic + ); + assert_eq!( + Eccentricity::new_unchecked(1.0 + 5e-11).classify(1e-10), + ConicRegime::Parabolic + ); + assert_eq!( + Eccentricity::new_unchecked(1.5).classify(1e-10), + ConicRegime::Hyperbolic + ); + } +} diff --git a/src/elements.rs b/src/elements.rs index 6683e23..60719ef 100644 --- a/src/elements.rs +++ b/src/elements.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Typed Keplerian orbital elements and Cartesian conversions. @@ -33,18 +33,9 @@ use crate::eccentricity::Eccentricity; use crate::state::CartesianState; use crate::vec3::{cross, dot, norm, scale, sub}; -const EPS: f64 = 1.0e-10; +pub use crate::eccentricity::ConicRegime; -/// Conic regime classified by eccentricity. -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum ConicRegime { - /// Bounded ellipse, `0 ≤ e < 1`. - Elliptic, - /// Parabola, `e = 1` within numerical tolerance. - Parabolic, - /// Unbounded hyperbola, `e > 1`. - Hyperbolic, -} +const EPS: f64 = 1.0e-10; /// Errors returned by element validation and Cartesian conversion. #[derive(Debug, Clone, Copy, PartialEq, thiserror::Error)] @@ -66,6 +57,14 @@ pub enum ConversionError { /// The conversion encountered a degenerate geometry. #[error("degenerate orbital geometry: {0}")] Degenerate(&'static str), + /// Semi-major axis sign is inconsistent with the eccentricity regime. + #[error("semi-major axis {sma} km is incoherent with eccentricity {ecc}")] + IncoherentRegime { + /// Semi-major axis value. + sma: f64, + /// Eccentricity value. + ecc: f64, + }, } /// Classical Keplerian elements with typed distance and angular quantities. @@ -135,6 +134,14 @@ impl KeplerianElements { if !(0.0..=PI).contains(&inclination.value()) { return Err(ConversionError::InvalidInclination(inclination.value())); } + let e = eccentricity.value(); + let a = semi_major_axis.value(); + if e < 1.0 - EPS && a <= 0.0 { + return Err(ConversionError::IncoherentRegime { sma: a, ecc: e }); + } + if e > 1.0 + EPS && a >= 0.0 { + return Err(ConversionError::IncoherentRegime { sma: a, ecc: e }); + } Ok(Self { semi_major_axis, eccentricity, @@ -160,6 +167,11 @@ impl KeplerianElements { /// Converts elements to a typed Cartesian state. /// + /// # Errors + /// + /// Returns [`ConversionError`] if `mu ≤ 0`, semi-latus rectum `p ≤ 0`, or + /// the orbit denominator `1 + e cos(ν) ≤ ε`. + /// /// # Examples /// /// ``` @@ -181,22 +193,35 @@ impl KeplerianElements { /// Radians::new(0.0), /// Radians::new(0.0), /// ).unwrap(); - /// let state = el.to_cartesian::(GravitationalParameter::new(398600.4418)); + /// let state = el.try_to_cartesian::(GravitationalParameter::new(398600.4418)).unwrap(); /// let r = state.position(); /// assert!((r.x().value().hypot(r.y().value()).hypot(r.z().value()) - 7000.0).abs() < 1.0); /// ``` - #[must_use] - pub fn to_cartesian>( + pub fn try_to_cartesian>( &self, mu: GravitationalParameter, - ) -> CartesianState { + ) -> Result, ConversionError> { + let mu_val = mu.value(); + if !mu_val.is_finite() || mu_val <= 0.0 { + return Err(ConversionError::Degenerate( + "non-positive or non-finite gravitational parameter", + )); + } let a = self.semi_major_axis.value(); let e = self.eccentricity.value(); let nu = self.true_anomaly.value(); let p = a * (1.0 - e * e); + if p <= 0.0 { + return Err(ConversionError::Degenerate( + "non-positive semi-latus rectum", + )); + } let denom = 1.0 + e * nu.cos(); + if denom <= EPS { + return Err(ConversionError::Degenerate("degenerate orbit denominator")); + } let r = p / denom; - let root = (mu.value() / p).sqrt(); + let root = (mu_val / p).sqrt(); let r_pqw = [r * nu.cos(), r * nu.sin(), 0.0]; let v_pqw = [-root * nu.sin(), root * (e + nu.cos()), 0.0]; let r_ijk = rotate_pqw( @@ -211,10 +236,10 @@ impl KeplerianElements { self.inclination.value(), self.arg_periapsis.value(), ); - CartesianState::new( + Ok(CartesianState::new( Position::::new(r_ijk[0], r_ijk[1], r_ijk[2]), Velocity::::new(v_ijk[0], v_ijk[1], v_ijk[2]), - ) + )) } /// Converts a typed Cartesian state into classical Keplerian elements. @@ -395,10 +420,218 @@ mod tests { Radians::new(1.1), ) .unwrap(); - let st = el.to_cartesian::
(GravitationalParameter::new(398600.4418)); + let st = el + .try_to_cartesian::
(GravitationalParameter::new(398600.4418)) + .unwrap(); let back = KeplerianElements::from_cartesian(&st, GravitationalParameter::new(398600.4418)) .unwrap(); assert!((back.semi_major_axis.value() - el.semi_major_axis.value()).abs() < 1e-8); assert!((back.eccentricity.value() - el.eccentricity.value()).abs() < 1e-12); } + + #[test] + fn new_rejects_negative_eccentricity() { + let err = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(-0.1), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ); + assert!(matches!(err, Err(ConversionError::InvalidEccentricity(_)))); + } + + #[test] + fn new_rejects_inclination_out_of_range() { + let err = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(0.1), + Radians::new(PI + 0.1), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ); + assert!(matches!(err, Err(ConversionError::InvalidInclination(_)))); + } + + #[test] + fn conic_kind_classifies_parabolic_and_hyperbolic() { + let parabolic = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(1.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ) + .unwrap(); + assert_eq!(parabolic.conic_kind(), ConicRegime::Parabolic); + + let hyperbolic = KeplerianElements::::new( + Kilometers::new(-40000.0), + Eccentricity::new_unchecked(1.5), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ) + .unwrap(); + assert_eq!(hyperbolic.conic_kind(), ConicRegime::Hyperbolic); + } + + #[test] + fn from_cartesian_rejects_degenerate_states() { + let mu = GravitationalParameter::new(398600.4418); + let zero_position = CartesianState::::new( + Position::::new(0.0, 0.0, 0.0), + Velocity::::new(0.0, 7.5, 0.0), + ); + assert!(KeplerianElements::::from_cartesian(&zero_position, mu).is_err()); + + let zero_angular_momentum = CartesianState::::new( + Position::::new(7000.0, 0.0, 0.0), + Velocity::::new(7.5, 0.0, 0.0), + ); + assert!(KeplerianElements::::from_cartesian(&zero_angular_momentum, mu).is_err()); + } + + #[test] + fn from_cartesian_equatorial_eccentric_orbit() { + let mu = GravitationalParameter::new(398600.4418); + let state = CartesianState::::new( + Position::::new(7000.0, 0.0, 0.0), + Velocity::::new(0.0, 6.5, 0.0), + ); + let el = KeplerianElements::::from_cartesian(&state, mu); + assert!(el.is_ok(), "{el:?}"); + } + + #[test] + fn new_rejects_incoherent_regime_signs() { + let elliptic = KeplerianElements::::new( + Kilometers::new(-7000.0), + Eccentricity::new_unchecked(0.5), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ); + assert!(matches!( + elliptic, + Err(ConversionError::IncoherentRegime { .. }) + )); + + let hyperbolic = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(1.5), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ); + assert!(matches!( + hyperbolic, + Err(ConversionError::IncoherentRegime { .. }) + )); + } + + #[test] + fn try_to_cartesian_rejects_invalid_mu_and_degenerate_geometry() { + let el = KeplerianElements::::new( + Kilometers::new(-10000.0), + Eccentricity::new_unchecked(1.5), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ) + .unwrap(); + assert!(matches!( + el.try_to_cartesian::
(GravitationalParameter::new(-1.0)), + Err(ConversionError::Degenerate(_)) + )); + + let parabolic = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(1.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + ) + .unwrap(); + assert!(matches!( + parabolic.try_to_cartesian::
(GravitationalParameter::new(398600.4418)), + Err(ConversionError::Degenerate(_)) + )); + } + + #[test] + fn try_to_cartesian_rejects_degenerate_orbit_denominator() { + let nu = (-1.0_f64 / 1.5).acos(); + let el = KeplerianElements::::new( + Kilometers::new(-10000.0), + Eccentricity::new_unchecked(1.5), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(0.0), + Radians::new(nu), + ) + .unwrap(); + assert!(el + .try_to_cartesian::
(GravitationalParameter::new(398600.4418)) + .is_err()); + } + + #[test] + fn from_cartesian_rejects_invalid_mu() { + let state = CartesianState::::new( + Position::::new(7000.0, 0.0, 0.0), + Velocity::::new(0.0, 7.5, 0.0), + ); + assert!(KeplerianElements::::from_cartesian( + &state, + GravitationalParameter::new(-1.0) + ) + .is_err()); + assert!(matches!( + KeplerianElements::::from_cartesian( + &state, + GravitationalParameter::new(f64::NAN) + ), + Err(ConversionError::NonFiniteValue { .. }) + )); + } + + #[test] + fn from_cartesian_circular_inclined_orbit() { + let mu = GravitationalParameter::new(398600.4418); + let r = 7000.0_f64; + let v_circ = (mu.value() / r).sqrt(); + let vy = v_circ * core::f64::consts::FRAC_1_SQRT_2; + let vz = v_circ * core::f64::consts::FRAC_1_SQRT_2; + let state = CartesianState::::new( + Position::::new(r, 0.0, 0.0), + Velocity::::new(0.0, vy, vz), + ); + let el = KeplerianElements::::from_cartesian(&state, mu); + assert!(el.is_ok(), "{el:?}"); + } + + #[test] + fn try_to_cartesian_produces_finite_state() { + let mu = GravitationalParameter::new(398600.4418); + let el = KeplerianElements::::new( + Kilometers::new(7000.0), + Eccentricity::new_unchecked(0.1), + Radians::new(0.5), + Radians::new(1.0), + Radians::new(0.3), + Radians::new(0.7), + ) + .unwrap(); + let state = el.try_to_cartesian::
(mu).unwrap(); + assert!(state.position().x().value().is_finite()); + } } diff --git a/src/error.rs b/src/error.rs index d6973bf..9723ac0 100644 --- a/src/error.rs +++ b/src/error.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Top-level error aggregation for crate workflows. @@ -57,3 +57,36 @@ impl From for KeplerError { Self::Lambert(value) } } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn kepler_error_from_anomaly() { + let inner = anomaly::AnomalyError::InvalidEccentricity(1.5); + let e = KeplerError::from(inner); + assert!(matches!(e, KeplerError::Anomaly(_))); + } + + #[test] + fn kepler_error_from_conversion() { + let inner = elements::ConversionError::InvalidEccentricity(-1.0); + let e = KeplerError::from(inner); + assert!(matches!(e, KeplerError::Conversion(_))); + } + + #[test] + fn kepler_error_from_propagation() { + let inner = problem::PropagationError::ParabolicUnsupported; + let e = KeplerError::from(inner); + assert!(matches!(e, KeplerError::Propagation(_))); + } + + #[test] + fn kepler_error_from_lambert() { + let inner = lambert::LambertError::ZeroPosition; + let e = KeplerError::from(inner); + assert!(matches!(e, KeplerError::Lambert(_))); + } +} diff --git a/src/lambert/error.rs b/src/lambert/error.rs index 3ca2630..dc4b0c0 100644 --- a/src/lambert/error.rs +++ b/src/lambert/error.rs @@ -12,7 +12,8 @@ //! - Izzo, D. (2014). *Revisiting Lambert's Problem*. /// Errors returned by the Lambert solver. -#[derive(Debug, thiserror::Error, PartialEq)] +#[derive(Debug, Clone, thiserror::Error, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub enum LambertError { /// Gravitational parameter must be strictly positive (km³/s²). #[error("non-positive gravitational parameter ({0})")] diff --git a/src/lambert/izzo.rs b/src/lambert/izzo.rs index eb9e558..cbe443c 100644 --- a/src/lambert/izzo.rs +++ b/src/lambert/izzo.rs @@ -25,6 +25,7 @@ use super::error::LambertError; /// Direction of revolution about the central body. #[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub enum LambertBranch { /// Prograde (counter-clockwise as seen from `+z`). Prograde, @@ -34,6 +35,7 @@ pub enum LambertBranch { /// Side selection for the multi-revolution branch. #[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub enum NRevBranch { /// Low-energy / short-period side (smaller semi-major axis). Left, @@ -42,7 +44,8 @@ pub enum NRevBranch { } /// Diagnostics produced by a single Lambert solve. -#[derive(Debug, Clone, Copy)] +#[derive(Debug, Clone, Copy, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct LambertDiagnostics { /// Householder iterations consumed before convergence. pub iterations: u32, diff --git a/src/lambert/mod.rs b/src/lambert/mod.rs index 9c262a2..50b0caa 100644 --- a/src/lambert/mod.rs +++ b/src/lambert/mod.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Typed Lambert boundary-value solving. diff --git a/src/lambert/typed.rs b/src/lambert/typed.rs index cfcdebe..268522e 100644 --- a/src/lambert/typed.rs +++ b/src/lambert/typed.rs @@ -26,7 +26,8 @@ use super::izzo::{ }; /// Typed Lambert solution — departure / arrival velocities plus diagnostics. -#[derive(Debug, Clone, Copy)] +#[derive(Debug, Clone, Copy, PartialEq)] +#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] pub struct TypedLambertSolution { /// Departure velocity at `r1`, in the same frame as the input positions. pub v1: Velocity, @@ -210,4 +211,31 @@ mod tests { other => panic!("unexpected error: {other:?}"), } } + + #[test] + fn n_rev_valid_geometry_paths_are_exercised() { + let r1 = Position::<(), ICRS, Kilometer>::new(7000.0, 0.0, 0.0); + let r2 = Position::<(), ICRS, Kilometer>::new(0.0, 7000.0, 0.0); + let tof = Second::new(10800.0); + let mu = GravitationalParameter::new(398600.4418); + + let _ = lambert_n_rev( + r1, + r2, + tof, + mu, + LambertBranch::Prograde, + 1, + NRevBranch::Left, + ); + let _ = lambert_n_rev( + r1, + r2, + tof, + mu, + LambertBranch::Retrograde, + 1, + NRevBranch::Right, + ); + } } diff --git a/src/lib.rs b/src/lib.rs index f75ca90..f52ecdb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Domain-agnostic Keplerian dynamics on typed quantities. @@ -18,6 +18,7 @@ #![forbid(unsafe_code)] #![deny(missing_docs)] +#![cfg_attr(not(feature = "std"), no_std)] #![doc = include_str!("../README.md")] #[cfg(feature = "alloc")] @@ -28,6 +29,7 @@ pub mod eccentricity; pub mod elements; pub mod error; pub mod lambert; +pub mod prelude; pub mod problem; pub mod state; pub mod transfer; diff --git a/src/prelude.rs b/src/prelude.rs new file mode 100644 index 0000000..220bdee --- /dev/null +++ b/src/prelude.rs @@ -0,0 +1,24 @@ +// SPDX-License-Identifier: AGPL-3.0-only +// Copyright (C) 2026 Vallés Puig, Ramon + +//! Convenience re-exports of the most commonly used public items. +//! +//! ``` +//! use keplerian::prelude::*; +//! ``` + +pub use crate::anomaly::{ + AnomalyError, AnomalyOptions, EccentricAnomaly, HyperbolicAnomaly, MeanAnomaly, + ParabolicAnomaly, TrueAnomaly, +}; +pub use crate::eccentricity::{ConicRegime, Eccentricity, EccentricityError}; +pub use crate::elements::{ConversionError, KeplerianElements}; +pub use crate::error::KeplerError; +pub use crate::lambert::{LambertBranch, LambertError, NRevBranch}; +pub use crate::problem::{KeplerProblem, PropagationError}; +pub use crate::state::CartesianState; +pub use crate::transfer::{ + escape_speed, hohmann_delta_v, orbital_period, specific_angular_momentum, + specific_orbital_energy, try_escape_speed, try_hohmann_delta_v, try_orbital_period, + try_vis_viva_speed, vis_viva_speed, HohmannResult, TransferError, +}; diff --git a/src/problem.rs b/src/problem.rs index fd8aa8c..d83d948 100644 --- a/src/problem.rs +++ b/src/problem.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Typed central-force two-body propagation context. @@ -26,7 +26,7 @@ use qtty::Second; use crate::anomaly::{ hyperbolic_from_mean, hyperbolic_from_true, mean_from_hyperbolic, mean_from_true, - true_from_hyperbolic, true_from_mean, AnomalyError, AnomalyOptions, + true_from_hyperbolic, true_from_mean, AnomalyError, AnomalyOptions, MeanAnomaly, TrueAnomaly, }; use crate::elements::{ConicRegime, ConversionError, KeplerianElements}; use crate::state::CartesianState; @@ -51,8 +51,7 @@ pub enum PropagationError { /// A central-force Kepler problem for a center/frame pair. #[derive(Debug, Clone, Copy)] pub struct KeplerProblem { - /// Standard gravitational parameter of the central potential. - pub mu: GravitationalParameter, + mu: GravitationalParameter, _marker: PhantomData<(C, F)>, } @@ -133,15 +132,15 @@ impl KeplerProblem { let true_anomaly = match el.conic_kind() { ConicRegime::Elliptic => { let n = (self.mu.value() / (a * a * a)).sqrt(); - let m = mean_from_true(el.true_anomaly, el.eccentricity) - + qtty::angular::Radians::new(n * dt.value()); + let m0 = mean_from_true(TrueAnomaly::new(el.true_anomaly), el.eccentricity); + let m = MeanAnomaly::from_value(m0.value() + n * dt.value()); true_from_mean(m, el.eccentricity, AnomalyOptions::default())? } ConicRegime::Hyperbolic => { let n = (-self.mu.value() / (a * a * a)).sqrt(); - let f0 = hyperbolic_from_true(el.true_anomaly, el.eccentricity); - let m = mean_from_hyperbolic(f0, el.eccentricity) - + qtty::angular::Radians::new(n * dt.value()); + let f0 = hyperbolic_from_true(TrueAnomaly::new(el.true_anomaly), el.eccentricity); + let m0 = mean_from_hyperbolic(f0, el.eccentricity); + let m = MeanAnomaly::from_value(m0.value() + n * dt.value()); true_from_hyperbolic( hyperbolic_from_mean(m, el.eccentricity, AnomalyOptions::default())?, el.eccentricity, @@ -155,9 +154,9 @@ impl KeplerProblem { el.inclination, el.raan, el.arg_periapsis, - true_anomaly, + true_anomaly.radians(), )?; - Ok(next.to_cartesian::(self.mu)) + Ok(next.try_to_cartesian::(self.mu)?) } } @@ -226,4 +225,38 @@ mod tests { < 1e-8 ); } + + #[test] + fn mu_accessor_returns_configured_mu() { + let mu = GravitationalParameter::new(398600.4418); + let p = KeplerProblem::::new(mu); + assert_eq!(p.mu().value(), 398600.4418); + } + + #[test] + fn hyperbolic_propagation_returns_finite_state() { + let mu = GravitationalParameter::new(398600.4418); + let r = 7000.0_f64; + let v_esc = (2.0 * mu.value() / r).sqrt(); + let state = CartesianState::::new( + Position::::new(r, 0.0, 0.0), + Velocity::::new(0.0, v_esc * 1.3, 0.0), + ); + let result = KeplerProblem::::new(mu).propagate(&state, Second::new(100.0)); + assert!(result.is_ok(), "hyperbolic propagation failed: {result:?}"); + } + + #[test] + fn parabolic_state_returns_error() { + let mu = GravitationalParameter::new(398600.4418); + let r = 7000.0_f64; + let v_par = (2.0 * mu.value() / r).sqrt(); + let state = CartesianState::::new( + Position::::new(r, 0.0, 0.0), + Velocity::::new(0.0, v_par, 0.0), + ); + assert!(KeplerProblem::::new(mu) + .propagate(&state, Second::new(100.0)) + .is_err()); + } } diff --git a/src/search.rs b/src/search.rs index 2aeb55b..30c31e8 100644 --- a/src/search.rs +++ b/src/search.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Caller-driven Lambert transfer-grid search. @@ -55,7 +55,7 @@ pub struct TransferCandidate { /// Arrival velocity from Lambert. pub v2: Velocity, /// Sum of Lambert endpoint speed magnitudes. - pub total_dv: KmPerSeconds, + pub endpoint_speed_sum: KmPerSeconds, } /// Outcome stored for each search grid cell. @@ -110,7 +110,7 @@ where flight_time: *tof, v1: sol.v1, v2: sol.v2, - total_dv: KmPerSeconds::new(total), + endpoint_speed_sum: KmPerSeconds::new(total), }) } }, @@ -188,4 +188,85 @@ mod tests { assert_eq!(out.cells[0].len(), 1); assert!(matches!(out.cells[0][0], CellOutcome::ProviderFailed(_))); } + + struct FixedProvider { + pos: [f64; 3], + } + impl TrajectoryProvider for FixedProvider { + type Error = &'static str; + fn position_at(&self, _: Second) -> Result, Self::Error> { + Ok(Position::::new( + self.pos[0], + self.pos[1], + self.pos[2], + )) + } + } + + #[test] + fn search_success_covers_speed_helper() { + let grid = SearchGrid { + departures: alloc::vec![Second::new(0.0)], + flight_times: alloc::vec![Second::new(4560.0)], + }; + let out = lambert_search( + &FixedProvider { + pos: [15945.34, 0.0, 0.0], + }, + &FixedProvider { + pos: [12214.84, 10249.47, 0.0], + }, + grid, + GravitationalParameter::new(398600.4418), + LambertBranch::Prograde, + ); + assert_eq!(out.cells.len(), 1); + assert!(matches!(out.cells[0][0], CellOutcome::Success(_))); + } + + struct TargetFails; + impl TrajectoryProvider for TargetFails { + type Error = &'static str; + fn position_at(&self, _: Second) -> Result, Self::Error> { + Err("target down") + } + } + + #[test] + fn search_target_provider_failure() { + let grid = SearchGrid { + departures: alloc::vec![Second::new(0.0)], + flight_times: alloc::vec![Second::new(4560.0)], + }; + let out = lambert_search( + &FixedProvider { + pos: [15945.34, 0.0, 0.0], + }, + &TargetFails, + grid, + GravitationalParameter::new(398600.4418), + LambertBranch::Prograde, + ); + assert!(matches!(out.cells[0][0], CellOutcome::ProviderFailed(_))); + } + + #[test] + fn search_lambert_failed_cell() { + let grid = SearchGrid { + departures: alloc::vec![Second::new(0.0)], + flight_times: alloc::vec![Second::new(4560.0)], + }; + let out = lambert_search( + &FixedProvider { + pos: [0.0, 0.0, 0.0], + }, + &FixedProvider { + pos: [0.0, 0.0, 0.0], + }, + grid, + GravitationalParameter::new(398600.4418), + LambertBranch::Prograde, + ); + assert!(matches!(out.cells[0][0], CellOutcome::LambertFailed(_))); + } } diff --git a/src/state.rs b/src/state.rs index a830b3c..b5c3649 100644 --- a/src/state.rs +++ b/src/state.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Typed Cartesian two-body state vectors. @@ -107,4 +107,25 @@ mod tests { assert_eq!(s.position().z().value(), 3.0); assert_eq!(s.velocity().x().value(), 4.0); } + + #[test] + fn clone_preserves_velocity() { + let s = CartesianState::::new( + Position::::new(7000.0, 0.0, 0.0), + Velocity::::new(0.0, 7.5, 0.0), + ); + #[allow(clippy::clone_on_copy)] + let s2 = Clone::clone(&s); + assert_eq!(s2.velocity().x().value(), 0.0); + assert_eq!(s2.velocity().y().value(), 7.5); + } + + #[test] + fn velocity_accessor_returns_velocity() { + let s = CartesianState::::new( + Position::::new(1.0, 2.0, 3.0), + Velocity::::new(4.0, 5.0, 6.0), + ); + assert_eq!(s.velocity().z().value(), 6.0); + } } diff --git a/src/transfer.rs b/src/transfer.rs index 1ddee32..7ce405b 100644 --- a/src/transfer.rs +++ b/src/transfer.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Analytic transfer and invariant helpers for two-body motion. @@ -29,6 +29,23 @@ use crate::problem::KeplerProblem; use crate::state::CartesianState; use crate::vec3::{cross, norm}; +/// Errors returned by fallible transfer helpers. +#[derive(Debug, Clone, Copy, PartialEq, thiserror::Error)] +pub enum TransferError { + /// Gravitational parameter is not strictly positive or not finite. + #[error("invalid gravitational parameter: {0}")] + InvalidGravitationalParameter(f64), + /// A radius argument is not strictly positive or not finite. + #[error("invalid {0} radius: {1}")] + InvalidRadius(&'static str, f64), + /// Semi-major axis is not strictly positive or not finite. + #[error("invalid semi-major axis: {0}")] + InvalidSemiMajorAxis(f64), + /// Computation produced a non-finite result. + #[error("non-finite result")] + NonFiniteResult, +} + /// Result of an ideal coplanar Hohmann transfer between circular orbits. #[derive(Debug, Clone, Copy, PartialEq)] pub struct HohmannResult { @@ -221,6 +238,172 @@ pub fn escape_speed(mu: GravitationalParameter, r: Kilometers) -> KmPerSeconds { KmPerSeconds::new((2.0 * mu.value() / r.value()).sqrt()) } +/// Fallible version of [`orbital_period`]. +/// +/// # Errors +/// +/// Returns [`TransferError`] if `mu ≤ 0`, `a ≤ 0`, or the result is not finite. +/// +/// # Examples +/// +/// ``` +/// use keplerian::problem::KeplerProblem; +/// use keplerian::transfer::try_orbital_period; +/// use qtty::dynamics::GravitationalParameter; +/// use qtty::length::Kilometers; +/// +/// #[derive(Debug, Clone, Copy)] struct C; impl affn::centers::ReferenceCenter for C { type Params = (); fn center_name() -> &'static str { "C" } } +/// #[derive(Debug, Clone, Copy)] struct F; impl affn::frames::ReferenceFrame for F { fn frame_name() -> &'static str { "F" } } +/// +/// let p = KeplerProblem::::new(GravitationalParameter::new(398600.4418)); +/// let t = try_orbital_period(&p, Kilometers::new(7000.0)).unwrap(); +/// assert!((t.value() - 5840.0).abs() < 200.0); +/// ``` +pub fn try_orbital_period( + problem: &KeplerProblem, + semi_major_axis: Kilometers, +) -> Result { + let mu = problem.mu().value(); + if !mu.is_finite() || mu <= 0.0 { + return Err(TransferError::InvalidGravitationalParameter(mu)); + } + let a = semi_major_axis.value(); + if !a.is_finite() || a <= 0.0 { + return Err(TransferError::InvalidSemiMajorAxis(a)); + } + let t = 2.0 * core::f64::consts::PI * (a * a * a / mu).sqrt(); + if !t.is_finite() { + return Err(TransferError::NonFiniteResult); + } + Ok(Second::new(t)) +} + +/// Fallible version of [`hohmann_delta_v`]. +/// +/// # Errors +/// +/// Returns [`TransferError`] if inputs are non-positive or the result is not finite. +/// +/// # Examples +/// +/// ``` +/// use keplerian::transfer::try_hohmann_delta_v; +/// use qtty::dynamics::GravitationalParameter; +/// use qtty::length::Kilometers; +/// +/// let h = try_hohmann_delta_v( +/// GravitationalParameter::new(398600.4418), +/// Kilometers::new(6678.0), +/// Kilometers::new(42164.0), +/// ).unwrap(); +/// assert!((h.total.value() - 3.91).abs() < 0.1); +/// ``` +pub fn try_hohmann_delta_v( + mu: GravitationalParameter, + r1: Kilometers, + r2: Kilometers, +) -> Result { + let mu_val = mu.value(); + if !mu_val.is_finite() || mu_val <= 0.0 { + return Err(TransferError::InvalidGravitationalParameter(mu_val)); + } + let r1_val = r1.value(); + if !r1_val.is_finite() || r1_val <= 0.0 { + return Err(TransferError::InvalidRadius("r1", r1_val)); + } + let r2_val = r2.value(); + if !r2_val.is_finite() || r2_val <= 0.0 { + return Err(TransferError::InvalidRadius("r2", r2_val)); + } + let result = hohmann_delta_v(mu, r1, r2); + if !result.total.value().is_finite() { + return Err(TransferError::NonFiniteResult); + } + Ok(result) +} + +/// Fallible version of [`vis_viva_speed`]. +/// +/// # Errors +/// +/// Returns [`TransferError`] if inputs are non-positive or the result is not finite. +/// +/// # Examples +/// +/// ``` +/// use keplerian::transfer::try_vis_viva_speed; +/// use qtty::dynamics::GravitationalParameter; +/// use qtty::length::Kilometers; +/// +/// let v = try_vis_viva_speed( +/// GravitationalParameter::new(398600.4418), +/// Kilometers::new(7000.0), +/// Kilometers::new(7000.0), +/// ).unwrap(); +/// assert!((v.value() - 7.546).abs() < 0.01); +/// ``` +pub fn try_vis_viva_speed( + mu: GravitationalParameter, + r: Kilometers, + a: Kilometers, +) -> Result { + let mu_val = mu.value(); + if !mu_val.is_finite() || mu_val <= 0.0 { + return Err(TransferError::InvalidGravitationalParameter(mu_val)); + } + let r_val = r.value(); + if !r_val.is_finite() || r_val <= 0.0 { + return Err(TransferError::InvalidRadius("r", r_val)); + } + let a_val = a.value(); + if !a_val.is_finite() || a_val == 0.0 { + return Err(TransferError::InvalidSemiMajorAxis(a_val)); + } + let v = vis_viva_speed(mu, r, a); + if !v.value().is_finite() { + return Err(TransferError::NonFiniteResult); + } + Ok(v) +} + +/// Fallible version of [`escape_speed`]. +/// +/// # Errors +/// +/// Returns [`TransferError`] if inputs are non-positive or the result is not finite. +/// +/// # Examples +/// +/// ``` +/// use keplerian::transfer::try_escape_speed; +/// use qtty::dynamics::GravitationalParameter; +/// use qtty::length::Kilometers; +/// +/// let v = try_escape_speed( +/// GravitationalParameter::new(398600.4418), +/// Kilometers::new(6378.0), +/// ).unwrap(); +/// assert!((v.value() - 11.18).abs() < 0.1); +/// ``` +pub fn try_escape_speed( + mu: GravitationalParameter, + r: Kilometers, +) -> Result { + let mu_val = mu.value(); + if !mu_val.is_finite() || mu_val <= 0.0 { + return Err(TransferError::InvalidGravitationalParameter(mu_val)); + } + let r_val = r.value(); + if !r_val.is_finite() || r_val <= 0.0 { + return Err(TransferError::InvalidRadius("r", r_val)); + } + let v = escape_speed(mu, r); + if !v.value().is_finite() { + return Err(TransferError::NonFiniteResult); + } + Ok(v) +} + #[cfg(test)] mod tests { use super::*; @@ -268,4 +451,125 @@ mod tests { ); assert!((specific_angular_momentum(&s).value() - 52500.0).abs() < 1e-12); } + + #[test] + fn specific_orbital_energy_and_period() { + let mu = GravitationalParameter::new(398600.4418); + let problem = KeplerProblem::::new(mu); + let state = CartesianState::::new( + Position::::new(7000.0, 0.0, 0.0), + Velocity::::new(0.0, 7.546, 0.0), + ); + + assert!(specific_orbital_energy(&state, mu).value() < 0.0); + assert!(specific_angular_momentum(&state).value() > 0.0); + + let t = orbital_period(&problem, Kilometers::new(7000.0)).unwrap(); + assert!((t.value() - 5840.0).abs() < 200.0); + } + + #[test] + fn orbital_period_negative_sma_returns_none() { + let problem = KeplerProblem::::new(GravitationalParameter::new(398600.4418)); + assert!(orbital_period(&problem, Kilometers::new(-7000.0)).is_none()); + } + + #[test] + fn try_orbital_period_validates_inputs() { + let problem = KeplerProblem::::new(GravitationalParameter::new(398600.4418)); + let t = try_orbital_period(&problem, Kilometers::new(7000.0)).unwrap(); + assert!((t.value() - 5840.0).abs() < 200.0); + assert!(matches!( + try_orbital_period(&problem, Kilometers::new(0.0)), + Err(TransferError::InvalidSemiMajorAxis(_)) + )); + assert!(matches!( + try_orbital_period(&problem, Kilometers::new(-7000.0)), + Err(TransferError::InvalidSemiMajorAxis(_)) + )); + } + + #[test] + fn try_hohmann_delta_v_validates_inputs() { + let mu = GravitationalParameter::new(398600.4418); + let result = + try_hohmann_delta_v(mu, Kilometers::new(6678.0), Kilometers::new(42164.0)).unwrap(); + assert!((result.total.value() - 3.91).abs() < 0.3); + + assert!(matches!( + try_hohmann_delta_v( + GravitationalParameter::new(-1.0), + Kilometers::new(6678.0), + Kilometers::new(42164.0) + ), + Err(TransferError::InvalidGravitationalParameter(_)) + )); + assert!(matches!( + try_hohmann_delta_v( + GravitationalParameter::new(0.0), + Kilometers::new(6678.0), + Kilometers::new(42164.0) + ), + Err(TransferError::InvalidGravitationalParameter(_)) + )); + assert!(matches!( + try_hohmann_delta_v(mu, Kilometers::new(0.0), Kilometers::new(42164.0)), + Err(TransferError::InvalidRadius(_, _)) + )); + assert!(matches!( + try_hohmann_delta_v(mu, Kilometers::new(-6678.0), Kilometers::new(42164.0)), + Err(TransferError::InvalidRadius(_, _)) + )); + assert!(matches!( + try_hohmann_delta_v(mu, Kilometers::new(6678.0), Kilometers::new(-1.0)), + Err(TransferError::InvalidRadius(_, _)) + )); + } + + #[test] + fn try_vis_viva_speed_validates_inputs() { + let mu = GravitationalParameter::new(398600.4418); + let v = try_vis_viva_speed(mu, Kilometers::new(7000.0), Kilometers::new(7000.0)).unwrap(); + assert!((v.value() - 7.546).abs() < 0.05); + + assert!(matches!( + try_vis_viva_speed( + GravitationalParameter::new(-1.0), + Kilometers::new(7000.0), + Kilometers::new(7000.0) + ), + Err(TransferError::InvalidGravitationalParameter(_)) + )); + assert!(matches!( + try_vis_viva_speed(mu, Kilometers::new(0.0), Kilometers::new(7000.0)), + Err(TransferError::InvalidRadius(_, _)) + )); + assert!(matches!( + try_vis_viva_speed(mu, Kilometers::new(7000.0), Kilometers::new(0.0)), + Err(TransferError::InvalidSemiMajorAxis(_)) + )); + } + + #[test] + fn try_escape_speed_validates_inputs() { + let mu = GravitationalParameter::new(398600.4418); + let v = try_escape_speed(mu, Kilometers::new(6378.0)).unwrap(); + assert!((v.value() - 11.18).abs() < 0.1); + + assert!(matches!( + try_escape_speed(GravitationalParameter::new(0.0), Kilometers::new(6378.0)), + Err(TransferError::InvalidGravitationalParameter(_)) + )); + assert!(matches!( + try_escape_speed(mu, Kilometers::new(-1.0)), + Err(TransferError::InvalidRadius(_, _)) + )); + assert!(matches!( + try_escape_speed( + GravitationalParameter::new(f64::NAN), + Kilometers::new(6378.0) + ), + Err(TransferError::InvalidGravitationalParameter(_)) + )); + } } diff --git a/src/vec3.rs b/src/vec3.rs index 9b90222..e2d288d 100644 --- a/src/vec3.rs +++ b/src/vec3.rs @@ -1,4 +1,4 @@ -// SPDX-License-Identifier: AGPL-3.0-or-later +// SPDX-License-Identifier: AGPL-3.0-only // Copyright (C) 2026 Vallés Puig, Ramon //! Private numeric helpers for `[f64; 3]` vector algebra. diff --git a/tests/anomaly_property.rs b/tests/anomaly_property.rs index 78f8cf2..6570647 100644 --- a/tests/anomaly_property.rs +++ b/tests/anomaly_property.rs @@ -1,15 +1,14 @@ //! Property tests for typed anomaly conversions. -use keplerian::anomaly::{eccentric_from_mean, mean_from_eccentric, AnomalyOptions}; +use keplerian::anomaly::{eccentric_from_mean, mean_from_eccentric, AnomalyOptions, MeanAnomaly}; use keplerian::Eccentricity; use proptest::prelude::*; -use qtty::angular::Radians; proptest! { #[test] fn elliptic_mean_round_trips(m in 0.0_f64..core::f64::consts::TAU, e in 0.0_f64..0.9) { let ecc = Eccentricity::new(e).unwrap(); - let ea = eccentric_from_mean(Radians::new(m), ecc, AnomalyOptions::default()).unwrap(); + let ea = eccentric_from_mean(MeanAnomaly::from_value(m), ecc, AnomalyOptions::default()).unwrap(); let m2 = mean_from_eccentric(ea, ecc); prop_assert!((m2.value() - m).abs() < 1e-10); } diff --git a/tests/coverage_extras.rs b/tests/coverage_extras.rs deleted file mode 100644 index 3616755..0000000 --- a/tests/coverage_extras.rs +++ /dev/null @@ -1,547 +0,0 @@ -//! Tests targeting coverage gaps: error conversions, state clone, transfer -//! helpers, conic kinds, anomaly edge cases, Lambert n-rev, and search grid. - -use affn::cartesian::{Position, Velocity}; -use affn::centers::ReferenceCenter; -use affn::frames::ReferenceFrame; -use affn::frames::ICRS; -use keplerian::anomaly::{ - eccentric_from_mean, hyperbolic_from_mean, hyperbolic_from_true, kepler_elliptic, - kepler_hyperbolic, kepler_parabolic, mean_from_hyperbolic, true_from_hyperbolic, - AnomalyOptions, -}; -use keplerian::eccentricity::Eccentricity; -use keplerian::elements::{ConversionError, KeplerianElements}; -use keplerian::error::KeplerError; -use keplerian::lambert::{lambert_n_rev, LambertBranch, LambertError, NRevBranch}; -use keplerian::problem::{KeplerProblem, PropagationError}; -use keplerian::state::CartesianState; -use keplerian::transfer::{orbital_period, specific_angular_momentum, specific_orbital_energy}; -use qtty::angular::Radians; -use qtty::dynamics::GravitationalParameter; -use qtty::length::{Kilometer, Kilometers}; -use qtty::Second; - -#[derive(Debug, Clone, Copy)] -struct C; -impl ReferenceCenter for C { - type Params = (); - fn center_name() -> &'static str { - "C" - } -} - -#[derive(Debug, Clone, Copy)] -struct F; -impl ReferenceFrame for F { - fn frame_name() -> &'static str { - "F" - } -} - -// ── error.rs ────────────────────────────────────────────────────────────────── - -#[test] -fn kepler_error_from_anomaly() { - let inner = keplerian::anomaly::AnomalyError::InvalidEccentricity(1.5); - let e = KeplerError::from(inner); - assert!(matches!(e, KeplerError::Anomaly(_))); -} - -#[test] -fn kepler_error_from_conversion() { - let inner = ConversionError::InvalidEccentricity(-1.0); - let e = KeplerError::from(inner); - assert!(matches!(e, KeplerError::Conversion(_))); -} - -#[test] -fn kepler_error_from_propagation() { - let inner = PropagationError::ParabolicUnsupported; - let e = KeplerError::from(inner); - assert!(matches!(e, KeplerError::Propagation(_))); -} - -#[test] -fn kepler_error_from_lambert() { - let inner = LambertError::ZeroPosition; - let e = KeplerError::from(inner); - assert!(matches!(e, KeplerError::Lambert(_))); -} - -// ── state.rs ────────────────────────────────────────────────────────────────── - -#[test] -fn cartesian_state_clone_and_velocity() { - let pos = Position::::new(7000.0, 0.0, 0.0); - let vel = Velocity::::new(0.0, 7.5, 0.0); - let s = CartesianState::::new(pos, vel); - // CartesianState is Copy when Position and Velocity are Copy; - // use explicit Clone::clone to exercise the Clone impl body. - #[allow(clippy::clone_on_copy)] - let s2 = Clone::clone(&s); - assert_eq!(s2.velocity().x().value(), 0.0); - assert_eq!(s2.velocity().y().value(), 7.5); -} - -// ── transfer.rs ─────────────────────────────────────────────────────────────── - -#[test] -fn specific_orbital_energy_and_period() { - let mu = GravitationalParameter::new(398600.4418); - let problem = KeplerProblem::::new(mu); - - let pos = Position::::new(7000.0, 0.0, 0.0); - let vel = Velocity::::new(0.0, 7.546, 0.0); - let state = CartesianState::::new(pos, vel); - - let eps = specific_orbital_energy(&state, mu); - assert!(eps.value() < 0.0, "bound orbit has negative energy"); - - let h = specific_angular_momentum(&state); - assert!(h.value() > 0.0); - - let t = orbital_period(&problem, Kilometers::new(7000.0)).unwrap(); - assert!((t.value() - 5840.0).abs() < 200.0); -} - -#[test] -fn orbital_period_negative_sma_returns_none() { - let mu = GravitationalParameter::new(398600.4418); - let problem = KeplerProblem::::new(mu); - assert!(orbital_period(&problem, Kilometers::new(-7000.0)).is_none()); -} - -// ── problem.rs ──────────────────────────────────────────────────────────────── - -#[test] -fn kepler_problem_mu_accessor() { - let mu = GravitationalParameter::new(398600.4418); - let p = KeplerProblem::::new(mu); - assert_eq!(p.mu().value(), 398600.4418); -} - -#[test] -fn kepler_problem_hyperbolic_propagation() { - // Build a clearly hyperbolic orbit (e ≈ 1.5) by giving v > escape speed. - let mu = GravitationalParameter::new(398600.4418); - let r = 7000.0_f64; - let v_esc = (2.0 * mu.value() / r).sqrt(); - let pos = Position::::new(r, 0.0, 0.0); - let vel = Velocity::::new(0.0, v_esc * 1.3, 0.0); - let state = CartesianState::::new(pos, vel); - let problem = KeplerProblem::::new(mu); - let result = problem.propagate(&state, Second::new(100.0)); - assert!(result.is_ok(), "hyperbolic propagation failed: {result:?}"); -} - -#[test] -fn kepler_problem_parabolic_returns_error() { - // Near-parabolic: energy ≈ 0 → from_cartesian returns Degenerate("parabolic orbit") - // which propagate() wraps as PropagationError::Conversion. - let mu = GravitationalParameter::new(398600.4418); - let r = 7000.0_f64; - let v_par = (2.0 * mu.value() / r).sqrt(); // exact escape speed - let pos = Position::::new(r, 0.0, 0.0); - let vel = Velocity::::new(0.0, v_par, 0.0); - let state = CartesianState::::new(pos, vel); - let problem = KeplerProblem::::new(mu); - assert!(problem.propagate(&state, Second::new(100.0)).is_err()); -} - -// ── elements.rs ─────────────────────────────────────────────────────────────── - -#[test] -fn elements_new_rejects_negative_eccentricity() { - let err = KeplerianElements::::new( - Kilometers::new(7000.0), - Eccentricity::new_unchecked(-0.1), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - ); - assert!(matches!(err, Err(ConversionError::InvalidEccentricity(_)))); -} - -#[test] -fn elements_new_rejects_inclination_out_of_range() { - use core::f64::consts::PI; - let err = KeplerianElements::::new( - Kilometers::new(7000.0), - Eccentricity::new_unchecked(0.1), - Radians::new(PI + 0.1), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - ); - assert!(matches!(err, Err(ConversionError::InvalidInclination(_)))); -} - -#[test] -fn conic_kind_parabolic() { - // Eccentricity::new_unchecked allows e = 1.0 which is parabolic. - use keplerian::elements::ConicRegime; - let el = KeplerianElements::::new( - Kilometers::new(7000.0), - Eccentricity::new_unchecked(1.0), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - ) - .unwrap(); - assert_eq!(el.conic_kind(), ConicRegime::Parabolic); -} - -#[test] -fn conic_kind_hyperbolic() { - use keplerian::elements::ConicRegime; - let el = KeplerianElements::::new( - Kilometers::new(-40000.0), - Eccentricity::new_unchecked(1.5), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - Radians::new(0.0), - ) - .unwrap(); - assert_eq!(el.conic_kind(), ConicRegime::Hyperbolic); -} - -#[test] -fn from_cartesian_degenerate_zero_position() { - let mu = GravitationalParameter::new(398600.4418); - let pos = Position::::new(0.0, 0.0, 0.0); - let vel = Velocity::::new(0.0, 7.5, 0.0); - let state = CartesianState::::new(pos, vel); - assert!(KeplerianElements::::from_cartesian(&state, mu).is_err()); -} - -#[test] -fn from_cartesian_degenerate_zero_angular_momentum() { - let mu = GravitationalParameter::new(398600.4418); - let pos = Position::::new(7000.0, 0.0, 0.0); - // radial velocity → h = r × v = 0 - let vel = Velocity::::new(7.5, 0.0, 0.0); - let state = CartesianState::::new(pos, vel); - assert!(KeplerianElements::::from_cartesian(&state, mu).is_err()); -} - -#[test] -fn from_cartesian_equatorial_eccentric_orbit() { - // Equatorial (inc ≈ 0) non-circular orbit: tests argp from equatorial branch. - let mu = GravitationalParameter::new(398600.4418); - let pos = Position::::new(7000.0, 0.0, 0.0); - // Slightly off-circular, equatorial: vy < v_circ, vz = 0. - let vel = Velocity::::new(0.0, 6.5, 0.0); - let state = CartesianState::::new(pos, vel); - let el = KeplerianElements::::from_cartesian(&state, mu); - assert!(el.is_ok(), "{el:?}"); -} - -// ── anomaly.rs ──────────────────────────────────────────────────────────────── - -#[test] -fn elliptic_circular_orbit_returns_mean_anomaly() { - // e = 0 → E = M (identity) - let zero = Eccentricity::new(0.0).unwrap(); - let m = Radians::new(1.23); - let e = eccentric_from_mean(m, zero, AnomalyOptions::default()).unwrap(); - assert!((e.value() - m.value()).abs() < 1e-14); -} - -#[test] -fn hyperbolic_mean_anomaly_zero_returns_zero() { - let ecc = Eccentricity::new_unchecked(1.5); - let f = hyperbolic_from_mean(Radians::new(0.0), ecc, AnomalyOptions::default()).unwrap(); - assert!((f).abs() < 1e-14); -} - -#[test] -fn hyperbolic_mean_anomaly_large_branch() { - // |M| > 50 * e triggers the log initial guess branch. - let ecc = Eccentricity::new_unchecked(1.2); - let m = Radians::new(100.0); - let f = hyperbolic_from_mean(m, ecc, AnomalyOptions::default()); - assert!(f.is_ok()); -} - -#[test] -fn kepler_parabolic_round_trips() { - let m = 0.5_f64; - let d = kepler_parabolic(m); - let m_back = d + d.powi(3) / 3.0; - assert!((m_back - m).abs() < 1e-12); -} - -#[test] -fn hyperbolic_anomaly_round_trips() { - let ecc = Eccentricity::new_unchecked(2.0); - let nu = Radians::new(0.8); - let f = hyperbolic_from_true(nu, ecc); - let nu2 = true_from_hyperbolic(f, ecc); - assert!((nu2.value() - nu.value()).abs() < 1e-12); - let m = mean_from_hyperbolic(f, ecc); - let f2 = hyperbolic_from_mean(m, ecc, AnomalyOptions::default()).unwrap(); - assert!((f2 - f).abs() < 1e-12); -} - -#[test] -fn hyperbolic_kepler_rejects_nan_mean_anomaly() { - let ecc = Eccentricity::new_unchecked(1.5); - let err = hyperbolic_from_mean(Radians::new(f64::NAN), ecc, AnomalyOptions::default()); - assert!(err.is_err()); -} - -#[test] -fn elliptic_kepler_rejects_nan_mean_anomaly() { - let ecc = Eccentricity::new_unchecked(0.5); - let err = eccentric_from_mean(Radians::new(f64::NAN), ecc, AnomalyOptions::default()); - assert!(err.is_err()); -} - -// ── lambert/typed.rs ────────────────────────────────────────────────────────── - -#[test] -fn lambert_n_rev_valid_case() { - // Two positions separated by ~90°; 2-hour TOF, 1 extra revolution. - let r1 = Position::<(), ICRS, Kilometer>::new(7000.0, 0.0, 0.0); - let r2 = Position::<(), ICRS, Kilometer>::new(0.0, 7000.0, 0.0); - let tof = Second::new(10800.0); - let mu = GravitationalParameter::new(398600.4418); - // May succeed or fail depending on geometry; just exercise the code path. - let _ = lambert_n_rev( - r1, - r2, - tof, - mu, - LambertBranch::Prograde, - 1, - NRevBranch::Left, - ); -} - -#[test] -fn lambert_n_rev_retrograde_case() { - let r1 = Position::<(), ICRS, Kilometer>::new(7000.0, 0.0, 0.0); - let r2 = Position::<(), ICRS, Kilometer>::new(0.0, 7000.0, 0.0); - let tof = Second::new(10800.0); - let mu = GravitationalParameter::new(398600.4418); - let _ = lambert_n_rev( - r1, - r2, - tof, - mu, - LambertBranch::Retrograde, - 1, - NRevBranch::Right, - ); -} - -// ── search.rs ───────────────────────────────────────────────────────────────── - -#[cfg(feature = "alloc")] -mod search_tests { - extern crate alloc; - - use super::*; - use keplerian::search::{lambert_search, CellOutcome, SearchGrid, TrajectoryProvider}; - - struct FixedProvider { - pos: [f64; 3], - } - impl TrajectoryProvider for FixedProvider { - type Error = &'static str; - fn position_at(&self, _: Second) -> Result, Self::Error> { - Ok(Position::::new( - self.pos[0], - self.pos[1], - self.pos[2], - )) - } - } - - #[test] - fn search_success_covers_speed_helper() { - let grid = SearchGrid { - departures: alloc::vec![Second::new(0.0)], - flight_times: alloc::vec![Second::new(4560.0)], - }; - let out = lambert_search( - &FixedProvider { - pos: [15945.34, 0.0, 0.0], - }, - &FixedProvider { - pos: [12214.84, 10249.47, 0.0], - }, - grid, - GravitationalParameter::new(398600.4418), - LambertBranch::Prograde, - ); - assert_eq!(out.cells.len(), 1); - assert!(matches!(out.cells[0][0], CellOutcome::Success(_))); - } - - struct TargetFails; - impl TrajectoryProvider for TargetFails { - type Error = &'static str; - fn position_at(&self, _: Second) -> Result, Self::Error> { - Err("target down") - } - } - - #[test] - fn search_target_provider_failure() { - let grid = SearchGrid { - departures: alloc::vec![Second::new(0.0)], - flight_times: alloc::vec![Second::new(4560.0)], - }; - let out = lambert_search( - &FixedProvider { - pos: [15945.34, 0.0, 0.0], - }, - &TargetFails, - grid, - GravitationalParameter::new(398600.4418), - LambertBranch::Prograde, - ); - assert!(matches!(out.cells[0][0], CellOutcome::ProviderFailed(_))); - } - - /// `r1 == r2` forces a `ZeroPosition` Lambert failure → covers `CellOutcome::LambertFailed`. - #[test] - fn search_lambert_failed_cell() { - let grid = SearchGrid { - departures: alloc::vec![Second::new(0.0)], - flight_times: alloc::vec![Second::new(4560.0)], - }; - let out = lambert_search( - &FixedProvider { - pos: [0.0, 0.0, 0.0], - }, - &FixedProvider { - pos: [0.0, 0.0, 0.0], - }, - grid, - GravitationalParameter::new(398600.4418), - LambertBranch::Prograde, - ); - assert!(matches!(out.cells[0][0], CellOutcome::LambertFailed(_))); - } -} - -// ── eccentricity.rs ─────────────────────────────────────────────────────────── - -#[test] -fn eccentricity_is_hyperbolic() { - assert!(Eccentricity::new_unchecked(1.5).is_hyperbolic()); - assert!(!Eccentricity::new_unchecked(0.5).is_hyperbolic()); -} - -// ── state.rs: velocity() return value ───────────────────────────────────────── - -#[test] -fn cartesian_state_velocity_accessor() { - let pos = Position::::new(1.0, 2.0, 3.0); - let vel = Velocity::::new(4.0, 5.0, 6.0); - let s = CartesianState::::new(pos, vel); - let vref = s.velocity(); - assert_eq!(vref.z().value(), 6.0); -} - -// ── elements.rs: non-positive mu and NaN ───────────────────────────────────── - -#[test] -fn from_cartesian_non_positive_mu() { - let pos = Position::::new(7000.0, 0.0, 0.0); - let vel = Velocity::::new(0.0, 7.5, 0.0); - let state = CartesianState::::new(pos, vel); - let mu = GravitationalParameter::new(-1.0); - let err = KeplerianElements::::from_cartesian(&state, mu); - assert!(err.is_err()); -} - -#[test] -fn from_cartesian_nan_mu() { - let pos = Position::::new(7000.0, 0.0, 0.0); - let vel = Velocity::::new(0.0, 7.5, 0.0); - let state = CartesianState::::new(pos, vel); - let mu = GravitationalParameter::new(f64::NAN); - let err = KeplerianElements::::from_cartesian(&state, mu); - assert!(matches!(err, Err(ConversionError::NonFiniteValue { .. }))); -} - -#[test] -fn from_cartesian_circular_inclined_orbit() { - // Inclined circular orbit: ecc ≈ 0 and nmag > 0 → exercises the - // `else if nmag > EPS` branch for nu computation. - let mu = GravitationalParameter::new(398600.4418); - let r = 7000.0_f64; - let v_circ = (mu.value() / r).sqrt(); - // Tilt into the i=45° plane: velocity in y-z plane. - let vy = v_circ * std::f64::consts::FRAC_1_SQRT_2; - let vz = v_circ * std::f64::consts::FRAC_1_SQRT_2; - let pos = Position::::new(r, 0.0, 0.0); - let vel = Velocity::::new(0.0, vy, vz); - let state = CartesianState::::new(pos, vel); - let el = KeplerianElements::::from_cartesian(&state, mu); - assert!(el.is_ok(), "{el:?}"); -} - -#[test] -fn elements_to_cartesian_and_back() { - // Explicitly calls to_cartesian, covering its rotate_pqw invocations. - let mu = GravitationalParameter::new(398600.4418); - let el = KeplerianElements::::new( - Kilometers::new(7000.0), - Eccentricity::new_unchecked(0.1), - Radians::new(0.5), - Radians::new(1.0), - Radians::new(0.3), - Radians::new(0.7), - ) - .unwrap(); - let state = el.to_cartesian::(mu); - assert!(state.position().x().value().is_finite()); -} - -// ── anomaly.rs: bisection fallback paths ───────────────────────────────────── - -#[test] -fn elliptic_bisection_fallback_with_zero_tol() { - // max_iter = 2, tol = 0.0: Newton fails to reach exact zero residual, - // triggers bisection, exercises lines 110-112 and 421-429. - let ecc = Eccentricity::new_unchecked(0.5); - let opts = AnomalyOptions { - max_iter: 2, - tol: 0.0, - }; - let result = kepler_elliptic(Radians::new(1.0), ecc, opts); - // Either converges or not, but the code paths are exercised. - let _ = result; -} - -#[test] -fn hyperbolic_bisection_fallback_with_zero_tol() { - // max_iter = 2, tol = 0.0: exercises hyperbolic bisection lines. - let ecc = Eccentricity::new_unchecked(1.5); - let opts = AnomalyOptions { - max_iter: 2, - tol: 0.0, - }; - let result = kepler_hyperbolic(Radians::new(1.0), ecc, opts); - let _ = result; -} - -#[test] -fn hyperbolic_bisection_large_mean_anomaly() { - // Large |M| triggers the upper-bracket expansion loop (lines 459-460). - let ecc = Eccentricity::new_unchecked(1.1); - let opts = AnomalyOptions { - max_iter: 2, - tol: 0.0, - }; - let result = kepler_hyperbolic(Radians::new(50.0), ecc, opts); - let _ = result; -}