From 7dff52734453a5f567d51911b6d3fe5b15badb07 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Thu, 15 Aug 2024 07:50:15 -0600 Subject: [PATCH] Still trying to pinpoint the modeling differences --- examples/04_lro_od/main.rs | 45 +++++++++++++++++++++++------- examples/04_lro_od/plot_od_rslt.py | 33 ++++++++++++---------- src/od/process/mod.rs | 6 ++-- 3 files changed, 57 insertions(+), 27 deletions(-) diff --git a/examples/04_lro_od/main.rs b/examples/04_lro_od/main.rs index d2c21da9..c289f1b4 100644 --- a/examples/04_lro_od/main.rs +++ b/examples/04_lro_od/main.rs @@ -26,7 +26,7 @@ use nyx::{ GroundStation, }, propagators::Propagator, - Orbit, Spacecraft, + Orbit, Spacecraft, State, }; use std::{collections::BTreeMap, error::Error, path::PathBuf, str::FromStr, sync::Arc}; @@ -61,6 +61,9 @@ fn main() -> Result<(), Box> { // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case. // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO. let lro_frame = Frame::from_ephem_j2000(-85); + + // We also use a different gravitational parameter for the Moon to have a better estimation + // To build the trajectory we need to provide a spacecraft template. let sc_template = Spacecraft::builder() .dry_mass_kg(1018.0) // Launch masses @@ -74,7 +77,7 @@ fn main() -> Result<(), Box> { .build(); // Now we can build the trajectory from the BSP file. // We'll arbitrarily set the tracking arc to 48 hours with a one minute time step. - let trajectory = Traj::from_bsp( + let traj_as_flown = Traj::from_bsp( lro_frame, MOON_J2000, almanac.clone(), @@ -86,7 +89,7 @@ fn main() -> Result<(), Box> { Some("LRO".to_string()), )?; - println!("{trajectory}"); + println!("{traj_as_flown}"); // Load the Deep Space Network ground stations. // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML. @@ -117,7 +120,7 @@ fn main() -> Result<(), Box> { // Build the tracking arc simulation to generate a "standard measurement". let mut trk = TrackingArcSim::::with_seed( devices, - trajectory.clone(), + traj_as_flown.clone(), configs, 12345, )?; @@ -164,8 +167,7 @@ fn main() -> Result<(), Box> { // We define the solar radiation pressure, using the default solar flux and accounting only // for the eclipsing caused by the Earth and Moon. // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity. - let srp_dyn = - SolarPressure::default_no_estimation(vec![EARTH_J2000, MOON_J2000], almanac.clone())?; + let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?; // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models. @@ -179,7 +181,7 @@ fn main() -> Result<(), Box> { // For an OD arc, we need to start with an initial estimate and covariance. // The ephem published by NASA does not include the covariance. Instead, we'll make one up! - let sc_seed = trajectory.first().with_stm(); + let sc_seed = traj_as_flown.first().with_stm(); let sc = SpacecraftUncertainty::builder() .nominal(sc_seed) @@ -200,11 +202,16 @@ fn main() -> Result<(), Box> { // Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame. let kf = KF::new( initial_estimate, - SNC3::from_diagonal(2 * Unit::Minute, &[1e-16, 1e-16, 1e-16]), + SNC3::from_diagonal(2 * Unit::Minute, &[1e-15, 1e-15, 1e-15]), ); // We'll set up the OD process to reject measurements whose residuals are mover than 4 sigmas away from what we expect. - let mut odp = ODProcess::ckf(setup.with(sc_seed, almanac.clone()), kf, None, almanac); + let mut odp = ODProcess::ckf( + setup.with(sc_seed, almanac.clone()), + kf, + None, + almanac.clone(), + ); odp.process_arc::(&arc)?; // Let's run a smoother just to see that the filter won't run it if the RSS error is small. @@ -218,10 +225,28 @@ fn main() -> Result<(), Box> { let od_trajectory = odp.to_traj()?; // Build the RIC difference. od_trajectory.ric_diff_to_parquet( - &trajectory, + &traj_as_flown, "./04_lro_od_truth_error.parquet", ExportCfg::default(), )?; + // For reference, let's build the trajectory with Nyx's models from that LRO state. + let (_, traj_as_sim) = setup + .with(*traj_as_flown.first(), almanac.clone()) + .until_epoch_with_traj(traj_as_flown.last().epoch())?; + + // Build the RIC differences. + od_trajectory.ric_diff_to_parquet( + &traj_as_sim, + "./04_lro_od_sim_error.parquet", + ExportCfg::default(), + )?; + + traj_as_sim.ric_diff_to_parquet( + &traj_as_flown, + "./04_lro_sim_truth_error.parquet", + ExportCfg::default(), + )?; + Ok(()) } diff --git a/examples/04_lro_od/plot_od_rslt.py b/examples/04_lro_od/plot_od_rslt.py index fdbff8b6..15c9b1bb 100644 --- a/examples/04_lro_od/plot_od_rslt.py +++ b/examples/04_lro_od/plot_od_rslt.py @@ -33,18 +33,21 @@ px.line(df, x="Epoch (UTC)", y=["cr", "Cr + Sigma", "Cr - Sigma"]).show() # Load the RIC diff. - df_ric = pl.read_parquet("./04_lro_od_truth_error.parquet") - df_ric = df_ric.with_columns( - pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f") - ).sort("Epoch (UTC)", descending=False) - # Plot the RIC difference - px.line( - df_ric, - x="Epoch (UTC)", - y=["Delta X (RIC) (km)", "Delta Y (RIC) (km)", "Delta Z (RIC) (km)"], - ).show() - px.line( - df_ric, - x="Epoch (UTC)", - y=["Delta Vx (RIC) (km/s)", "Delta Vy (RIC) (km/s)", "Delta Vz (RIC) (km/s)"], - ).show() + for fname, errname in [("04_lro_od_truth_error", "OD vs Flown"), ("04_lro_od_sim_error", "OD vs Sim"), ("04_lro_sim_truth_error", "Sim vs Flown")]: + df_ric = pl.read_parquet(f"./{fname}.parquet") + df_ric = df_ric.with_columns( + pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f") + ).sort("Epoch (UTC)", descending=False) + # Plot the RIC difference + px.line( + df_ric, + x="Epoch (UTC)", + y=["Delta X (RIC) (km)", "Delta Y (RIC) (km)", "Delta Z (RIC) (km)"], + title=f"Position error with {errname} ({fname})" + ).show() + px.line( + df_ric, + x="Epoch (UTC)", + y=["Delta Vx (RIC) (km/s)", "Delta Vy (RIC) (km/s)", "Delta Vz (RIC) (km/s)"], + title=f"Velocity error with {errname} ({fname})" + ).show() diff --git a/src/od/process/mod.rs b/src/od/process/mod.rs index 75558684..511ff896 100644 --- a/src/od/process/mod.rs +++ b/src/od/process/mod.rs @@ -507,6 +507,7 @@ where let mut traj: Traj = Traj::new(); let mut msr_accepted_cnt: usize = 0; + let tick = Epoch::now().unwrap(); for (msr_cnt, (device_name, msr)) in measurements.iter().enumerate() { let next_msr_epoch = msr.epoch(); @@ -639,7 +640,7 @@ where if !reported[msr_prct] { let num_rejected = msr_cnt - msr_accepted_cnt.saturating_sub(1); let msg = format!( - "{:>3}% done ({msr_accepted_cnt:.0} measurements accepted, {:.0} rejected)", + "{:>3}% done - {msr_accepted_cnt:.0} measurements accepted, {:.0} rejected", 10 * msr_prct, num_rejected ); if msr_accepted_cnt < num_rejected { @@ -671,8 +672,9 @@ where // Always report the 100% mark if !reported[10] { + let tock_time = Epoch::now().unwrap() - tick; info!( - "100% done ({msr_accepted_cnt:.0} measurements accepted, {:.0} rejected)", + "100% done - {msr_accepted_cnt:.0} measurements accepted, {:.0} rejected (done in {tock_time})", num_msrs - msr_accepted_cnt ); }