Struct ODSolution

Source
pub struct ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,
{ pub estimates: Vec<EstType>, pub residuals: Vec<Option<Residual<MsrSize>>>, pub gains: Vec<Option<OMatrix<f64, <StateType as State>::Size, MsrSize>>>, pub filter_smoother_ratios: Vec<Option<OVector<f64, <StateType as State>::Size>>>, pub devices: BTreeMap<String, Trk>, pub measurement_types: IndexSet<MeasurementType>, }
Expand description

The ODSolution structure is designed to manage and analyze the results of an OD process, including smoothing. It provides various functionalities such as splitting solutions by tracker or measurement type, joining solutions, and performing statistical analyses.

Note: Many methods in this structure assume that the solution has been split into subsets using the split() method. Calling these methods without first splitting will make analysis of operations results less obvious.

§Fields

  • estimates: A vector of state estimates generated during the OD process.
  • residuals: A vector of residuals corresponding to the state estimates.
  • gains: Filter gains used for measurement updates. These are set to None after running the smoother.
  • filter_smoother_ratios: Filter-smoother consistency ratios. These are set to None before running the smoother.
  • devices: A map of tracking devices used in the OD process.
  • measurement_types: A set of unique measurement types used in the OD process.

Implementation detail: these are not stored in vectors to allow for multiple estimates at the same time, e.g. when there are simultaneous measurements of angles and the filter processes each as a scalar.

Fields§

§estimates: Vec<EstType>

Vector of estimates available after a pass

§residuals: Vec<Option<Residual<MsrSize>>>

Vector of residuals available after a pass

§gains: Vec<Option<OMatrix<f64, <StateType as State>::Size, MsrSize>>>

Vector of filter gains used for each measurement update, all None after running the smoother.

§filter_smoother_ratios: Vec<Option<OVector<f64, <StateType as State>::Size>>>

Filter-smoother consistency ratios, all None before running the smoother.

§devices: BTreeMap<String, Trk>

Tracking devices

§measurement_types: IndexSet<MeasurementType>

Implementations§

Source§

impl<MsrSize: DimName, Trk: TrackerSensitivity<Spacecraft, Spacecraft>> ODSolution<Spacecraft, KfEstimate<Spacecraft>, MsrSize, Trk>
where DefaultAllocator: Allocator<MsrSize> + Allocator<MsrSize, <Spacecraft as State>::Size> + Allocator<Const<1>, MsrSize> + Allocator<<Spacecraft as State>::Size> + Allocator<<Spacecraft as State>::Size, <Spacecraft as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<Spacecraft as State>::Size, MsrSize> + Allocator<<Spacecraft as State>::VecLength>,

Source

pub fn to_parquet<P: AsRef<Path>>( &self, path: P, cfg: ExportCfg, ) -> Result<PathBuf, ODError>

Store the estimates and residuals in a parquet file

