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 toNone
after running the smoother.filter_smoother_ratios
: Filter-smoother consistency ratios. These are set toNone
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>,
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>,
Sourcepub fn to_parquet<P: AsRef<Path>>(
&self,
path: P,
cfg: ExportCfg,
) -> Result<PathBuf, ODError>
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?
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
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>,
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>,
Sourcepub fn unique(&self) -> IndexSet<(String, MeasurementType)>
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)}
.
Sourcepub fn drop_time_updates(self) -> Self
pub fn drop_time_updates(self) -> Self
Returns this OD solution without any time update
Sourcepub fn filter_by_msr_type(self, msr_type: MeasurementType) -> Self
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.
Sourcepub fn filter_by_tracker(self, tracker: String) -> Self
pub fn filter_by_tracker(self, tracker: String) -> Self
Returns this OD solution with only data from the desired tracker, dropping all time updates.
Sourcepub fn exclude_tracker(self, excluded_tracker: String) -> Self
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
Sourcepub fn split(self) -> Vec<Self>
pub fn split(self) -> Vec<Self>
Split this OD solution per tracker and per measurement type, dropping all time updates.
Sourcepub fn merge(self, other: Self) -> Self
pub fn merge(self, other: Self) -> Self
Merge this OD solution with another one, returning a new OD solution.
pub fn at(&self, epoch: Epoch) -> Option<ODRecord<StateType, EstType, MsrSize>>
Source§impl<MsrSize, Trk> ODSolution<Spacecraft, KfEstimate<Spacecraft>, MsrSize, Trk>where
MsrSize: DimName + Debug + Clone,
Trk: TrackerSensitivity<Spacecraft, Spacecraft> + Clone,
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>,
<DefaultAllocator as Allocator<<Spacecraft as State>::VecLength>>::Buffer<f64>: Send,
<DefaultAllocator as Allocator<<Spacecraft as State>::Size>>::Buffer<f64>: Copy,
<DefaultAllocator as Allocator<<Spacecraft as State>::Size, <Spacecraft as State>::Size>>::Buffer<f64>: Copy,
impl<MsrSize, Trk> ODSolution<Spacecraft, KfEstimate<Spacecraft>, MsrSize, Trk>where
MsrSize: DimName + Debug + Clone,
Trk: TrackerSensitivity<Spacecraft, Spacecraft> + Clone,
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>,
<DefaultAllocator as Allocator<<Spacecraft as State>::VecLength>>::Buffer<f64>: Send,
<DefaultAllocator as Allocator<<Spacecraft as State>::Size>>::Buffer<f64>: Copy,
<DefaultAllocator as Allocator<<Spacecraft as State>::Size, <Spacecraft as State>::Size>>::Buffer<f64>: Copy,
Sourcepub fn from_parquet<P: AsRef<Path>>(
path: P,
devices: BTreeMap<String, Trk>,
) -> Result<Self, InputOutputError>
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>,
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>,
Sourcepub fn smooth(self, almanac: Arc<Almanac>) -> Result<Self, ODError>
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:
- Gains will be scrubbed because the smoother process does not recompute the gain.
- Prefit residuals, ratios, and measurement covariances are not updated, as these depend on the filtering process.
- 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>,
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>,
Sourcepub fn rms_prefit_residuals(&self) -> f64
pub fn rms_prefit_residuals(&self) -> f64
Returns the root mean square of the prefit residuals
Sourcepub fn rms_postfit_residuals(&self) -> f64
pub fn rms_postfit_residuals(&self) -> f64
Returns the root mean square of the postfit residuals
Sourcepub fn rms_residual_ratios(&self) -> f64
pub fn rms_residual_ratios(&self) -> f64
Returns the root mean square of the prefit residual ratios
Sourcepub fn residual_ratio_within_threshold(
&self,
threshold: f64,
) -> Result<f64, ODError>
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?
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}
Sourcepub fn ks_test_normality(&self) -> Result<f64, ODError>
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.
Sourcepub fn is_normal(&self, alpha: Option<f64>) -> Result<bool, ODError>
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?
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>,
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>,
pub fn new( devices: BTreeMap<String, Trk>, measurement_types: IndexSet<MeasurementType>, ) -> Self
Sourcepub fn results(
&self,
) -> Zip<Iter<'_, EstType>, Iter<'_, Option<Residual<MsrSize>>>>
pub fn results( &self, ) -> Zip<Iter<'_, EstType>, Iter<'_, Option<Residual<MsrSize>>>>
Returns a zipper iterator on the estimates and the associated residuals.
Sourcepub fn is_filter_run(&self) -> bool
pub fn is_filter_run(&self) -> bool
Returns True if this is the result of a filter run
Sourcepub fn is_smoother_run(&self) -> bool
pub fn is_smoother_run(&self) -> bool
Returns True if this is the result of a smoother run
Sourcepub fn to_traj(&self) -> Result<Traj<StateType>, NyxError>
pub fn to_traj(&self) -> Result<Traj<StateType>, NyxError>
Builds the navigation trajectory for the estimated state only
Examples found in repository?
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}
Sourcepub fn accepted_residuals(&self) -> Vec<Residual<MsrSize>>
pub fn accepted_residuals(&self) -> Vec<Residual<MsrSize>>
Returns the accepted residuals.
Sourcepub fn rejected_residuals(&self) -> Vec<Residual<MsrSize>>
pub fn rejected_residuals(&self) -> Vec<Residual<MsrSize>>
Returns the rejected residuals.
Examples found in repository?
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>,
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>
fn clone(&self) -> ODSolution<StateType, EstType, MsrSize, Trk>
1.0.0 · Source§fn clone_from(&mut self, source: &Self)
fn clone_from(&mut self, source: &Self)
source
. Read moreSource§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>,
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§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>,
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§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>,
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>,
Auto Trait Implementations§
impl<StateType, EstType, MsrSize, Trk> Freeze for ODSolution<StateType, EstType, MsrSize, Trk>where
DefaultAllocator: Sized,
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> BorrowMut<T> for Twhere
T: ?Sized,
impl<T> BorrowMut<T> for Twhere
T: ?Sized,
Source§fn borrow_mut(&mut self) -> &mut T
fn borrow_mut(&mut self) -> &mut T
Source§impl<T> CloneToUninit for Twhere
T: Clone,
impl<T> CloneToUninit for Twhere
T: Clone,
§impl<T> Instrument for T
impl<T> Instrument for T
§fn instrument(self, span: Span) -> Instrumented<Self>
fn instrument(self, span: Span) -> Instrumented<Self>
§fn in_current_span(self) -> Instrumented<Self>
fn in_current_span(self) -> Instrumented<Self>
Source§impl<T> IntoEither for T
impl<T> IntoEither for T
Source§fn into_either(self, into_left: bool) -> Either<Self, Self>
fn into_either(self, into_left: bool) -> Either<Self, Self>
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 moreSource§fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
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
impl<T> Pointable for T
§impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
§fn to_subset(&self) -> Option<SS>
fn to_subset(&self) -> Option<SS>
self
from the equivalent element of its
superset. Read more§fn is_in_subset(&self) -> bool
fn is_in_subset(&self) -> bool
self
is actually part of its subset T
(and can be converted to it).§fn to_subset_unchecked(&self) -> SS
fn to_subset_unchecked(&self) -> SS
self.to_subset
but without any property checks. Always succeeds.§fn from_subset(element: &SS) -> SP
fn from_subset(element: &SS) -> SP
self
to the equivalent element of its superset.