diff --git a/Cargo.toml b/Cargo.toml index a1673a7..8fb287f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -84,8 +84,8 @@ name = "fit_full_iod_parallel" path = "examples/run_full_iod_parallel.rs" required-features = ["parallel"] -#[profile.test] -#debug = false +[profile.dev] +debug = false [profile.bench] opt-level = 3 diff --git a/src/differential_orbit_correction/diff_cor.rs b/src/differential_orbit_correction/diff_cor.rs index 98c2afe..688604a 100644 --- a/src/differential_orbit_correction/diff_cor.rs +++ b/src/differential_orbit_correction/diff_cor.rs @@ -55,9 +55,12 @@ use crate::{ outlier_rejection::{update_observation_selection, OutlierRejectionConfig}, single_iteration::single_iteration, }, - orbit_type::equinoctial_element::EquinoctialLimits, + orbit_type::{ + equinoctial_element::EquinoctialLimits, + uncertainty::{EquinoctialUncertainty, OrbitalCovariance}, + }, propagator::PropagatorKind, - EquinoctialElements, JPLEphem, OutfitError, + EquinoctialElements, JPLEphem, OrbitalElements, OutfitError, }; use super::least_square::OrbitalUncertainty; @@ -221,6 +224,24 @@ pub struct DifferentialCorrectionOutput { pub num_measurements: usize, } +/// Conversion from [`DifferentialCorrectionOutput`] to the more general +/// [`OrbitalElements`] type. +/// +/// The covariance is extracted from the `uncertainty` field and included in the +/// `OrbitalElements::Equinoctial` variant. +impl From for OrbitalElements { + fn from(output: DifferentialCorrectionOutput) -> Self { + let orb_covariance = OrbitalCovariance { + matrix: output.uncertainty.covariance, + }; + OrbitalElements::Equinoctial { + elements: output.elements, + uncertainty: Some(EquinoctialUncertainty::from_covariance(&orb_covariance)), + covariance: Some(orb_covariance), + } + } +} + // ───────────────────────────────────────────────────────────────────────────── // Main function // ───────────────────────────────────────────────────────────────────────────── diff --git a/src/differential_orbit_correction/obs_dataset_api.rs b/src/differential_orbit_correction/obs_dataset_api.rs index d095edd..b24c199 100644 --- a/src/differential_orbit_correction/obs_dataset_api.rs +++ b/src/differential_orbit_correction/obs_dataset_api.rs @@ -26,7 +26,7 @@ use crate::{ obs_fit_data::ObsFitData, }, initial_orbit_determination::{obs_dataset_api::run_iod_on_observations, IODParams}, - FullOrbitResult, JPLEphem, OrbitalElements, OutfitError, + FullOrbitResult, JPLEphem, OutfitError, }; use hifitime::ut1::Ut1Provider; use photom::{ @@ -210,7 +210,13 @@ fn run_differential_correction_for_trajectory( let initial_equinoctial = match initial_orbits.and_then(|map| map.get(traj_id)) { Some(Ok(orbital_elements)) => { // Caller provided an initial orbit — convert to equinoctial. - orbital_elements.orbital_elements().to_equinoctial()? + orbital_elements + .orbital_elements() + .to_equinoctial()? + .as_equinoctial() + .ok_or(OutfitError::InvalidConversion( + "Conversion to equinoctial elements failed".to_string(), + ))? } Some(Err(_)) | None => { // No initial orbit — run IOD directly on the already-corrected @@ -218,7 +224,13 @@ fn run_differential_correction_for_trajectory( // dataset. This avoids reconstructing an ObsDataset and // rebuilding the cache. let iod_result = run_iod_on_observations(&observations, cache, jpl, iod_params, rng)?; - iod_result.orbital_elements().to_equinoctial()? + iod_result + .orbital_elements() + .to_equinoctial()? + .as_equinoctial() + .ok_or(OutfitError::InvalidConversion( + "Conversion to equinoctial elements failed".to_string(), + ))? } }; @@ -239,9 +251,9 @@ fn run_differential_correction_for_trajectory( )?; // ── Package the result ──────────────────────────────────────────────────── - let orbital_elements = OrbitalElements::Equinoctial(dc_output.elements); + let normalised_rms = dc_output.normalised_rms; Ok(FitOrbitResult::DifferentialCorrection(( - orbital_elements, - dc_output.normalised_rms, + dc_output.into(), + normalised_rms, ))) } diff --git a/src/ephemeris/mod.rs b/src/ephemeris/mod.rs index e5f7491..306b122 100644 --- a/src/ephemeris/mod.rs +++ b/src/ephemeris/mod.rs @@ -111,8 +111,8 @@ use hifitime::ut1::Ut1Provider; use crate::{ cache::observer_fixed_cache::ObserverFixedCache, - ephemeris::observation_ephemeris::check_elliptical_orbit, propagator::PropagatorKind, JPLEphem, - OrbitalElements, OutfitError, + ephemeris::observation_ephemeris::check_elliptical_orbit, propagator::PropagatorKind, + EquinoctialElements, JPLEphem, OrbitalElements, OutfitError, }; // --------------------------------------------------------------------------- @@ -291,7 +291,11 @@ impl OrbitalElements { /// Convert `self` to [`crate::EquinoctialElements`] for use in ephemeris /// computation. - fn to_equinoctial_for_ephemeris(&self) -> Result { - self.to_equinoctial() + fn to_equinoctial_for_ephemeris(&self) -> Result { + self.to_equinoctial()? + .as_equinoctial() + .ok_or(OutfitError::InvalidConversion( + "Conversion to equinoctial elements failed".to_string(), + )) } } diff --git a/src/initial_orbit_determination/gauss.rs b/src/initial_orbit_determination/gauss.rs index b3c4c49..0056ffa 100644 --- a/src/initial_orbit_determination/gauss.rs +++ b/src/initial_orbit_determination/gauss.rs @@ -76,13 +76,13 @@ //! GaussResult::PrelimOrbit(oe) //! | GaussResult::CorrectedOrbit(oe) => { //! match oe { -//! OrbitalElements::Keplerian(kepl) => { +//! OrbitalElements::Keplerian { elements: kepl, .. } => { //! println!("a [AU] = {}", kepl.semi_major_axis); //! } -//! OrbitalElements::Equinoctial(eq) => { +//! OrbitalElements::Equinoctial { elements: eq, .. } => { //! println!("lambda [rad] = {}", eq.mean_longitude); //! } -//! OrbitalElements::Cometary(com) => { +//! OrbitalElements::Cometary { elements: com, .. } => { //! println!("q [AU] = {}", com.perihelion_distance); //! } //! } @@ -1737,15 +1737,19 @@ pub(crate) mod gauss_test { // The values are very close to the ones obtained from the Rust implementation // The floating point differences is very close to one ulp // (unit in the last place) of the floating point representation - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 57_049.229_045_244_22, - semi_major_axis: 1.8014943988486352, - eccentricity: 0.283_514_142_249_080_7, - inclination: 0.20264170920820326, - ascending_node_longitude: 8.118_562_444_269_591E-3, - periapsis_argument: 1.244_795_311_814_302, - mean_anomaly: 0.44065425435816186, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 57_049.229_045_244_22, + semi_major_axis: 1.8014943988486352, + eccentricity: 0.283_514_142_249_080_7, + inclination: 0.20264170920820326, + ascending_node_longitude: 8.118_562_444_269_591E-3, + periapsis_argument: 1.244_795_311_814_302, + mean_anomaly: 0.44065425435816186, + }, + uncertainty: None, + covariance: None, + }; // Compare the prelim_orbit with the expected orbit assert!(approx_equal(prelim_orbit, &expected_orbit, tol)); @@ -1766,15 +1770,19 @@ pub(crate) mod gauss_test { let binding = a.prelim_orbit(&IODParams::default()).unwrap(); let prelim_orbit_a = binding.get_orbit(); - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 57049.22904525282, - semi_major_axis: 1.801490008178814, - eccentricity: 0.28350961635625993, - inclination: 0.20264261257939395, - ascending_node_longitude: 0.008105552171682476, - periapsis_argument: 1.244832121745955, - mean_anomaly: 0.4406444535028061, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 57049.22904525282, + semi_major_axis: 1.801490008178814, + eccentricity: 0.28350961635625993, + inclination: 0.20264261257939395, + ascending_node_longitude: 0.008105552171682476, + periapsis_argument: 1.244832121745955, + mean_anomaly: 0.4406444535028061, + }, + uncertainty: None, + covariance: None, + }; // Compare the prelim_orbit_a with the expected orbit assert!(approx_equal(prelim_orbit_a, &expected_orbit, tol)); @@ -1814,15 +1822,19 @@ pub(crate) mod gauss_test { // This is the expected orbit based on the Orbfit software // The values are obtained from the Orbfit output for the same observations - let expected_orbfit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 57_049.229_045_608_86, - semi_major_axis: 1.8013098187420686, - eccentricity: 0.28347096712267805, - inclination: 0.202_617_665_872_441_2, - ascending_node_longitude: 8.194_805_420_465_082E-3, - periapsis_argument: 1.2446747244785052, - mean_anomaly: 0.44073731381184733, - }); + let expected_orbfit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 57_049.229_045_608_86, + semi_major_axis: 1.8013098187420686, + eccentricity: 0.28347096712267805, + inclination: 0.202_617_665_872_441_2, + ascending_node_longitude: 8.194_805_420_465_082E-3, + periapsis_argument: 1.2446747244785052, + mean_anomaly: 0.44073731381184733, + }, + uncertainty: None, + covariance: None, + }; // Compare the prelim_orbit_b with the expected orbit // The epsilon value is set to 1e-13 for a very close match between the OrbFit and the Rust implementation diff --git a/src/initial_orbit_determination/gauss_result.rs b/src/initial_orbit_determination/gauss_result.rs index dd0bf9c..c57604e 100644 --- a/src/initial_orbit_determination/gauss_result.rs +++ b/src/initial_orbit_determination/gauss_result.rs @@ -223,15 +223,19 @@ mod gauss_results_tests { use std::format; fn dummy_orbit() -> OrbitalElements { - OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 59000.0, - semi_major_axis: 2.5, - eccentricity: 0.12, - inclination: 0.1, - ascending_node_longitude: 1.2, - periapsis_argument: 0.8, - mean_anomaly: 0.3, - }) + OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 59000.0, + semi_major_axis: 2.5, + eccentricity: 0.12, + inclination: 0.1, + ascending_node_longitude: 1.2, + periapsis_argument: 0.8, + mean_anomaly: 0.3, + }, + uncertainty: None, + covariance: None, + } } #[test] @@ -264,15 +268,19 @@ mod gauss_results_tests { #[test] fn test_display_format_summary() { - let orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 59001.5, - semi_major_axis: 1.234567, - eccentricity: 0.1, - inclination: 0.2, - ascending_node_longitude: 0.3, - periapsis_argument: 0.4, - mean_anomaly: 0.5, - }); + let orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 59001.5, + semi_major_axis: 1.234567, + eccentricity: 0.1, + inclination: 0.2, + ascending_node_longitude: 0.3, + periapsis_argument: 0.4, + mean_anomaly: 0.5, + }, + uncertainty: None, + covariance: None, + }; let result = GaussResult::CorrectedOrbit(orbit); @@ -284,6 +292,6 @@ mod gauss_results_tests { assert!(output.contains("a (semi-major axis) = 1.234567 AU")); assert!(output.contains("e (eccentricity) = 0.100000")); assert!(output.contains("i (inclination) = 0.200000 rad (11.459156°)")); - assert!(output.contains("Keplerian Elements @ epoch (MJD): 59001.500000")); + assert!(output.contains("Elements @ epoch (MJD): 59001.500000")); } } diff --git a/src/lib.rs b/src/lib.rs index 2a1ec28..f1944c8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -28,7 +28,12 @@ //! - Classical **Gauss method** for three observations. //! - **Orbital elements**: //! - Classical **Keplerian elements**, -//! - **Equinoctial elements** with conversions and two-body solver. +//! - **Equinoctial elements** with conversions and two-body solver, +//! - **Cometary elements** for parabolic and hyperbolic trajectories. +//! - **Uncertainty propagation**: +//! - Full **covariance matrix propagation** between orbital element representations, +//! - Standard deviations and correlations tracked via **Jacobian transformations**, +//! - Rigorous linear uncertainty propagation for orbit conversions. //! - **Reference frames & preprocessing**: //! - Precession, nutation (IAU 1980), aberration, and light-time correction, //! - Ecliptic ↔ equatorial conversions, RA/DEC parsing, time systems. @@ -47,6 +52,22 @@ //! - **Vaisalä method** for short arcs, //! - Full support for **hyperbolic trajectories**. //! +//! ## Why Track Uncertainties? +//! +//! Orbital element uncertainties quantify the confidence in derived orbits and are essential for: +//! +//! - **Orbit quality assessment** — distinguishing reliable solutions from poorly-constrained fits. +//! - **Prediction uncertainty** — forecasting position errors for future ephemerides. +//! - **Data association** — deciding whether observations belong to the same object. +//! - **Differential correction** — weighting observations and detecting outliers. +//! - **Mission planning** — evaluating pointing and timing requirements for follow-up observations. +//! +//! Outfit propagates uncertainties rigorously through coordinate transformations using **Jacobian +//! matrices** and **linear covariance propagation**, ensuring that correlations between orbital +//! parameters are preserved when converting between Keplerian, equinoctial, and cometary representations. +//! +//! For mathematical details and usage examples, see the [`orbit_type`] module documentation. +//! //! ## Workflow at a Glance //! //! 1. **Load observations** into an [`ObsDataset`](photom::observation_dataset::ObsDataset) @@ -177,7 +198,7 @@ //! //! - [`initial_orbit_determination`] — Gauss IOD algorithm, triplet generation, IOD parameters, public api to perform iod on an [`ObsDataset`](photom::observation_dataset::ObsDataset). //! - [`jpl_ephem`] — Ephemerides backends (Horizons/NAIF DE440). -//! - [`orbit_type`] — Orbital element representations (Keplerian, Equinoctial, Cometary). +//! - [`orbit_type`] — Orbital element representations (Keplerian, Equinoctial, Cometary) with full covariance propagation and uncertainty tracking for conversions between representations. //! - [`ref_system`] — Reference frame transformations. //! - [`constants`] — Physical constants and unit conversions. //! - [`conversion`] — RA/DEC parsing and coordinate utilities. diff --git a/src/orb_elem.rs b/src/orb_elem.rs index 54bbf0d..404f1f1 100644 --- a/src/orb_elem.rs +++ b/src/orb_elem.rs @@ -129,15 +129,19 @@ pub(crate) fn ccek1( let true_anomaly = sin_true_anomaly.atan2(cos_true_anomaly); let periapsis_argument = periapsis_argument_from_true_anomaly(true_anomaly); - OrbitalElements::Cometary(CometaryElements { - reference_epoch, - perihelion_distance, - eccentricity, - inclination, - ascending_node_longitude, - periapsis_argument, - true_anomaly, - }) + OrbitalElements::Cometary { + elements: CometaryElements { + reference_epoch, + perihelion_distance, + eccentricity, + inclination, + ascending_node_longitude, + periapsis_argument, + true_anomaly, + }, + uncertainty: None, + covariance: None, + } }; // ---- 3) Orbit classification -------------------------------------------- @@ -172,15 +176,19 @@ pub(crate) fn ccek1( let cos_periapsis = x1n * position_orbital.x + x2n * position_orbital.y; let periapsis_argument = wrap_0_2pi(sin_periapsis.atan2(cos_periapsis)); - OrbitalElements::Keplerian(KeplerianElements { - reference_epoch, - semi_major_axis, - eccentricity, - inclination, - ascending_node_longitude, - periapsis_argument, - mean_anomaly, - }) + OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch, + semi_major_axis, + eccentricity, + inclination, + ascending_node_longitude, + periapsis_argument, + mean_anomaly, + }, + uncertainty: None, + covariance: None, + } } else if reciprocal_semi_major_axis.abs() <= EPS_PARAB { // -------- Parabolic orbit parabolic_solution() @@ -201,15 +209,19 @@ pub(crate) fn ccek1( let perihelion_distance = semi_latus_rectum / (1.0 + eccentricity); let periapsis_argument = periapsis_argument_from_true_anomaly(true_anomaly); - OrbitalElements::Cometary(CometaryElements { - reference_epoch, - perihelion_distance, - eccentricity, - inclination, - ascending_node_longitude, - periapsis_argument, - true_anomaly, - }) + OrbitalElements::Cometary { + elements: CometaryElements { + reference_epoch, + perihelion_distance, + eccentricity, + inclination, + ascending_node_longitude, + periapsis_argument, + true_anomaly, + }, + uncertainty: None, + covariance: None, + } } } @@ -408,7 +420,7 @@ mod orb_elem_test { let orbit = ccek1(&position, &velocity, 0.0); match orbit { - OrbitalElements::Cometary(ce) => { + OrbitalElements::Cometary { elements: ce, .. } => { assert!( ce.perihelion_distance > 0.0, "Perihelion distance must be positive, got {q}", @@ -421,7 +433,7 @@ mod orb_elem_test { ); assert!(ce.true_anomaly.is_finite()); } - OrbitalElements::Keplerian(_) | OrbitalElements::Equinoctial(_) => { + OrbitalElements::Keplerian { .. } | OrbitalElements::Equinoctial { .. } => { panic!("Expected Cometary elements, got {orbit:#?}") } } @@ -453,7 +465,7 @@ mod orb_elem_test { let orbit = ccek1(&position, &velocity, 0.0); match orbit { - OrbitalElements::Keplerian(ke) => { + OrbitalElements::Keplerian { elements: ke, .. } => { // Eccentricity should be very close to 1 assert!( (ke.eccentricity - 1.0).abs() < 1e-5, @@ -462,7 +474,7 @@ mod orb_elem_test { ); } - OrbitalElements::Cometary(ce) => { + OrbitalElements::Cometary { elements: ce, .. } => { assert!( (ce.eccentricity - 1.0).abs() < 1e-5, "Near-parabolic eccentricity should be ~1, got {e}", @@ -470,7 +482,7 @@ mod orb_elem_test { ); } - OrbitalElements::Equinoctial(_) => { + OrbitalElements::Equinoctial { .. } => { panic!("Expected Keplerian or Cometary elements, got {orbit:#?}") } } @@ -552,34 +564,34 @@ mod orb_elem_test { // Extract invariants a, e, i (Keplerian) or q, e, i (Cometary). let (a1, e1, i1, om1) = match &orbit { - OrbitalElements::Keplerian(k) => ( + OrbitalElements::Keplerian { elements: k, .. } => ( k.semi_major_axis, k.eccentricity, k.inclination, k.ascending_node_longitude, ), - OrbitalElements::Cometary(c) => ( + OrbitalElements::Cometary { elements: c, .. } => ( c.perihelion_distance, c.eccentricity, c.inclination, c.ascending_node_longitude, ), - OrbitalElements::Equinoctial(_) => panic!("Unexpected equinoctial here"), + OrbitalElements::Equinoctial { .. } => panic!("Unexpected equinoctial here"), }; let (a2, e2, i2, om2) = match &orbit2 { - OrbitalElements::Keplerian(k) => ( + OrbitalElements::Keplerian { elements: k, .. } => ( k.semi_major_axis, k.eccentricity, k.inclination, k.ascending_node_longitude, ), - OrbitalElements::Cometary(c) => ( + OrbitalElements::Cometary { elements: c, .. } => ( c.perihelion_distance, c.eccentricity, c.inclination, c.ascending_node_longitude, ), - OrbitalElements::Equinoctial(_) => panic!("Unexpected equinoctial here"), + OrbitalElements::Equinoctial { .. } => panic!("Unexpected equinoctial here"), }; // Invariants under active Z-rotation of (r, v): @@ -636,7 +648,7 @@ mod orb_elem_test { let wrap_0_2pi = |ang: f64| ang.rem_euclid(TWO_PI); match orbit { - OrbitalElements::Keplerian(ke) => { + OrbitalElements::Keplerian { elements: ke, .. } => { // Eccentricity must match the control within tolerance prop_assert!((e_ctrl - ke.eccentricity).abs() < 1e-6); @@ -668,7 +680,7 @@ mod orb_elem_test { prop_assert!((-1e-12..=std::f64::consts::PI + 1e-12).contains(&inc)); } - OrbitalElements::Cometary(ce) => { + OrbitalElements::Cometary { elements: ce, .. } => { // Eccentricity must match the control within tolerance prop_assert!((e_ctrl - ce.eccentricity).abs() < 1e-6); @@ -692,7 +704,7 @@ mod orb_elem_test { prop_assert!((-1e-12..=std::f64::consts::PI + 1e-12).contains(&inc)); } - OrbitalElements::Equinoctial(_) => { + OrbitalElements::Equinoctial { .. } => { // Not tested here prop_assert!(false, "Equinoctial elements not tested"); } diff --git a/src/orbit_type/cometary_element.rs b/src/orbit_type/cometary_element.rs index 686beaa..569927a 100644 --- a/src/orbit_type/cometary_element.rs +++ b/src/orbit_type/cometary_element.rs @@ -1,25 +1,181 @@ +//! Cometary (perihelion-based) orbital elements for parabolic and hyperbolic trajectories +//! +//! This module implements the **cometary orbital element representation**, which uses +//! perihelion distance $q$ and true anomaly $\nu$ instead of semi-major axis $a$ and +//! mean anomaly $M$. This parameterization is **natural and numerically stable** for +//! cometary and interstellar objects on parabolic ($e = 1$) or hyperbolic ($e > 1$) orbits. +//! +//! ## Motivation +//! +//! For unbound or marginally-bound trajectories, classical Keplerian elements become +//! problematic: +//! +//! - **Parabolic orbits** ($e = 1$): semi-major axis $a = \infty$ is undefined +//! - **Hyperbolic orbits** ($e > 1$): $a < 0$ is mathematically valid but less intuitive +//! +//! In contrast, the **perihelion distance** $q = a(1 - e)$ remains **finite and positive** +//! for all eccentricities, making it a robust observable directly linked to astrometric +//! measurements near perihelion passage. +//! +//! ## Element definition +//! +//! Cometary elements consist of six parameters: +//! +//! | Symbol | Name | Units | Notes | +//! |--------|------|-------|-------| +//! | $q$ | Perihelion distance | AU | Always positive | +//! | $e$ | Eccentricity | — | $e \geq 1$ for cometary/interstellar objects | +//! | $i$ | Inclination | radians | $0 \leq i \leq \pi$ | +//! | $\Omega$ | Longitude of ascending node | radians | $0 \leq \Omega < 2\pi$ | +//! | $\omega$ | Argument of periapsis | radians | $0 \leq \omega < 2\pi$ | +//! | $\nu$ | True anomaly | radians | Angle from periapsis at epoch | +//! +//! ## Conversions +//! +//! This module provides conversions to standard representations via `TryFrom` traits: +//! +//! - **Cometary → Keplerian**: Compute $a = q/(1-e)$, derive $M$ from $\nu$ via +//! eccentric/hyperbolic anomaly (see [`CometaryElements::hyperbolic_mean_anomaly`]) +//! +//! - **Cometary → Equinoctial**: Chain through Keplerian as intermediate step +//! +//! Conversions **fail** for parabolic orbits ($|e - 1| < 10^{-12}$) since $a$ becomes +//! numerically unstable. +//! +//! ## Uncertainty propagation +//! +//! Uncertainties are propagated through coordinate transformations using **analytical Jacobians**: +//! +//! - [`jacobian_to_keplerian`](CometaryElements::jacobian_to_keplerian) — +//! $J = \partial(a,e,i,\Omega,\omega,M) / \partial(q,e,i,\Omega,\omega,\nu)$ +//! +//! - [`jacobian_to_equinoctial`](CometaryElements::jacobian_to_equinoctial) — +//! Computed via chain rule: $J_{\text{Com→Eq}} = J_{\text{Kep→Eq}} \cdot J_{\text{Com→Kep}}$ +//! +//! The Jacobian matrices enable **linear covariance propagation** (see +//! [`OrbitalCovariance::propagate`](crate::orbit_type::uncertainty::OrbitalCovariance::propagate)) +//! when converting between element sets. +//! +//! ## Usage +//! +//! ```rust +//! use outfit::orbit_type::cometary_element::CometaryElements; +//! use outfit::orbit_type::keplerian_element::KeplerianElements; +//! +//! // Define hyperbolic orbit elements (e.g., 'Oumuamua-like) +//! let cometary = CometaryElements { +//! reference_epoch: 58000.0, // MJD +//! perihelion_distance: 0.25, // AU +//! eccentricity: 1.2, // hyperbolic +//! inclination: 0.5, // radians +//! ascending_node_longitude: 1.0, +//! periapsis_argument: 0.3, +//! true_anomaly: 0.1, +//! }; +//! +//! // Convert to Keplerian (a < 0 for hyperbolic) +//! let keplerian = KeplerianElements::try_from(&cometary).unwrap(); +//! assert!(keplerian.semi_major_axis < 0.0); +//! +//! // Compute Jacobian for uncertainty propagation +//! let jacobian = cometary.jacobian_to_keplerian(); +//! ``` +//! +//! ## Implementation notes +//! +//! - **Hyperbolic mean anomaly** is computed using $\text{atanh}$ with numerical clamping +//! to avoid infinities when approaching asymptotic directions (see +//! [`hyperbolic_mean_anomaly`](CometaryElements::hyperbolic_mean_anomaly)) +//! +//! - **Parabolic tolerance**: orbits with $|e - 1| < 10^{-12}$ are rejected to avoid +//! catastrophic cancellation in $a = q/(1-e)$ +//! +//! - **Test coverage**: unit tests verify round-trip conversions, Jacobian accuracy via +//! finite differences, and edge cases (parabolic rejection, asymptotic behavior) +//! - Milani & Gronchi, *Theory of Orbit Determination* (2010). + +use nalgebra::{Matrix6, Vector6}; + use crate::{ orbit_type::{equinoctial_element::EquinoctialElements, keplerian_element::KeplerianElements}, outfit_errors::OutfitError, }; -/// # Cometary orbital elements +/// Cometary orbital elements with uncertainty propagation +/// +/// Cometary (perihelion-based) elements use perihelion distance $q$ and true anomaly $\nu$ +/// instead of semi-major axis $a$ and mean anomaly $M$. This representation is **natural for +/// parabolic and hyperbolic orbits**, where $a$ is infinite ($e = 1$) or negative ($e > 1$). +/// +/// ## Element definition +/// +/// The six cometary elements are: +/// +/// 1. **q** — Perihelion distance (AU): $q = a(1 - e)$ +/// 2. **e** — Eccentricity (unitless): $e \geq 1$ for cometary trajectories +/// 3. **i** — Inclination (radians) +/// 4. **Ω** — Longitude of ascending node (radians) +/// 5. **ω** — Argument of periapsis (radians) +/// 6. **ν** — True anomaly at epoch (radians) +/// +/// ## When to use cometary elements +/// +/// - **Parabolic orbits ($e = 1$)**: $a = \infty$ is undefined in Keplerian/equinoctial forms; +/// cometary elements remain well-defined +/// - **Hyperbolic orbits ($e > 1$)**: $a < 0$ is mathematically valid but less intuitive than $q > 0$ +/// - **Comets and interstellar objects**: observed near perihelion where $q$ and $\nu$ are +/// well-determined from astrometry +/// +/// For elliptic orbits ($e < 1$), Keplerian or equinoctial elements are generally preferred. +/// +/// ## Uncertainty representation and propagation +/// +/// Uncertainties in cometary elements are represented through: +/// +/// - **Standard deviations** on each element $(q, e, i, \Omega, \omega, \nu)$ via +/// [`CometaryUncertainty`](crate::orbit_type::uncertainty::CometaryUncertainty) +/// +/// - **Covariance matrices** capturing correlations via +/// [`OrbitalCovariance`](crate::orbit_type::uncertainty::OrbitalCovariance) +/// +/// When converting to Keplerian or equinoctial elements, uncertainties are propagated using +/// analytical Jacobian matrices: /// -/// Cometary (perihelion-based) elements are convenient for **parabolic and -/// hyperbolic** solutions, where the semi-major axis is not finite (parabola) -/// or negative (hyperbola). +/// - [`jacobian_to_keplerian`](CometaryElements::jacobian_to_keplerian) — +/// $J_{\text{Com} \to \text{Kep}} = \partial(a,e,i,\Omega,\omega,M) / \partial(q,e,i,\Omega,\omega,\nu)$ /// -/// Units & conventions -/// -------------------- -/// - Distances in **AU**; angles in **radians**; epochs in **MJD (TDB)**. -/// - State is assumed heliocentric, equatorial mean J2000. -/// - For hyperbolic motion: `a < 0`, `e > 1`; for parabolic: `e = 1`. +/// - [`jacobian_to_equinoctial`](CometaryElements::jacobian_to_equinoctial) — +/// Computed via chain rule: $J_{\text{Com} \to \text{Eq}} = J_{\text{Kep} \to \text{Eq}} \cdot J_{\text{Com} \to \text{Kep}}$ /// -/// See also -/// ------------ -/// * [`KeplerianElements`] – Classical elements `(a, e, i, Ω, ω, M)` (supports elliptic & hyperbolic). -/// * [`EquinoctialElements`] – Non-singular elements `(a, h, k, p, q, λ)`. -/// * [`CometaryElements::hyperbolic_mean_anomaly`] – Returns the hyperbolic mean anomaly from `(e, ν)` +/// ## Mathematical relations +/// +/// The semi-major axis relates to perihelion distance by: +/// +/// $$ +/// a = \frac{q}{1 - e} +/// $$ +/// +/// - For $e < 1$ (ellipse): $a > q > 0$ +/// - For $e = 1$ (parabola): $a = \infty$ +/// - For $e > 1$ (hyperbola): $a < 0$ and $q > 0$ +/// +/// The mean anomaly $M$ is computed from true anomaly $\nu$ via: +/// +/// - **Elliptic** ($e < 1$): through eccentric anomaly $E$ +/// - **Hyperbolic** ($e > 1$): through hyperbolic anomaly $H$ (see [`hyperbolic_mean_anomaly`](CometaryElements::hyperbolic_mean_anomaly)) +/// +/// ## Units and conventions +/// +/// - Distances: **AU** +/// - Angles: **radians** +/// - Epochs: **MJD (TDB)** +/// - Reference frame: heliocentric ecliptic J2000 +/// +/// ## See also +/// +/// * [`KeplerianElements`] — Classical elements $(a, e, i, \Omega, \omega, M)$ +/// * [`EquinoctialElements`] — Non-singular elements $(a, h, k, p, q, \lambda)$ +/// * [`hyperbolic_mean_anomaly`](CometaryElements::hyperbolic_mean_anomaly) — Computes $M$ from $(e, \nu)$ for $e > 1$ #[derive(Debug, Clone, PartialEq)] pub struct CometaryElements { /// Reference epoch of the element set (MJD, TDB). @@ -131,6 +287,140 @@ impl CometaryElements { mean_anomaly: m, // interpreted as hyperbolic mean anomaly M }) } + + /// Compute the Jacobian of the cometary-to-Keplerian transformation. + /// + /// Given the cometary element vector + /// $\mathbf{x} = [q, e, i, \Omega, \omega, \nu]^\top$, + /// this method returns the $6 \times 6$ matrix + /// + /// $$J = \frac{\partial \mathbf{y}}{\partial \mathbf{x}}, \quad + /// \mathbf{y} = [a, e, i, \Omega, \omega, M]^\top$$ + /// + /// Derivation + /// ---------- + /// The semi-major axis is recovered via: + /// + /// $$a = \frac{q}{1 - e}$$ + /// + /// The mean anomaly $M$ is obtained from the true anomaly $\nu$ through + /// the eccentric anomaly $E$ (Kepler's equation). The closed-form + /// partial derivatives are: + /// + /// $$\frac{\partial a}{\partial q} = \frac{1}{1-e}, \quad + /// \frac{\partial a}{\partial e} = \frac{q}{(1-e)^2}$$ + /// + /// $$\frac{\partial M}{\partial \nu} = \frac{(1-e^2)^{3/2}}{(1 + e\cos\nu)^2}$$ + /// + /// $$\frac{\partial M}{\partial e} = -\frac{\sin\nu\,(2 + e\cos\nu)}{(1 + e\cos\nu)^2}$$ + /// + /// Arguments + /// --------- + /// * `&self` – Cometary elements $(q, e, i, \Omega, \omega, \nu)$. + /// + /// Return + /// ------ + /// * `Matrix6` – The $6 \times 6$ Jacobian $\partial\mathbf{y}/\partial\mathbf{x}$. + /// + /// Notes + /// ----- + /// - Valid for elliptic orbits only ($0 \leq e < 1$). + /// Hyperbolic and parabolic cases require a separate treatment. + /// + /// See also + /// -------- + /// * [`CometaryElements::jacobian_to_equinoctial`] – Jacobian to equinoctial via chain rule. + pub fn jacobian_to_keplerian(&self) -> Matrix6 { + let q = self.perihelion_distance; + let e = self.eccentricity; + let nu = self.true_anomaly; + + let one_minus_e = 1.0 - e; + let cos_nu = nu.cos(); + let sin_nu = nu.sin(); + let denom = 1.0 + e * cos_nu; + + // a = q / (1 - e) holds for both elliptic (e<1, a>0) and hyperbolic (e>1, a<0). + let da_dq = 1.0 / one_minus_e; + let da_de = q / one_minus_e.powi(2); + + // Mean-anomaly partial derivatives differ between the elliptic and hyperbolic regimes. + let (dm_de, dm_dnu) = if e < 1.0 { + // Elliptic: M = E - e·sin(E); ∂M/∂ν = (1-e²)^(3/2)/(1+e·cos ν)² + // ∂M/∂e = -sin(ν)·(2+e·cos ν)/(1+e·cos ν)² + let dm_de = -sin_nu * (2.0 + e * cos_nu) / denom.powi(2); + let dm_dnu = (1.0 - e * e).powf(1.5) / denom.powi(2); + (dm_de, dm_dnu) + } else { + // Hyperbolic: M = e·sinh H - H; ∂M/∂ν = (e²-1)^(3/2)/(1+e·cos ν)² + // ∂M/∂e = sin(ν)·√(e²-1)·(2+e·cos ν)/(1+e·cos ν)² + let e2m1_sqrt = (e * e - 1.0).sqrt(); + let dm_de = sin_nu * e2m1_sqrt * (2.0 + e * cos_nu) / denom.powi(2); + let dm_dnu = e2m1_sqrt.powi(3) / denom.powi(2); + (dm_de, dm_dnu) + }; + + // Row ordering (target): [a, e, i, Ω, ω, M] + // Column ordering (source): [q, e, i, Ω, ω, ν] + + let col_q = Vector6::new(da_dq, 0.0, 0.0, 0.0, 0.0, 0.0); + + let col_e = Vector6::new( + da_de, // ∂a/∂e + 1.0, // ∂e/∂e + 0.0, // ∂i/∂e + 0.0, // ∂Ω/∂e + 0.0, // ∂ω/∂e + dm_de, // ∂M/∂e + ); + + let col_i = Vector6::new(0.0, 0.0, 1.0, 0.0, 0.0, 0.0); + let col_big_omega = Vector6::new(0.0, 0.0, 0.0, 1.0, 0.0, 0.0); + let col_omega = Vector6::new(0.0, 0.0, 0.0, 0.0, 1.0, 0.0); + + let col_nu = Vector6::new( + 0.0, // ∂a/∂ν + 0.0, // ∂e/∂ν + 0.0, // ∂i/∂ν + 0.0, // ∂Ω/∂ν + 0.0, // ∂ω/∂ν + dm_dnu, // ∂M/∂ν + ); + + Matrix6::from_columns(&[col_q, col_e, col_i, col_big_omega, col_omega, col_nu]) + } + + /// Compute the Jacobian of the cometary-to-equinoctial transformation. + /// + /// Uses the chain rule through the Keplerian representation: + /// + /// $$J_{\text{Com}\to\text{Eq}} = J_{\text{Kep}\to\text{Eq}} \cdot J_{\text{Com}\to\text{Kep}}$$ + /// + /// where $J_{\text{Kep}\to\text{Eq}}$ is evaluated at the Keplerian elements + /// obtained by converting `self`. + /// + /// Arguments + /// --------- + /// * `&self` – Cometary elements $(q, e, i, \Omega, \omega, t_p)$. + /// * `mu` – Gravitational parameter $\mu = GM_\odot$ (AU³/day²). + /// + /// Return + /// ------ + /// * `Ok(Matrix6)` – The $6 \times 6$ Jacobian matrix + /// $ \frac{\partial\mathbf{y}_\text{eq}}{ \partial \mathbf{x} _\text{com}} $. + /// * `Err(OutfitError)` – If the cometary-to-Keplerian conversion fails + /// (e.g. parabolic orbit $e = 1$). + /// + /// See also + /// -------- + /// * [`CometaryElements::jacobian_to_keplerian`] – First factor of the chain rule. + /// * [`KeplerianElements::jacobian_to_equinoctial`] – Second factor of the chain rule. + pub fn jacobian_to_equinoctial(&self) -> Result, OutfitError> { + let kep = KeplerianElements::try_from(self)?; + let j_kep_to_eq = kep.jacobian_to_equinoctial(); + let j_com_to_kep = self.jacobian_to_keplerian(); + Ok(j_kep_to_eq * j_com_to_kep) + } } impl TryFrom for KeplerianElements { @@ -224,11 +514,7 @@ impl fmt::Display for CometaryElements { /// * [`KeplerianElements`] – Classical representation used in conversions. fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let rad_to_deg = 180.0 / std::f64::consts::PI; - writeln!( - f, - "Cometary Elements @ epoch (MJD): {:.6}", - self.reference_epoch - )?; + writeln!(f, "Elements @ epoch (MJD): {:.6}", self.reference_epoch)?; writeln!(f, "------------------------------------------------")?; writeln!( f, @@ -504,10 +790,7 @@ mod cometary_element_tests { }; let s = format!("{ce}"); - assert!( - s.contains("Cometary Elements @ epoch (MJD)"), - "header missing" - ); + assert!(s.contains("Elements @ epoch (MJD)"), "header missing"); assert!(s.contains("q (perihelion distance)"), "q line missing"); assert!(s.contains("e (eccentricity)"), "e line missing"); assert!( @@ -740,3 +1023,137 @@ mod cometary_element_proptests { } } } + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + use nalgebra::Matrix6; + + fn make_cometary() -> CometaryElements { + CometaryElements { + reference_epoch: 60000.0, + perihelion_distance: 1.5, + eccentricity: 1.5, // must be > 1 for hyperbolic CometaryElements + inclination: 0.5, + ascending_node_longitude: 1.2, + periapsis_argument: 0.8, + true_anomaly: 0.6, + } + } + + fn cometary_fields(c: &CometaryElements) -> [f64; 6] { + [ + c.perihelion_distance, + c.eccentricity, + c.inclination, + c.ascending_node_longitude, + c.periapsis_argument, + c.true_anomaly, + ] + } + + fn make_cometary_from_fields(base: &CometaryElements, f: [f64; 6]) -> CometaryElements { + CometaryElements { + reference_epoch: base.reference_epoch, + perihelion_distance: f[0], + eccentricity: f[1], + inclination: f[2], + ascending_node_longitude: f[3], + periapsis_argument: f[4], + true_anomaly: f[5], + } + } + + fn cometary_to_kep_arr(c: &CometaryElements) -> [f64; 6] { + let k = KeplerianElements::try_from(c).unwrap(); + [ + k.semi_major_axis, + k.eccentricity, + k.inclination, + k.ascending_node_longitude, + k.periapsis_argument, + k.mean_anomaly, + ] + } + + fn cometary_to_eq_arr(c: &CometaryElements) -> [f64; 6] { + let eq = EquinoctialElements::from(KeplerianElements::try_from(c).unwrap()); + [ + eq.semi_major_axis, + eq.eccentricity_sin_lon, + eq.eccentricity_cos_lon, + eq.tan_half_incl_sin_node, + eq.tan_half_incl_cos_node, + eq.mean_longitude, + ] + } + + fn fd_jacobian [f64; 6]>( + c: &CometaryElements, + eps: f64, + output: F, + ) -> Matrix6 { + let fields = cometary_fields(c); + Matrix6::from_fn(|row, col| { + let mut fwd = fields; + let mut bwd = fields; + fwd[col] += eps; + bwd[col] -= eps; + let y_fwd = output(&make_cometary_from_fields(c, fwd)); + let y_bwd = output(&make_cometary_from_fields(c, bwd)); + (y_fwd[row] - y_bwd[row]) / (2.0 * eps) + }) + } + + #[test] + fn test_jacobian_cometary_to_keplerian_against_fd() { + let c = make_cometary(); + let analytical = c.jacobian_to_keplerian(); + let numerical = fd_jacobian(&c, 1e-6, cometary_to_kep_arr); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!( + analytical[(row, col)], + numerical[(row, col)], + epsilon = 1e-7 + ); + } + } + } + + #[test] + fn test_jacobian_cometary_to_equinoctial_against_fd() { + let c = make_cometary(); + let analytical = c.jacobian_to_equinoctial().unwrap(); + let numerical = fd_jacobian(&c, 1e-6, cometary_to_eq_arr); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!( + analytical[(row, col)], + numerical[(row, col)], + epsilon = 1e-7 + ); + } + } + } + + #[test] + fn test_jacobian_cometary_chain_rule_consistency() { + let c = make_cometary(); + let kep = KeplerianElements::try_from(&c).unwrap(); + + let j_com_to_kep = c.jacobian_to_keplerian(); + let j_kep_to_eq = kep.jacobian_to_equinoctial(); + let expected = j_kep_to_eq * j_com_to_kep; + let result = c.jacobian_to_equinoctial().unwrap(); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!(result[(row, col)], expected[(row, col)], epsilon = 1e-14); + } + } + } +} diff --git a/src/orbit_type/equinoctial_element.rs b/src/orbit_type/equinoctial_element.rs index 2e63de5..415ca02 100644 --- a/src/orbit_type/equinoctial_element.rs +++ b/src/orbit_type/equinoctial_element.rs @@ -59,9 +59,41 @@ //! //! ## Advantages of equinoctial elements //! -//! * Avoid singularities when `e → 0` or `i → 0` -//! * Smooth derivatives, ideal for gradient-based fitting -//! * Directly compatible with least-squares adjustment and orbit uncertainty analysis +//! Equinoctial elements provide several critical advantages over classical Keplerian elements: +//! +//! * **Non-singular for $e < 1$ and $0 \leq i < \pi$** — no undefined elements for circular +//! or equatorial orbits, unlike Keplerian $\omega$ and $\Omega$ +//! +//! * **Smooth, well-conditioned derivatives** — Jacobian matrices remain well-behaved throughout +//! the valid domain, enabling robust numerical optimization and uncertainty propagation +//! +//! * **Ideal for orbit determination** — gradient-based least-squares fitting converges reliably +//! without special handling of degenerate cases +//! +//! * **Superior uncertainty propagation** — covariance matrices transform smoothly without +//! singular amplification near $e = 0$ or $i = 0$ +//! +//! * **Compatible with two-body propagation** — Keplerian dynamics equations extend naturally +//! to equinoctial form with regular behavior +//! +//! ## Uncertainty representation +//! +//! Uncertainties in equinoctial elements are represented through: +//! +//! - **Standard deviations** on each element $(a, h, k, p, q, \lambda)$ via +//! [`EquinoctialUncertainty`](crate::orbit_type::uncertainty::EquinoctialUncertainty) +//! +//! - **Covariance matrices** capturing correlations between elements via +//! [`OrbitalCovariance`](crate::orbit_type::uncertainty::OrbitalCovariance) +//! +//! When converting to/from Keplerian elements, uncertainties are propagated using analytical +//! Jacobian matrices: +//! +//! - [`jacobian_to_keplerian`](EquinoctialElements::jacobian_to_keplerian) — +//! $J_{\text{Eq} \to \text{Kep}} = \partial(a,e,i,\Omega,\omega,M) / \partial(a,h,k,p,q,\lambda)$ +//! +//! These Jacobians handle singular points by zeroing derivatives when denominators vanish +//! (threshold $\epsilon = 10^{-12}$), ensuring numerical stability at all orbital configurations. //! //! ## Example //! @@ -91,10 +123,11 @@ //! //! - [`KeplerianElements`](crate::orbit_type::keplerian_element::KeplerianElements) //! - Milani & Gronchi, *Theory of Orbit Determination* (2010). + use core::f64; use std::{f64::consts::PI, fmt}; -use nalgebra::{Matrix3x6, Matrix6, Matrix6x3, Vector3}; +use nalgebra::{Matrix3x6, Matrix6, Matrix6x3, Vector3, Vector6}; use roots::{find_root_newton_raphson, SimpleConvergency}; use crate::{ @@ -933,6 +966,177 @@ impl EquinoctialElements { dvel_delem: dvel_delem_at_t1, }) } + + /// Compute the Jacobian of the equinoctial-to-Keplerian transformation. + /// + /// Given the equinoctial element vector + /// $\mathbf{x} = [a, h, k, p, q, \lambda]^\top$, + /// this method returns the $6 \times 6$ matrix + /// + /// $$J = \frac{\partial \mathbf{y}}{\partial \mathbf{x}}, \quad + /// \mathbf{y} = [a, e, i, \Omega, \omega, M]^\top$$ + /// + /// where the column ordering matches the source vector $\mathbf{x}$ and + /// the row ordering matches the target vector $\mathbf{y}$. + /// + /// Derivation + /// ---------- + /// The mapping relies on the intermediate quantities: + /// + /// $$e = \sqrt{h^2 + k^2}, \quad \varpi = \text{atan2}(h, k)$$ + /// $$t = \sqrt{p^2 + q^2}, \quad \Omega = \text{atan2}(p, q)$$ + /// $$i = 2\arctan(t), \quad \omega = \varpi - \Omega, \quad M = \lambda - \varpi$$ + /// + /// The partial derivatives of $\varpi$ with respect to $h$ and $k$ are: + /// + /// $$\frac{\partial \varpi}{\partial h} = \frac{k}{e^2}, \quad + /// \frac{\partial \varpi}{\partial k} = -\frac{h}{e^2}$$ + /// + /// The partial derivatives of $i$ with respect to $p$ and $q$ are: + /// + /// $$\frac{\partial i}{\partial p} = \frac{2p}{t(1+t^2)}, \quad + /// \frac{\partial i}{\partial q} = \frac{2q}{t(1+t^2)}$$ + /// + /// Degenerate cases and regularization + /// ------------------------------------ + /// + /// The transformation has singularities at certain orbital configurations: + /// + /// * **Near-circular orbits ($e \approx 0$)**: $\varpi = \text{atan2}(h, k)$ is undefined + /// when $h^2 + k^2 \to 0$. Derivatives $\partial \varpi / \partial h$ and + /// $\partial \varpi / \partial k$ involve $e^{-2}$, which diverges. + /// + /// **Regularization**: When $e < \epsilon$ (with $\epsilon = 10^{-12}$), these + /// derivatives are set to zero. Physically, the periapsis direction has no meaning + /// for a circular orbit. + /// + /// * **Near-equatorial orbits ($i \approx 0$, i.e., $t = \tan(i/2) \approx 0$)**: + /// $\Omega = \text{atan2}(p, q)$ is undefined when $p^2 + q^2 \to 0$. Derivatives + /// $\partial \Omega / \partial p$ and $\partial \Omega / \partial q$ involve $t^{-2}$, + /// which diverges. + /// + /// **Regularization**: When $t < \epsilon$, these derivatives are set to zero. + /// The ascending node direction is undefined for an equatorial orbit. + /// + /// These regularizations prevent numerical overflow but may **underestimate uncertainties** + /// in $\omega$ and $\Omega$ for degenerate orbits. For such cases, equinoctial elements + /// should remain the primary representation. + /// + /// ## Usage in uncertainty propagation + /// + /// This Jacobian transforms covariance matrices from equinoctial to Keplerian: + /// + /// $$ + /// \Sigma_{\text{Kep}} = J \, \Sigma_{\text{Eq}} \, J^\top + /// $$ + /// + /// See [`OrbitalCovariance::propagate`](crate::orbit_type::uncertainty::OrbitalCovariance::propagate) + /// and [`OrbitalElements::to_keplerian`](crate::orbit_type::OrbitalElements::to_keplerian). + /// + /// ## Arguments + /// + /// * `&self` – Equinoctial elements $(a, h, k, p, q, \lambda)$ + /// + /// ## Return + /// + /// * `Matrix6` – The $6 \times 6$ Jacobian matrix $\partial\mathbf{y}/\partial\mathbf{x}$ + /// where rows correspond to $[a, e, i, \Omega, \omega, M]$ and columns to $[a, h, k, p, q, \lambda]$ + /// + /// ## See also + /// + /// * [`KeplerianElements::jacobian_to_equinoctial`](crate::orbit_type::keplerian_element::KeplerianElements::jacobian_to_equinoctial) — Inverse Jacobian + /// * [`OrbitalElements::to_keplerian`](crate::orbit_type::OrbitalElements::to_keplerian) — High-level conversion with uncertainty + pub fn jacobian_to_keplerian(&self) -> Matrix6 { + let h = self.eccentricity_sin_lon; + let k = self.eccentricity_cos_lon; + let p = self.tan_half_incl_sin_node; + let q = self.tan_half_incl_cos_node; + + let eps = 1.0e-12; + + // --- eccentricity and longitude of periapsis --- + let e = (h * h + k * k).sqrt(); + let e_sq = e * e; + + // ∂varpi/∂h = k/e², ∂varpi/∂k = -h/e² (zero when e ≈ 0) + let (d_varpi_dh, d_varpi_dk) = if e < eps { + (0.0, 0.0) + } else { + (k / e_sq, -h / e_sq) + }; + + // --- inclination node --- + let t = (p * p + q * q).sqrt(); + let t_sq = t * t; + + // ∂i/∂p = 2p / (t(1+t²)), ∂i/∂q = 2q / (t(1+t²)) (zero when t ≈ 0) + let (d_i_dp, d_i_dq) = if t < eps { + (0.0, 0.0) + } else { + let denom = t * (1.0 + t_sq); + (2.0 * p / denom, 2.0 * q / denom) + }; + + // ∂Ω/∂p = q/t², ∂Ω/∂q = -p/t² (zero when t ≈ 0) + let (d_omega_node_dp, d_omega_node_dq) = if t < eps { + (0.0, 0.0) + } else { + (q / t_sq, -p / t_sq) + }; + + // Each Vector6 is one column: ∂y/∂x_j for source variable x_j. + // Row ordering (target): [a, e, i, Ω, ω, M] + // Column ordering (source): [a, h, k, p, q, λ] + + let col_a = Vector6::new(1.0, 0.0, 0.0, 0.0, 0.0, 0.0); + + let col_h = Vector6::new( + 0.0, // ∂a/∂h + h / e.max(eps), // ∂e/∂h = h/e + 0.0, // ∂i/∂h + 0.0, // ∂Ω/∂h + d_varpi_dh, // ∂ω/∂h = ∂varpi/∂h (∂Ω/∂h = 0) + -d_varpi_dh, // ∂M/∂h = -∂varpi/∂h + ); + + let col_k = Vector6::new( + 0.0, // ∂a/∂k + k / e.max(eps), // ∂e/∂k = k/e + 0.0, // ∂i/∂k + 0.0, // ∂Ω/∂k + d_varpi_dk, // ∂ω/∂k = ∂varpi/∂k + -d_varpi_dk, // ∂M/∂k = -∂varpi/∂k + ); + + let col_p = Vector6::new( + 0.0, // ∂a/∂p + 0.0, // ∂e/∂p + d_i_dp, // ∂i/∂p + d_omega_node_dp, // ∂Ω/∂p + -d_omega_node_dp, // ∂ω/∂p = -∂Ω/∂p (varpi independent of p) + 0.0, // ∂M/∂p + ); + + let col_q = Vector6::new( + 0.0, // ∂a/∂q + 0.0, // ∂e/∂q + d_i_dq, // ∂i/∂q + d_omega_node_dq, // ∂Ω/∂q + -d_omega_node_dq, // ∂ω/∂q = -∂Ω/∂q + 0.0, // ∂M/∂q + ); + + let col_lambda = Vector6::new( + 0.0, // ∂a/∂λ + 0.0, // ∂e/∂λ + 0.0, // ∂i/∂λ + 0.0, // ∂Ω/∂λ + 0.0, // ∂ω/∂λ + 1.0, // ∂M/∂λ + ); + + Matrix6::from_columns(&[col_a, col_h, col_k, col_p, col_q, col_lambda]) + } } impl From for KeplerianElements { @@ -966,11 +1170,7 @@ impl From<&EquinoctialElements> for KeplerianElements { impl fmt::Display for EquinoctialElements { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let rad_to_deg = 180.0 / std::f64::consts::PI; - writeln!( - f, - "Equinoctial Elements @ epoch (MJD): {:.6}", - self.reference_epoch - )?; + writeln!(f, "Elements @ epoch (MJD): {:.6}", self.reference_epoch)?; writeln!(f, "------------------------------------------------")?; writeln!( f, @@ -1385,3 +1585,161 @@ mod test_equinoctial_element { } } } + +#[cfg(test)] +mod tests_jacobian_equinoctial_to_keplerian { + use super::*; + use approx::assert_abs_diff_eq; + + fn fd_jacobian_eq_to_kep(eq: &EquinoctialElements, eps: f64) -> Matrix6 { + let fields = [ + eq.semi_major_axis, + eq.eccentricity_sin_lon, + eq.eccentricity_cos_lon, + eq.tan_half_incl_sin_node, + eq.tan_half_incl_cos_node, + eq.mean_longitude, + ]; + + let kep_to_arr = |k: KeplerianElements| { + [ + k.semi_major_axis, + k.eccentricity, + k.inclination, + k.ascending_node_longitude, + k.periapsis_argument, + k.mean_anomaly, + ] + }; + + let make_eq = |f: [f64; 6]| EquinoctialElements { + reference_epoch: eq.reference_epoch, + semi_major_axis: f[0], + eccentricity_sin_lon: f[1], + eccentricity_cos_lon: f[2], + tan_half_incl_sin_node: f[3], + tan_half_incl_cos_node: f[4], + mean_longitude: f[5], + }; + + let mut columns = [[0.0f64; 6]; 6]; + + for col in 0..6 { + let mut fwd = fields; + let mut bwd = fields; + fwd[col] += eps; + bwd[col] -= eps; + + let kep_fwd = kep_to_arr(KeplerianElements::from(&make_eq(fwd))); + let kep_bwd = kep_to_arr(KeplerianElements::from(&make_eq(bwd))); + + for row in 0..6 { + columns[col][row] = (kep_fwd[row] - kep_bwd[row]) / (2.0 * eps); + } + } + + Matrix6::from_fn(|row, col| columns[col][row]) + } + + #[test] + fn test_jacobian_to_keplerian_against_finite_differences() { + let eq = EquinoctialElements { + reference_epoch: 60000.0, + semi_major_axis: 2.5, + eccentricity_sin_lon: 0.15, + eccentricity_cos_lon: 0.25, + tan_half_incl_sin_node: 0.10, + tan_half_incl_cos_node: 0.20, + mean_longitude: 1.8, + }; + + let analytical = eq.jacobian_to_keplerian(); + let numerical = fd_jacobian_eq_to_kep(&eq, 1e-6); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!( + analytical[(row, col)], + numerical[(row, col)], + epsilon = 1e-7 + ); + } + } + } + + /// Near-circular orbit: e ≈ 0, degenerate 1/e² derivatives should be zeroed. + #[test] + fn test_jacobian_to_keplerian_near_circular() { + // Use h=k small enough that e = sqrt(h²+k²) < eps (1e-12) in the Jacobian, + // triggering the degenerate-zero branch for ∂ω/∂h, ∂ω/∂k, ∂M/∂h, ∂M/∂k. + // Note: ∂e/∂h = h/e → 1/√2 for h=k, which is well-defined and NOT zero. + let eq = EquinoctialElements { + reference_epoch: 60000.0, + semi_major_axis: 1.0, + eccentricity_sin_lon: 1e-13, + eccentricity_cos_lon: 1e-13, + tan_half_incl_sin_node: 0.10, + tan_half_incl_cos_node: 0.20, + mean_longitude: 1.0, + }; + + let j = eq.jacobian_to_keplerian(); + + // Only the 1/e²-diverging derivatives are zeroed by the regularisation + assert_abs_diff_eq!(j[(4, 1)], 0.0, epsilon = 1e-10); // ∂ω/∂h + assert_abs_diff_eq!(j[(4, 2)], 0.0, epsilon = 1e-10); // ∂ω/∂k + assert_abs_diff_eq!(j[(5, 1)], 0.0, epsilon = 1e-10); // ∂M/∂h + assert_abs_diff_eq!(j[(5, 2)], 0.0, epsilon = 1e-10); // ∂M/∂k + } + + /// Near-equatorial orbit: t ≈ 0, degenerate 1/t²-node derivatives should be zeroed. + #[test] + fn test_jacobian_to_keplerian_near_equatorial() { + // Use p=q small enough that t = sqrt(p²+q²) < eps (1e-12), triggering the + // degenerate-zero branch for ∂i/∂p, ∂Ω/∂p, ∂ω/∂p (and their q-counterparts). + let eq = EquinoctialElements { + reference_epoch: 60000.0, + semi_major_axis: 1.5, + eccentricity_sin_lon: 0.10, + eccentricity_cos_lon: 0.20, + tan_half_incl_sin_node: 1e-13, + tan_half_incl_cos_node: 1e-13, + mean_longitude: 2.0, + }; + + let j = eq.jacobian_to_keplerian(); + + // Both the bounded (∂i/∂p) and diverging (∂Ω/∂p, ∂ω/∂p) entries are zeroed + assert_abs_diff_eq!(j[(2, 3)], 0.0, epsilon = 1e-10); // ∂i/∂p + assert_abs_diff_eq!(j[(3, 3)], 0.0, epsilon = 1e-10); // ∂Ω/∂p + assert_abs_diff_eq!(j[(4, 3)], 0.0, epsilon = 1e-10); // ∂ω/∂p + } + + /// Round-trip consistency: J_kep_to_eq · J_eq_to_kep ≈ I₆. + #[test] + fn test_jacobian_round_trip_identity() { + let eq = EquinoctialElements { + reference_epoch: 60000.0, + semi_major_axis: 2.5, + eccentricity_sin_lon: 0.15, + eccentricity_cos_lon: 0.25, + tan_half_incl_sin_node: 0.10, + tan_half_incl_cos_node: 0.20, + mean_longitude: 1.8, + }; + + let kep = KeplerianElements::from(&eq); + + let j_eq_to_kep = eq.jacobian_to_keplerian(); + let j_kep_to_eq = kep.jacobian_to_equinoctial(); + + let product = j_eq_to_kep * j_kep_to_eq; + let identity = Matrix6::identity(); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!(product[(row, col)], identity[(row, col)], epsilon = 1e-10); + } + } + } +} diff --git a/src/orbit_type/keplerian_element.rs b/src/orbit_type/keplerian_element.rs index 88dda3b..5b2229c 100644 --- a/src/orbit_type/keplerian_element.rs +++ b/src/orbit_type/keplerian_element.rs @@ -34,17 +34,46 @@ //! - Angles: **radians** //! - Time: **days** (epoch usually in **MJD**, TDB/TT scale) //! -//! ## Degeneracies +//! ## Degeneracies and singularities //! -//! Classical Keplerian elements suffer from singularities: +//! Classical Keplerian elements suffer from singularities in certain orbital configurations: //! -//! - **Circular orbits (`e → 0`)**: periapsis argument ω becomes undefined. -//! → conventionally set to `0.0` during conversion. +//! - **Circular orbits ($e \to 0$)**: the periapsis argument $\omega$ becomes undefined because +//! there is no distinct periapsis point. Conventionally set to $\omega = 0$ during conversion, +//! but derivatives $\partial \omega / \partial h$ and $\partial \omega / \partial k$ are +//! singular ($\sim e^{-2}$). //! -//! - **Equatorial orbits (`i → 0`)**: ascending node Ω becomes undefined. -//! → conventionally set to `0.0` during conversion. +//! - **Equatorial orbits ($i \to 0$)**: the ascending node longitude $\Omega$ becomes undefined +//! because the orbital plane coincides with the reference plane. Conventionally set to $\Omega = 0$, +//! with singular derivatives $\partial \Omega / \partial p$ and $\partial \Omega / \partial q$ +//! ($\sim i^{-1}$). //! -//! For robust numerical work, the [`EquinoctialElements`](crate::orbit_type::equinoctial_element::EquinoctialElements) representation is recommended. +//! - **Circular equatorial orbits ($e \to 0$ and $i \to 0$)**: both $\omega$ and $\Omega$ are +//! undefined. Only the mean longitude $\lambda = \Omega + \omega + M$ remains well-defined. +//! +//! These singularities impact **uncertainty propagation**: near-circular or near-equatorial orbits +//! will have artificially inflated or ill-defined uncertainties in $\omega$ or $\Omega$ when +//! converted from equinoctial elements. +//! +//! **Recommendation**: For numerical work and uncertainty analysis on circular or equatorial orbits, +//! use [`EquinoctialElements`](crate::orbit_type::equinoctial_element::EquinoctialElements) as the +//! primary representation to avoid singularities. +//! +//! ## Uncertainty propagation from equinoctial elements +//! +//! When converting from equinoctial $(a, h, k, p, q, \lambda)$ to Keplerian $(a, e, i, \Omega, \omega, M)$, +//! uncertainties are propagated using the Jacobian matrix computed by [`jacobian_to_equinoctial`](KeplerianElements::jacobian_to_equinoctial): +//! +//! $$ +//! \Sigma_{\text{Kep}} = J_{\text{Eq} \to \text{Kep}} \, \Sigma_{\text{Eq}} \, J_{\text{Eq} \to \text{Kep}}^\top +//! $$ +//! +//! where $J_{\text{Eq} \to \text{Kep}}$ is the inverse Jacobian +//! (computed in [`EquinoctialElements::jacobian_to_keplerian`](crate::orbit_type::equinoctial_element::EquinoctialElements::jacobian_to_keplerian)). +//! +//! Near singular points, derivatives are regularized by setting them to zero when denominators +//! become small (threshold $\epsilon = 10^{-12}$). This prevents numerical overflow but may +//! underestimate uncertainties in degenerate cases. //! //! ## Example //! @@ -84,6 +113,8 @@ //! - [`principal_angle`](crate::kepler::principal_angle) – helper to normalize angular elements. //! - Milani & Gronchi, *Theory of Orbit Determination* (2010). +use nalgebra::{Matrix6, Vector6}; + use crate::{kepler::principal_angle, orbit_type::equinoctial_element::EquinoctialElements}; use std::fmt; @@ -200,6 +231,152 @@ impl KeplerianElements { mean_anomaly, } } + + /// Compute the Jacobian of the Keplerian-to-equinoctial transformation + /// + /// This method returns the $6 \times 6$ Jacobian matrix of the transformation from + /// Keplerian elements $(a, e, i, \Omega, \omega, M)$ to equinoctial elements $(a, h, k, p, q, \lambda)$. + /// + /// ## Mathematical formulation + /// + /// The equinoctial elements are defined by: + /// + /// $$ + /// \begin{aligned} + /// a &= a \\ + /// h &= e \sin(\varpi) \\ + /// k &= e \cos(\varpi) \\ + /// p &= \tan(i/2) \sin(\Omega) \\ + /// q &= \tan(i/2) \cos(\Omega) \\ + /// \lambda &= \varpi + M = \Omega + \omega + M + /// \end{aligned} + /// $$ + /// + /// where $\varpi = \Omega + \omega$ is the longitude of periapsis. + /// + /// ## Jacobian structure + /// + /// The Jacobian $J = \partial \mathbf{y} / \partial \mathbf{x}$ has the following structure: + /// + /// - **Column 1** ($\partial / \partial a$): Only $\partial a / \partial a = 1$ is nonzero + /// + /// - **Column 2** ($\partial / \partial e$): + /// $$\frac{\partial h}{\partial e} = \sin(\varpi), \quad \frac{\partial k}{\partial e} = \cos(\varpi)$$ + /// + /// - **Column 3** ($\partial / \partial i$): + /// $$\frac{\partial p}{\partial i} = \frac{1}{2 \cos^2(i/2)} \sin(\Omega), \quad + /// \frac{\partial q}{\partial i} = \frac{1}{2 \cos^2(i/2)} \cos(\Omega)$$ + /// + /// - **Column 4** ($\partial / \partial \Omega$): + /// $$\frac{\partial h}{\partial \Omega} = e \cos(\varpi), \quad + /// \frac{\partial k}{\partial \Omega} = -e \sin(\varpi)$$ + /// $$\frac{\partial p}{\partial \Omega} = \tan(i/2) \cos(\Omega), \quad + /// \frac{\partial q}{\partial \Omega} = -\tan(i/2) \sin(\Omega)$$ + /// $$\frac{\partial \lambda}{\partial \Omega} = 1$$ + /// + /// - **Column 5** ($\partial / \partial \omega$): + /// Same as $\partial / \partial \Omega$ for $h, k, \lambda$ but zero for $p, q$ + /// + /// - **Column 6** ($\partial / \partial M$): + /// Only $\partial \lambda / \partial M = 1$ is nonzero + /// + /// ## Handling singularities + /// + /// No singularities occur in this forward transformation (Keplerian → Equinoctial) because: + /// - $h, k$ are well-defined for all $e \geq 0$ + /// - $p, q$ are well-defined for all $i \in [0, \pi]$ + /// + /// This makes the Jacobian numerically stable for all orbital configurations. + /// + /// ## Usage in uncertainty propagation + /// + /// This Jacobian is used to transform covariance matrices from Keplerian to equinoctial: + /// + /// $$ + /// \Sigma_{\text{Eq}} = J \, \Sigma_{\text{Kep}} \, J^\top + /// $$ + /// + /// See [`OrbitalCovariance::propagate`](crate::orbit_type::uncertainty::OrbitalCovariance::propagate) + /// for the propagation implementation. + /// + /// ## Return + /// + /// A $6 \times 6$ matrix where: + /// - **Rows** correspond to target elements: $[a, h, k, p, q, \lambda]$ + /// - **Columns** correspond to source elements: $[a, e, i, \Omega, \omega, M]$ + /// + /// ## See also + /// + /// - [`EquinoctialElements::jacobian_to_keplerian`](crate::orbit_type::equinoctial_element::EquinoctialElements::jacobian_to_keplerian) — Inverse Jacobian + /// - [`OrbitalElements::to_equinoctial`](crate::orbit_type::OrbitalElements::to_equinoctial) — High-level conversion with uncertainty propagation + pub fn jacobian_to_equinoctial(&self) -> Matrix6 { + let e = self.eccentricity; + let i = self.inclination; + let varpi = self.ascending_node_longitude + self.periapsis_argument; + let big_omega = self.ascending_node_longitude; + + let sin_varpi = varpi.sin(); + let cos_varpi = varpi.cos(); + let sin_omega = big_omega.sin(); + let cos_omega = big_omega.cos(); + + let half_i = i / 2.0; + let tan_half_i = half_i.tan(); + let d_tan_half_i_d_i = 0.5 / half_i.cos().powi(2); + + // Each Vector6 is one column: ∂y/∂x_j for source variable x_j. + // Row ordering (target): [a, h, k, p, q, λ] + // Column ordering (source): [a, e, i, Ω, ω, M] + + let col_a = Vector6::new(1.0, 0.0, 0.0, 0.0, 0.0, 0.0); + + let col_e = Vector6::new( + 0.0, // ∂a/∂e + sin_varpi, // ∂h/∂e + cos_varpi, // ∂k/∂e + 0.0, // ∂p/∂e + 0.0, // ∂q/∂e + 0.0, // ∂λ/∂e + ); + + let col_i = Vector6::new( + 0.0, // ∂a/∂i + 0.0, // ∂h/∂i + 0.0, // ∂k/∂i + d_tan_half_i_d_i * sin_omega, // ∂p/∂i + d_tan_half_i_d_i * cos_omega, // ∂q/∂i + 0.0, // ∂λ/∂i + ); + + let col_big_omega = Vector6::new( + 0.0, // ∂a/∂Ω + e * cos_varpi, // ∂h/∂Ω + -e * sin_varpi, // ∂k/∂Ω + tan_half_i * cos_omega, // ∂p/∂Ω + -tan_half_i * sin_omega, // ∂q/∂Ω + 1.0, // ∂λ/∂Ω + ); + + let col_omega = Vector6::new( + 0.0, // ∂a/∂ω + e * cos_varpi, // ∂h/∂ω + -e * sin_varpi, // ∂k/∂ω + 0.0, // ∂p/∂ω + 0.0, // ∂q/∂ω + 1.0, // ∂λ/∂ω + ); + + let col_m = Vector6::new( + 0.0, // ∂a/∂M + 0.0, // ∂h/∂M + 0.0, // ∂k/∂M + 0.0, // ∂p/∂M + 0.0, // ∂q/∂M + 1.0, // ∂λ/∂M + ); + + Matrix6::from_columns(&[col_a, col_e, col_i, col_big_omega, col_omega, col_m]) + } } impl From for EquinoctialElements { @@ -252,11 +429,7 @@ impl From<&KeplerianElements> for EquinoctialElements { impl fmt::Display for KeplerianElements { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { let rad_to_deg = 180.0 / std::f64::consts::PI; - writeln!( - f, - "Keplerian Elements @ epoch (MJD): {:.6}", - self.reference_epoch - )?; + writeln!(f, "Elements @ epoch (MJD): {:.6}", self.reference_epoch)?; writeln!(f, "-------------------------------------------")?; writeln!( f, @@ -328,3 +501,109 @@ pub(crate) mod test_keplerian_element { ); } } + +#[cfg(test)] +mod tests_jacobian_keplerian { + use super::*; + use approx::assert_abs_diff_eq; + + /// Finite-difference Jacobian for reference. + fn fd_jacobian_kep_to_eq(kep: &KeplerianElements, eps: f64) -> Matrix6 { + let fields: [f64; 6] = [ + kep.semi_major_axis, + kep.eccentricity, + kep.inclination, + kep.ascending_node_longitude, + kep.periapsis_argument, + kep.mean_anomaly, + ]; + + let eq_to_arr = |e: &EquinoctialElements| -> [f64; 6] { + [ + e.semi_major_axis, + e.eccentricity_sin_lon, + e.eccentricity_cos_lon, + e.tan_half_incl_sin_node, + e.tan_half_incl_cos_node, + e.mean_longitude, + ] + }; + + let mut columns = [[0.0f64; 6]; 6]; + + for col in 0..6 { + let mut fwd = fields; + let mut bwd = fields; + fwd[col] += eps; + bwd[col] -= eps; + + let make_kep = |f: [f64; 6]| KeplerianElements { + reference_epoch: kep.reference_epoch, + semi_major_axis: f[0], + eccentricity: f[1], + inclination: f[2], + ascending_node_longitude: f[3], + periapsis_argument: f[4], + mean_anomaly: f[5], + }; + + let eq_fwd = eq_to_arr(&EquinoctialElements::from(&make_kep(fwd))); + let eq_bwd = eq_to_arr(&EquinoctialElements::from(&make_kep(bwd))); + + for row in 0..6 { + columns[col][row] = (eq_fwd[row] - eq_bwd[row]) / (2.0 * eps); + } + } + + Matrix6::from_fn(|row, col| columns[col][row]) + } + + #[test] + fn test_jacobian_to_equinoctial_against_finite_differences() { + let kep = KeplerianElements { + reference_epoch: 60000.0, + semi_major_axis: 2.5, + eccentricity: 0.3, + inclination: 0.5, + ascending_node_longitude: 1.2, + periapsis_argument: 0.8, + mean_anomaly: 2.1, + }; + + let analytical = kep.jacobian_to_equinoctial(); + let numerical = fd_jacobian_kep_to_eq(&kep, 1e-6); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!( + analytical[(row, col)], + numerical[(row, col)], + epsilon = 1e-7 + ); + } + } + } + + /// Near-circular orbit: e ≈ 0, the h/k row should remain well-defined. + #[test] + fn test_jacobian_to_equinoctial_near_circular() { + let kep = KeplerianElements { + reference_epoch: 60000.0, + semi_major_axis: 1.0, + eccentricity: 1e-8, + inclination: 0.3, + ascending_node_longitude: 0.5, + periapsis_argument: 0.2, + mean_anomaly: 1.0, + }; + + let j = kep.jacobian_to_equinoctial(); + let numerical = fd_jacobian_kep_to_eq(&kep, 1e-6); + + for row in 0..6 { + for col in 0..6 { + assert_abs_diff_eq!(j[(row, col)], numerical[(row, col)], epsilon = 1e-6); + } + } + } +} diff --git a/src/orbit_type/mod.rs b/src/orbit_type/mod.rs index ff816a7..4317696 100644 --- a/src/orbit_type/mod.rs +++ b/src/orbit_type/mod.rs @@ -1,19 +1,140 @@ -//! # Orbital element representations +//! Orbital element representations and uncertainty propagation //! -//! This module defines multiple **canonical orbital element sets** and the -//! associated conversions between them: +//! This module defines multiple **canonical orbital element sets** with full support for +//! **uncertainty propagation** during conversions between representations. It provides +//! three distinct parameterizations of Keplerian orbits, each with advantages for different +//! orbital regimes: //! -//! - [`keplerian_element`](crate::orbit_type::keplerian_element) — Classical Keplerian elements `(a, e, i, Ω, ω, M)`, -//! valid for elliptic and hyperbolic orbits. -//! - [`equinoctial_element`](crate::orbit_type::equinoctial_element) — Equinoctial elements `(a, h, k, p, q, λ)`, -//! a **non-singular formulation** well suited for orbit determination near +//! - [`keplerian_element`](crate::orbit_type::keplerian_element) — Classical Keplerian elements $(a, e, i, \Omega, \omega, M)$, +//! widely used for elliptic and hyperbolic orbits but singular for circular and equatorial cases. +//! - [`equinoctial_element`](crate::orbit_type::equinoctial_element) — Equinoctial elements $(a, h, k, p, q, \lambda)$, +//! a **non-singular formulation** ideal for orbit determination and propagation near //! zero eccentricity or inclination. -//! - [`cometary_element`](crate::orbit_type::cometary_element) — Perihelion-based representation `(q, e, i, Ω, ω, ν)`, -//! convenient for parabolic and hyperbolic orbits. +//! - [`cometary_element`](crate::orbit_type::cometary_element) — Perihelion-based representation $(q, e, i, \Omega, \omega, \nu)$, +//! natural for parabolic and hyperbolic orbits where perihelion distance $q$ is better +//! defined than semi-major axis $a$. //! -//! The [`OrbitalElements`](crate::orbit_type::OrbitalElements) enum acts as a **type-erased wrapper** that can hold -//! any of these three representations, while providing uniform constructors and -//! conversion methods. +//! The [`OrbitalElements`] enum acts as a **type-erased wrapper** that can hold any of these +//! three representations along with optional **uncertainty** and **covariance** information, +//! providing uniform constructors and rigorous conversion methods with uncertainty propagation. +//! +//! ## Choosing an orbital element representation +//! +//! Each representation has optimal use cases: +//! +//! **Keplerian elements** $(a, e, i, \Omega, \omega, M)$: +//! - Intuitive physical interpretation: size, shape, orientation, and position in orbit +//! - Standard in classical celestial mechanics and literature +//! - Singular when $e \to 0$ (circular) or $i \to 0$ (equatorial): derivatives become undefined +//! - Best for: communication, visualization, moderate eccentricity/inclination orbits +//! +//! **Equinoctial elements** $(a, h, k, p, q, \lambda)$: +//! - Non-singular for all $e < 1$ and $0 \leq i < \pi$ +//! - Smooth, well-defined derivatives enable robust orbit fitting and least-squares optimization +//! - Ideal for numerical propagation and uncertainty analysis +//! - Best for: orbit determination, propagation, covariance propagation, near-circular or near-equatorial orbits +//! +//! **Cometary elements** $(q, e, i, \Omega, \omega, \nu)$: +//! - Uses perihelion distance $q = a(1 - e)$ instead of semi-major axis +//! - Well-behaved for $e \geq 1$ (parabolic and hyperbolic orbits) +//! - Natural for cometary and hyperbolic trajectories +//! - Best for: comets, interstellar objects, hyperbolic encounters +//! +//! ## Uncertainty representation and propagation +//! +//! Orbital uncertainties arise from measurement errors, numerical approximations, and model limitations. +//! This module represents uncertainties in two complementary ways: +//! +//! 1. **Standard deviations**: individual $1\sigma$ uncertainties on each element +//! (see [`KeplerianUncertainty`](crate::orbit_type::uncertainty::KeplerianUncertainty), [`EquinoctialUncertainty`](crate::orbit_type::uncertainty::EquinoctialUncertainty), [`CometaryUncertainty`](crate::orbit_type::uncertainty::CometaryUncertainty)) +//! +//! 2. **Covariance matrices**: full $6 \times 6$ symmetric positive semi-definite matrices capturing +//! element correlations (see [`OrbitalCovariance`](crate::orbit_type::uncertainty::OrbitalCovariance)) +//! +//! When converting between representations, **covariance matrices are transformed** using the Jacobian +//! of the conversion: +//! +//! $$ +//! \Sigma_y = J \, \Sigma_x \, J^\top +//! $$ +//! +//! where $\Sigma_x$ is the covariance in the source representation, $\Sigma_y$ is the covariance in the +//! target representation, and $J = \partial \mathbf{y} / \partial \mathbf{x}$ is the Jacobian matrix of +//! partial derivatives evaluated at the nominal element values. +//! +//! This **linear covariance propagation** preserves the statistical properties of uncertainties under +//! first-order approximation, which is accurate for small uncertainties relative to the element values. +//! +//! ## Mathematical background: Jacobian matrices +//! +//! Each element representation provides Jacobian methods for conversions: +//! +//! ### Keplerian ↔ Equinoctial +//! +//! The transformation between Keplerian $(a, e, i, \Omega, \omega, M)$ and equinoctial +//! $(a, h, k, p, q, \lambda)$ is defined by: +//! +//! $$ +//! \begin{aligned} +//! h &= e \sin(\Omega + \omega) \\ +//! k &= e \cos(\Omega + \omega) \\ +//! p &= \tan(i/2) \sin(\Omega) \\ +//! q &= \tan(i/2) \cos(\Omega) \\ +//! \lambda &= \Omega + \omega + M +//! \end{aligned} +//! $$ +//! +//! The Jacobian $J_{\text{Kep} \to \text{Eq}}$ is computed analytically in +//! [`KeplerianElements::jacobian_to_equinoctial`], and the inverse in +//! [`EquinoctialElements::jacobian_to_keplerian`]. +//! +//! **Singular cases**: When $e \approx 0$, the longitude of periapsis $\varpi = \Omega + \omega$ +//! is undefined, and derivatives involving $\partial \varpi / \partial h$ and +//! $\partial \varpi / \partial k$ become singular. These are handled by setting the derivatives +//! to zero when $e < \epsilon$ (typically $\epsilon = 10^{-12}$). Similarly, when $i \approx 0$, +//! $\Omega$ is undefined and related derivatives are zeroed. +//! +//! ### Cometary ↔ Keplerian +//! +//! The cometary representation $(q, e, i, \Omega, \omega, \nu)$ relates to Keplerian elements through: +//! +//! $$ +//! a = \frac{q}{1 - e}, \quad M = M(e, \nu) +//! $$ +//! +//! where $M(e, \nu)$ is the mean anomaly computed from the true anomaly $\nu$ via the eccentric +//! anomaly $E$ (for $e < 1$) or hyperbolic anomaly $H$ (for $e > 1$). The Jacobian +//! $J_{\text{Com} \to \text{Kep}}$ is computed in [`CometaryElements::jacobian_to_keplerian`]. +//! +//! For conversions to equinoctial, the **chain rule** is applied: +//! +//! $$ +//! J_{\text{Com} \to \text{Eq}} = J_{\text{Kep} \to \text{Eq}} \cdot J_{\text{Com} \to \text{Kep}} +//! $$ +//! +//! This is implemented in [`CometaryElements::jacobian_to_equinoctial`]. +//! +//! ## Conversion methods with uncertainty propagation +//! +//! The [`OrbitalElements`] enum provides conversion methods that automatically propagate covariance: +//! +//! - [`OrbitalElements::to_keplerian`] — Convert to Keplerian, propagating covariance if present +//! - [`OrbitalElements::to_equinoctial`] — Convert to equinoctial, propagating covariance if present +//! +//! If a covariance matrix is attached to the source elements, it is transformed using the appropriate +//! Jacobian and attached to the result. The 1-σ uncertainties are then recomputed from the diagonal +//! of the transformed covariance. +//! +//! **Error handling**: Conversions that are mathematically undefined (e.g., parabolic cometary → Keplerian, +//! where $a = \infty$) return `Err(OutfitError::InvalidConversion)`. +//! +//! ## Units and conventions +//! +//! - **Lengths**: astronomical units (AU) +//! - **Angles**: radians +//! - **Time**: Modified Julian Date (MJD) in TDB scale for epochs +//! - **Velocities**: AU/day +//! - **Reference frame**: heliocentric ecliptic J2000 //! //! ## Typical workflow //! @@ -28,18 +149,39 @@ //! // Build canonical orbital elements from state //! let elems = OrbitalElements::from_orbital_state(&r, &v, 2460000.5); //! -//! // Convert to Keplerian form if possible +//! // Convert to Keplerian form if possible (uncertainty propagated if present) //! if let Ok(kep) = elems.to_keplerian() { -//! println!("semi-major axis = {}", kep.semi_major_axis); +//! if let OrbitalElements::Keplerian { elements, uncertainty, .. } = kep { +//! println!("semi-major axis = {}", elements.semi_major_axis); +//! if let Some(unc) = uncertainty { +//! println!(" ± {} AU", unc.semi_major_axis); +//! } +//! } //! } //! ``` +//! +//! ## See also +//! +//! - [`KeplerianElements`](crate::orbit_type::keplerian_element::KeplerianElements) — Classical elements with Jacobian methods +//! - [`EquinoctialElements`](crate::orbit_type::equinoctial_element::EquinoctialElements) — Non-singular elements with propagation +//! - [`CometaryElements`](crate::orbit_type::cometary_element::CometaryElements) — Perihelion-based elements for hyperbolic orbits +//! - [`uncertainty`](crate::orbit_type::uncertainty) — Uncertainty structures and covariance propagation +//! +//! ## References +//! +//! - Milani & Gronchi, *Theory of Orbit Determination* (2010), Chapter 5 +//! - Walker et al., "A Set of Modified Equinoctial Orbital Elements", *Celestial Mechanics* 36 (1985) use nalgebra::Vector3; use crate::{ orb_elem::ccek1, orbit_type::{ - cometary_element::CometaryElements, equinoctial_element::EquinoctialElements, + cometary_element::CometaryElements, + equinoctial_element::EquinoctialElements, keplerian_element::KeplerianElements, + uncertainty::{ + CometaryUncertainty, EquinoctialUncertainty, KeplerianUncertainty, OrbitalCovariance, + }, }, OutfitError, }; @@ -53,6 +195,9 @@ pub mod keplerian_element; /// Cometary (parabolic/hyperbolic) orbital elements and related conversions. pub mod cometary_element; +/// Uncertainty structures for orbital elements. +pub mod uncertainty; + /// Canonical orbital elements in multiple representations. /// /// This enum acts as a sum type over several orbital-element parameterizations. @@ -74,9 +219,21 @@ pub mod cometary_element; /// * [`crate::orbit_type::OrbitalElements::to_keplerian`] – Conversion with domain checks. #[derive(Debug, Clone, PartialEq)] pub enum OrbitalElements { - Keplerian(KeplerianElements), - Equinoctial(EquinoctialElements), - Cometary(CometaryElements), + Keplerian { + elements: KeplerianElements, + uncertainty: Option, + covariance: Option, + }, + Equinoctial { + elements: EquinoctialElements, + uncertainty: Option, + covariance: Option, + }, + Cometary { + elements: CometaryElements, + uncertainty: Option, + covariance: Option, + }, } impl OrbitalElements { @@ -121,66 +278,249 @@ impl OrbitalElements { ccek1(position, velocity, reference_epoch) } - /// Convert to Keplerian elements, if possible. + /// Convert to Keplerian representation with full uncertainty propagation /// - /// This conversion may fail if the current representation is not suitable - /// for Keplerian elements. In particular, `Cometary` elements with - /// eccentricity e < 1 cannot be converted to Keplerian form. + /// This method converts the orbital elements to Keplerian form $(a, e, i, \Omega, \omega, M)$ + /// and, if a covariance matrix is present, propagates it through the transformation using + /// the appropriate Jacobian matrix. /// - /// Errors - /// ------ - /// Returns an `OutfitError::InvalidOrbit` if the conversion is not possible. - pub fn to_keplerian(&self) -> Result { + /// ## Uncertainty propagation + /// + /// When the source elements have an attached covariance matrix $\Sigma_x$, this method: + /// + /// 1. Computes the Jacobian $J = \partial \mathbf{y}_{\text{Kep}} / \partial \mathbf{x}_{\text{src}}$ + /// where $\mathbf{x}_{\text{src}}$ is the source parameterization and + /// $\mathbf{y}_{\text{Kep}} = [a, e, i, \Omega, \omega, M]^\top$ + /// + /// 2. Transforms the covariance using **linear covariance propagation**: + /// $$\Sigma_{\text{Kep}} = J \, \Sigma_x \, J^\top$$ + /// + /// 3. Extracts 1-σ uncertainties from the diagonal: $\sigma_i = \sqrt{\Sigma_{ii}}$ + /// + /// This transformation preserves statistical properties under first-order approximation, + /// which is accurate when uncertainties are small relative to element values. + /// + /// ## Conversions + /// + /// - **From Keplerian**: returns a clone (no transformation needed) + /// - **From Equinoctial**: uses [`EquinoctialElements::jacobian_to_keplerian`] + /// - **From Cometary**: uses [`CometaryElements::jacobian_to_keplerian`] + /// + /// ## Return + /// + /// * `Ok(OrbitalElements::Keplerian)` – Converted elements with propagated + /// uncertainty and covariance when available + /// * `Err(OutfitError)` – If the conversion is mathematically undefined + /// (e.g., parabolic cometary elements with $e = 1$, where $a = \infty$) + pub fn to_keplerian(&self) -> Result { match self { - OrbitalElements::Keplerian(ke) => Ok(ke.clone()), - OrbitalElements::Equinoctial(ee) => Ok(KeplerianElements::from(ee)), - OrbitalElements::Cometary(ce) => KeplerianElements::try_from(ce), + OrbitalElements::Keplerian { .. } => Ok(self.clone()), + + OrbitalElements::Equinoctial { + elements, + covariance, + .. + } => { + let kep = KeplerianElements::from(elements); + let jacobian = elements.jacobian_to_keplerian(); + + let new_cov = covariance.as_ref().map(|c| c.propagate(&jacobian)); + let new_unc = new_cov.as_ref().map(KeplerianUncertainty::from_covariance); + + Ok(OrbitalElements::Keplerian { + elements: kep, + uncertainty: new_unc, + covariance: new_cov, + }) + } + + OrbitalElements::Cometary { + elements, + covariance, + .. + } => { + let kep = KeplerianElements::try_from(elements)?; + let jacobian = elements.jacobian_to_keplerian(); + + let new_cov = covariance.as_ref().map(|c| c.propagate(&jacobian)); + let new_unc = new_cov.as_ref().map(KeplerianUncertainty::from_covariance); + + Ok(OrbitalElements::Keplerian { + elements: kep, + uncertainty: new_unc, + covariance: new_cov, + }) + } } } - /// Convert to Equinoctial elements, if possible. + /// Convert to equinoctial representation with full uncertainty propagation /// - /// This conversion may fail if the current representation is not suitable - /// for Equinoctial elements. In particular, `Cometary` elements with - /// eccentricity e < 1 cannot be converted to Equinoctial form. + /// This method converts the orbital elements to equinoctial form $(a, h, k, p, q, \lambda)$ + /// and, if a covariance matrix is present, propagates it through the transformation using + /// the appropriate Jacobian matrix. /// - /// Errors - /// ------ - /// Returns an `OutfitError::InvalidOrbit` if the conversion is not possible. - pub fn to_equinoctial(&self) -> Result { + /// ## Uncertainty propagation + /// + /// When the source elements have an attached covariance matrix $\Sigma_x$, this method: + /// + /// 1. Computes the Jacobian $J = \partial \mathbf{y}_{\text{Eq}} / \partial \mathbf{x}_{\text{src}}$ + /// where $\mathbf{x}_{\text{src}}$ is the source parameterization and + /// $\mathbf{y}_{\text{Eq}} = [a, h, k, p, q, \lambda]^\top$ + /// + /// 2. Transforms the covariance using **linear covariance propagation**: + /// $$\Sigma_{\text{Eq}} = J \, \Sigma_x \, J^\top$$ + /// + /// 3. Extracts 1-σ uncertainties from the diagonal: $\sigma_i = \sqrt{\Sigma_{ii}}$ + /// + /// The equinoctial representation is **non-singular** for $e < 1$ and $0 \leq i < \pi$, + /// making it ideal for uncertainty analysis and propagation when Keplerian elements + /// would be near-singular. + /// + /// ## Conversions + /// + /// - **From Equinoctial**: returns a clone (no transformation needed) + /// - **From Keplerian**: uses [`KeplerianElements::jacobian_to_equinoctial`] + /// - **From Cometary**: uses [`CometaryElements::jacobian_to_equinoctial`] (chain rule through Keplerian) + /// + /// ## Return + /// + /// * `Ok(OrbitalElements::Equinoctial)` – Converted elements with propagated + /// uncertainty and covariance when available + /// * `Err(OutfitError)` – If the conversion is mathematically undefined + pub fn to_equinoctial(&self) -> Result { match self { - OrbitalElements::Keplerian(ke) => Ok(EquinoctialElements::from(ke)), - OrbitalElements::Equinoctial(ee) => Ok(ee.clone()), - OrbitalElements::Cometary(ce) => EquinoctialElements::try_from(ce), + OrbitalElements::Equinoctial { .. } => Ok(self.clone()), + + OrbitalElements::Keplerian { + elements, + covariance, + .. + } => { + let eq = EquinoctialElements::from(elements); + let jacobian = elements.jacobian_to_equinoctial(); + + let new_cov = covariance.as_ref().map(|c| c.propagate(&jacobian)); + let new_unc = new_cov + .as_ref() + .map(EquinoctialUncertainty::from_covariance); + + Ok(OrbitalElements::Equinoctial { + elements: eq, + uncertainty: new_unc, + covariance: new_cov, + }) + } + + OrbitalElements::Cometary { + elements, + covariance, + .. + } => { + let eq = EquinoctialElements::try_from(elements)?; + let jacobian = elements.jacobian_to_equinoctial()?; + + let new_cov = covariance.as_ref().map(|c| c.propagate(&jacobian)); + let new_unc = new_cov + .as_ref() + .map(EquinoctialUncertainty::from_covariance); + + Ok(OrbitalElements::Equinoctial { + elements: eq, + uncertainty: new_unc, + covariance: new_cov, + }) + } } } /// Get a reference to the underlying [`KeplerianElements`] if this is `Keplerian`. - pub fn as_keplerian(&self) -> Option<&KeplerianElements> { - if let OrbitalElements::Keplerian(ref k) = self { - Some(k) + pub fn as_keplerian_ref(&self) -> Option<&KeplerianElements> { + if let OrbitalElements::Keplerian { elements, .. } = self { + Some(elements) + } else { + None + } + } + + /// Get the owned underlying [`KeplerianElements`] if this is `Keplerian`. + pub fn as_keplerian(self) -> Option { + if let OrbitalElements::Keplerian { elements, .. } = self { + Some(elements) } else { None } } /// Get a reference to the underlying [`EquinoctialElements`] if this is `Equinoctial`. - pub fn as_equinoctial(&self) -> Option<&EquinoctialElements> { - if let OrbitalElements::Equinoctial(ref e) = self { - Some(e) + pub fn as_equinoctial_ref(&self) -> Option<&EquinoctialElements> { + if let OrbitalElements::Equinoctial { elements, .. } = self { + Some(elements) + } else { + None + } + } + + /// Get the owned underlying [`EquinoctialElements`] if this is `Equinoctial`. + pub fn as_equinoctial(self) -> Option { + if let OrbitalElements::Equinoctial { elements, .. } = self { + Some(elements) } else { None } } /// Get a reference to the underlying [`CometaryElements`] if this is `Cometary`. - pub fn as_cometary(&self) -> Option<&CometaryElements> { - if let OrbitalElements::Cometary(ref c) = self { - Some(c) + pub fn as_cometary_ref(&self) -> Option<&CometaryElements> { + if let OrbitalElements::Cometary { elements, .. } = self { + Some(elements) + } else { + None + } + } + + /// Get the owned underlying [`CometaryElements`] if this is `Cometary`. + pub fn as_cometary(self) -> Option { + if let OrbitalElements::Cometary { elements, .. } = self { + Some(elements) } else { None } } + + /// Convert to [`KeplerianElements`], propagating covariance if present. + /// + /// Shorthand for `.to_keplerian()?.as_keplerian()`. + /// + /// Return + /// ------ + /// * `Ok(KeplerianElements)` – Converted elements. + /// * `Err(OutfitError)` – If the conversion is not defined for the current + /// element set (e.g. parabolic cometary elements). + pub fn into_keplerian(self) -> Result { + self.to_keplerian()? + .as_keplerian() + .ok_or(OutfitError::InvalidConversion( + "Conversion to Keplerian elements failed".to_string(), + )) + } + + /// Convert to [`EquinoctialElements`], propagating covariance if present. + /// + /// Shorthand for `.to_equinoctial()?.as_equinoctial()`. + /// + /// Return + /// ------ + /// * `Ok(EquinoctialElements)` – Converted elements. + /// * `Err(OutfitError)` – If the conversion is not defined for the current + /// element set (e.g. hyperbolic cometary elements). + pub fn into_equinoctial(self) -> Result { + self.to_equinoctial()? + .as_equinoctial() + .ok_or(OutfitError::InvalidConversion( + "Conversion to equinoctial elements failed".to_string(), + )) + } } use std::fmt; @@ -188,17 +528,17 @@ use std::fmt; impl fmt::Display for OrbitalElements { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { match self { - OrbitalElements::Keplerian(k) => { - writeln!(f, "[Keplerian representation]")?; - write!(f, "{k}") + OrbitalElements::Keplerian { elements, .. } => { + writeln!(f, "[Keplerian]")?; + write!(f, "{elements}") } - OrbitalElements::Equinoctial(e) => { - writeln!(f, "[Equinoctial representation]")?; - write!(f, "{e}") + OrbitalElements::Equinoctial { elements, .. } => { + writeln!(f, "[Equinoctial]")?; + write!(f, "{elements}") } - OrbitalElements::Cometary(c) => { - writeln!(f, "[Cometary representation]")?; - write!(f, "{c}") + OrbitalElements::Cometary { elements, .. } => { + writeln!(f, "[Cometary]")?; + write!(f, "{elements}") } } } @@ -217,7 +557,10 @@ pub(crate) mod orbit_type_test { tol: f64, ) -> bool { match (current, other) { - (OrbitalElements::Keplerian(ke1), OrbitalElements::Keplerian(ke2)) => { + ( + OrbitalElements::Keplerian { elements: ke1, .. }, + OrbitalElements::Keplerian { elements: ke2, .. }, + ) => { abs_diff_eq!(ke1.semi_major_axis, ke2.semi_major_axis, epsilon = tol) && abs_diff_eq!(ke1.eccentricity, ke2.eccentricity, epsilon = tol) && abs_diff_eq!(ke1.inclination, ke2.inclination, epsilon = tol) @@ -233,7 +576,10 @@ pub(crate) mod orbit_type_test { ) && abs_diff_eq!(ke1.mean_anomaly, ke2.mean_anomaly, epsilon = tol) } - (OrbitalElements::Equinoctial(ee1), OrbitalElements::Equinoctial(ee2)) => { + ( + OrbitalElements::Equinoctial { elements: ee1, .. }, + OrbitalElements::Equinoctial { elements: ee2, .. }, + ) => { abs_diff_eq!(ee1.semi_major_axis, 0.0, epsilon = tol) && abs_diff_eq!( ee1.eccentricity_sin_lon, @@ -257,7 +603,10 @@ pub(crate) mod orbit_type_test { ) && abs_diff_eq!(ee1.mean_longitude, ee2.mean_longitude, epsilon = tol) } - (OrbitalElements::Cometary(ce1), OrbitalElements::Cometary(ce2)) => { + ( + OrbitalElements::Cometary { elements: ce1, .. }, + OrbitalElements::Cometary { elements: ce2, .. }, + ) => { abs_diff_eq!( ce1.perihelion_distance, ce2.perihelion_distance, @@ -317,7 +666,7 @@ pub(crate) mod orbit_type_test { let elems = OrbitalElements::from_orbital_state(&r, &v, jd_tdb); match elems { - OrbitalElements::Keplerian(ke) => { + OrbitalElements::Keplerian { elements: ke, .. } => { assert!(ke.semi_major_axis > 0.0); // Loose target for a sanity check, not a golden number. assert_abs_diff_eq!(ke.semi_major_axis, 1.8155, epsilon = 5e-3); @@ -343,7 +692,7 @@ pub(crate) mod orbit_type_test { let elems = OrbitalElements::from_orbital_state(&r, &v, jd_tdb); match elems { - OrbitalElements::Cometary(ce) => { + OrbitalElements::Cometary { elements: ce, .. } => { assert!(ce.eccentricity >= 1.0); assert!(ce.perihelion_distance > 0.0); } @@ -364,11 +713,17 @@ pub(crate) mod orbit_type_test { periapsis_argument: deg(30.0), mean_anomaly: deg(40.0), }; - let oe = OrbitalElements::Keplerian(ke.clone()); + let oe = OrbitalElements::Keplerian { + elements: ke.clone(), + uncertainty: None, + covariance: None, + }; let back = oe .to_keplerian() - .expect("Keplerian -> Keplerian should succeed"); + .expect("Keplerian -> Keplerian should succeed") + .as_keplerian() + .expect("Failed to convert to Keplerian"); assert_abs_diff_eq!(back.semi_major_axis, ke.semi_major_axis, epsilon = 1e-14); assert_abs_diff_eq!(back.eccentricity, ke.eccentricity, epsilon = 1e-14); @@ -397,11 +752,17 @@ pub(crate) mod orbit_type_test { periapsis_argument: deg(35.0), mean_anomaly: deg(45.0), }; - let oe = OrbitalElements::Keplerian(ke.clone()); + let oe = OrbitalElements::Keplerian { + elements: ke.clone(), + uncertainty: None, + covariance: None, + }; let eq = oe .to_equinoctial() - .expect("Keplerian -> Equinoctial should succeed"); + .expect("Keplerian -> Equinoctial should succeed") + .as_equinoctial() + .expect("Failed to convert to Equinoctial"); let ke_back = KeplerianElements::from(&eq); // Use a mix of absolute and relative checks to be robust to scaling. @@ -433,9 +794,17 @@ pub(crate) mod orbit_type_test { mean_anomaly: deg(10.0), }; let eq = EquinoctialElements::from(&ke); - let oe = OrbitalElements::Equinoctial(eq.clone()); + let oe = OrbitalElements::Equinoctial { + elements: eq.clone(), + uncertainty: None, + covariance: None, + }; - let back_via_enum = oe.to_keplerian().expect("Eq -> Kep should succeed"); + let back_via_enum = oe + .to_keplerian() + .expect("Eq -> Kep should succeed") + .as_keplerian() + .expect("Failed to convert to Keplerian"); let back_via_direct = KeplerianElements::from(&eq); assert_abs_diff_eq!( @@ -484,11 +853,17 @@ pub(crate) mod orbit_type_test { periapsis_argument: deg(45.0), true_anomaly: deg(5.0), }; - let oe = OrbitalElements::Cometary(ce); + let oe = OrbitalElements::Cometary { + elements: ce, + uncertainty: None, + covariance: None, + }; let ke = oe .to_keplerian() - .expect("Cometary(e>1) -> Keplerian should succeed for hyperbolic orbits"); + .expect("Cometary(e>1) -> Keplerian should succeed for hyperbolic orbits") + .as_keplerian() + .expect("Failed to convert to Keplerian"); assert!(ke.eccentricity >= 1.0, "expected e >= 1 for hyperbola"); assert!(ke.semi_major_axis < 0.0, "expected a < 0 for hyperbola"); @@ -508,9 +883,16 @@ pub(crate) mod orbit_type_test { periapsis_argument: deg(30.0), true_anomaly: deg(15.0), }; - let oe = OrbitalElements::Cometary(ce); + let oe = OrbitalElements::Cometary { + elements: ce, + uncertainty: None, + covariance: None, + }; if let Ok(eq) = oe.to_equinoctial() { + let eq = eq + .as_equinoctial() + .expect("Failed to convert to Equinoctial after successful conversion"); assert!( eq.semi_major_axis < 0.0, "equinoctial a should be < 0 for hyperbolic orbits" @@ -542,20 +924,32 @@ pub(crate) mod orbit_type_test { true_anomaly: 0.0, }; - let oe_k = OrbitalElements::Keplerian(ke.clone()); - assert!(oe_k.as_keplerian().is_some()); - assert!(oe_k.as_equinoctial().is_none()); - assert!(oe_k.as_cometary().is_none()); - - let oe_e = OrbitalElements::Equinoctial(eq); - assert!(oe_e.as_keplerian().is_none()); - assert!(oe_e.as_equinoctial().is_some()); - assert!(oe_e.as_cometary().is_none()); - - let oe_c = OrbitalElements::Cometary(ce); - assert!(oe_c.as_keplerian().is_none()); - assert!(oe_c.as_equinoctial().is_none()); - assert!(oe_c.as_cometary().is_some()); + let oe_k = OrbitalElements::Keplerian { + elements: ke.clone(), + uncertainty: None, + covariance: None, + }; + assert!(oe_k.as_keplerian_ref().is_some()); + assert!(oe_k.as_equinoctial_ref().is_none()); + assert!(oe_k.as_cometary_ref().is_none()); + + let oe_e = OrbitalElements::Equinoctial { + elements: eq, + uncertainty: None, + covariance: None, + }; + assert!(oe_e.as_keplerian_ref().is_none()); + assert!(oe_e.as_equinoctial_ref().is_some()); + assert!(oe_e.as_cometary_ref().is_none()); + + let oe_c = OrbitalElements::Cometary { + elements: ce, + uncertainty: None, + covariance: None, + }; + assert!(oe_c.as_keplerian_ref().is_none()); + assert!(oe_c.as_equinoctial_ref().is_none()); + assert!(oe_c.as_cometary_ref().is_some()); } // ---------- Display formatting ---------- @@ -582,13 +976,34 @@ pub(crate) mod orbit_type_test { true_anomaly: deg(0.0), }; - let s_k = format!("{}", OrbitalElements::Keplerian(ke)); - assert!(s_k.starts_with("[Keplerian representation]")); - - let s_e = format!("{}", OrbitalElements::Equinoctial(eq)); - assert!(s_e.starts_with("[Equinoctial representation]")); - - let s_c = format!("{}", OrbitalElements::Cometary(ce)); - assert!(s_c.starts_with("[Cometary representation]")); + let s_k = format!( + "{}", + OrbitalElements::Keplerian { + elements: ke, + uncertainty: None, + covariance: None, + } + ); + assert!(s_k.starts_with("[Keplerian]")); + + let s_e = format!( + "{}", + OrbitalElements::Equinoctial { + elements: eq, + uncertainty: None, + covariance: None, + } + ); + assert!(s_e.starts_with("[Equinoctial]")); + + let s_c = format!( + "{}", + OrbitalElements::Cometary { + elements: ce, + uncertainty: None, + covariance: None, + } + ); + assert!(s_c.starts_with("[Cometary]")); } } diff --git a/src/orbit_type/uncertainty.rs b/src/orbit_type/uncertainty.rs new file mode 100644 index 0000000..fad4362 --- /dev/null +++ b/src/orbit_type/uncertainty.rs @@ -0,0 +1,417 @@ +//! Orbital uncertainty representation and covariance propagation +//! +//! This module defines uncertainty structures for orbital elements and provides mechanisms +//! for propagating uncertainties through coordinate transformations. +//! +//! ## Overview +//! +//! Orbital uncertainties arise from multiple sources: +//! - **Measurement errors** in astrometric observations (position, timing) +//! - **Numerical approximations** in orbit fitting and propagation +//! - **Model limitations** (unmodeled perturbations, approximations in dynamics) +//! +//! These uncertainties are represented in two complementary forms: +//! +//! 1. **Standard deviations** ($1\sigma$ uncertainties) — individual uncertainties on each element, +//! extracted from the diagonal of the covariance matrix +//! 2. **Covariance matrices** — full $6 \times 6$ symmetric matrices capturing correlations between elements +//! +//! ## Uncertainty structures +//! +//! Three uncertainty structures correspond to the three orbital element representations: +//! +//! - [`KeplerianUncertainty`] — for classical Keplerian elements $(a, e, i, \Omega, \omega, M)$ +//! - [`EquinoctialUncertainty`] — for equinoctial elements $(a, h, k, p, q, \lambda)$ +//! - [`CometaryUncertainty`] — for cometary elements $(q, e, i, \Omega, \omega, \nu)$ +//! +//! Each structure stores the standard deviation $\sigma_i$ for each element $x_i$, extracted from +//! the diagonal of the covariance matrix: $\sigma_i = \sqrt{\Sigma_{ii}}$. +//! +//! All units match those of the parent element struct (AU for distances, radians for angles). +//! +//! ## Covariance matrix representation +//! +//! The [`OrbitalCovariance`] structure holds a full $6 \times 6$ **symmetric positive semi-definite** +//! covariance matrix $\Sigma$. The matrix element $\Sigma_{ij}$ represents the covariance between +//! orbital elements $x_i$ and $x_j$: +//! +//! $$ +//! \Sigma_{ij} = \text{Cov}(x_i, x_j) = E[(x_i - \mu_i)(x_j - \mu_j)] +//! $$ +//! +//! where $\mu_i$ is the expected value (nominal element value), with $\mu_i = \mathbb{E}(x_i)$ and $\mathbb{E}(\cdot)$ denoting expectation. +//! +//! **Diagonal entries** $\Sigma_{ii} = \sigma_i^2$ are the variances of each element. +//! +//! **Off-diagonal entries** $\Sigma_{ij}$ ($i \neq j$) capture correlations. The correlation coefficient is: +//! +//! $$ +//! \rho_{ij} = \frac{\Sigma_{ij}}{\sigma_i \sigma_j} \in [-1, 1] +//! $$ +//! +//! ## Linear covariance propagation +//! +//! When transforming orbital elements from one representation to another, the covariance matrix +//! is propagated using **first-order linear approximation**: +//! +//! $$ +//! \Sigma_y = J \, \Sigma_x \, J^\top +//! $$ +//! +//! where: +//! - $\Sigma_x$ is the covariance in the source representation $\mathbf{x}$ +//! - $\Sigma_y$ is the covariance in the target representation $\mathbf{y}$ +//! - $J = \frac{\partial \mathbf{y}}{\partial \mathbf{x}}$ is the $6 \times 6$ Jacobian matrix +//! evaluated at the nominal element values +//! +//! This transformation preserves the **statistical properties** of the uncertainty distribution +//! under the assumption that: +//! 1. The transformation $\mathbf{y} = f(\mathbf{x})$ is smooth and differentiable +//! 2. Uncertainties are small relative to element values (linear approximation valid) +//! 3. The distribution is approximately Gaussian (covariance fully characterizes uncertainty) +//! +//! ## Jacobian computation +//! +//! Each orbital element representation provides methods to compute analytical Jacobians: +//! +//! - [`KeplerianElements::jacobian_to_equinoctial`](crate::orbit_type::keplerian_element::KeplerianElements::jacobian_to_equinoctial) +//! — $J_{\text{Kep} \to \text{Eq}} = \partial(a,h,k,p,q,\lambda) / \partial(a,e,i,\Omega,\omega,M)$ +//! +//! - [`EquinoctialElements::jacobian_to_keplerian`](crate::orbit_type::equinoctial_element::EquinoctialElements::jacobian_to_keplerian) +//! — $J_{\text{Eq} \to \text{Kep}} = \partial(a,e,i,\Omega,\omega,M) / \partial(a,h,k,p,q,\lambda)$ +//! +//! - [`CometaryElements::jacobian_to_keplerian`](crate::orbit_type::cometary_element::CometaryElements::jacobian_to_keplerian) +//! — $J_{\text{Com} \to \text{Kep}} = \partial(a,e,i,\Omega,\omega,M) / \partial(q,e,i,\Omega,\omega,\nu)$ +//! +//! - [`CometaryElements::jacobian_to_equinoctial`](crate::orbit_type::cometary_element::CometaryElements::jacobian_to_equinoctial) +//! — $J_{\text{Com} \to \text{Eq}}$ computed via chain rule +//! +//! These Jacobians are computed analytically for numerical accuracy and efficiency. +//! +//! ## Handling singularities +//! +//! Some orbital element representations have **singularities** where derivatives are undefined: +//! +//! - **Keplerian elements**: singular when $e \to 0$ (circular) or $i \to 0$ (equatorial) +//! - **Equinoctial elements**: non-singular for $e < 1$ and $0 \leq i < \pi$ +//! +//! When computing Jacobians near singular points, derivatives involving undefined quantities +//! (e.g., $\partial \varpi / \partial h$ when $e = 0$) are set to zero. This ensures numerical +//! stability but may underestimate uncertainties in degenerate configurations. +//! +//! For orbits near circular or equatorial, **equinoctial elements** should be preferred as the +//! primary representation to avoid singularities in both the transformation and its Jacobian. +//! +//! ## Usage +//! +//! Uncertainties are typically created by orbit determination codes and attached to +//! [`OrbitalElements`](crate::orbit_type::OrbitalElements) variants. When converting between +//! representations, the covariance is automatically propagated: +//! +//! ```rust +//! # use outfit::orbit_type::{OrbitalElements, keplerian_element::KeplerianElements}; +//! # use outfit::orbit_type::uncertainty::{KeplerianUncertainty, OrbitalCovariance}; +//! # use nalgebra::Matrix6; +//! // Assume we have Keplerian elements with uncertainty +//! let kep = KeplerianElements { +//! reference_epoch: 60000.0, +//! semi_major_axis: 2.5, +//! eccentricity: 0.1, +//! inclination: 0.2, +//! ascending_node_longitude: 1.0, +//! periapsis_argument: 0.5, +//! mean_anomaly: 2.0, +//! }; +//! +//! // With an associated covariance (example: diagonal matrix) +//! let cov = OrbitalCovariance { +//! matrix: Matrix6::from_diagonal_element(1e-6), +//! }; +//! +//! let oe = OrbitalElements::Keplerian { +//! elements: kep, +//! uncertainty: Some(KeplerianUncertainty::from_covariance(&cov)), +//! covariance: Some(cov), +//! }; +//! +//! // Convert to equinoctial — covariance automatically propagated +//! let oe_eq = oe.to_equinoctial().unwrap(); +//! +//! // Uncertainty in equinoctial representation now available +//! if let OrbitalElements::Equinoctial { uncertainty, .. } = oe_eq { +//! if let Some(unc) = uncertainty { +//! println!("Uncertainty in h: {}", unc.eccentricity_sin_lon); +//! } +//! } +//! ``` +//! +//! ## See also +//! +//! - [`OrbitalElements`](crate::orbit_type::OrbitalElements) — Container with uncertainty +//! - [`OrbitalCovariance::propagate`] — Covariance transformation method +//! - Module-level documentation in [`crate::orbit_type`] for conversion details +//! +//! ## References +//! +//! - Tapley, Schutz, & Born, *Statistical Orbit Determination* (2004), Chapter 4 +//! - Milani & Gronchi, *Theory of Orbit Determination* (2010), Chapter 5 +use nalgebra::Matrix6; +use photom::Radians; + +/// One-sigma uncertainties on Keplerian elements. +/// +/// Units +/// ----- +/// * `semi_major_axis`: AU. +/// * `eccentricity`: unitless. +/// * `inclination`: radians. +/// * `ascending_node_longitude`: radians. +/// * `periapsis_argument`: radians. +/// * `mean_anomaly`: radians. +/// +/// Notes +/// ----- +/// The reference epoch is treated as exact (no uncertainty propagated here). +/// +/// See also +/// -------- +/// * [`crate::orbit_type::keplerian_element::KeplerianElements`] – Associated element struct. +/// * [`OrbitalCovariance`] – Full 6×6 covariance when available. +#[derive(Debug, Clone, PartialEq)] +pub struct KeplerianUncertainty { + pub semi_major_axis: f64, + pub eccentricity: f64, + pub inclination: Radians, + pub ascending_node_longitude: Radians, + pub periapsis_argument: Radians, + pub mean_anomaly: Radians, +} + +/// One-sigma uncertainties on equinoctial elements. +/// +/// Units +/// ----- +/// * `semi_major_axis`: AU. +/// * `eccentricity_sin_lon`: unitless. +/// * `eccentricity_cos_lon`: unitless. +/// * `tan_half_incl_sin_node`: unitless. +/// * `tan_half_incl_cos_node`: unitless. +/// * `mean_longitude`: radians. +/// +/// See also +/// -------- +/// * [`crate::orbit_type::equinoctial_element::EquinoctialElements`] – Associated element struct. +/// * [`OrbitalCovariance`] – Full 6×6 covariance when available. +#[derive(Debug, Clone, PartialEq)] +pub struct EquinoctialUncertainty { + pub semi_major_axis: f64, + pub eccentricity_sin_lon: f64, + pub eccentricity_cos_lon: f64, + pub tan_half_incl_sin_node: f64, + pub tan_half_incl_cos_node: f64, + pub mean_longitude: Radians, +} + +/// One-sigma uncertainties on cometary elements. +/// +/// Units +/// ----- +/// * `perihelion_distance`: AU. +/// * `eccentricity`: unitless. +/// * `inclination`: radians. +/// * `ascending_node_longitude`: radians. +/// * `periapsis_argument`: radians. +/// * `true_anomaly`: radians. +/// +/// See also +/// -------- +/// * [`crate::orbit_type::cometary_element::CometaryElements`] – Associated element struct. +/// * [`OrbitalCovariance`] – Full 6×6 covariance when available. +#[derive(Debug, Clone, PartialEq)] +pub struct CometaryUncertainty { + pub perihelion_distance: f64, + pub eccentricity: f64, + pub inclination: Radians, + pub ascending_node_longitude: Radians, + pub periapsis_argument: Radians, + pub true_anomaly: Radians, +} + +impl KeplerianUncertainty { + /// Extract 1-σ standard deviations from the diagonal of a covariance matrix. + /// + /// Element ordering: $[a, e, i, \Omega, \omega, M]$. + pub fn from_covariance(cov: &OrbitalCovariance) -> Self { + let v = cov.variances(); + KeplerianUncertainty { + semi_major_axis: v[0].sqrt(), + eccentricity: v[1].sqrt(), + inclination: v[2].sqrt(), + ascending_node_longitude: v[3].sqrt(), + periapsis_argument: v[4].sqrt(), + mean_anomaly: v[5].sqrt(), + } + } +} + +impl EquinoctialUncertainty { + /// Extract 1-σ standard deviations from the diagonal of a covariance matrix. + /// + /// Element ordering: $[a, h, k, p, q, \lambda]$. + pub fn from_covariance(cov: &OrbitalCovariance) -> Self { + let v = cov.variances(); + EquinoctialUncertainty { + semi_major_axis: v[0].sqrt(), + eccentricity_sin_lon: v[1].sqrt(), + eccentricity_cos_lon: v[2].sqrt(), + tan_half_incl_sin_node: v[3].sqrt(), + tan_half_incl_cos_node: v[4].sqrt(), + mean_longitude: v[5].sqrt(), + } + } +} + +impl CometaryUncertainty { + /// Extract 1-σ standard deviations from the diagonal of a covariance matrix. + /// + /// Element ordering: $[q, e, i, \Omega, \omega, \nu]$. + pub fn from_covariance(cov: &OrbitalCovariance) -> Self { + let v = cov.variances(); + CometaryUncertainty { + perihelion_distance: v[0].sqrt(), + eccentricity: v[1].sqrt(), + inclination: v[2].sqrt(), + ascending_node_longitude: v[3].sqrt(), + periapsis_argument: v[4].sqrt(), + true_anomaly: v[5].sqrt(), + } + } +} + +/// Full 6×6 covariance matrix for a set of orbital elements. +/// +/// The matrix is stored as a flat array of 36 elements in **row-major** order. +/// The parameterization (Keplerian, equinoctial, or cometary) is determined by +/// the variant it accompanies inside [`crate::orbit_type::OrbitalElements`]. +/// +/// The matrix is assumed to be **symmetric positive semi-definite**. +/// Only the upper triangle is guaranteed to be written by solvers; the lower +/// triangle is kept consistent by convention. +/// +/// Element ordering +/// ---------------- +/// For Keplerian elements: `[a, e, i, Ω, ω, M]`. +/// For equinoctial elements: `[a, h, k, p, q, λ]`. +/// For cometary elements: `[q, e, i, Ω, ω, ν]`. +/// +/// See also +/// -------- +/// * [`crate::orbit_type::OrbitalElements`] – Container that holds this matrix. +/// * [`KeplerianUncertainty`], [`EquinoctialUncertainty`], [`CometaryUncertainty`] – Diagonal 1-σ summaries. +#[derive(Debug, Clone, PartialEq)] +pub struct OrbitalCovariance { + /// Row-major 6×6 covariance matrix. + pub matrix: Matrix6, +} + +impl OrbitalCovariance { + /// Returns the diagonal entries, i.e., the variances of each element. + /// + /// The square roots of the returned values equal the 1-σ standard + /// deviations stored in the matching uncertainty struct + /// ([`KeplerianUncertainty`], [`EquinoctialUncertainty`], or [`CometaryUncertainty`]). + #[inline] + pub fn variances(&self) -> [f64; 6] { + [ + self.matrix[(0, 0)], + self.matrix[(1, 1)], + self.matrix[(2, 2)], + self.matrix[(3, 3)], + self.matrix[(4, 4)], + self.matrix[(5, 5)], + ] + } + + /// Propagate this covariance matrix through a coordinate transformation + /// + /// This method implements **linear covariance propagation**, transforming uncertainties + /// from one orbital element representation to another using the Jacobian matrix of the + /// transformation. + /// + /// ## Mathematical formulation + /// + /// Given a transformation $\mathbf{y} = f(\mathbf{x})$ where $\mathbf{x}$ and $\mathbf{y}$ + /// are 6-element orbital parameter vectors, the covariance in the target space is: + /// + /// $$ + /// \Sigma_y = J \, \Sigma_x \, J^\top + /// $$ + /// + /// where: + /// - $\Sigma_x$ is the $6 \times 6$ covariance matrix in the source representation (this matrix) + /// - $\Sigma_y$ is the $6 \times 6$ covariance matrix in the target representation (returned) + /// - $J = \frac{\partial \mathbf{y}}{\partial \mathbf{x}} \bigg|_{\mathbf{x}_0}$ + /// is the Jacobian matrix evaluated at the nominal element values $\mathbf{x}_0$ + /// + /// ## Derivation + /// + /// Under a first-order Taylor expansion around the nominal point $\mathbf{x}_0$: + /// + /// $$ + /// \mathbf{y} \approx f(\mathbf{x}_0) + J(\mathbf{x} - \mathbf{x}_0) + /// $$ + /// + /// The covariance of $\mathbf{y}$ is: + /// + /// $$ + /// \begin{aligned} + /// \Sigma_y &= E[(\mathbf{y} - E[\mathbf{y}])(\mathbf{y} - E[\mathbf{y}])^\top] \\ + /// &= E[J(\mathbf{x} - \mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0)^\top J^\top] \\ + /// &= J \, E[(\mathbf{x} - \mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0)^\top] \, J^\top \\ + /// &= J \, \Sigma_x \, J^\top + /// ) -> OrbitalCovariance { + OrbitalCovariance { + matrix: jacobian * self.matrix * jacobian.transpose(), + } + } +} diff --git a/src/trajectory.rs b/src/trajectory.rs index 3757eed..54151dc 100644 --- a/src/trajectory.rs +++ b/src/trajectory.rs @@ -488,7 +488,13 @@ impl TrajectoryFit for Vec<&Observation> { }; // 4.b) Convert to the element set required by the scorer. - let equinoctial_elements = gauss_res.get_orbit().to_equinoctial()?; + let equinoctial_elements = gauss_res + .get_orbit() + .to_equinoctial()? + .as_equinoctial() + .ok_or(OutfitError::InvalidConversion( + "Conversion to equinoctial elements failed".to_string(), + ))?; // 4.c) Score orbit vs. full observation set (RMS residual). let rms = match self.rms_orbit_error( @@ -817,15 +823,19 @@ mod test_obs_ext { let orbit = best_orbit.orbital_elements(); - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 57049.22904475403, - semi_major_axis: 1.8017609974509807, - eccentricity: 0.2835733667643381, - inclination: 0.20267686119302475, - ascending_node_longitude: 0.00799201841873464, - periapsis_argument: 1.245034216916367, - mean_anomaly: 0.4405089048961484, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 57049.22904475403, + semi_major_axis: 1.8017609974509807, + eccentricity: 0.2835733667643381, + inclination: 0.20267686119302475, + ascending_node_longitude: 0.00799201841873464, + periapsis_argument: 1.245034216916367, + mean_anomaly: 0.4405089048961484, + }, + uncertainty: None, + covariance: None, + }; assert!(approx_equal(orbit, &expected_orbit, 1e-14)); assert_relative_eq!( diff --git a/tests/common/mod.rs b/tests/common/mod.rs index 2d75f21..96b5918 100644 --- a/tests/common/mod.rs +++ b/tests/common/mod.rs @@ -1,11 +1,22 @@ use approx::abs_diff_eq; -use outfit::orbit_type::OrbitalElements; +use outfit::{orbit_type::uncertainty::OrbitalCovariance, OrbitalElements}; #[allow(dead_code)] pub fn approx_equal(current: &OrbitalElements, other: &OrbitalElements, tol: f64) -> bool { match (current, other) { - (OrbitalElements::Keplerian(ke1), OrbitalElements::Keplerian(ke2)) => { - abs_diff_eq!(ke1.semi_major_axis, ke2.semi_major_axis, epsilon = tol) + ( + OrbitalElements::Keplerian { + elements: ke1, + uncertainty: unc1, + covariance: cov1, + }, + OrbitalElements::Keplerian { + elements: ke2, + uncertainty: unc2, + covariance: cov2, + }, + ) => { + let elements_eq = abs_diff_eq!(ke1.semi_major_axis, ke2.semi_major_axis, epsilon = tol) && abs_diff_eq!(ke1.eccentricity, ke2.eccentricity, epsilon = tol) && abs_diff_eq!(ke1.inclination, ke2.inclination, epsilon = tol) && abs_diff_eq!( @@ -18,10 +29,43 @@ pub fn approx_equal(current: &OrbitalElements, other: &OrbitalElements, tol: f64 ke2.periapsis_argument, epsilon = tol ) - && abs_diff_eq!(ke1.mean_anomaly, ke2.mean_anomaly, epsilon = tol) + && abs_diff_eq!(ke1.mean_anomaly, ke2.mean_anomaly, epsilon = tol); + + let uncertainty_eq = match (unc1, unc2) { + (None, None) => true, + (Some(u1), Some(u2)) => { + abs_diff_eq!(u1.semi_major_axis, u2.semi_major_axis, epsilon = tol) + && abs_diff_eq!(u1.eccentricity, u2.eccentricity, epsilon = tol) + && abs_diff_eq!(u1.inclination, u2.inclination, epsilon = tol) + && abs_diff_eq!( + u1.ascending_node_longitude, + u2.ascending_node_longitude, + epsilon = tol + ) + && abs_diff_eq!(u1.periapsis_argument, u2.periapsis_argument, epsilon = tol) + && abs_diff_eq!(u1.mean_anomaly, u2.mean_anomaly, epsilon = tol) + } + _ => false, + }; + + let covariance_eq = approx_equal_covariance(cov1, cov2, tol); + + elements_eq && uncertainty_eq && covariance_eq } - (OrbitalElements::Equinoctial(ee1), OrbitalElements::Equinoctial(ee2)) => { - abs_diff_eq!(ee1.semi_major_axis, ee2.semi_major_axis, epsilon = tol) + + ( + OrbitalElements::Equinoctial { + elements: ee1, + uncertainty: unc1, + covariance: cov1, + }, + OrbitalElements::Equinoctial { + elements: ee2, + uncertainty: unc2, + covariance: cov2, + }, + ) => { + let elements_eq = abs_diff_eq!(ee1.semi_major_axis, ee2.semi_major_axis, epsilon = tol) && abs_diff_eq!( ee1.eccentricity_sin_lon, ee2.eccentricity_sin_lon, @@ -42,10 +86,55 @@ pub fn approx_equal(current: &OrbitalElements, other: &OrbitalElements, tol: f64 ee2.tan_half_incl_cos_node, epsilon = tol ) - && abs_diff_eq!(ee1.mean_longitude, ee2.mean_longitude, epsilon = tol) + && abs_diff_eq!(ee1.mean_longitude, ee2.mean_longitude, epsilon = tol); + + let uncertainty_eq = match (unc1, unc2) { + (None, None) => true, + (Some(u1), Some(u2)) => { + abs_diff_eq!(u1.semi_major_axis, u2.semi_major_axis, epsilon = tol) + && abs_diff_eq!( + u1.eccentricity_sin_lon, + u2.eccentricity_sin_lon, + epsilon = tol + ) + && abs_diff_eq!( + u1.eccentricity_cos_lon, + u2.eccentricity_cos_lon, + epsilon = tol + ) + && abs_diff_eq!( + u1.tan_half_incl_sin_node, + u2.tan_half_incl_sin_node, + epsilon = tol + ) + && abs_diff_eq!( + u1.tan_half_incl_cos_node, + u2.tan_half_incl_cos_node, + epsilon = tol + ) + && abs_diff_eq!(u1.mean_longitude, u2.mean_longitude, epsilon = tol) + } + _ => false, + }; + + let covariance_eq = approx_equal_covariance(cov1, cov2, tol); + + elements_eq && uncertainty_eq && covariance_eq } - (OrbitalElements::Cometary(ce1), OrbitalElements::Cometary(ce2)) => { - abs_diff_eq!( + + ( + OrbitalElements::Cometary { + elements: ce1, + uncertainty: unc1, + covariance: cov1, + }, + OrbitalElements::Cometary { + elements: ce2, + uncertainty: unc2, + covariance: cov2, + }, + ) => { + let elements_eq = abs_diff_eq!( ce1.perihelion_distance, ce2.perihelion_distance, epsilon = tol @@ -61,8 +150,55 @@ pub fn approx_equal(current: &OrbitalElements, other: &OrbitalElements, tol: f64 ce2.periapsis_argument, epsilon = tol ) - && abs_diff_eq!(ce1.true_anomaly, ce2.true_anomaly, epsilon = tol) + && abs_diff_eq!(ce1.true_anomaly, ce2.true_anomaly, epsilon = tol); + + let uncertainty_eq = match (unc1, unc2) { + (None, None) => true, + (Some(u1), Some(u2)) => { + abs_diff_eq!( + u1.perihelion_distance, + u2.perihelion_distance, + epsilon = tol + ) && abs_diff_eq!(u1.eccentricity, u2.eccentricity, epsilon = tol) + && abs_diff_eq!(u1.inclination, u2.inclination, epsilon = tol) + && abs_diff_eq!( + u1.ascending_node_longitude, + u2.ascending_node_longitude, + epsilon = tol + ) + && abs_diff_eq!(u1.periapsis_argument, u2.periapsis_argument, epsilon = tol) + && abs_diff_eq!(u1.true_anomaly, u2.true_anomaly, epsilon = tol) + } + _ => false, + }; + + let covariance_eq = approx_equal_covariance(cov1, cov2, tol); + + elements_eq && uncertainty_eq && covariance_eq } - _ => false, // Different types cannot be equal + + _ => false, + } +} + +/// Compare two optional [`OrbitalCovariance`] matrices entry-wise within `tol`. +/// +/// * `(None, None)` → `true` +/// * `(Some, Some)` → all 36 entries within `tol` +/// * mixed → `false` +#[allow(dead_code)] +fn approx_equal_covariance( + cov1: &Option, + cov2: &Option, + tol: f64, +) -> bool { + match (cov1, cov2) { + (None, None) => true, + (Some(c1), Some(c2)) => c1 + .matrix + .iter() + .zip(c2.matrix.iter()) + .all(|(a, b)| abs_diff_eq!(a, b, epsilon = tol)), + _ => false, } } diff --git a/tests/test_diff_cor.rs b/tests/test_diff_cor.rs index be26c06..d1fdc9c 100644 --- a/tests/test_diff_cor.rs +++ b/tests/test_diff_cor.rs @@ -6,6 +6,7 @@ use hifitime::ut1::Ut1Provider; use outfit::jpl_ephem::naif::naif_ids::{ planet_bary::PlanetaryBary, solar_system_bary::SolarSystemBary, NaifIds, }; +use outfit::orbit_type::uncertainty::{EquinoctialUncertainty, OrbitalCovariance}; use outfit::{ orbit_type::{equinoctial_element::EquinoctialElements, OrbitalElements}, propagator::{NBodyConfig, PropagatorKind}, @@ -96,15 +97,78 @@ fn test_diff_cor() { .as_ref() .expect("K09R05F (2015AB) should converge"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 57049.2684537375, - semi_major_axis: 1.801837227645679, - eccentricity_sin_lon: 0.26941036025991355, - eccentricity_cos_lon: 0.08909600747061494, - tan_half_incl_sin_node: 0.0008708024189761142, - tan_half_incl_cos_node: 0.10166598640878513, - mean_longitude: 1.6929834276945714, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 57049.2684537375, + semi_major_axis: 1.801837227645679, + eccentricity_sin_lon: 0.26941036025991355, + eccentricity_cos_lon: 0.08909600747061494, + tan_half_incl_sin_node: 0.0008708024189761142, + tan_half_incl_cos_node: 0.10166598640878513, + mean_longitude: 1.6929834276945714, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 1.3935756201273647e-6, + eccentricity_sin_lon: 2.399103573371585e-6, + eccentricity_cos_lon: 9.380584628466963e-6, + tan_half_incl_sin_node: 4.2486965596206456e-7, + tan_half_incl_cos_node: 9.938054593077774e-7, + mean_longitude: 1.5699462542222023e-5, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 1.942053009013369e-12, + -3.7365542822268565e-13, + 1.250111987715944e-11, + -3.8069560012308287e-13, + 5.495356218939393e-13, + -2.1061628726935973e-11, + ], + [ + -3.736554282226888e-13, + 5.7556979557643085e-12, + -8.919579576942644e-12, + 6.829258011452513e-13, + -2.190283688325579e-12, + 1.4156679672214094e-11, + ], + [ + 1.2501119877159442e-11, + -8.919579576942621e-12, + 8.799536797183067e-11, + -3.157563107997367e-12, + 5.930188854586023e-12, + -1.472073140503015e-10, + ], + [ + -3.806956001230829e-13, + 6.829258011452509e-13, + -3.157563107997368e-12, + 1.8051422455732311e-13, + -3.5751562142662264e-13, + 5.229181995216352e-12, + ], + [ + 5.495356218939391e-13, + -2.1902836883255787e-12, + 5.930188854586025e-12, + -3.5751562142662264e-13, + 9.876492909499423e-13, + -9.67328953098736e-12, + ], + [ + -2.1061628726935976e-11, + 1.4156679672214063e-11, + -1.472073140503015e-10, + 5.229181995216351e-12, + -9.673289530987361e-12, + 2.464731241146324e-10, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), nr_tol), @@ -123,15 +187,78 @@ fn test_diff_cor() { .as_ref() .expect("33803 should converge"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 60465.26777915681, - semi_major_axis: 2.190614169340076, - eccentricity_sin_lon: -0.13393967896355405, - eccentricity_cos_lon: 0.1533932583177835, - tan_half_incl_sin_node: 0.002997272576917091, - tan_half_incl_cos_node: -0.05948928702443621, - mean_longitude: 4.224671691074116, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 60465.26777915681, + semi_major_axis: 2.190614169340076, + eccentricity_sin_lon: -0.13393967896355405, + eccentricity_cos_lon: 0.1533932583177835, + tan_half_incl_sin_node: 0.002997272576917091, + tan_half_incl_cos_node: -0.05948928702443621, + mean_longitude: 4.224671691074116, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 2.1400421559849134e-5, + eccentricity_sin_lon: 1.364670439647764e-5, + eccentricity_cos_lon: 5.318530114145479e-6, + tan_half_incl_sin_node: 3.44968775225327e-7, + tan_half_incl_cos_node: 8.503880052285401e-7, + mean_longitude: 2.664301205078454e-5, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 4.5797804293925557e-10, + -2.443785426064791e-10, + 7.203221689097433e-11, + -1.883169629832777e-12, + -6.3279112379918766e-12, + 4.3441160814862357e-10, + ], + [ + -2.443785426064796e-10, + 1.8623254088484216e-10, + -6.032986816763725e-11, + 7.999773867024745e-15, + -6.598752075412107e-13, + -3.5829528431457476e-10, + ], + [ + 7.203221689097439e-11, + -6.032986816763721e-11, + 2.8286762575072326e-11, + 2.0398130597296797e-14, + 1.4218640626998597e-13, + 1.2758725519460455e-10, + ], + [ + -1.883169629832779e-12, + 7.99977386702494e-15, + 2.0398130597296844e-14, + 1.190034558804622e-13, + 2.64333826423024e-13, + 3.756599803475119e-13, + ], + [ + -6.327911237991877e-12, + -6.598752075412104e-13, + 1.4218640626998607e-13, + 2.64333826423024e-13, + 7.231597594365756e-13, + 2.605687909220327e-12, + ], + [ + 4.3441160814862383e-10, + -3.582952843145747e-10, + 1.2758725519460457e-10, + 3.7565998034751195e-13, + 2.6056879092203274e-12, + 7.098500911382502e-10, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), nr_tol), @@ -150,15 +277,78 @@ fn test_diff_cor() { .as_ref() .expect("8467 should converge with rms_divergence_ratio=10"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 60672.2443617134, - semi_major_axis: 3.2073734821020743, - eccentricity_sin_lon: 0.053597752212361474, - eccentricity_cos_lon: -0.023229330026225303, - tan_half_incl_sin_node: 0.0028890355813102732, - tan_half_incl_cos_node: 0.09179492536540514, - mean_longitude: 0.626741395885302, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 60672.2443617134, + semi_major_axis: 3.2073734821020743, + eccentricity_sin_lon: 0.053597752212361474, + eccentricity_cos_lon: -0.023229330026225303, + tan_half_incl_sin_node: 0.0028890355813102732, + tan_half_incl_cos_node: 0.09179492536540514, + mean_longitude: 0.626741395885302, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 0.00758317975106881, + eccentricity_sin_lon: 0.002478406542589576, + eccentricity_cos_lon: 0.0007443879537814839, + tan_half_incl_sin_node: 4.277383244080703e-5, + tan_half_incl_cos_node: 5.706392699913953e-5, + mean_longitude: 0.00333399562783862, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 5.750461513702002e-5, + 1.8729896457450725e-5, + 5.604248768814215e-6, + -3.2370073744381016e-7, + -4.297318085854602e-7, + 2.504633450274609e-5, + ], + [ + 1.8729896457450735e-5, + 6.1424989903508165e-6, + 1.8071841318216132e-6, + -1.0560687892019813e-7, + -1.409247502206143e-7, + 8.250952263039232e-6, + ], + [ + 5.604248768814217e-6, + 1.807184131821612e-6, + 5.541134257349846e-7, + -3.14728840772654e-8, + -4.14717463955493e-8, + 2.4005716002617356e-6, + ], + [ + -3.237007374438101e-7, + -1.0560687892019811e-7, + -3.147288407726542e-8, + 1.8296007416742358e-9, + 2.435346888714026e-9, + -1.4137265325860534e-7, + ], + [ + -4.2973180858546056e-7, + -1.4092475022061433e-7, + -4.1471746395549346e-8, + 2.4353468887140264e-9, + 3.2562917645631254e-9, + -1.8928599918199224e-7, + ], + [ + 2.50463345027461e-5, + 8.250952263039232e-6, + 2.400571600261738e-6, + -1.4137265325860537e-7, + -1.8928599918199224e-7, + 1.1115526846447033e-5, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), nr_tol), @@ -236,7 +426,7 @@ fn test_diff_cor_nbody() { .expect("K09R05F should converge under N-body propagation"); let elem = match orbit.orbital_elements() { - OrbitalElements::Equinoctial(e) => e, + OrbitalElements::Equinoctial { elements, .. } => elements, _ => panic!("Expected equinoctial elements for K09R05F"), }; @@ -271,7 +461,7 @@ fn test_diff_cor_nbody() { .expect("33803 should converge under N-body propagation"); let elem = match orbit.orbital_elements() { - OrbitalElements::Equinoctial(e) => e, + OrbitalElements::Equinoctial { elements, .. } => elements, _ => panic!("Expected equinoctial elements for 33803"), }; @@ -303,7 +493,7 @@ fn test_diff_cor_nbody() { .expect("8467 should converge under N-body propagation"); let elem = match orbit.orbital_elements() { - OrbitalElements::Equinoctial(e) => e, + OrbitalElements::Equinoctial { elements, .. } => elements, _ => panic!("Expected equinoctial elements for 8467"), }; @@ -398,15 +588,78 @@ fn test_diff_cor_nbody_nonregression() { .as_ref() .expect("8467 should converge"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 60672.2443617134, - semi_major_axis: 3.2064058028481552, - eccentricity_sin_lon: 0.05300520970081429, - eccentricity_cos_lon: -0.023197692700634636, - tan_half_incl_sin_node: 0.002896813138792025, - tan_half_incl_cos_node: 0.09181010554057693, - mean_longitude: 0.6256995904459722, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 60672.2443617134, + semi_major_axis: 3.2064058028481552, + eccentricity_sin_lon: 0.05300520970081429, + eccentricity_cos_lon: -0.023197692700634636, + tan_half_incl_sin_node: 0.002896813138792025, + tan_half_incl_cos_node: 0.09181010554057693, + mean_longitude: 0.6256995904459722, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 0.007572375820104381, + eccentricity_sin_lon: 0.0024777464468933156, + eccentricity_cos_lon: 0.0007445419051153811, + tan_half_incl_sin_node: 4.2789628256661375e-5, + tan_half_incl_cos_node: 5.7090614265788426e-5, + mean_longitude: 0.003334899745150928, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 5.73408755609015e-5, + 1.869803895909995e-5, + 5.597518961032454e-6, + -3.23358524526529e-7, + -4.293184367946207e-7, + 2.5017097632757542e-5, + ], + [ + 1.8698038959099974e-5, + 6.139227455092451e-6, + 1.8070844259588129e-6, + -1.0561768228833672e-7, + -1.409534091855587e-7, + 8.251003158870094e-6, + ], + [ + 5.597518961032457e-6, + 1.8070844259588129e-6, + 5.543426484728413e-7, + -3.149123145262512e-8, + -4.1500376105573355e-8, + 2.4017490220016504e-6, + ], + [ + -3.233585245265291e-7, + -1.0561768228833674e-7, + -3.149123145262512e-8, + 1.8309522863432738e-9, + 2.4373900210699776e-9, + -1.4146340017585484e-7, + ], + [ + -4.2931843679462117e-7, + -1.4095340918555872e-7, + -4.1500376105573355e-8, + 2.4373900210699772e-9, + 3.259338237245045e-9, + -1.894260958072586e-7, + ], + [ + 2.501709763275752e-5, + 8.251003158870094e-6, + 2.401749022001649e-6, + -1.414634001758548e-7, + -1.894260958072586e-7, + 1.1121556310207724e-5, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), tol), @@ -425,15 +678,78 @@ fn test_diff_cor_nbody_nonregression() { .as_ref() .expect("33803 should converge"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 60465.26777915681, - semi_major_axis: 2.190348311458185, - eccentricity_sin_lon: -0.13373910921857446, - eccentricity_cos_lon: 0.15339157238172804, - tan_half_incl_sin_node: 0.0029876412023201416, - tan_half_incl_cos_node: -0.05950692044872062, - mean_longitude: 4.224365041422834, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 60465.26777915681, + semi_major_axis: 2.190348311458185, + eccentricity_sin_lon: -0.13373910921857446, + eccentricity_cos_lon: 0.15339157238172804, + tan_half_incl_sin_node: 0.0029876412023201416, + tan_half_incl_cos_node: -0.05950692044872062, + mean_longitude: 4.224365041422834, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 2.1385808329040844e-5, + eccentricity_sin_lon: 1.3645976741407893e-5, + eccentricity_cos_lon: 5.318680335330679e-6, + tan_half_incl_sin_node: 3.4498285885877785e-7, + tan_half_incl_cos_node: 8.504424577287931e-7, + mean_longitude: 2.6647348205193914e-5, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 4.573527978864727e-10, + -2.441395550149477e-10, + 7.195967928874385e-11, + -1.8832832793515527e-12, + -6.328517179823655e-12, + 4.340505064567979e-10, + ], + [ + -2.4413955501494734e-10, + 1.862126812270452e-10, + -6.032972260112714e-11, + 8.172697046412589e-15, + -6.596088825542575e-13, + -3.58336068263501e-10, + ], + [ + 7.195967928874368e-11, + -6.032972260112711e-11, + 2.8288360509433262e-11, + 2.035203222614985e-14, + 1.4220205527899883e-13, + 1.2761184787385386e-10, + ], + [ + -1.8832832793515527e-12, + 8.17269704641279e-15, + 2.0352032226149824e-14, + 1.1901317290637546e-13, + 2.6436375372263e-13, + 3.754248073793805e-13, + ], + [ + -6.328517179823653e-12, + -6.596088825542574e-13, + 1.4220205527899883e-13, + 2.6436375372262996e-13, + 7.2325237390779e-13, + 2.605903338872906e-12, + ], + [ + 4.340505064567979e-10, + -3.5833606826350087e-10, + 1.2761184787385386e-10, + 3.754248073793805e-13, + 2.6059033388729057e-12, + 7.100811663688513e-10, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), tol), @@ -452,15 +768,78 @@ fn test_diff_cor_nbody_nonregression() { .as_ref() .expect("K09R05F should converge"); - let expected = OrbitalElements::Equinoctial(EquinoctialElements { - reference_epoch: 57049.2684537375, - semi_major_axis: 1.8021517900042052, - eccentricity_sin_lon: 0.2694922786015968, - eccentricity_cos_lon: 0.08955282358108035, - tan_half_incl_sin_node: 0.0008974287327937245, - tan_half_incl_cos_node: 0.10167548786557225, - mean_longitude: 1.6921653421358704, - }); + let expected = OrbitalElements::Equinoctial { + elements: EquinoctialElements { + reference_epoch: 57049.2684537375, + semi_major_axis: 1.8021517900042052, + eccentricity_sin_lon: 0.2694922786015968, + eccentricity_cos_lon: 0.08955282358108035, + tan_half_incl_sin_node: 0.0008974287327937245, + tan_half_incl_cos_node: 0.10167548786557225, + mean_longitude: 1.6921653421358704, + }, + uncertainty: Some(EquinoctialUncertainty { + semi_major_axis: 1.910876358918557e-6, + eccentricity_sin_lon: 2.7271919973585478e-6, + eccentricity_cos_lon: 1.2559941333300101e-5, + tan_half_incl_sin_node: 6.143234310625764e-7, + tan_half_incl_cos_node: 1.1476173256703189e-6, + mean_longitude: 2.1064465635865037e-5, + }), + covariance: Some(OrbitalCovariance { + matrix: [ + [ + 3.651448459073842e-12, + -4.87907485491453e-13, + 2.321298362132558e-11, + -3.7695250201166625e-13, + 8.511532638002078e-13, + -3.91138523482157e-11, + ], + [ + -4.879074854914533e-13, + 7.437576190456506e-12, + -1.1647669978804286e-11, + 9.359797430147383e-13, + -2.8577594338429333e-12, + 1.853502993770551e-11, + ], + [ + 2.3212983621325566e-11, + -1.164766997880434e-11, + 1.577521262959403e-10, + -3.47676746499932e-12, + 8.610023673871895e-12, + -2.644913915663376e-10, + ], + [ + -3.7695250201166625e-13, + 9.359797430147385e-13, + -3.4767674649993202e-12, + 3.7739327795249603e-13, + -5.048815271306508e-13, + 5.7505636344116006e-12, + ], + [ + 8.511532638002078e-13, + -2.857759433842935e-12, + 8.610023673871898e-12, + -5.048815271306507e-13, + 1.3170255261786945e-12, + -1.4110008489365913e-11, + ], + [ + -3.911385234821569e-11, + 1.8535029937705585e-11, + -2.6449139156633765e-10, + 5.750563634411601e-12, + -1.4110008489365913e-11, + 4.437117125245391e-10, + ], + ] + .into(), + }), + }; assert!( approx_equal(&expected, orbit.orbital_elements(), tol), diff --git a/tests/test_gauss_iod.rs b/tests/test_gauss_iod.rs index 4302b05..23c4f1c 100644 --- a/tests/test_gauss_iod.rs +++ b/tests/test_gauss_iod.rs @@ -22,39 +22,51 @@ struct ExpectedResult { fn expected_results() -> Vec { vec![ ExpectedResult { - orbit: OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 57049.2684537375, - semi_major_axis: 1.801740835743616, - eccentricity: 0.28356259478492557, - inclination: 0.2026828189979528, - ascending_node_longitude: 0.007951791820548622, - periapsis_argument: 1.2450647642587158, - mean_anomaly: 0.4408048786626789, - }), + orbit: OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 57049.2684537375, + semi_major_axis: 1.801740835743616, + eccentricity: 0.28356259478492557, + inclination: 0.2026828189979528, + ascending_node_longitude: 0.007951791820548622, + periapsis_argument: 1.2450647642587158, + mean_anomaly: 0.4408048786626789, + }, + uncertainty: None, + covariance: None, + }, rms: 66.97479288637471, }, ExpectedResult { - orbit: OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 60672.2443617134, - semi_major_axis: 3.2199380906809876, - eccentricity: 0.0624192099888107, - inclination: 0.1829771029880289, - ascending_node_longitude: 0.030775930195064964, - periapsis_argument: 1.9053705720223801, - mean_anomaly: 4.980622835177979, - }), + orbit: OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 60672.2443617134, + semi_major_axis: 3.2199380906809876, + eccentricity: 0.0624192099888107, + inclination: 0.1829771029880289, + ascending_node_longitude: 0.030775930195064964, + periapsis_argument: 1.9053705720223801, + mean_anomaly: 4.980622835177979, + }, + uncertainty: None, + covariance: None, + }, rms: 0.5739558189489471, }, ExpectedResult { - orbit: OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 60465.26777915681, - semi_major_axis: 2.1874983804796972, - eccentricity: 0.20256414489486008, - inclination: 0.11906245183260411, - ascending_node_longitude: 3.0918063960305293, - periapsis_argument: 2.4793248309745692, - mean_anomaly: 4.934465465531324, - }), + orbit: OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 60465.26777915681, + semi_major_axis: 2.1874983804796972, + eccentricity: 0.20256414489486008, + inclination: 0.11906245183260411, + ascending_node_longitude: 3.0918063960305293, + periapsis_argument: 2.4793248309745692, + mean_anomaly: 4.934465465531324, + }, + uncertainty: None, + covariance: None, + }, rms: 18.963755528781288, }, ] diff --git a/tests/test_iod_from_polars.rs b/tests/test_iod_from_polars.rs index 3b82ee1..353d01c 100644 --- a/tests/test_iod_from_polars.rs +++ b/tests/test_iod_from_polars.rs @@ -77,15 +77,19 @@ fn test_iod_from_polars() { let best_orbit = full_orbit.remove(&"14226".into()).unwrap().unwrap(); let orbit = best_orbit.orbital_elements(); - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 60894.372896385554, - semi_major_axis: 0.5415009930884174, - eccentricity: 0.9027228307040831, - inclination: 0.31200939353818746, - ascending_node_longitude: 5.550735343593096, - periapsis_argument: 3.1638244350882596, - mean_anomaly: 2.7888128618151495, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 60894.372896385554, + semi_major_axis: 0.5415009930884174, + eccentricity: 0.9027228307040831, + inclination: 0.31200939353818746, + ascending_node_longitude: 5.550735343593096, + periapsis_argument: 3.1638244350882596, + mean_anomaly: 2.7888128618151495, + }, + uncertainty: None, + covariance: None, + }; assert!(approx_equal(&expected_orbit, orbit, test_epsilon)); assert_relative_eq!( @@ -99,15 +103,19 @@ fn test_iod_from_polars() { let best_orbit = full_orbit.remove(&"29757".into()).unwrap().unwrap(); let orbit = best_orbit.orbital_elements(); - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 60835.25573266984, - semi_major_axis: 0.635955220245824, - eccentricity: 0.5904849550180079, - inclination: 0.2263529126300279, - ascending_node_longitude: 4.366539949885583, - periapsis_argument: 3.3107966723035602, - mean_anomaly: 3.0157533331616966, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 60835.25573266984, + semi_major_axis: 0.635955220245824, + eccentricity: 0.5904849550180079, + inclination: 0.2263529126300279, + ascending_node_longitude: 4.366539949885583, + periapsis_argument: 3.3107966723035602, + mean_anomaly: 3.0157533331616966, + }, + uncertainty: None, + covariance: None, + }; assert!(approx_equal(&expected_orbit, orbit, test_epsilon)); assert_relative_eq!( @@ -121,15 +129,19 @@ fn test_iod_from_polars() { let best_orbit = full_orbit.remove(&"95777".into()).unwrap().unwrap(); let orbit = best_orbit.orbital_elements(); - let expected_orbit = OrbitalElements::Keplerian(KeplerianElements { - reference_epoch: 60894.252965553926, - semi_major_axis: 1.24701989952089, - eccentricity: 0.2082069422196415, - inclination: 0.08116316040114972, - ascending_node_longitude: 2.49554922649176, - periapsis_argument: 2.5470318525197477, - mean_anomaly: 0.2983936748249412, - }); + let expected_orbit = OrbitalElements::Keplerian { + elements: KeplerianElements { + reference_epoch: 60894.252965553926, + semi_major_axis: 1.24701989952089, + eccentricity: 0.2082069422196415, + inclination: 0.08116316040114972, + ascending_node_longitude: 2.49554922649176, + periapsis_argument: 2.5470318525197477, + mean_anomaly: 0.2983936748249412, + }, + uncertainty: None, + covariance: None, + }; assert!(approx_equal(&expected_orbit, orbit, test_epsilon)); assert_relative_eq!(