Examples found in repository?
examples/02_jwst_covar_monte_carlo/main.rs (line 128)
26fn main() -> Result<(), Box<dyn Error>> {
27    pel::init();
28    // Dynamics models require planetary constants and ephemerides to be defined.
29    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
31
32    // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
33    // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
34    let mut latest_jwst_ephem = MetaFile {
35        uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
36        crc32: None,
37    };
38    latest_jwst_ephem.process(true)?;
39
40    // Load this ephem in the general Almanac we're using for this analysis.
41    let almanac = Arc::new(
42        MetaAlmanac::latest()
43            .map_err(Box::new)?
44            .load_from_metafile(latest_jwst_ephem, true)?,
45    );
46
47    // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
48    // in the BSP. We need this ID in order to query the ephemeris.
49    const JWST_NAIF_ID: i32 = -170;
50    // Let's build a frame in the J2000 orientation centered on the JWST.
51    const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
52
53    // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
54    let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
55    println!("JWST defined from {earliest_epoch} to {latest_epoch}");
56    // Fetch the state, printing it in the Earth J2000 frame.
57    let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
58    println!("{jwst_orbit:x}");
59
60    // Build the spacecraft
61    // SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
62    // SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
63    let jwst = Spacecraft::builder()
64        .orbit(jwst_orbit)
65        .srp(SRPData {
66            area_m2: 21.197 * 14.162,
67            coeff_reflectivity: 1.56,
68        })
69        .mass(Mass::from_dry_mass(6200.0))
70        .build();
71
72    // Build up the spacecraft uncertainty builder.
73    // We can use the spacecraft uncertainty structure to build this up.
74    // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
75    // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
76    // Nyx can also estimate the deviation of the spacecraft parameters.
77    let jwst_uncertainty = SpacecraftUncertainty::builder()
78        .nominal(jwst)
79        .frame(LocalFrame::RIC)
80        .x_km(0.5)
81        .y_km(0.3)
82        .z_km(1.5)
83        .vx_km_s(1e-4)
84        .vy_km_s(0.6e-3)
85        .vz_km_s(3e-3)
86        .build();
87
88    println!("{jwst_uncertainty}");
89
90    // Build the Kalman filter estimate.
91    // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
92    // but this approach requires quite a bit more boilerplate code.
93    let jwst_estimate = jwst_uncertainty.to_estimate()?;
94
95    // Set up the spacecraft dynamics.
96    // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
97    // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
98
99    let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
100    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102    // Finalize setting up the dynamics.
103    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
104
105    // Build the propagator set up to use for the whole analysis.
106    let setup = Propagator::default(dynamics);
107
108    // All of the analysis will use this duration.
109    let prediction_duration = 6.5 * Unit::Day;
110
111    // === Covariance mapping ===
112    // For the covariance mapping / prediction, we'll use the common orbit determination approach.
113    // This is done by setting up a spacecraft Kalman filter OD process, and predicting for the analysis duration.
114
115    // Build the propagation instance for the OD process.
116    let odp = SpacecraftKalmanOD::new(
117        setup.clone(),
118        KalmanVariant::DeviationTracking,
119        None,
120        BTreeMap::new(),
121        almanac.clone(),
122    );
123
124    // The prediction step is 1 minute by default, configured in the OD process, i.e. how often we want to know the covariance.
125    assert_eq!(odp.max_step, 1_i64.minutes());
126    // Finally, predict, and export the trajectory with covariance to a parquet file.
127    let od_sol = odp.predict_for(jwst_estimate, prediction_duration)?;
128    od_sol.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
129
130    // === Monte Carlo framework ===
131    // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
132
133    let my_mc = MonteCarlo::new(
134        jwst, // Nominal state
135        jwst_estimate.to_random_variable()?,
136        "02_jwst".to_string(), // Scenario name
137        None, // No specific seed specified, so one will be drawn from the computer's entropy.
138    );
139
140    let num_runs = 5_000;
141    let rslts = my_mc.run_until_epoch(
142        setup,
143        almanac.clone(),
144        jwst.epoch() + prediction_duration,
145        num_runs,
146    );
147
148    assert_eq!(rslts.runs.len(), num_runs);
149    // Finally, export these results, computing the eclipse percentage for all of these results.
150
151    // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
152    let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
153    let umbra_event = eclipse_loc.to_umbra_event();
154    let penumbra_event = eclipse_loc.to_penumbra_event();
155
156    rslts.to_parquet(
157        "02_jwst_monte_carlo.parquet",
158        Some(vec![&umbra_event, &penumbra_event]),
159        ExportCfg::default(),
160        almanac,
161    )?;
162
163    Ok(())
164}
More examples
Hide additional examples
examples/04_lro_od/main.rs (line 290)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
261    let odp = SpacecraftKalmanOD::new(
262        setup,
263        KalmanVariant::ReferenceUpdate,
264        Some(ResidRejectCrit::default()),
265        devices,
266        almanac.clone(),
267    )
268    .with_process_noise(process_noise);
269
270    let od_sol = odp.process_arc(initial_estimate, &arc)?;
271
272    let ric_err = traj_as_flown
273        .at(od_sol.estimates.last().unwrap().epoch())?
274        .orbit
275        .ric_difference(&od_sol.estimates.last().unwrap().orbital_state())?;
276    println!("== RIC at end ==");
277    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
278    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
279
280    println!(
281        "Num residuals rejected: #{}",
282        od_sol.rejected_residuals().len()
283    );
284    println!(
285        "Percentage within +/-3: {}",
286        od_sol.residual_ratio_within_threshold(3.0).unwrap()
287    );
288    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289
290    od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
291
292    // In our case, we have the truth trajectory from NASA.
293    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
294    // Export the OD trajectory first.
295    let od_trajectory = od_sol.to_traj()?;
296    // Build the RIC difference.
297    od_trajectory.ric_diff_to_parquet(
298        &traj_as_flown,
299        "./04_lro_od_truth_error.parquet",
300        ExportCfg::default(),
301    )?;
302
303    Ok(())
304}
Source§

impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source

pub fn unique(&self) -> IndexSet<(String, MeasurementType)>

Returns a set of tuples of tracker and measurement types in this OD solution, e.g. {(Canberra, Range), (Canberra, Doppler)}.

Source

pub fn drop_time_updates(self) -> Self

Returns this OD solution without any time update

Source

pub fn filter_by_msr_type(self, msr_type: MeasurementType) -> Self

Returns this OD solution with only data from the desired measurement type, dropping all time updates.

Source

pub fn filter_by_tracker(self, tracker: String) -> Self

Returns this OD solution with only data from the desired tracker, dropping all time updates.

Source

pub fn exclude_tracker(self, excluded_tracker: String) -> Self

Returns this OD solution with all data except from the desired tracker, including all time updates

Source

pub fn split(self) -> Vec<Self>

Split this OD solution per tracker and per measurement type, dropping all time updates.

Source

pub fn merge(self, other: Self) -> Self

Merge this OD solution with another one, returning a new OD solution.

Source

pub fn at(&self, epoch: Epoch) -> Option<ODRecord<StateType, EstType, MsrSize>>

Source§

impl<MsrSize, Trk> ODSolution<Spacecraft, KfEstimate<Spacecraft>, MsrSize, Trk>

Source

pub fn from_parquet<P: AsRef<Path>>( path: P, devices: BTreeMap<String, Trk>, ) -> Result<Self, InputOutputError>

Loads an OD solution from a Parquet file created by ODSolution::to_parquet.

The devices map must be provided by the caller as it contains potentially complex state (like Almanac references) that isn’t serialized in the Parquet file.

Note: This function currently assumes the StateType is Spacecraft and the estimate type is KfEstimate<Spacecraft>.

Source§

impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source

pub fn smooth(self, almanac: Arc<Almanac>) -> Result<Self, ODError>

Smoothes this OD solution, returning a new OD solution and the filter-smoother consistency ratios, with updated postfit residuals, and where the ratio now represents the filter-smoother consistency ratio.

Notes:

  1. Gains will be scrubbed because the smoother process does not recompute the gain.
  2. Prefit residuals, ratios, and measurement covariances are not updated, as these depend on the filtering process.
  3. Note: this function consumes the current OD solution to prevent reusing the wrong one.

To assess whether the smoothing process improved the solution, compare the RMS of the postfit residuals from the filter and the smoother process.

§Filter-Smoother consistency ratio

The filter-smoother consistency ratio is used to evaluate the consistency between the state estimates produced by a filter (e.g., Kalman filter) and a smoother. This ratio is called “filter smoother consistency test” in the ODTK MathSpec.

It is computed as follows:

§1. Define the State Estimates

Filter state estimate: $ \hat{X}_{f,k} $ This is the state estimate at time step $ k $ from the filter.

Smoother state estimate: $ \hat{X}_{s,k} $ This is the state estimate at time step $ k $ from the smoother.

§2. Define the Covariances

Filter covariance: $ P_{f,k} $ This is the covariance of the state estimate at time step $ k $ from the filter.

Smoother covariance: $ P_{s,k} $ This is the covariance of the state estimate at time step $ k $ from the smoother.

§3. Compute the Differences

State difference: $ \Delta X_k = \hat{X}{s,k} - \hat{X}{f,k} $

Covariance difference: $ \Delta P_k = P_{s,k} - P_{f,k} $

§4. Calculate the Consistency Ratio

For each element $ i $ of the state vector, compute the ratio:

$$ R_{i,k} = \frac{\Delta X_{i,k}}{\sqrt{\Delta P_{i,k}}} $$

§5. Evaluate Consistency
  • If $ |R_{i,k}| \leq 3 $ for all $ i $ and $ k $, the filter-smoother consistency test is satisfied, indicating good consistency.
  • If $ |R_{i,k}| > 3 $ for any $ i $ or $ k $, the test fails, suggesting potential modeling inconsistencies or issues with the estimation process.
Source§

impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source

pub fn rms_prefit_residuals(&self) -> f64

Returns the root mean square of the prefit residuals

Source

pub fn rms_postfit_residuals(&self) -> f64

Returns the root mean square of the postfit residuals

Source

pub fn rms_residual_ratios(&self) -> f64

Returns the root mean square of the prefit residual ratios

Source

