add n-body propagation and associated tests

Cargo.toml

@@ -31,7 +31,7 @@ hifitime = { version = "4.2.0", default-features = false, features = [
     "ut1",
     "std",
 ] }
-nalgebra = { version = "0.33.2" }
+nalgebra = { version = "0.34.2" }
 ordered-float = { version = "5.0.0", default-features = false }
 smallvec = { version = "1.14.0", default-features = false }
 thiserror = { version = "2.0.12", default-features = false }
@@ -58,6 +58,7 @@ tokio = { version = "1.44.1", default-features = false, features = [
 ] }
 tokio-stream = { version = "0.1.17", default-features = false }
 rayon = { version = "1.12.0", default-features = false, optional = true }
+differential-equations = { version = "0.6.1", default-features = false }
 
 [dev-dependencies]
 approx = { version = "0.5.1", default-features = false }

examples/run_full_iod.rs

@@ -57,6 +57,8 @@ fn outfit_error_label(err: &OutfitError) -> &'static str {
         OutfitError::BizarreOrbit => "BizarreOrbit",
         OutfitError::DifferentialCorrectionDiverged => "DifferentialCorrectionDiverged",
         OutfitError::DifferentialCorrectionFailed(_) => "DifferentialCorrectionFailed",
+        OutfitError::EphemerisBodyNotSupported(_) => "EphemerisBodyNotSupported",
+        OutfitError::NBodyPropagationFailed(_) => "NBodyPropagationFailed",
     }
 }
 

examples/run_full_iod_parallel.rs

@@ -59,6 +59,8 @@ fn outfit_error_label(err: &OutfitError) -> &'static str {
         OutfitError::BizarreOrbit => "DiffCorBizarreOrbit",
         OutfitError::DifferentialCorrectionDiverged => "DiffCorDiverged",
         OutfitError::DifferentialCorrectionFailed(_) => "DiffCorFailed",
+        OutfitError::EphemerisBodyNotSupported(_) => "EphemerisBodyNotSupported",
+        OutfitError::NBodyPropagationFailed(_) => "NBodyPropagationFailed",
     }
 }
 

src/differential_orbit_correction/diff_cor.rs

@@ -56,6 +56,7 @@ use crate::{
         single_iteration::single_iteration,
     },
     orbit_type::equinoctial_element::EquinoctialLimits,
+    propagator::PropagatorKind,
     EquinoctialElements, JPLEphem, OutfitError,
 };
 
@@ -159,6 +160,13 @@ pub struct DifferentialCorrectionConfig {
     ///
     /// Default: `[true; 6]` (all elements free).
     pub free_elements: [bool; 6],
+
+    /// Propagator to use for computing predicted observations and partials.
+    ///
+    /// - [`PropagatorKind::TwoBody`] (default): analytic Keplerian propagation.
+    /// - [`PropagatorKind::NBody`]: numerical DOP853 N-body integration with
+    ///   user-specified perturbing bodies.
+    pub propagator: PropagatorKind,
 }
 
 impl Default for DifferentialCorrectionConfig {
@@ -175,6 +183,7 @@ impl Default for DifferentialCorrectionConfig {
             outlier_rejection_config: OutlierRejectionConfig::default(),
             orbital_limits: EquinoctialLimits::default(),
             free_elements: [true; 6],
+            propagator: PropagatorKind::TwoBody,
         }
     }
 }
@@ -300,19 +309,9 @@ pub fn run_differential_correction(
                 cache,
                 jpl,
                 true,
+                &config.propagator,
             )?;
 
-            println!(
-                "Corrected elements after iteration {total_newton_iterations}: {:?}",
-                iter_result.corrected_elements
-            );
-            println!(
-                "IterResult correction norm: {}, normalised RMS: {}, num active measurements: {}\n\n",
-                iter_result.correction_norm,
-                iter_result.normalised_rms,
-                iter_result.num_measurements
-            );
-
             // ── Check inversion ──────────────────────────────────────────────
             if !iter_result.uncertainty.inversion_succeeded {
                 return Err(OutfitError::DifferentialCorrectionFailed(
@@ -395,12 +394,6 @@ pub fn run_differential_correction(
         if num_selection_changes == 0 {
             break;
         }
-
-        println!("last normalised RMS after outer pass {outer_pass}: {last_normalised_rms}");
-        println!("number of selection changes in outlier rejection: {num_selection_changes}");
-        println!("number of active measurements after outlier rejection: {last_num_measurements}");
-        println!("stagnation count after outer pass {outer_pass}: {stagnation_count}");
-        println!("\n ====== End of outer pass {outer_pass} ====== \n");
     }
 
     // ── Final covariance rescaling ────────────────────────────────────────────

src/differential_orbit_correction/obs_dataset_api.rs

@@ -222,8 +222,6 @@ fn run_differential_correction_for_trajectory(
         }
     };
 
-    println!("\n\nTrajectory {traj_id}: initial equinoctial elements: {initial_equinoctial:?}\n\n");
-
     // ── Build per-observation fit data from the error-model uncertainties ─────
     let obs_fit_data: Vec<ObsFitData> = observations
         .iter()
@@ -240,8 +238,6 @@ fn run_differential_correction_for_trajectory(
         diff_cor_config,
     )?;
 
-    println!("\ndc_output for trajectory {traj_id}: {dc_output:?}\n");
-
     // ── Package the result ────────────────────────────────────────────────────
     let orbital_elements = OrbitalElements::Equinoctial(dc_output.elements);
     Ok(FitOrbitResult::DifferentialCorrection((