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 247-252)
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}