pub fn residual_ratio_within_threshold( &self, threshold: f64, ) -> Result<f64, ODError>

Computes the fraction of residual ratios that lie within ±threshold.

Examples found in repository?
examples/04_lro_od/main.rs (line 286)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
261    let odp = SpacecraftKalmanOD::new(
262        setup,
263        KalmanVariant::ReferenceUpdate,
264        Some(ResidRejectCrit::default()),
265        devices,
266        almanac.clone(),
267    )
268    .with_process_noise(process_noise);
269
270    let od_sol = odp.process_arc(initial_estimate, &arc)?;
271
272    let ric_err = traj_as_flown
273        .at(od_sol.estimates.last().unwrap().epoch())?
274        .orbit
275        .ric_difference(&od_sol.estimates.last().unwrap().orbital_state())?;
276    println!("== RIC at end ==");
277    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
278    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
279
280    println!(
281        "Num residuals rejected: #{}",
282        od_sol.rejected_residuals().len()
283    );
284    println!(
285        "Percentage within +/-3: {}",
286        od_sol.residual_ratio_within_threshold(3.0).unwrap()
287    );
288    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289
290    od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
291
292    // In our case, we have the truth trajectory from NASA.
293    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
294    // Export the OD trajectory first.
295    let od_trajectory = od_sol.to_traj()?;
296    // Build the RIC difference.
297    od_trajectory.ric_diff_to_parquet(
298        &traj_as_flown,
299        "./04_lro_od_truth_error.parquet",
300        ExportCfg::default(),
301    )?;
302
303    Ok(())
304}
Source

pub fn ks_test_normality(&self) -> Result<f64, ODError>

Computes the Kolmogorov–Smirnov statistic for the aggregated residual ratios, by tracker and measurement type.

Returns Ok(ks_statistic) if residuals are available.

Source

pub fn is_normal(&self, alpha: Option<f64>) -> Result<bool, ODError>

Checks whether the residual ratios pass a normality test at a given significance level alpha, default to 0.05.

This uses a simplified KS-test threshold: D_alpha = c(α) / √n. For example, for α = 0.05, c(α) is approximately 1.36. α=0.05 means a 5% probability of Type I error (incorrectly rejecting the null hypothesis when it is true). This threshold is standard in many statistical tests to balance sensitivity and false positives

The computation of the c(alpha) is from https://en.wikipedia.org/w/index.php?title=Kolmogorov%E2%80%93Smirnov_test&oldid=1280260701#Two-sample_Kolmogorov%E2%80%93Smirnov_test

Returns Ok(true) if the residuals are consistent with a normal distribution, Ok(false) if not, or None if no residuals are available.

Examples found in repository?
examples/04_lro_od/main.rs (line 288)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
261    let odp = SpacecraftKalmanOD::new(
262        setup,
263        KalmanVariant::ReferenceUpdate,
264        Some(ResidRejectCrit::default()),
265        devices,
266        almanac.clone(),
267    )
268    .with_process_noise(process_noise);
269
270    let od_sol = odp.process_arc(initial_estimate, &arc)?;
271
272    let ric_err = traj_as_flown
273        .at(od_sol.estimates.last().unwrap().epoch())?
274        .orbit
275        .ric_difference(&od_sol.estimates.last().unwrap().orbital_state())?;
276    println!("== RIC at end ==");
277    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
278    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
279
280    println!(
281        "Num residuals rejected: #{}",
282        od_sol.rejected_residuals().len()
283    );
284    println!(
285        "Percentage within +/-3: {}",
286        od_sol.residual_ratio_within_threshold(3.0).unwrap()
287    );
288    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289
290    od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
291
292    // In our case, we have the truth trajectory from NASA.
293    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
294    // Export the OD trajectory first.
295    let od_trajectory = od_sol.to_traj()?;
296    // Build the RIC difference.
297    od_trajectory.ric_diff_to_parquet(
298        &traj_as_flown,
299        "./04_lro_od_truth_error.parquet",
300        ExportCfg::default(),
301    )?;
302
303    Ok(())
304}
Source§

impl<StateType, EstType, MsrSize, Trk> ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source

pub fn new( devices: BTreeMap<String, Trk>, measurement_types: IndexSet<MeasurementType>, ) -> Self

Source

pub fn results( &self, ) -> Zip<Iter<'_, EstType>, Iter<'_, Option<Residual<MsrSize>>>>

Returns a zipper iterator on the estimates and the associated residuals.

Source

