Struct Frame
pub struct Frame {
pub ephemeris_id: i32,
pub orientation_id: i32,
pub mu_km3_s2: Option<f64>,
pub shape: Option<Ellipsoid>,
}Expand description
A Frame uniquely defined by its ephemeris center and orientation. Refer to FrameDetail for frames combined with parameters.
:type ephemeris_id: int :type orientation_id: int :type mu_km3_s2: float, optional :type shape: Ellipsoid, optional :rtype: Frame
Fields§
§ephemeris_id: i32§orientation_id: i32§mu_km3_s2: Option<f64>Gravity parameter of this frame, only defined on celestial frames
shape: Option<Ellipsoid>Shape of the geoid of this frame, only defined on geodetic frames
Implementations§
§impl Frame
impl Frame
pub const fn new(ephemeris_id: i32, orientation_id: i32) -> Frame
pub const fn new(ephemeris_id: i32, orientation_id: i32) -> Frame
Constructs a new frame given its ephemeris and orientations IDs, without defining anything else (so this is not a valid celestial frame, although the data could be populated later).
pub const fn from_ephem_j2000(ephemeris_id: i32) -> Frame
pub const fn from_ephem_j2000(ephemeris_id: i32) -> Frame
Examples found in repository?
More examples
35fn main() -> Result<(), Box<dyn Error>> {
36 pel::init();
37
38 // ====================== //
39 // === ALMANAC SET UP === //
40 // ====================== //
41
42 // Dynamics models require planetary constants and ephemerides to be defined.
43 // Let's start by grabbing those by using ANISE's MetaAlmanac.
44
45 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
46 .iter()
47 .collect();
48
49 let meta = data_folder.join("lro-dynamics.dhall");
50
51 // Load this ephem in the general Almanac we're using for this analysis.
52 let mut almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
53 .map_err(Box::new)?
54 .process(true)
55 .map_err(Box::new)?;
56
57 let mut moon_pc = almanac.get_planetary_data_from_id(MOON).unwrap();
58 moon_pc.mu_km3_s2 = 4902.74987;
59 almanac.set_planetary_data_from_id(MOON, moon_pc).unwrap();
60
61 let mut earth = almanac.get_planetary_data_from_id(EARTH).unwrap();
62 earth.mu_km3_s2 = 398600.436;
63 almanac.set_planetary_data_from_id(EARTH, earth).unwrap();
64
65 // Save this new kernel for reuse.
66 // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
67 almanac
68 .planetary_data
69 .values()
70 .next()
71 .unwrap()
72 .save_as(&data_folder.join("lro-specific.pca"), true)?;
73
74 // Lock the almanac (an Arc is a read only structure).
75 let almanac = Arc::new(almanac);
76
77 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
78 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
79 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
80 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
81 let lro_frame = Frame::from_ephem_j2000(-85);
82
83 // To build the trajectory we need to provide a spacecraft template.
84 let sc_template = Spacecraft::builder()
85 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
86 .srp(SRPData {
87 // SRP configuration is arbitrary, but we will be estimating it anyway.
88 area_m2: 3.9 * 2.7,
89 coeff_reflectivity: 0.96,
90 })
91 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
92 .build();
93 // Now we can build the trajectory from the BSP file.
94 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
95 let traj_as_flown = Traj::from_bsp(
96 lro_frame,
97 MOON_J2000,
98 almanac.clone(),
99 sc_template,
100 5.seconds(),
101 Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
102 Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
103 Aberration::LT,
104 Some("LRO".to_string()),
105 )?;
106
107 println!("{traj_as_flown}");
108
109 // ====================== //
110 // === MODEL MATCHING === //
111 // ====================== //
112
113 // Set up the spacecraft dynamics.
114
115 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
116 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
117 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
118
119 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
120 // We're using the GRAIL JGGRX model.
121 let mut jggrx_meta = MetaFile {
122 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
123 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
124 };
125 // And let's download it if we don't have it yet.
126 jggrx_meta.process(true)?;
127
128 // Build the spherical harmonics.
129 // The harmonics must be computed in the body fixed frame.
130 // We're using the long term prediction of the Moon principal axes frame.
131 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
132 let sph_harmonics = Harmonics::from_stor(
133 almanac.frame_info(moon_pa_frame)?,
134 HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
135 );
136
137 // Include the spherical harmonics into the orbital dynamics.
138 orbital_dyn.accel_models.push(sph_harmonics);
139
140 // We define the solar radiation pressure, using the default solar flux and accounting only
141 // for the eclipsing caused by the Earth and Moon.
142 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
143 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
144
145 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
146 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
147 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
148
149 println!("{dynamics}");
150
151 // Now we can build the propagator.
152 let setup = Propagator::default_dp78(dynamics.clone());
153
154 // For reference, let's build the trajectory with Nyx's models from that LRO state.
155 let (sim_final, traj_as_sim) = setup
156 .with(*traj_as_flown.first(), almanac.clone())
157 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
158
159 println!("SIM INIT: {:x}", traj_as_flown.first());
160 println!("SIM FINAL: {sim_final:x}");
161 // Compute RIC difference between SIM and LRO ephem
162 let sim_lro_delta = sim_final
163 .orbit
164 .ric_difference(&traj_as_flown.last().orbit)?;
165 println!("{traj_as_sim}");
166 println!(
167 "SIM v LRO - RIC Position (m): {:.3}",
168 sim_lro_delta.radius_km * 1e3
169 );
170 println!(
171 "SIM v LRO - RIC Velocity (m/s): {:.3}",
172 sim_lro_delta.velocity_km_s * 1e3
173 );
174
175 traj_as_sim.ric_diff_to_parquet(
176 &traj_as_flown,
177 "./data/04_output/04_lro_sim_truth_error.parquet",
178 ExportCfg::default(),
179 )?;
180
181 // ==================== //
182 // === OD SIMULATOR === //
183 // ==================== //
184
185 // 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
186 // and the truth LRO state.
187
188 // Therefore, we will actually run an estimation from a dispersed LRO state.
189 // The sc_seed is the true LRO state from the BSP.
190 let sc_seed = *traj_as_flown.first();
191
192 // Load the Deep Space Network ground stations.
193 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
194 let ground_station_file: PathBuf = [
195 env!("CARGO_MANIFEST_DIR"),
196 "examples",
197 "04_lro_od",
198 "dsn-network.yaml",
199 ]
200 .iter()
201 .collect();
202
203 let devices = GroundStation::load_named(ground_station_file)?;
204
205 let mut proc_devices = devices.clone();
206
207 // Increase the noise in the devices to accept more measurements.
208 for gs in proc_devices.values_mut() {
209 if let Some(noise) = &mut gs
210 .stochastic_noises
211 .as_mut()
212 .unwrap()
213 .get_mut(&MeasurementType::Range)
214 {
215 *noise.white_noise.as_mut().unwrap() *= 3.0;
216 }
217 }
218
219 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
220 // Nyx can build a tracking schedule for you based on the first station with access.
221 let trkconfg_yaml: PathBuf = [
222 env!("CARGO_MANIFEST_DIR"),
223 "examples",
224 "04_lro_od",
225 "tracking-cfg.yaml",
226 ]
227 .iter()
228 .collect();
229
230 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
231
232 // Build the tracking arc simulation to generate a "standard measurement".
233 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
234 devices.clone(),
235 traj_as_flown.clone(),
236 configs,
237 123, // Set a seed for reproducibility
238 )?;
239
240 trk.build_schedule(almanac.clone())?;
241 let arc = trk.generate_measurements(almanac.clone())?;
242 // Save the simulated tracking data
243 arc.to_parquet_simple("./data/04_output/04_lro_simulated_tracking.parquet")?;
244
245 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
246 println!("{arc}");
247
248 // Now that we have simulated measurements, we'll run the orbit determination.
249
250 // ===================== //
251 // === OD ESTIMATION === //
252 // ===================== //
253
254 let sc = SpacecraftUncertainty::builder()
255 .nominal(sc_seed)
256 .frame(LocalFrame::RIC)
257 .x_km(0.5)
258 .y_km(0.5)
259 .z_km(0.5)
260 .vx_km_s(5e-3)
261 .vy_km_s(5e-3)
262 .vz_km_s(5e-3)
263 .build();
264
265 // Build the filter initial estimate, which we will reuse in the filter.
266 let mut initial_estimate = sc.to_estimate()?;
267 initial_estimate.covar *= 3.0;
268
269 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
270
271 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
272 let process_noise = ProcessNoise3D::from_velocity_km_s(
273 &[1e-10, 1e-10, 1e-10],
274 1 * Unit::Hour,
275 10 * Unit::Minute,
276 None,
277 );
278
279 println!("{process_noise}");
280
281 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
282 let odp = SpacecraftKalmanOD::new(
283 setup,
284 KalmanVariant::ReferenceUpdate,
285 Some(ResidRejectCrit::default()),
286 proc_devices,
287 almanac.clone(),
288 )
289 .with_process_noise(process_noise);
290
291 let od_sol = odp.process_arc(initial_estimate, &arc)?;
292
293 let final_est = od_sol.estimates.last().unwrap();
294
295 println!("{final_est}");
296
297 let ric_err = traj_as_flown
298 .at(final_est.epoch())?
299 .orbit
300 .ric_difference(&final_est.orbital_state())?;
301 println!("== RIC at end ==");
302 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
303 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
304
305 println!(
306 "Num residuals rejected: #{}",
307 od_sol.rejected_residuals().len()
308 );
309 println!(
310 "Percentage within +/-3: {}",
311 od_sol.residual_ratio_within_threshold(3.0).unwrap()
312 );
313 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
314
315 od_sol.to_parquet(
316 "./data/04_output/04_lro_od_results.parquet",
317 ExportCfg::default(),
318 )?;
319
320 // Create the ephemeris
321 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
322 let ephem_start = ephem.start_epoch().unwrap();
323 let ephem_end = ephem.end_epoch().unwrap();
324 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
325 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
326 ephem
327 .covar_at(
328 epoch,
329 anise::ephemerides::ephemeris::LocalFrame::RIC,
330 &almanac,
331 )
332 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
333 }
334 // Export as BSP!
335 ephem
336 .write_spice_bsp(-85, "./data/04_output/04_lro_rebuilt.bsp", None)
337 .expect("could not built BSP");
338 let new_almanac = Almanac::default()
339 .load("./data/04_output/04_lro_rebuilt.bsp")
340 .unwrap();
341 new_almanac.describe(None, None, None, None, None, None, None, None);
342 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
343
344 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
345 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
346
347 // In our case, we have the truth trajectory from NASA.
348 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
349 // Export the OD trajectory first.
350 let od_trajectory = od_sol.to_traj()?;
351 // Build the RIC difference.
352 od_trajectory.ric_diff_to_parquet(
353 &traj_as_flown,
354 "./data/04_output/04_lro_od_truth_error.parquet",
355 ExportCfg::default(),
356 )?;
357
358 Ok(())
359}pub const fn from_orient_ssb(orientation_id: i32) -> Frame
pub fn from_name(center: &str, ref_frame: &str) -> Result<Frame, AlmanacError>
pub fn from_name(center: &str, ref_frame: &str) -> Result<Frame, AlmanacError>
Attempts to create a new frame from its center and reference frame name. This function is compatible with the CCSDS OEM names.
pub fn with_ellipsoid(self, shape: Ellipsoid) -> Frame
pub fn with_ellipsoid(self, shape: Ellipsoid) -> Frame
Define Ellipsoid shape and return a new Frame
§impl Frame
impl Frame
pub const fn with_ephem(&self, new_ephem_id: i32) -> Frame
pub const fn with_ephem(&self, new_ephem_id: i32) -> Frame
Returns a copy of this Frame whose ephemeris ID is set to the provided ID
:type new_ephem_id: int :rtype: Frame
pub const fn with_orient(&self, new_orient_id: i32) -> Frame
pub const fn with_orient(&self, new_orient_id: i32) -> Frame
Returns a copy of this Frame whose orientation ID is set to the provided ID
:type new_orient_id: int :rtype: Frame
Examples found in repository?
35fn main() -> Result<(), Box<dyn Error>> {
36 pel::init();
37
38 // ====================== //
39 // === ALMANAC SET UP === //
40 // ====================== //
41
42 // Dynamics models require planetary constants and ephemerides to be defined.
43 // Let's start by grabbing those by using ANISE's MetaAlmanac.
44
45 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
46 .iter()
47 .collect();
48
49 let meta = data_folder.join("lro-dynamics.dhall");
50
51 // Load this ephem in the general Almanac we're using for this analysis.
52 let mut almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
53 .map_err(Box::new)?
54 .process(true)
55 .map_err(Box::new)?;
56
57 let mut moon_pc = almanac.get_planetary_data_from_id(MOON).unwrap();
58 moon_pc.mu_km3_s2 = 4902.74987;
59 almanac.set_planetary_data_from_id(MOON, moon_pc).unwrap();
60
61 let mut earth = almanac.get_planetary_data_from_id(EARTH).unwrap();
62 earth.mu_km3_s2 = 398600.436;
63 almanac.set_planetary_data_from_id(EARTH, earth).unwrap();
64
65 // Save this new kernel for reuse.
66 // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
67 almanac
68 .planetary_data
69 .values()
70 .next()
71 .unwrap()
72 .save_as(&data_folder.join("lro-specific.pca"), true)?;
73
74 // Lock the almanac (an Arc is a read only structure).
75 let almanac = Arc::new(almanac);
76
77 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
78 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
79 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
80 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
81 let lro_frame = Frame::from_ephem_j2000(-85);
82
83 // To build the trajectory we need to provide a spacecraft template.
84 let sc_template = Spacecraft::builder()
85 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
86 .srp(SRPData {
87 // SRP configuration is arbitrary, but we will be estimating it anyway.
88 area_m2: 3.9 * 2.7,
89 coeff_reflectivity: 0.96,
90 })
91 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
92 .build();
93 // Now we can build the trajectory from the BSP file.
94 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
95 let traj_as_flown = Traj::from_bsp(
96 lro_frame,
97 MOON_J2000,
98 almanac.clone(),
99 sc_template,
100 5.seconds(),
101 Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
102 Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
103 Aberration::LT,
104 Some("LRO".to_string()),
105 )?;
106
107 println!("{traj_as_flown}");
108
109 // ====================== //
110 // === MODEL MATCHING === //
111 // ====================== //
112
113 // Set up the spacecraft dynamics.
114
115 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
116 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
117 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
118
119 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
120 // We're using the GRAIL JGGRX model.
121 let mut jggrx_meta = MetaFile {
122 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
123 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
124 };
125 // And let's download it if we don't have it yet.
126 jggrx_meta.process(true)?;
127
128 // Build the spherical harmonics.
129 // The harmonics must be computed in the body fixed frame.
130 // We're using the long term prediction of the Moon principal axes frame.
131 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
132 let sph_harmonics = Harmonics::from_stor(
133 almanac.frame_info(moon_pa_frame)?,
134 HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
135 );
136
137 // Include the spherical harmonics into the orbital dynamics.
138 orbital_dyn.accel_models.push(sph_harmonics);
139
140 // We define the solar radiation pressure, using the default solar flux and accounting only
141 // for the eclipsing caused by the Earth and Moon.
142 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
143 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
144
145 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
146 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
147 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
148
149 println!("{dynamics}");
150
151 // Now we can build the propagator.
152 let setup = Propagator::default_dp78(dynamics.clone());
153
154 // For reference, let's build the trajectory with Nyx's models from that LRO state.
155 let (sim_final, traj_as_sim) = setup
156 .with(*traj_as_flown.first(), almanac.clone())
157 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
158
159 println!("SIM INIT: {:x}", traj_as_flown.first());
160 println!("SIM FINAL: {sim_final:x}");
161 // Compute RIC difference between SIM and LRO ephem
162 let sim_lro_delta = sim_final
163 .orbit
164 .ric_difference(&traj_as_flown.last().orbit)?;
165 println!("{traj_as_sim}");
166 println!(
167 "SIM v LRO - RIC Position (m): {:.3}",
168 sim_lro_delta.radius_km * 1e3
169 );
170 println!(
171 "SIM v LRO - RIC Velocity (m/s): {:.3}",
172 sim_lro_delta.velocity_km_s * 1e3
173 );
174
175 traj_as_sim.ric_diff_to_parquet(
176 &traj_as_flown,
177 "./data/04_output/04_lro_sim_truth_error.parquet",
178 ExportCfg::default(),
179 )?;
180
181 // ==================== //
182 // === OD SIMULATOR === //
183 // ==================== //
184
185 // 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
186 // and the truth LRO state.
187
188 // Therefore, we will actually run an estimation from a dispersed LRO state.
189 // The sc_seed is the true LRO state from the BSP.
190 let sc_seed = *traj_as_flown.first();
191
192 // Load the Deep Space Network ground stations.
193 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
194 let ground_station_file: PathBuf = [
195 env!("CARGO_MANIFEST_DIR"),
196 "examples",
197 "04_lro_od",
198 "dsn-network.yaml",
199 ]
200 .iter()
201 .collect();
202
203 let devices = GroundStation::load_named(ground_station_file)?;
204
205 let mut proc_devices = devices.clone();
206
207 // Increase the noise in the devices to accept more measurements.
208 for gs in proc_devices.values_mut() {
209 if let Some(noise) = &mut gs
210 .stochastic_noises
211 .as_mut()
212 .unwrap()
213 .get_mut(&MeasurementType::Range)
214 {
215 *noise.white_noise.as_mut().unwrap() *= 3.0;
216 }
217 }
218
219 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
220 // Nyx can build a tracking schedule for you based on the first station with access.
221 let trkconfg_yaml: PathBuf = [
222 env!("CARGO_MANIFEST_DIR"),
223 "examples",
224 "04_lro_od",
225 "tracking-cfg.yaml",
226 ]
227 .iter()
228 .collect();
229
230 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
231
232 // Build the tracking arc simulation to generate a "standard measurement".
233 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
234 devices.clone(),
235 traj_as_flown.clone(),
236 configs,
237 123, // Set a seed for reproducibility
238 )?;
239
240 trk.build_schedule(almanac.clone())?;
241 let arc = trk.generate_measurements(almanac.clone())?;
242 // Save the simulated tracking data
243 arc.to_parquet_simple("./data/04_output/04_lro_simulated_tracking.parquet")?;
244
245 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
246 println!("{arc}");
247
248 // Now that we have simulated measurements, we'll run the orbit determination.
249
250 // ===================== //
251 // === OD ESTIMATION === //
252 // ===================== //
253
254 let sc = SpacecraftUncertainty::builder()
255 .nominal(sc_seed)
256 .frame(LocalFrame::RIC)
257 .x_km(0.5)
258 .y_km(0.5)
259 .z_km(0.5)
260 .vx_km_s(5e-3)
261 .vy_km_s(5e-3)
262 .vz_km_s(5e-3)
263 .build();
264
265 // Build the filter initial estimate, which we will reuse in the filter.
266 let mut initial_estimate = sc.to_estimate()?;
267 initial_estimate.covar *= 3.0;
268
269 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
270
271 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
272 let process_noise = ProcessNoise3D::from_velocity_km_s(
273 &[1e-10, 1e-10, 1e-10],
274 1 * Unit::Hour,
275 10 * Unit::Minute,
276 None,
277 );
278
279 println!("{process_noise}");
280
281 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
282 let odp = SpacecraftKalmanOD::new(
283 setup,
284 KalmanVariant::ReferenceUpdate,
285 Some(ResidRejectCrit::default()),
286 proc_devices,
287 almanac.clone(),
288 )
289 .with_process_noise(process_noise);
290
291 let od_sol = odp.process_arc(initial_estimate, &arc)?;
292
293 let final_est = od_sol.estimates.last().unwrap();
294
295 println!("{final_est}");
296
297 let ric_err = traj_as_flown
298 .at(final_est.epoch())?
299 .orbit
300 .ric_difference(&final_est.orbital_state())?;
301 println!("== RIC at end ==");
302 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
303 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
304
305 println!(
306 "Num residuals rejected: #{}",
307 od_sol.rejected_residuals().len()
308 );
309 println!(
310 "Percentage within +/-3: {}",
311 od_sol.residual_ratio_within_threshold(3.0).unwrap()
312 );
313 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
314
315 od_sol.to_parquet(
316 "./data/04_output/04_lro_od_results.parquet",
317 ExportCfg::default(),
318 )?;
319
320 // Create the ephemeris
321 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
322 let ephem_start = ephem.start_epoch().unwrap();
323 let ephem_end = ephem.end_epoch().unwrap();
324 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
325 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
326 ephem
327 .covar_at(
328 epoch,
329 anise::ephemerides::ephemeris::LocalFrame::RIC,
330 &almanac,
331 )
332 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
333 }
334 // Export as BSP!
335 ephem
336 .write_spice_bsp(-85, "./data/04_output/04_lro_rebuilt.bsp", None)
337 .expect("could not built BSP");
338 let new_almanac = Almanac::default()
339 .load("./data/04_output/04_lro_rebuilt.bsp")
340 .unwrap();
341 new_almanac.describe(None, None, None, None, None, None, None, None);
342 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
343
344 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
345 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
346
347 // In our case, we have the truth trajectory from NASA.
348 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
349 // Export the OD trajectory first.
350 let od_trajectory = od_sol.to_traj()?;
351 // Build the RIC difference.
352 od_trajectory.ric_diff_to_parquet(
353 &traj_as_flown,
354 "./data/04_output/04_lro_od_truth_error.parquet",
355 ExportCfg::default(),
356 )?;
357
358 Ok(())
359}pub const fn is_celestial(&self) -> bool
pub const fn is_celestial(&self) -> bool
Returns whether this is a celestial frame
:rtype: bool
pub const fn is_geodetic(&self) -> bool
pub const fn is_geodetic(&self) -> bool
Returns whether this is a geodetic frame
:rtype: bool
pub const fn ephem_origin_id_match(&self, other_id: i32) -> bool
pub const fn ephem_origin_id_match(&self, other_id: i32) -> bool
Returns true if the ephemeris origin is equal to the provided ID
:type other_id: int :rtype: bool
pub const fn orient_origin_id_match(&self, other_id: i32) -> bool
pub const fn orient_origin_id_match(&self, other_id: i32) -> bool
Returns true if the orientation origin is equal to the provided ID
:type other_id: int :rtype: bool
pub const fn ephem_origin_match(&self, other: Frame) -> bool
pub const fn ephem_origin_match(&self, other: Frame) -> bool
Returns true if the ephemeris origin is equal to the provided frame
:type other: Frame :rtype: bool
pub const fn orient_origin_match(&self, other: Frame) -> bool
pub const fn orient_origin_match(&self, other: Frame) -> bool
Returns true if the orientation origin is equal to the provided frame
:type other: Frame :rtype: bool
pub fn strip(&mut self)
pub fn strip(&mut self)
Removes the graviational parameter and the shape information from this frame. Use this to prevent astrodynamical computations.
:rtype: None
pub fn mu_km3_s2(&self) -> Result<f64, PhysicsError>
pub fn mu_km3_s2(&self) -> Result<f64, PhysicsError>
Returns the gravitational parameters of this frame, if defined
:rtype: float
pub fn with_mu_km3_s2(&self, mu_km3_s2: f64) -> Frame
pub fn with_mu_km3_s2(&self, mu_km3_s2: f64) -> Frame
Returns a copy of this frame with the graviational parameter set to the new value.
:type mu_km3_s2: float :rtype: Frame
pub fn mean_equatorial_radius_km(&self) -> Result<f64, PhysicsError>
pub fn mean_equatorial_radius_km(&self) -> Result<f64, PhysicsError>
Returns the mean equatorial radius in km, if defined
:rtype: float
pub fn semi_major_radius_km(&self) -> Result<f64, PhysicsError>
pub fn semi_major_radius_km(&self) -> Result<f64, PhysicsError>
Returns the semi major radius of the tri-axial ellipoid shape of this frame, if defined
:rtype: float
pub fn flattening(&self) -> Result<f64, PhysicsError>
pub fn flattening(&self) -> Result<f64, PhysicsError>
Returns the flattening ratio (unitless)
:rtype: float
pub fn polar_radius_km(&self) -> Result<f64, PhysicsError>
pub fn polar_radius_km(&self) -> Result<f64, PhysicsError>
Returns the polar radius in km, if defined
:rtype: float
Trait Implementations§
§impl<'de> Deserialize<'de> for Frame
impl<'de> Deserialize<'de> for Frame
§fn deserialize<__D>(
__deserializer: __D,
) -> Result<Frame, <__D as Deserializer<'de>>::Error>where
__D: Deserializer<'de>,
fn deserialize<__D>(
__deserializer: __D,
) -> Result<Frame, <__D as Deserializer<'de>>::Error>where
__D: Deserializer<'de>,
§impl Serialize for Frame
impl Serialize for Frame
§fn serialize<__S>(
&self,
__serializer: __S,
) -> Result<<__S as Serializer>::Ok, <__S as Serializer>::Error>where
__S: Serializer,
fn serialize<__S>(
&self,
__serializer: __S,
) -> Result<<__S as Serializer>::Ok, <__S as Serializer>::Error>where
__S: Serializer,
§impl StaticType for Frame
impl StaticType for Frame
§fn static_type() -> SimpleType
fn static_type() -> SimpleType
impl Copy for Frame
impl StructuralPartialEq for Frame
Auto Trait Implementations§
impl Freeze for Frame
impl RefUnwindSafe for Frame
impl Send for Frame
impl Sync for Frame
impl Unpin for Frame
impl UnwindSafe for Frame
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,
Source§impl<T> FromDhall for Twhere
T: DeserializeOwned,
impl<T> FromDhall for Twhere
T: DeserializeOwned,
fn from_dhall(v: &Value) -> Result<T, Error>
§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.