Struct Frame
pub struct Frame {
pub ephemeris_id: i32,
pub orientation_id: i32,
pub force_inertial: bool,
pub mu_km3_s2: Option<f64>,
pub shape: Option<Ellipsoid>,
pub frozen_epoch: Option<Epoch>,
}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 force_inertial: bool :type mu_km3_s2: float, optional :type shape: Ellipsoid, optional :type frozen_epoch: Epoch, optional :rtype: Frame
Fields§
§ephemeris_id: i32§orientation_id: i32§force_inertial: boolIf true, the DCM of this frame will always force its time derivative to zero, inertially fixing the frame.
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
frozen_epoch: Option<Epoch>If set, the DCM will always be evaluated at the provided epoch, freezing it in time.
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
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 // Dynamics models require planetary constants and ephemerides to be defined.
42 // Let's start by grabbing those by using ANISE's MetaAlmanac.
43
44 let output_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "../data", "04_output"]
45 .iter()
46 .collect();
47
48 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
49 .iter()
50 .collect();
51
52 let meta = data_folder.join("lro-dynamics.dhall");
53
54 // Load this ephem in the general Almanac we're using for this analysis.
55 let almanac = Arc::new(
56 MetaAlmanac::new(meta.to_string_lossy().as_ref())
57 .map_err(Box::new)?
58 .process(true)
59 .map_err(Box::new)?,
60 );
61
62 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
63 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
64 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
65 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
66 let lro_frame = Frame::from_ephem_j2000(-85);
67
68 // To build the trajectory we need to provide a spacecraft template.
69 let sc_template = Spacecraft::builder()
70 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
71 .srp(SRPData {
72 // SRP configuration is arbitrary, but we will be estimating it anyway.
73 area_m2: 3.9 * 2.7,
74 coeff_reflectivity: 0.96,
75 })
76 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
77 .build();
78 // Now we can build the trajectory from the BSP file.
79 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
80 let traj_as_flown = Traj::from_bsp(
81 lro_frame,
82 MOON_J2000,
83 &almanac,
84 sc_template,
85 5.seconds(),
86 Some(Epoch::from_str("2024-01-01 01:00:00 UTC")?),
87 Some(Epoch::from_str("2024-01-02 01:00:00 UTC")?),
88 None,
89 Some("LRO".to_string()),
90 )?;
91
92 println!("{traj_as_flown}");
93
94 // ====================== //
95 // === MODEL MATCHING === //
96 // ====================== //
97
98 // Set up the spacecraft dynamics.
99
100 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
101 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
102 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
103
104 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
105 // We're using the GRAIL JGGRX model.
106 let mut jggrx_meta = MetaFile {
107 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
108 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
109 };
110 // And let's download it if we don't have it yet.
111 jggrx_meta.process(true)?;
112
113 // Build the spherical harmonics.
114 // The harmonics must be computed in the body fixed frame.
115 // We're using the long term prediction of the Moon principal axes frame.
116 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
117 let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
118 &jggrx_meta.uri,
119 81,
120 81,
121 almanac.frame_info(moon_pa_frame)?,
122 )?);
123
124 // Include the spherical harmonics into the orbital dynamics.
125 orbital_dyn.accel_models.push(sph_harmonics);
126
127 // We define the solar radiation pressure, using the default solar flux and accounting only
128 // for the eclipsing caused by the Earth and Moon.
129 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
130 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;
131
132 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
133 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
134 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
135
136 println!("{dynamics}");
137
138 // Now we can build the propagator.
139 let setup = Propagator::default_dp78(dynamics.clone());
140
141 // For reference, let's build the trajectory with Nyx's models from that LRO state.
142 let (sim_final, traj_as_sim) = setup
143 .with(*traj_as_flown.first(), almanac.clone())
144 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
145
146 println!("SIM INIT: {:x}", traj_as_flown.first());
147 println!("SIM FINAL: {sim_final:x}");
148 // Compute RIC difference between SIM and LRO ephem
149 let sim_lro_delta = sim_final
150 .orbit
151 .ric_difference(&traj_as_flown.last().orbit)?;
152 println!("{traj_as_sim}");
153 println!(
154 "SIM v LRO - RIC Position (m): {:.3}",
155 sim_lro_delta.radius_km * 1e3
156 );
157 println!(
158 "SIM v LRO - RIC Velocity (m/s): {:.3}",
159 sim_lro_delta.velocity_km_s * 1e3
160 );
161
162 traj_as_sim.ric_diff_to_parquet(
163 &traj_as_flown,
164 output_folder.join("./04_lro_sim_truth_error.parquet"),
165 ExportCfg::default(),
166 )?;
167
168 // ==================== //
169 // === OD SIMULATOR === //
170 // ==================== //
171
172 // 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
173 // and the truth LRO state.
174
175 // Therefore, we will actually run an estimation from a dispersed LRO state.
176 // The sc_seed is the true LRO state from the BSP.
177 let sc_seed = *traj_as_flown.first();
178
179 // Load the Deep Space Network ground stations.
180 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
181 let ground_station_file: PathBuf = [
182 env!("CARGO_MANIFEST_DIR"),
183 "examples",
184 "04_lro_od",
185 "dsn-network.yaml",
186 ]
187 .iter()
188 .collect();
189
190 let devices = GroundStation::load_named(ground_station_file)?;
191
192 let proc_devices = devices.clone();
193
194 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
195 // Nyx can build a tracking schedule for you based on the first station with access.
196 let trkconfg_yaml: PathBuf = [
197 env!("CARGO_MANIFEST_DIR"),
198 "examples",
199 "04_lro_od",
200 "tracking-cfg.yaml",
201 ]
202 .iter()
203 .collect();
204
205 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
206
207 // Build the tracking arc simulation to generate a "standard measurement".
208 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
209 devices.clone(),
210 traj_as_flown.clone(),
211 configs,
212 123, // Set a seed for reproducibility
213 )?;
214
215 trk.build_schedule(&almanac)?;
216 let arc = trk.generate_measurements(&almanac)?;
217 // Save the simulated tracking data
218 arc.to_parquet_simple(output_folder.join("04_lro_simulated_tracking.parquet"))?;
219
220 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
221 println!("{arc}");
222
223 // Now that we have simulated measurements, we'll run the orbit determination.
224
225 // ===================== //
226 // === OD ESTIMATION === //
227 // ===================== //
228
229 let sc = SpacecraftUncertainty::builder()
230 .nominal(sc_seed)
231 .frame(LocalFrame::RIC)
232 .x_km(0.5)
233 .y_km(0.5)
234 .z_km(0.5)
235 .vx_km_s(5e-3)
236 .vy_km_s(5e-3)
237 .vz_km_s(5e-3)
238 .build();
239
240 // Build the filter initial estimate, which we will reuse in the filter.
241 let mut initial_estimate = sc.to_estimate()?;
242 initial_estimate.covar *= 1.5;
243
244 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
245
246 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
247 let process_noise = ProcessNoise3D::from_velocity_km_s(
248 &[5e-13, 5e-13, 5e-13],
249 1 * Unit::Hour,
250 10 * Unit::Minute,
251 None,
252 );
253
254 println!("{process_noise}");
255
256 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
257 let odp = SpacecraftKalmanOD::new(
258 setup,
259 KalmanVariant::ReferenceUpdate,
260 Some(SigmaRejection::default()),
261 proc_devices,
262 almanac.clone(),
263 )
264 .with_process_noise(process_noise);
265
266 let od_sol = odp.process_arc(initial_estimate, &arc)?;
267
268 let final_est = od_sol.estimates.last().unwrap();
269
270 println!("{final_est}");
271
272 let ric_err = traj_as_flown
273 .at(final_est.epoch())?
274 .orbit
275 .ric_difference(&final_est.orbital_state())?;
276 println!("== RIC at end ==");
277 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
278 println!("RIC Velocity (m/s): {:.3}", 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 od_sol.nis_consistency(None)?.log();
290
291 od_sol.to_parquet(
292 output_folder.join("04_lro_od_results.parquet"),
293 ExportCfg::default(),
294 )?;
295
296 // Create the ephemeris
297 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
298 let ephem_start = ephem.start_epoch().unwrap();
299 let ephem_end = ephem.end_epoch().unwrap();
300 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
301 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
302 ephem
303 .covar_at(
304 epoch,
305 anise::ephemerides::ephemeris::LocalFrame::RIC,
306 &almanac,
307 )
308 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
309 }
310 // Export as BSP!
311 ephem
312 .write_spice_bsp(
313 -85,
314 output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
315 None,
316 )
317 .expect("could not built BSP");
318 let new_almanac = Almanac::default()
319 .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
320 .unwrap();
321 new_almanac.describe(None, None, None, None, None, None, None, None);
322 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
323
324 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
325 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
326
327 // In our case, we have the truth trajectory from NASA.
328 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
329 // Export the OD trajectory first.
330 let od_trajectory = od_sol.to_traj()?;
331 // Build the RIC difference.
332 od_trajectory.ric_diff_to_parquet(
333 &traj_as_flown,
334 output_folder.join("04_lro_od_truth_error.parquet"),
335 ExportCfg::default(),
336 )?;
337
338 Ok(())
339}pub const fn from_orient_ssb(orientation_id: i32) -> Frame
pub const fn new_inertial(ephemeris_id: i32, 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 fn py_new(
ephemeris_id: i32,
orientation_id: i32,
mu_km3_s2: Option<f64>,
shape: Option<Ellipsoid>,
) -> Frame
pub fn py_new( ephemeris_id: i32, orientation_id: i32, mu_km3_s2: Option<f64>, shape: Option<Ellipsoid>, ) -> Frame
Initializes a new Frame provided its ephemeris and orientation identifiers, and optionally its gravitational parameter (in km^3/s^2) and optionally its shape (cf. [Ellipsoid]).
§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 = [
46 env!("CARGO_MANIFEST_DIR"),
47 "examples",
48 "06_lunar_orbit_determination",
49 ]
50 .iter()
51 .collect();
52
53 let meta = data_folder.join("metaalmanac.dhall");
54
55 // Load this ephem in the general Almanac we're using for this analysis.
56 let almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
57 .map_err(Box::new)?
58 .process(true)
59 .map_err(Box::new)?;
60
61 // Lock the almanac (an Arc is a read only structure).
62 let almanac = Arc::new(almanac);
63
64 // Build a nominal trajectory
65 // TODO: Switch this to a sequence once the OD over a spacecraft sequence is implemented.
66
67 let epoch = Epoch::from_gregorian_utc_at_noon(2024, 2, 29);
68 let moon_j2000 = almanac.frame_info(MOON_J2000)?;
69
70 // To build the trajectory we need to provide a spacecraft template.
71 let orbiter = Spacecraft::builder()
72 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0))
73 .srp(SRPData {
74 area_m2: 3.9 * 2.7,
75 coeff_reflectivity: 0.96,
76 })
77 .orbit(Orbit::try_keplerian_altitude(
78 150.0, 0.00212, 33.6, 45.0, 45.0, 0.0, epoch, moon_j2000,
79 )?) // Setting a zero orbit here because it's just a template
80 .build();
81
82 // ========================== //
83 // === BUILD NOMINAL TRAJ === //
84 // ========================== //
85
86 // Set up the spacecraft dynamics.
87
88 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
89 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
90 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
91
92 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
93 // We're using the GRAIL JGGRX model.
94 let mut jggrx_meta = MetaFile {
95 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
96 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
97 };
98 // And let's download it if we don't have it yet.
99 jggrx_meta.process(true)?;
100
101 // Build the spherical harmonics.
102 // The harmonics must be computed in the body fixed frame.
103 // We're using the long term prediction of the Moon principal axes frame.
104 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
105 let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
106 &jggrx_meta.uri,
107 80,
108 80,
109 almanac.frame_info(moon_pa_frame)?,
110 )?);
111
112 // Include the spherical harmonics into the orbital dynamics.
113 orbital_dyn.accel_models.push(sph_harmonics);
114
115 // We define the solar radiation pressure, using the default solar flux and accounting only
116 // for the eclipsing caused by the Earth and Moon.
117 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
118 let srp_dyn = SolarPressure::new(vec![MOON_J2000], &almanac)?;
119
120 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
121 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
122 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
123
124 println!("{dynamics}");
125
126 let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::default());
127
128 let truth_traj = setup
129 .with(orbiter, almanac.clone())
130 .for_duration_with_traj(Unit::Day * 2)?
131 .1;
132
133 // ==================== //
134 // === OD SIMULATOR === //
135 // ==================== //
136
137 // Load the Deep Space Network ground stations.
138 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
139 let ground_station_file = data_folder.join("dsn-network.yaml");
140 let devices = GroundStation::load_named(ground_station_file)?;
141
142 let proc_devices = devices.clone();
143
144 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
145 // Nyx can build a tracking schedule for you based on the first station with access.
146 let configs: BTreeMap<String, TrkConfig> =
147 TrkConfig::load_named(data_folder.join("tracking-cfg.yaml"))?;
148
149 // Build the tracking arc simulation to generate a "standard measurement".
150 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
151 devices.clone(),
152 truth_traj.clone(),
153 configs,
154 123, // Set a seed for reproducibility
155 )?;
156
157 trk.build_schedule(&almanac)?;
158 let arc = trk.generate_measurements(&almanac)?;
159 // Save the simulated tracking data
160 arc.to_parquet_simple("./data/04_output/06_lunar_simulated_tracking.parquet")?;
161
162 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
163 println!("{arc}");
164
165 // Now that we have simulated measurements, we'll run the orbit determination.
166
167 // ===================== //
168 // === OD ESTIMATION === //
169 // ===================== //
170
171 let sc = SpacecraftUncertainty::builder()
172 .nominal(orbiter)
173 .frame(LocalFrame::RIC)
174 .x_km(0.5)
175 .y_km(0.5)
176 .z_km(0.5)
177 .vx_km_s(5e-3)
178 .vy_km_s(5e-3)
179 .vz_km_s(5e-3)
180 .build();
181
182 // Build the filter initial estimate, which we will reuse in the filter.
183 let initial_estimate = sc.to_estimate()?;
184
185 println!("== FILTER STATE ==\n{orbiter:x}\n{initial_estimate}");
186
187 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
188 let process_noise = ProcessNoise3D::from_velocity_km_s(
189 &[1e-14, 1e-14, 1e-14],
190 1 * Unit::Hour,
191 10 * Unit::Minute,
192 None,
193 );
194
195 println!("{process_noise}");
196
197 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
198 let odp = SpacecraftKalmanScalarOD::new(
199 setup,
200 KalmanVariant::ReferenceUpdate,
201 Some(SigmaRejection::default()),
202 proc_devices,
203 almanac.clone(),
204 )
205 .with_process_noise(process_noise);
206
207 let od_sol = odp.process_arc(initial_estimate, &arc)?;
208
209 let final_est = od_sol.estimates.last().unwrap();
210
211 println!("{final_est}");
212
213 let ric_err = truth_traj
214 .at(final_est.epoch())?
215 .orbit
216 .ric_difference(&final_est.orbital_state())?;
217 println!("== RIC at end ==");
218 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
219 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
220
221 println!(
222 "Num residuals rejected: #{}",
223 od_sol.rejected_residuals().len()
224 );
225 println!(
226 "Percentage within +/-3: {}",
227 od_sol.residual_ratio_within_threshold(3.0).unwrap()
228 );
229 println!("Whitened residuals normal? {}", od_sol.is_normal(None)?);
230 println!("NIS consistency: {}", od_sol.nis_consistency(None)?);
231
232 od_sol.to_parquet(
233 "./data/04_output/06_lunar_od_results.parquet",
234 ExportCfg::default(),
235 )?;
236
237 let od_trajectory = od_sol.to_traj()?;
238 // Build the RIC difference.
239 od_trajectory.ric_diff_to_parquet(
240 &truth_traj,
241 "./data/04_output/06_lunar_od_truth_error.parquet",
242 ExportCfg::default(),
243 )?;
244
245 Ok(())
246}More examples
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 // Dynamics models require planetary constants and ephemerides to be defined.
42 // Let's start by grabbing those by using ANISE's MetaAlmanac.
43
44 let output_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "../data", "04_output"]
45 .iter()
46 .collect();
47
48 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
49 .iter()
50 .collect();
51
52 let meta = data_folder.join("lro-dynamics.dhall");
53
54 // Load this ephem in the general Almanac we're using for this analysis.
55 let almanac = Arc::new(
56 MetaAlmanac::new(meta.to_string_lossy().as_ref())
57 .map_err(Box::new)?
58 .process(true)
59 .map_err(Box::new)?,
60 );
61
62 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
63 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
64 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
65 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
66 let lro_frame = Frame::from_ephem_j2000(-85);
67
68 // To build the trajectory we need to provide a spacecraft template.
69 let sc_template = Spacecraft::builder()
70 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
71 .srp(SRPData {
72 // SRP configuration is arbitrary, but we will be estimating it anyway.
73 area_m2: 3.9 * 2.7,
74 coeff_reflectivity: 0.96,
75 })
76 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
77 .build();
78 // Now we can build the trajectory from the BSP file.
79 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
80 let traj_as_flown = Traj::from_bsp(
81 lro_frame,
82 MOON_J2000,
83 &almanac,
84 sc_template,
85 5.seconds(),
86 Some(Epoch::from_str("2024-01-01 01:00:00 UTC")?),
87 Some(Epoch::from_str("2024-01-02 01:00:00 UTC")?),
88 None,
89 Some("LRO".to_string()),
90 )?;
91
92 println!("{traj_as_flown}");
93
94 // ====================== //
95 // === MODEL MATCHING === //
96 // ====================== //
97
98 // Set up the spacecraft dynamics.
99
100 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
101 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
102 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
103
104 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
105 // We're using the GRAIL JGGRX model.
106 let mut jggrx_meta = MetaFile {
107 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
108 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
109 };
110 // And let's download it if we don't have it yet.
111 jggrx_meta.process(true)?;
112
113 // Build the spherical harmonics.
114 // The harmonics must be computed in the body fixed frame.
115 // We're using the long term prediction of the Moon principal axes frame.
116 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
117 let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
118 &jggrx_meta.uri,
119 81,
120 81,
121 almanac.frame_info(moon_pa_frame)?,
122 )?);
123
124 // Include the spherical harmonics into the orbital dynamics.
125 orbital_dyn.accel_models.push(sph_harmonics);
126
127 // We define the solar radiation pressure, using the default solar flux and accounting only
128 // for the eclipsing caused by the Earth and Moon.
129 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
130 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;
131
132 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
133 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
134 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
135
136 println!("{dynamics}");
137
138 // Now we can build the propagator.
139 let setup = Propagator::default_dp78(dynamics.clone());
140
141 // For reference, let's build the trajectory with Nyx's models from that LRO state.
142 let (sim_final, traj_as_sim) = setup
143 .with(*traj_as_flown.first(), almanac.clone())
144 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
145
146 println!("SIM INIT: {:x}", traj_as_flown.first());
147 println!("SIM FINAL: {sim_final:x}");
148 // Compute RIC difference between SIM and LRO ephem
149 let sim_lro_delta = sim_final
150 .orbit
151 .ric_difference(&traj_as_flown.last().orbit)?;
152 println!("{traj_as_sim}");
153 println!(
154 "SIM v LRO - RIC Position (m): {:.3}",
155 sim_lro_delta.radius_km * 1e3
156 );
157 println!(
158 "SIM v LRO - RIC Velocity (m/s): {:.3}",
159 sim_lro_delta.velocity_km_s * 1e3
160 );
161
162 traj_as_sim.ric_diff_to_parquet(
163 &traj_as_flown,
164 output_folder.join("./04_lro_sim_truth_error.parquet"),
165 ExportCfg::default(),
166 )?;
167
168 // ==================== //
169 // === OD SIMULATOR === //
170 // ==================== //
171
172 // 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
173 // and the truth LRO state.
174
175 // Therefore, we will actually run an estimation from a dispersed LRO state.
176 // The sc_seed is the true LRO state from the BSP.
177 let sc_seed = *traj_as_flown.first();
178
179 // Load the Deep Space Network ground stations.
180 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
181 let ground_station_file: PathBuf = [
182 env!("CARGO_MANIFEST_DIR"),
183 "examples",
184 "04_lro_od",
185 "dsn-network.yaml",
186 ]
187 .iter()
188 .collect();
189
190 let devices = GroundStation::load_named(ground_station_file)?;
191
192 let proc_devices = devices.clone();
193
194 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
195 // Nyx can build a tracking schedule for you based on the first station with access.
196 let trkconfg_yaml: PathBuf = [
197 env!("CARGO_MANIFEST_DIR"),
198 "examples",
199 "04_lro_od",
200 "tracking-cfg.yaml",
201 ]
202 .iter()
203 .collect();
204
205 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
206
207 // Build the tracking arc simulation to generate a "standard measurement".
208 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
209 devices.clone(),
210 traj_as_flown.clone(),
211 configs,
212 123, // Set a seed for reproducibility
213 )?;
214
215 trk.build_schedule(&almanac)?;
216 let arc = trk.generate_measurements(&almanac)?;
217 // Save the simulated tracking data
218 arc.to_parquet_simple(output_folder.join("04_lro_simulated_tracking.parquet"))?;
219
220 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
221 println!("{arc}");
222
223 // Now that we have simulated measurements, we'll run the orbit determination.
224
225 // ===================== //
226 // === OD ESTIMATION === //
227 // ===================== //
228
229 let sc = SpacecraftUncertainty::builder()
230 .nominal(sc_seed)
231 .frame(LocalFrame::RIC)
232 .x_km(0.5)
233 .y_km(0.5)
234 .z_km(0.5)
235 .vx_km_s(5e-3)
236 .vy_km_s(5e-3)
237 .vz_km_s(5e-3)
238 .build();
239
240 // Build the filter initial estimate, which we will reuse in the filter.
241 let mut initial_estimate = sc.to_estimate()?;
242 initial_estimate.covar *= 1.5;
243
244 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
245
246 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
247 let process_noise = ProcessNoise3D::from_velocity_km_s(
248 &[5e-13, 5e-13, 5e-13],
249 1 * Unit::Hour,
250 10 * Unit::Minute,
251 None,
252 );
253
254 println!("{process_noise}");
255
256 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
257 let odp = SpacecraftKalmanOD::new(
258 setup,
259 KalmanVariant::ReferenceUpdate,
260 Some(SigmaRejection::default()),
261 proc_devices,
262 almanac.clone(),
263 )
264 .with_process_noise(process_noise);
265
266 let od_sol = odp.process_arc(initial_estimate, &arc)?;
267
268 let final_est = od_sol.estimates.last().unwrap();
269
270 println!("{final_est}");
271
272 let ric_err = traj_as_flown
273 .at(final_est.epoch())?
274 .orbit
275 .ric_difference(&final_est.orbital_state())?;
276 println!("== RIC at end ==");
277 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
278 println!("RIC Velocity (m/s): {:.3}", 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 od_sol.nis_consistency(None)?.log();
290
291 od_sol.to_parquet(
292 output_folder.join("04_lro_od_results.parquet"),
293 ExportCfg::default(),
294 )?;
295
296 // Create the ephemeris
297 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
298 let ephem_start = ephem.start_epoch().unwrap();
299 let ephem_end = ephem.end_epoch().unwrap();
300 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
301 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
302 ephem
303 .covar_at(
304 epoch,
305 anise::ephemerides::ephemeris::LocalFrame::RIC,
306 &almanac,
307 )
308 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
309 }
310 // Export as BSP!
311 ephem
312 .write_spice_bsp(
313 -85,
314 output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
315 None,
316 )
317 .expect("could not built BSP");
318 let new_almanac = Almanac::default()
319 .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
320 .unwrap();
321 new_almanac.describe(None, None, None, None, None, None, None, None);
322 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
323
324 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
325 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
326
327 // In our case, we have the truth trajectory from NASA.
328 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
329 // Export the OD trajectory first.
330 let od_trajectory = od_sol.to_traj()?;
331 // Build the RIC difference.
332 od_trajectory.ric_diff_to_parquet(
333 &traj_as_flown,
334 output_folder.join("04_lro_od_truth_error.parquet"),
335 ExportCfg::default(),
336 )?;
337
338 Ok(())
339}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
pub fn is_dynamic(&self) -> bool
pub fn is_dynamic(&self) -> bool
Returns true if this is a dynamic frame, e.g. Mean/True of Date/Epoch
:rtype: bool
Trait Implementations§
impl Copy for Frame
impl DerefToPyAny for Frame
§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 Encode for Frame
impl Encode for Frame
§fn encoded_len(&self) -> Result<Length, Error>
fn encoded_len(&self) -> Result<Length, Error>
§fn encode(&self, encoder: &mut impl Writer) -> Result<(), Error>
fn encode(&self, encoder: &mut impl Writer) -> Result<(), Error>
Writer].§fn encode_to_slice<'a>(&self, buf: &'a mut [u8]) -> Result<&'a [u8], Error>
fn encode_to_slice<'a>(&self, buf: &'a mut [u8]) -> Result<&'a [u8], Error>
§fn encode_to_vec(&self, buf: &mut Vec<u8>) -> Result<Length, Error>
fn encode_to_vec(&self, buf: &mut Vec<u8>) -> Result<Length, Error>
§impl<'py> IntoPyObject<'py> for Frame
impl<'py> IntoPyObject<'py> for Frame
§impl PyTypeInfo for Frame
impl PyTypeInfo for Frame
§const NAME: &'static str = <Self as ::pyo3::PyClass>::NAME
const NAME: &'static str = <Self as ::pyo3::PyClass>::NAME
prefer using ::type_object(py).name() to get the correct runtime value
§const MODULE: Option<&'static str> = <Self as ::pyo3::impl_::pyclass::PyClassImpl>::MODULE
const MODULE: Option<&'static str> = <Self as ::pyo3::impl_::pyclass::PyClassImpl>::MODULE
prefer using ::type_object(py).module() to get the correct runtime value
§fn type_object_raw(py: Python<'_>) -> *mut PyTypeObject
fn type_object_raw(py: Python<'_>) -> *mut PyTypeObject
§fn type_object(py: Python<'_>) -> Bound<'_, PyType>
fn type_object(py: Python<'_>) -> Bound<'_, PyType>
§fn is_type_of(object: &Bound<'_, PyAny>) -> bool
fn is_type_of(object: &Bound<'_, PyAny>) -> bool
object is an instance of this type or a subclass of this type.§fn is_exact_type_of(object: &Bound<'_, PyAny>) -> bool
fn is_exact_type_of(object: &Bound<'_, PyAny>) -> bool
object is an instance of this type.§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
Available on crate feature metaload only.
impl StaticType for Frame
metaload only.§fn static_type() -> SimpleType
fn static_type() -> SimpleType
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 UnsafeUnpin for Frame
impl UnwindSafe for Frame
Blanket Implementations§
impl<T> Allocation for T
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
impl<ST, DT> CastableFrom<ST, Initialized, Initialized> for DT
impl<ST, DT> CastableFrom<ST, Uninit, Uninit> for DT
Source§impl<T> CloneToUninit for Twhere
T: Clone,
impl<T> CloneToUninit for Twhere
T: Clone,
impl<T> DecodeOwned for Twhere
T: for<'a> Decode<'a>,
impl<T> DeserializeOwned for Twhere
T: for<'de> Deserialize<'de>,
Source§impl<T> FromDhall for Twhere
T: DeserializeOwned,
impl<T> FromDhall for Twhere
T: DeserializeOwned,
fn from_dhall(v: &Value) -> Result<T, Error>
impl<'py, T> FromPyObjectOwned<'py> for Twhere
T: for<'a> FromPyObject<'a, 'py>,
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<'py, T> IntoPyObjectExt<'py> for Twhere
T: IntoPyObject<'py>,
impl<'py, T> IntoPyObjectExt<'py> for Twhere
T: IntoPyObject<'py>,
§fn into_bound_py_any(self, py: Python<'py>) -> Result<Bound<'py, PyAny>, PyErr>
fn into_bound_py_any(self, py: Python<'py>) -> Result<Bound<'py, PyAny>, PyErr>
self into an owned Python object, dropping type information.§fn into_py_any(self, py: Python<'py>) -> Result<Py<PyAny>, PyErr>
fn into_py_any(self, py: Python<'py>) -> Result<Py<PyAny>, PyErr>
self into an owned Python object, dropping type information and unbinding it
from the 'py lifetime.§fn into_pyobject_or_pyerr(self, py: Python<'py>) -> Result<Self::Output, PyErr>
fn into_pyobject_or_pyerr(self, py: Python<'py>) -> Result<Self::Output, PyErr>
self into a Python object. Read more§impl<T> Pointable for T
impl<T> Pointable for T
§impl<T> PyErrArguments for T
impl<T> PyErrArguments for T
§impl<T> PyTypeCheck for Twhere
T: PyTypeInfo,
impl<T> PyTypeCheck for Twhere
T: PyTypeInfo,
§fn type_check(object: &Bound<'_, PyAny>) -> bool
fn type_check(object: &Bound<'_, PyAny>) -> bool
§fn classinfo_object(py: Python<'_>) -> Bound<'_, PyAny>
fn classinfo_object(py: Python<'_>) -> Bound<'_, PyAny>
isinstance and issubclass function. Read moreimpl<T> Read<Exclusive, BecauseExclusive> for Twhere
T: ?Sized,
impl<T> Scalar 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.