pub fn is_filter_run(&self) -> bool

Returns True if this is the result of a filter run

Source

pub fn is_smoother_run(&self) -> bool

Returns True if this is the result of a smoother run

Source

pub fn to_traj(&self) -> Result<Traj<StateType>, NyxError>

Builds the navigation trajectory for the estimated state only

Examples found in repository?
examples/04_lro_od/main.rs (line 295)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
261    let odp = SpacecraftKalmanOD::new(
262        setup,
263        KalmanVariant::ReferenceUpdate,
264        Some(ResidRejectCrit::default()),
265        devices,
266        almanac.clone(),
267    )
268    .with_process_noise(process_noise);
269
270    let od_sol = odp.process_arc(initial_estimate, &arc)?;
271
272    let ric_err = traj_as_flown
273        .at(od_sol.estimates.last().unwrap().epoch())?
274        .orbit
275        .ric_difference(&od_sol.estimates.last().unwrap().orbital_state())?;
276    println!("== RIC at end ==");
277    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
278    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
279
280    println!(
281        "Num residuals rejected: #{}",
282        od_sol.rejected_residuals().len()
283    );
284    println!(
285        "Percentage within +/-3: {}",
286        od_sol.residual_ratio_within_threshold(3.0).unwrap()
287    );
288    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289
290    od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
291
292    // In our case, we have the truth trajectory from NASA.
293    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
294    // Export the OD trajectory first.
295    let od_trajectory = od_sol.to_traj()?;
296    // Build the RIC difference.
297    od_trajectory.ric_diff_to_parquet(
298        &traj_as_flown,
299        "./04_lro_od_truth_error.parquet",
300        ExportCfg::default(),
301    )?;
302
303    Ok(())
304}
Source

pub fn accepted_residuals(&self) -> Vec<Residual<MsrSize>>

Returns the accepted residuals.

Source

pub fn rejected_residuals(&self) -> Vec<Residual<MsrSize>>

Returns the rejected residuals.

