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 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>§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 189-194)
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 true,
110 almanac.frame_info(moon_pa_frame)?,
111 )?);
112
113 // Include the spherical harmonics into the orbital dynamics.
114 orbital_dyn.accel_models.push(sph_harmonics);
115
116 // We define the solar radiation pressure, using the default solar flux and accounting only
117 // for the eclipsing caused by the Earth and Moon.
118 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
119 let srp_dyn = SolarPressure::new(vec![MOON_J2000], &almanac)?;
120
121 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
122 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
123 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
124
125 println!("{dynamics}");
126
127 let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::default());
128
129 let truth_traj = setup
130 .with(orbiter, almanac.clone())
131 .for_duration_with_traj(Unit::Day * 2)?
132 .1;
133
134 // ==================== //
135 // === OD SIMULATOR === //
136 // ==================== //
137
138 // Load the Deep Space Network ground stations.
139 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
140 let ground_station_file = data_folder.join("dsn-network.yaml");
141 let devices = GroundStation::load_named(ground_station_file)?;
142
143 let proc_devices = devices.clone();
144
145 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
146 // Nyx can build a tracking schedule for you based on the first station with access.
147 let configs: BTreeMap<String, TrkConfig> =
148 TrkConfig::load_named(data_folder.join("tracking-cfg.yaml"))?;
149
150 // Build the tracking arc simulation to generate a "standard measurement".
151 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
152 devices.clone(),
153 truth_traj.clone(),
154 configs,
155 123, // Set a seed for reproducibility
156 )?;
157
158 trk.build_schedule(almanac.clone())?;
159 let arc = trk.generate_measurements(almanac.clone())?;
160 // Save the simulated tracking data
161 arc.to_parquet_simple("./data/04_output/06_lunar_simulated_tracking.parquet")?;
162
163 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
164 println!("{arc}");
165
166 // Now that we have simulated measurements, we'll run the orbit determination.
167
168 // ===================== //
169 // === OD ESTIMATION === //
170 // ===================== //
171
172 let sc = SpacecraftUncertainty::builder()
173 .nominal(orbiter)
174 .frame(LocalFrame::RIC)
175 .x_km(0.5)
176 .y_km(0.5)
177 .z_km(0.5)
178 .vx_km_s(5e-3)
179 .vy_km_s(5e-3)
180 .vz_km_s(5e-3)
181 .build();
182
183 // Build the filter initial estimate, which we will reuse in the filter.
184 let initial_estimate = sc.to_estimate()?;
185
186 println!("== FILTER STATE ==\n{orbiter:x}\n{initial_estimate}");
187
188 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
189 let process_noise = ProcessNoise3D::from_velocity_km_s(
190 &[1e-14, 1e-14, 1e-14],
191 1 * Unit::Hour,
192 10 * Unit::Minute,
193 None,
194 );
195
196 println!("{process_noise}");
197
198 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
199 let odp = SpacecraftKalmanScalarOD::new(
200 setup,
201 KalmanVariant::ReferenceUpdate,
202 Some(ResidRejectCrit::default()),
203 proc_devices,
204 almanac.clone(),
205 )
206 .with_process_noise(process_noise);
207
208 let od_sol = odp.process_arc(initial_estimate, &arc)?;
209
210 let final_est = od_sol.estimates.last().unwrap();
211
212 println!("{final_est}");
213
214 let ric_err = truth_traj
215 .at(final_est.epoch())?
216 .orbit
217 .ric_difference(&final_est.orbital_state())?;
218 println!("== RIC at end ==");
219 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
220 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
221
222 println!(
223 "Num residuals rejected: #{}",
224 od_sol.rejected_residuals().len()
225 );
226 println!(
227 "Percentage within +/-3: {}",
228 od_sol.residual_ratio_within_threshold(3.0).unwrap()
229 );
230 println!("Whitened residuals normal? {}", od_sol.is_normal(None)?);
231 println!("NIS test success? {}", od_sol.is_nis_consistent(None)?);
232
233 od_sol.to_parquet(
234 "./data/04_output/06_lunar_od_results.parquet",
235 ExportCfg::default(),
236 )?;
237
238 let od_trajectory = od_sol.to_traj()?;
239 // Build the RIC difference.
240 od_trajectory.ric_diff_to_parquet(
241 &truth_traj,
242 "./data/04_output/06_lunar_od_truth_error.parquet",
243 ExportCfg::default(),
244 )?;
245
246 Ok(())
247}More examples
nyx-core/examples/04_lro_od/main.rs (lines 279-284)
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.clone(),
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 true,
141 almanac.frame_info(moon_pa_frame)?,
142 )?);
143
144 // Include the spherical harmonics into the orbital dynamics.
145 orbital_dyn.accel_models.push(sph_harmonics);
146
147 // We define the solar radiation pressure, using the default solar flux and accounting only
148 // for the eclipsing caused by the Earth and Moon.
149 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
150 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;
151
152 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
153 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
154 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
155
156 println!("{dynamics}");
157
158 // Now we can build the propagator.
159 let setup = Propagator::default_dp78(dynamics.clone());
160
161 // For reference, let's build the trajectory with Nyx's models from that LRO state.
162 let (sim_final, traj_as_sim) = setup
163 .with(*traj_as_flown.first(), almanac.clone())
164 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
165
166 println!("SIM INIT: {:x}", traj_as_flown.first());
167 println!("SIM FINAL: {sim_final:x}");
168 // Compute RIC difference between SIM and LRO ephem
169 let sim_lro_delta = sim_final
170 .orbit
171 .ric_difference(&traj_as_flown.last().orbit)?;
172 println!("{traj_as_sim}");
173 println!(
174 "SIM v LRO - RIC Position (m): {:.3}",
175 sim_lro_delta.radius_km * 1e3
176 );
177 println!(
178 "SIM v LRO - RIC Velocity (m/s): {:.3}",
179 sim_lro_delta.velocity_km_s * 1e3
180 );
181
182 traj_as_sim.ric_diff_to_parquet(
183 &traj_as_flown,
184 output_folder.join("./04_lro_sim_truth_error.parquet"),
185 ExportCfg::default(),
186 )?;
187
188 // ==================== //
189 // === OD SIMULATOR === //
190 // ==================== //
191
192 // 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
193 // and the truth LRO state.
194
195 // Therefore, we will actually run an estimation from a dispersed LRO state.
196 // The sc_seed is the true LRO state from the BSP.
197 let sc_seed = *traj_as_flown.first();
198
199 // Load the Deep Space Network ground stations.
200 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
201 let ground_station_file: PathBuf = [
202 env!("CARGO_MANIFEST_DIR"),
203 "examples",
204 "04_lro_od",
205 "dsn-network.yaml",
206 ]
207 .iter()
208 .collect();
209
210 let devices = GroundStation::load_named(ground_station_file)?;
211
212 let mut proc_devices = devices.clone();
213
214 // Increase the noise in the devices to accept more measurements.
215 for gs in proc_devices.values_mut() {
216 if let Some(noise) = &mut gs
217 .stochastic_noises
218 .as_mut()
219 .unwrap()
220 .get_mut(&MeasurementType::Range)
221 {
222 *noise.white_noise.as_mut().unwrap() *= 3.0;
223 }
224 }
225
226 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
227 // Nyx can build a tracking schedule for you based on the first station with access.
228 let trkconfg_yaml: PathBuf = [
229 env!("CARGO_MANIFEST_DIR"),
230 "examples",
231 "04_lro_od",
232 "tracking-cfg.yaml",
233 ]
234 .iter()
235 .collect();
236
237 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
238
239 // Build the tracking arc simulation to generate a "standard measurement".
240 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
241 devices.clone(),
242 traj_as_flown.clone(),
243 configs,
244 123, // Set a seed for reproducibility
245 )?;
246
247 trk.build_schedule(almanac.clone())?;
248 let arc = trk.generate_measurements(almanac.clone())?;
249 // Save the simulated tracking data
250 arc.to_parquet_simple(output_folder.join("04_lro_simulated_tracking.parquet"))?;
251
252 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
253 println!("{arc}");
254
255 // Now that we have simulated measurements, we'll run the orbit determination.
256
257 // ===================== //
258 // === OD ESTIMATION === //
259 // ===================== //
260
261 let sc = SpacecraftUncertainty::builder()
262 .nominal(sc_seed)
263 .frame(LocalFrame::RIC)
264 .x_km(0.5)
265 .y_km(0.5)
266 .z_km(0.5)
267 .vx_km_s(5e-3)
268 .vy_km_s(5e-3)
269 .vz_km_s(5e-3)
270 .build();
271
272 // Build the filter initial estimate, which we will reuse in the filter.
273 let mut initial_estimate = sc.to_estimate()?;
274 initial_estimate.covar *= 3.0;
275
276 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
277
278 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
279 let process_noise = ProcessNoise3D::from_velocity_km_s(
280 &[1e-12, 1e-12, 1e-12],
281 1 * Unit::Hour,
282 10 * Unit::Minute,
283 None,
284 );
285
286 println!("{process_noise}");
287
288 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
289 let odp = SpacecraftKalmanOD::new(
290 setup,
291 KalmanVariant::ReferenceUpdate,
292 Some(ResidRejectCrit::default()),
293 proc_devices,
294 almanac.clone(),
295 )
296 .with_process_noise(process_noise);
297
298 let od_sol = odp.process_arc(initial_estimate, &arc)?;
299
300 let final_est = od_sol.estimates.last().unwrap();
301
302 println!("{final_est}");
303
304 let ric_err = traj_as_flown
305 .at(final_est.epoch())?
306 .orbit
307 .ric_difference(&final_est.orbital_state())?;
308 println!("== RIC at end ==");
309 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
310 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
311
312 println!(
313 "Num residuals rejected: #{}",
314 od_sol.rejected_residuals().len()
315 );
316 println!(
317 "Percentage within +/-3: {}",
318 od_sol.residual_ratio_within_threshold(3.0).unwrap()
319 );
320 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
321
322 od_sol.to_parquet(
323 output_folder.join("04_lro_od_results.parquet"),
324 ExportCfg::default(),
325 )?;
326
327 // Create the ephemeris
328 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
329 let ephem_start = ephem.start_epoch().unwrap();
330 let ephem_end = ephem.end_epoch().unwrap();
331 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
332 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
333 ephem
334 .covar_at(
335 epoch,
336 anise::ephemerides::ephemeris::LocalFrame::RIC,
337 &almanac,
338 )
339 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
340 }
341 // Export as BSP!
342 ephem
343 .write_spice_bsp(
344 -85,
345 output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
346 None,
347 )
348 .expect("could not built BSP");
349 let new_almanac = Almanac::default()
350 .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
351 .unwrap();
352 new_almanac.describe(None, None, None, None, None, None, None, None);
353 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
354
355 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
356 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
357
358 // In our case, we have the truth trajectory from NASA.
359 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
360 // Export the OD trajectory first.
361 let od_trajectory = od_sol.to_traj()?;
362 // Build the RIC difference.
363 od_trajectory.ric_diff_to_parquet(
364 &traj_as_flown,
365 output_folder.join("04_lro_od_truth_error.parquet"),
366 ExportCfg::default(),
367 )?;
368
369 Ok(())
370}