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