Examples found in repository?
examples/04_lro_od/main.rs (line 282)
33fn main() -> Result<(), Box<dyn Error>> {
34    pel::init();
35
36    // ====================== //
37    // === ALMANAC SET UP === //
38    // ====================== //
39
40    // Dynamics models require planetary constants and ephemerides to be defined.
41    // Let's start by grabbing those by using ANISE's MetaAlmanac.
42
43    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
44        .iter()
45        .collect();
46
47    let meta = data_folder.join("lro-dynamics.dhall");
48
49    // Load this ephem in the general Almanac we're using for this analysis.
50    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
51        .map_err(Box::new)?
52        .process(true)
53        .map_err(Box::new)?;
54
55    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
56    moon_pc.mu_km3_s2 = 4902.74987;
57    almanac.planetary_data.set_by_id(MOON, moon_pc)?;
58
59    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
60    earth_pc.mu_km3_s2 = 398600.436;
61    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
62
63    // Save this new kernel for reuse.
64    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
65    almanac
66        .planetary_data
67        .save_as(&data_folder.join("lro-specific.pca"), true)?;
68
69    // Lock the almanac (an Arc is a read only structure).
70    let almanac = Arc::new(almanac);
71
72    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
73    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
74    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
75    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
76    let lro_frame = Frame::from_ephem_j2000(-85);
77
78    // To build the trajectory we need to provide a spacecraft template.
79    let sc_template = Spacecraft::builder()
80        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
81        .srp(SRPData {
82            // SRP configuration is arbitrary, but we will be estimating it anyway.
83            area_m2: 3.9 * 2.7,
84            coeff_reflectivity: 0.96,
85        })
86        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
87        .build();
88    // Now we can build the trajectory from the BSP file.
89    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
90    let traj_as_flown = Traj::from_bsp(
91        lro_frame,
92        MOON_J2000,
93        almanac.clone(),
94        sc_template,
95        5.seconds(),
96        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
97        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
98        Aberration::LT,
99        Some("LRO".to_string()),
100    )?;
101
102    println!("{traj_as_flown}");
103
104    // ====================== //
105    // === MODEL MATCHING === //
106    // ====================== //
107
108    // Set up the spacecraft dynamics.
109
110    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
111    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
112    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
113
114    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
115    // We're using the GRAIL JGGRX model.
116    let mut jggrx_meta = MetaFile {
117        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
118        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
119    };
120    // And let's download it if we don't have it yet.
121    jggrx_meta.process(true)?;
122
123    // Build the spherical harmonics.
124    // The harmonics must be computed in the body fixed frame.
125    // We're using the long term prediction of the Moon principal axes frame.
126    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
127    let sph_harmonics = Harmonics::from_stor(
128        almanac.frame_from_uid(moon_pa_frame)?,
129        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
130    );
131
132    // Include the spherical harmonics into the orbital dynamics.
133    orbital_dyn.accel_models.push(sph_harmonics);
134
135    // We define the solar radiation pressure, using the default solar flux and accounting only
136    // for the eclipsing caused by the Earth and Moon.
137    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
138    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
139
140    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
141    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
142    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
143
144    println!("{dynamics}");
145
146    // Now we can build the propagator.
147    let setup = Propagator::default_dp78(dynamics.clone());
148
149    // For reference, let's build the trajectory with Nyx's models from that LRO state.
150    let (sim_final, traj_as_sim) = setup
151        .with(*traj_as_flown.first(), almanac.clone())
152        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
153
154    println!("SIM INIT:  {:x}", traj_as_flown.first());
155    println!("SIM FINAL: {sim_final:x}");
156    // Compute RIC difference between SIM and LRO ephem
157    let sim_lro_delta = sim_final
158        .orbit
159        .ric_difference(&traj_as_flown.last().orbit)?;
160    println!("{traj_as_sim}");
161    println!(
162        "SIM v LRO - RIC Position (m): {:.3}",
163        sim_lro_delta.radius_km * 1e3
164    );
165    println!(
166        "SIM v LRO - RIC Velocity (m/s): {:.3}",
167        sim_lro_delta.velocity_km_s * 1e3
168    );
169
170    traj_as_sim.ric_diff_to_parquet(
171        &traj_as_flown,
172        "./04_lro_sim_truth_error.parquet",
173        ExportCfg::default(),
174    )?;
175
176    // ==================== //
177    // === OD SIMULATOR === //
178    // ==================== //
179
180    // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
181    // and the truth LRO state.
182
183    // Therefore, we will actually run an estimation from a dispersed LRO state.
184    // The sc_seed is the true LRO state from the BSP.
185    let sc_seed = *traj_as_flown.first();
186
187    // Load the Deep Space Network ground stations.
188    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
189    let ground_station_file: PathBuf = [
190        env!("CARGO_MANIFEST_DIR"),
191        "examples",
192        "04_lro_od",
193        "dsn-network.yaml",
194    ]
195    .iter()
196    .collect();
197
198    let devices = GroundStation::load_named(ground_station_file)?;
199
200    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
201    // Nyx can build a tracking schedule for you based on the first station with access.
202    let trkconfg_yaml: PathBuf = [
203        env!("CARGO_MANIFEST_DIR"),
204        "examples",
205        "04_lro_od",
206        "tracking-cfg.yaml",
207    ]
208    .iter()
209    .collect();
210
211    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
212
213    // Build the tracking arc simulation to generate a "standard measurement".
214    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
215        devices.clone(),
216        traj_as_flown.clone(),
217        configs,
218    )?;
219
220    trk.build_schedule(almanac.clone())?;
221    let arc = trk.generate_measurements(almanac.clone())?;
222    // Save the simulated tracking data
223    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
224
225    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
226    println!("{arc}");
227
228    // Now that we have simulated measurements, we'll run the orbit determination.
229
230    // ===================== //
231    // === OD ESTIMATION === //
232    // ===================== //
233
234    let sc = SpacecraftUncertainty::builder()
235        .nominal(sc_seed)
236        .frame(LocalFrame::RIC)
237        .x_km(0.5)
238        .y_km(0.5)
239        .z_km(0.5)
240        .vx_km_s(5e-3)
241        .vy_km_s(5e-3)
242        .vz_km_s(5e-3)
243        .build();
244
245    // Build the filter initial estimate, which we will reuse in the filter.
246    let initial_estimate = sc.to_estimate()?;
247
248    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
249
250    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
251    let process_noise = ProcessNoise3D::from_velocity_km_s(
252        &[1.8e-9, 1.8e-9, 1.8e-9],
253        1 * Unit::Hour,
254        10 * Unit::Minute,
255        None,
256    );
257
258    println!("{process_noise}");
259
260    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
261    let odp = SpacecraftKalmanOD::new(
262        setup,
263        KalmanVariant::ReferenceUpdate,
264        Some(ResidRejectCrit::default()),
265        devices,
266        almanac.clone(),
267    )
268    .with_process_noise(process_noise);
269
270    let od_sol = odp.process_arc(initial_estimate, &arc)?;
271
272    let ric_err = traj_as_flown
273        .at(od_sol.estimates.last().unwrap().epoch())?
274        .orbit
275        .ric_difference(&od_sol.estimates.last().unwrap().orbital_state())?;
276    println!("== RIC at end ==");
277    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
278    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
279
280    println!(
281        "Num residuals rejected: #{}",
282        od_sol.rejected_residuals().len()
283    );
284    println!(
285        "Percentage within +/-3: {}",
286        od_sol.residual_ratio_within_threshold(3.0).unwrap()
287    );
288    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289
290    od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
291
292    // In our case, we have the truth trajectory from NASA.
293    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
294    // Export the OD trajectory first.
295    let od_trajectory = od_sol.to_traj()?;
296    // Build the RIC difference.
297    od_trajectory.ric_diff_to_parquet(
298        &traj_as_flown,
299        "./04_lro_od_truth_error.parquet",
300        ExportCfg::default(),
301    )?;
302
303    Ok(())
304}

