Differential correction stats errors on estimated orbits

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

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<DifferentialCorrectionOutput> 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
 // ─────────────────────────────────────────────────────────────────────────────

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,15 +210,27 @@ 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
             // observations, reusing the cache that was built for the full
             // 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,
     )))
 }

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<crate::EquinoctialElements, OutfitError> {
-        self.to_equinoctial()
+    fn to_equinoctial_for_ephemeris(&self) -> Result<EquinoctialElements, OutfitError> {
+        self.to_equinoctial()?
+            .as_equinoctial()
+            .ok_or(OutfitError::InvalidConversion(
+                "Conversion to equinoctial elements failed".to_string(),
+            ))
     }
 }

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

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"));
     }
 }

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.