Trait Implementations§

Source§

impl<StateType, EstType, MsrSize, Trk> Clone for ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType> + Clone, EstType: Estimate<StateType> + Clone, MsrSize: DimName + Clone, Trk: TrackerSensitivity<StateType, StateType> + Clone, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source§

fn clone(&self) -> ODSolution<StateType, EstType, MsrSize, Trk>

Returns a copy of the value. Read more
1.0.0 · Source§

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source. Read more
Source§

impl<StateType, EstType, MsrSize, Trk> Debug for ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType> + Debug, EstType: Estimate<StateType> + Debug, MsrSize: DimName + Debug, Trk: TrackerSensitivity<StateType, StateType> + Debug, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source§

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter. Read more
Source§

impl<StateType, EstType, MsrSize, Trk> Display for ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType>, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source§

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter. Read more
Source§

impl<StateType, EstType, MsrSize, Trk> PartialEq for ODSolution<StateType, EstType, MsrSize, Trk>
where StateType: Interpolatable + Add<OVector<f64, <StateType as State>::Size>, Output = StateType>, EstType: Estimate<StateType>, MsrSize: DimName, Trk: TrackerSensitivity<StateType, StateType> + PartialEq, <DefaultAllocator as Allocator<<StateType as State>::VecLength>>::Buffer<f64>: Send, DefaultAllocator: Allocator<<StateType as State>::Size> + Allocator<<StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<StateType as State>::Size, <StateType as State>::Size> + Allocator<<StateType as State>::Size, MsrSize>,

Source§

fn eq(&self, other: &Self) -> bool

Checks that the covariances are within 1e-8 in norm, the state vectors within 1e-6, the residual ratios within 1e-4, the gains and flight-smoother consistencies within 1e-8.

1.0.0 · Source§

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

Auto Trait Implementations§

§

impl<StateType, EstType, MsrSize, Trk> Freeze for ODSolution<StateType, EstType, MsrSize, Trk>

§

impl<StateType, EstType, MsrSize, Trk> !RefUnwindSafe for ODSolution<StateType, EstType, MsrSize, Trk>

§

impl<StateType, EstType, MsrSize, Trk> !Send for ODSolution<StateType, EstType, MsrSize, Trk>

§

impl<StateType, EstType, MsrSize, Trk> !Sync for ODSolution<StateType, EstType, MsrSize, Trk>

§

impl<StateType, EstType, MsrSize, Trk> !Unpin for ODSolution<StateType, EstType, MsrSize, Trk>

§

impl<StateType, EstType, MsrSize, Trk> !UnwindSafe for ODSolution<StateType, EstType, MsrSize, Trk>

Blanket Implementations§

Source§

impl<T> Any for T
where T: 'static + ?Sized,

Source§

fn type_id(&self) -> TypeId

Gets the TypeId of self. Read more
Source§

impl<T> Borrow<T> for T
where T: ?Sized,

Source§

fn borrow(&self) -> &T

Immutably borrows from an owned value. Read more
Source§

impl<T> BorrowMut<T> for T
where T: ?Sized,

Source§

fn borrow_mut(&mut self) -> &mut T

Mutably borrows from an owned value. Read more
Source§

impl<T> CloneToUninit for T
where T: Clone,

Source§

unsafe fn clone_to_uninit(&self, dest: *mut u8)

🔬This is a nightly-only experimental API. (clone_to_uninit)
Performs copy-assignment from self to dest. Read more
Source§

impl<T> From<T> for T

Source§

fn from(t: T) -> T

Returns the argument unchanged.

§

impl<T> Instrument for T

§

fn instrument(self, span: Span) -> Instrumented<Self>

Instruments this type with the provided [Span], returning an Instrumented wrapper. Read more
§

fn in_current_span(self) -> Instrumented<Self>

Instruments this type with the current Span, returning an Instrumented wrapper. Read more
Source§

impl<T, U> Into<U> for T
where U: From<T>,

Source§

fn into(self) -> U

Calls U::from(self).

That is, this conversion is whatever the implementation of From<T> for U chooses to do.

Source§

impl<T> IntoEither for T

Source§

fn into_either(self, into_left: bool) -> Either<Self, Self>

Converts self into a Left variant of Either<Self, Self> if into_left is true. Converts self into a Right variant of Either<Self, Self> otherwise. Read more
Source§

fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
where F: FnOnce(&Self) -> bool,

Converts self into a Left variant of Either<Self, Self> if into_left(&self) returns true. Converts self into a Right variant of Either<Self, Self> otherwise. Read more
§

impl<T> Pointable for T

§

const ALIGN: usize

The alignment of pointer.
§

type Init = T

The type for initializers.
§

unsafe fn init(init: <T as Pointable>::Init) -> usize

Initializes a with the given initializer. Read more
§

unsafe fn deref<'a>(ptr: usize) -> &'a T

Dereferences the given pointer. Read more
§

unsafe fn deref_mut<'a>(ptr: usize) -> &'a mut T

Mutably dereferences the given pointer. Read more
§

unsafe fn drop(ptr: usize)

Drops the object pointed to by the given pointer. Read more
Source§

impl<T> Same for T

Source§

type Output = T

Should always be Self
§

impl<SS, SP> SupersetOf<SS> for SP
where SS: SubsetOf<SP>,

§

fn to_subset(&self) -> Option<SS>

The inverse inclusion map: attempts to construct self from the equivalent element of its superset. Read more
§

fn is_in_subset(&self) -> bool

Checks if self is actually part of its subset T (and can be converted to it).
§

fn to_subset_unchecked(&self) -> SS

Use with care! Same as self.to_subset but without any property checks. Always succeeds.
§

fn from_subset(element: &SS) -> SP

The inclusion map: converts self to the equivalent element of its superset.
Source§

impl<T> ToOwned for T
where T: Clone,

Source§

type Owned = T

The resulting type after obtaining ownership.
Source§

fn to_owned(&self) -> T

Creates owned data from borrowed data, usually by cloning. Read more
Source§

fn clone_into(&self, target: &mut T)

Uses borrowed data to replace owned data, usually by cloning. Read more
Source§

impl<T> ToString for T
where T: Display + ?Sized,

Source§

fn to_string(&self) -> String

Converts the given value to a String. Read more
Source§

impl<T, U> TryFrom<U> for T
where U: Into<T>,

Source§

type Error = Infallible

The type returned in the event of a conversion error.
Source§

fn try_from(value: U) -> Result<T, <T as TryFrom<U>>::Error>

Performs the conversion.
Source§

impl<T, U> TryInto<U> for T
where U: TryFrom<T>,

Source§

type Error = <U as TryFrom<T>>::Error

The type returned in the event of a conversion error.
Source§

fn try_into(self) -> Result<U, <U as TryFrom<T>>::Error>

Performs the conversion.
§

impl<V, T> VZip<V> for T
where V: MultiLane<T>,

§

fn vzip(self) -> V

§

impl<T> WithSubscriber for T

§

fn with_subscriber<S>(self, subscriber: S) -> WithDispatch<Self>
where S: Into<Dispatch>,

Attaches the provided Subscriber to this type, returning a [WithDispatch] wrapper. Read more
§

fn with_current_subscriber(self) -> WithDispatch<Self>

Attaches the current default Subscriber to this type, returning a [WithDispatch] wrapper. Read more
§

impl<T> Allocation for T
where T: RefUnwindSafe + Send + Sync,

§

impl<T> ErasedDestructor for T
where T: 'static,

§

impl<T> MaybeSendSync for T

Source§

impl<T> Scalar for T
where T: 'static + Clone + PartialEq + Debug,