pub struct SpacecraftUncertainty {
pub nominal: Spacecraft,
pub frame: Option<LocalFrame>,
pub x_km: f64,
pub y_km: f64,
pub z_km: f64,
pub vx_km_s: f64,
pub vy_km_s: f64,
pub vz_km_s: f64,
pub coeff_reflectivity: f64,
pub coeff_drag: f64,
pub mass_kg: f64,
}
Expand description
Builds a spacecraft uncertainty in different local frames, dispersing any of the parameters of the spacecraft state.
§Usage
Use the TypeBuilder
trait, e.g SpacecraftUncertainty::builder().nominal(spacecraft).frame(LocalFrame::RIC).x_km(0.5).y_km(0.5).z_km(0.5).build()
to build an uncertainty on position in the RIC frame of 500 meters on R, I, and C, and zero on all other parameters (velocity components, Cr, Cd, mass).
Fields§
§nominal: Spacecraft
§frame: Option<LocalFrame>
§x_km: f64
§y_km: f64
§z_km: f64
§vx_km_s: f64
§vy_km_s: f64
§vz_km_s: f64
§coeff_reflectivity: f64
§coeff_drag: f64
§mass_kg: f64
Implementations§
Source§impl SpacecraftUncertainty
impl SpacecraftUncertainty
Sourcepub fn builder() -> SpacecraftUncertaintyBuilder<((), (), (), (), (), (), (), (), (), (), ())>
pub fn builder() -> SpacecraftUncertaintyBuilder<((), (), (), (), (), (), (), (), (), (), ())>
Create a builder for building SpacecraftUncertainty
.
On the builder, call .nominal(...)
, .frame(...)
(optional), .x_km(...)
(optional), .y_km(...)
(optional), .z_km(...)
(optional), .vx_km_s(...)
(optional), .vy_km_s(...)
(optional), .vz_km_s(...)
(optional), .coeff_reflectivity(...)
(optional), .coeff_drag(...)
(optional), .mass_kg(...)
(optional) to set the values of the fields.
Finally, call .build()
to create the instance of SpacecraftUncertainty
.
Examples found in repository?
26fn main() -> Result<(), Box<dyn Error>> {
27 pel::init();
28 // Dynamics models require planetary constants and ephemerides to be defined.
29 // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30 // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
31
32 // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
33 // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
34 let mut latest_jwst_ephem = MetaFile {
35 uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
36 crc32: None,
37 };
38 latest_jwst_ephem.process(true)?;
39
40 // Load this ephem in the general Almanac we're using for this analysis.
41 let almanac = Arc::new(
42 MetaAlmanac::latest()
43 .map_err(Box::new)?
44 .load_from_metafile(latest_jwst_ephem, true)?,
45 );
46
47 // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
48 // in the BSP. We need this ID in order to query the ephemeris.
49 const JWST_NAIF_ID: i32 = -170;
50 // Let's build a frame in the J2000 orientation centered on the JWST.
51 const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
52
53 // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
54 let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
55 println!("JWST defined from {earliest_epoch} to {latest_epoch}");
56 // Fetch the state, printing it in the Earth J2000 frame.
57 let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
58 println!("{jwst_orbit:x}");
59
60 // Build the spacecraft
61 // SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
62 // SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
63 let jwst = Spacecraft::builder()
64 .orbit(jwst_orbit)
65 .srp(SRPData {
66 area_m2: 21.197 * 14.162,
67 coeff_reflectivity: 1.56,
68 })
69 .mass(Mass::from_dry_mass(6200.0))
70 .build();
71
72 // Build up the spacecraft uncertainty builder.
73 // We can use the spacecraft uncertainty structure to build this up.
74 // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
75 // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
76 // Nyx can also estimate the deviation of the spacecraft parameters.
77 let jwst_uncertainty = SpacecraftUncertainty::builder()
78 .nominal(jwst)
79 .frame(LocalFrame::RIC)
80 .x_km(0.5)
81 .y_km(0.3)
82 .z_km(1.5)
83 .vx_km_s(1e-4)
84 .vy_km_s(0.6e-3)
85 .vz_km_s(3e-3)
86 .build();
87
88 println!("{jwst_uncertainty}");
89
90 // Build the Kalman filter estimate.
91 // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
92 // but this approach requires quite a bit more boilerplate code.
93 let jwst_estimate = jwst_uncertainty.to_estimate()?;
94
95 // Set up the spacecraft dynamics.
96 // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
97 // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
98
99 let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
100 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102 // Finalize setting up the dynamics.
103 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
104
105 // Build the propagator set up to use for the whole analysis.
106 let setup = Propagator::default(dynamics);
107
108 // All of the analysis will use this duration.
109 let prediction_duration = 6.5 * Unit::Day;
110
111 // === Covariance mapping ===
112 // For the covariance mapping / prediction, we'll use the common orbit determination approach.
113 // This is done by setting up a spacecraft Kalman filter OD process, and predicting for the analysis duration.
114
115 // Build the propagation instance for the OD process.
116 let odp = SpacecraftKalmanOD::new(
117 setup.clone(),
118 KalmanVariant::DeviationTracking,
119 None,
120 BTreeMap::new(),
121 almanac.clone(),
122 );
123
124 // The prediction step is 1 minute by default, configured in the OD process, i.e. how often we want to know the covariance.
125 assert_eq!(odp.max_step, 1_i64.minutes());
126 // Finally, predict, and export the trajectory with covariance to a parquet file.
127 let od_sol = odp.predict_for(jwst_estimate, prediction_duration)?;
128 od_sol.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
129
130 // === Monte Carlo framework ===
131 // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
132
133 let my_mc = MonteCarlo::new(
134 jwst, // Nominal state
135 jwst_estimate.to_random_variable()?,
136 "02_jwst".to_string(), // Scenario name
137 None, // No specific seed specified, so one will be drawn from the computer's entropy.
138 );
139
140 let num_runs = 5_000;
141 let rslts = my_mc.run_until_epoch(
142 setup,
143 almanac.clone(),
144 jwst.epoch() + prediction_duration,
145 num_runs,
146 );
147
148 assert_eq!(rslts.runs.len(), num_runs);
149 // Finally, export these results, computing the eclipse percentage for all of these results.
150
151 // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
152 let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
153 let umbra_event = eclipse_loc.to_umbra_event();
154 let penumbra_event = eclipse_loc.to_penumbra_event();
155
156 rslts.to_parquet(
157 "02_jwst_monte_carlo.parquet",
158 Some(vec![&umbra_event, &penumbra_event]),
159 ExportCfg::default(),
160 almanac,
161 )?;
162
163 Ok(())
164}
More examples
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 let manifest_dir =
42 PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
43
44 let out = manifest_dir.join("data/04_output/");
45
46 let almanac = Arc::new(
47 Almanac::new(
48 &manifest_dir
49 .join("data/01_planetary/pck08.pca")
50 .to_string_lossy(),
51 )
52 .unwrap()
53 .load(
54 &manifest_dir
55 .join("data/01_planetary/de440s.bsp")
56 .to_string_lossy(),
57 )
58 .unwrap(),
59 );
60
61 let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
62 let moon_iau = almanac.frame_from_uid(IAU_MOON_FRAME).unwrap();
63
64 let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
65 let nrho = Orbit::cartesian(
66 166_473.631_302_239_7,
67 -274_715.487_253_382_7,
68 -211_233.210_176_686_7,
69 0.933_451_604_520_018_4,
70 0.436_775_046_841_900_9,
71 -0.082_211_021_250_348_95,
72 epoch,
73 eme2k,
74 );
75
76 let tx_nrho_sc = Spacecraft::from(nrho);
77
78 let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
79 println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
80
81 let bodies = vec![EARTH, SUN];
82 let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
83
84 let setup = Propagator::rk89(
85 dynamics,
86 IntegratorOptions::builder().max_step(0.5.minutes()).build(),
87 );
88
89 /* == Propagate the NRHO vehicle == */
90 let prop_time = 1.1 * state_luna.period().unwrap();
91
92 let (nrho_final, mut tx_traj) = setup
93 .with(tx_nrho_sc, almanac.clone())
94 .for_duration_with_traj(prop_time)
95 .unwrap();
96
97 tx_traj.name = Some("NRHO Tx SC".to_string());
98
99 println!("{tx_traj}");
100
101 /* == Propagate an LLO vehicle == */
102 let llo_orbit =
103 Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
104
105 let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
106
107 let (_, llo_traj) = setup
108 .with(llo_sc, almanac.clone())
109 .until_epoch_with_traj(nrho_final.epoch())
110 .unwrap();
111
112 // Export the subset of the first two hours.
113 llo_traj
114 .clone()
115 .filter_by_offset(..2.hours())
116 .to_parquet_simple(out.join("05_caps_llo_truth.pq"), almanac.clone())?;
117
118 /* == Setup the interlink == */
119
120 let mut measurement_types = IndexSet::new();
121 measurement_types.insert(MeasurementType::Range);
122 measurement_types.insert(MeasurementType::Doppler);
123
124 let mut stochastics = IndexMap::new();
125
126 let sa45_csac_allan_dev = 1e-11;
127
128 stochastics.insert(
129 MeasurementType::Range,
130 StochasticNoise::from_hardware_range_km(
131 sa45_csac_allan_dev,
132 10.0.seconds(),
133 link_specific::ChipRate::StandardT4B,
134 link_specific::SN0::Average,
135 ),
136 );
137
138 stochastics.insert(
139 MeasurementType::Doppler,
140 StochasticNoise::from_hardware_doppler_km_s(
141 sa45_csac_allan_dev,
142 10.0.seconds(),
143 link_specific::CarrierFreq::SBand,
144 link_specific::CN0::Average,
145 ),
146 );
147
148 let interlink = InterlinkTxSpacecraft {
149 traj: tx_traj,
150 measurement_types,
151 integration_time: None,
152 timestamp_noise_s: None,
153 ab_corr: Aberration::LT,
154 stochastic_noises: Some(stochastics),
155 };
156
157 // Devices are the transmitter, which is our NRHO vehicle.
158 let mut devices = BTreeMap::new();
159 devices.insert("NRHO Tx SC".to_string(), interlink);
160
161 let mut configs = BTreeMap::new();
162 configs.insert(
163 "NRHO Tx SC".to_string(),
164 TrkConfig::builder()
165 .strands(vec![Strand {
166 start: epoch,
167 end: nrho_final.epoch(),
168 }])
169 .build(),
170 );
171
172 let mut trk_sim =
173 TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
174 println!("{trk_sim}");
175
176 let trk_data = trk_sim.generate_measurements(almanac.clone()).unwrap();
177 println!("{trk_data}");
178
179 trk_data
180 .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
181 .unwrap();
182
183 // Run a truth OD where we estimate the LLO position
184 let llo_uncertainty = SpacecraftUncertainty::builder()
185 .nominal(llo_sc)
186 .x_km(1.0)
187 .y_km(1.0)
188 .z_km(1.0)
189 .vx_km_s(1e-3)
190 .vy_km_s(1e-3)
191 .vz_km_s(1e-3)
192 .build();
193
194 let mut proc_devices = devices.clone();
195
196 // Define the initial estimate, randomized, seed for reproducibility
197 let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
198 // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
199 initial_estimate.covar *= 2.5;
200
201 // Increase the noise in the devices to accept more measurements.
202
203 for link in proc_devices.values_mut() {
204 for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
205 *noise.white_noise.as_mut().unwrap() *= 3.0;
206 }
207 }
208
209 let init_err = initial_estimate
210 .orbital_state()
211 .ric_difference(&llo_orbit)
212 .unwrap();
213
214 println!("initial estimate:\n{initial_estimate}");
215 println!("RIC errors = {init_err}",);
216
217 let odp = InterlinkKalmanOD::new(
218 setup.clone(),
219 KalmanVariant::ReferenceUpdate,
220 Some(ResidRejectCrit::default()),
221 proc_devices,
222 almanac.clone(),
223 );
224
225 // Shrink the data to process.
226 let arc = trk_data.filter_by_offset(..2.hours());
227
228 let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
229
230 println!("{od_sol}");
231
232 od_sol
233 .to_parquet(
234 out.join(format!("05_caps_interlink_od_sol.pq")),
235 ExportCfg::default(),
236 )
237 .unwrap();
238
239 let od_traj = od_sol.to_traj().unwrap();
240
241 od_traj
242 .ric_diff_to_parquet(
243 &llo_traj,
244 out.join(format!("05_caps_interlink_llo_est_error.pq")),
245 ExportCfg::default(),
246 )
247 .unwrap();
248
249 let final_est = od_sol.estimates.last().unwrap();
250 assert!(final_est.within_3sigma(), "should be within 3 sigma");
251
252 println!("ESTIMATE\n{final_est:x}\n");
253 let truth = llo_traj.at(final_est.epoch()).unwrap();
254 println!("TRUTH\n{truth:x}");
255
256 let final_err = truth
257 .orbit
258 .ric_difference(&final_est.orbital_state())
259 .unwrap();
260 println!("ERROR {final_err}");
261
262 // Build the residuals versus reference plot.
263 let rvr_sol = odp
264 .process_arc(initial_estimate, &arc.resid_vs_ref_check())
265 .unwrap();
266
267 rvr_sol
268 .to_parquet(
269 out.join(format!("05_caps_interlink_resid_v_ref.pq")),
270 ExportCfg::default(),
271 )
272 .unwrap();
273
274 let final_rvr = rvr_sol.estimates.last().unwrap();
275
276 println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
277 println!(
278 "Pure prop error {:.3} m",
279 final_rvr
280 .orbital_state()
281 .ric_difference(&final_est.orbital_state())
282 .unwrap()
283 .rmag_km()
284 * 1e3
285 );
286
287 Ok(())
288}
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 data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
45 .iter()
46 .collect();
47
48 let meta = data_folder.join("lro-dynamics.dhall");
49
50 // Load this ephem in the general Almanac we're using for this analysis.
51 let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
52 .map_err(Box::new)?
53 .process(true)
54 .map_err(Box::new)?;
55
56 let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
57 moon_pc.mu_km3_s2 = 4902.74987;
58 almanac.planetary_data.set_by_id(MOON, moon_pc)?;
59
60 let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
61 earth_pc.mu_km3_s2 = 398600.436;
62 almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
63
64 // Save this new kernel for reuse.
65 // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
66 almanac
67 .planetary_data
68 .save_as(&data_folder.join("lro-specific.pca"), true)?;
69
70 // Lock the almanac (an Arc is a read only structure).
71 let almanac = Arc::new(almanac);
72
73 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
74 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
75 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
76 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
77 let lro_frame = Frame::from_ephem_j2000(-85);
78
79 // To build the trajectory we need to provide a spacecraft template.
80 let sc_template = Spacecraft::builder()
81 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
82 .srp(SRPData {
83 // SRP configuration is arbitrary, but we will be estimating it anyway.
84 area_m2: 3.9 * 2.7,
85 coeff_reflectivity: 0.96,
86 })
87 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
88 .build();
89 // Now we can build the trajectory from the BSP file.
90 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
91 let traj_as_flown = Traj::from_bsp(
92 lro_frame,
93 MOON_J2000,
94 almanac.clone(),
95 sc_template,
96 5.seconds(),
97 Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
98 Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
99 Aberration::LT,
100 Some("LRO".to_string()),
101 )?;
102
103 println!("{traj_as_flown}");
104
105 // ====================== //
106 // === MODEL MATCHING === //
107 // ====================== //
108
109 // Set up the spacecraft dynamics.
110
111 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
112 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
113 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
114
115 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
116 // We're using the GRAIL JGGRX model.
117 let mut jggrx_meta = MetaFile {
118 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
119 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
120 };
121 // And let's download it if we don't have it yet.
122 jggrx_meta.process(true)?;
123
124 // Build the spherical harmonics.
125 // The harmonics must be computed in the body fixed frame.
126 // We're using the long term prediction of the Moon principal axes frame.
127 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
128 let sph_harmonics = Harmonics::from_stor(
129 almanac.frame_from_uid(moon_pa_frame)?,
130 HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
131 );
132
133 // Include the spherical harmonics into the orbital dynamics.
134 orbital_dyn.accel_models.push(sph_harmonics);
135
136 // We define the solar radiation pressure, using the default solar flux and accounting only
137 // for the eclipsing caused by the Earth and Moon.
138 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
139 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
140
141 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
142 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
143 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
144
145 println!("{dynamics}");
146
147 // Now we can build the propagator.
148 let setup = Propagator::default_dp78(dynamics.clone());
149
150 // For reference, let's build the trajectory with Nyx's models from that LRO state.
151 let (sim_final, traj_as_sim) = setup
152 .with(*traj_as_flown.first(), almanac.clone())
153 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
154
155 println!("SIM INIT: {:x}", traj_as_flown.first());
156 println!("SIM FINAL: {sim_final:x}");
157 // Compute RIC difference between SIM and LRO ephem
158 let sim_lro_delta = sim_final
159 .orbit
160 .ric_difference(&traj_as_flown.last().orbit)?;
161 println!("{traj_as_sim}");
162 println!(
163 "SIM v LRO - RIC Position (m): {:.3}",
164 sim_lro_delta.radius_km * 1e3
165 );
166 println!(
167 "SIM v LRO - RIC Velocity (m/s): {:.3}",
168 sim_lro_delta.velocity_km_s * 1e3
169 );
170
171 traj_as_sim.ric_diff_to_parquet(
172 &traj_as_flown,
173 "./04_lro_sim_truth_error.parquet",
174 ExportCfg::default(),
175 )?;
176
177 // ==================== //
178 // === OD SIMULATOR === //
179 // ==================== //
180
181 // 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
182 // and the truth LRO state.
183
184 // Therefore, we will actually run an estimation from a dispersed LRO state.
185 // The sc_seed is the true LRO state from the BSP.
186 let sc_seed = *traj_as_flown.first();
187
188 // Load the Deep Space Network ground stations.
189 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
190 let ground_station_file: PathBuf = [
191 env!("CARGO_MANIFEST_DIR"),
192 "examples",
193 "04_lro_od",
194 "dsn-network.yaml",
195 ]
196 .iter()
197 .collect();
198
199 let devices = GroundStation::load_named(ground_station_file)?;
200
201 let mut proc_devices = devices.clone();
202
203 // Increase the noise in the devices to accept more measurements.
204 for gs in proc_devices.values_mut() {
205 if let Some(noise) = &mut gs
206 .stochastic_noises
207 .as_mut()
208 .unwrap()
209 .get_mut(&MeasurementType::Range)
210 {
211 *noise.white_noise.as_mut().unwrap() *= 3.0;
212 }
213 }
214
215 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
216 // Nyx can build a tracking schedule for you based on the first station with access.
217 let trkconfg_yaml: PathBuf = [
218 env!("CARGO_MANIFEST_DIR"),
219 "examples",
220 "04_lro_od",
221 "tracking-cfg.yaml",
222 ]
223 .iter()
224 .collect();
225
226 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
227
228 // Build the tracking arc simulation to generate a "standard measurement".
229 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
230 devices.clone(),
231 traj_as_flown.clone(),
232 configs,
233 123, // Set a seed for reproducibility
234 )?;
235
236 trk.build_schedule(almanac.clone())?;
237 let arc = trk.generate_measurements(almanac.clone())?;
238 // Save the simulated tracking data
239 arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
240
241 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
242 println!("{arc}");
243
244 // Now that we have simulated measurements, we'll run the orbit determination.
245
246 // ===================== //
247 // === OD ESTIMATION === //
248 // ===================== //
249
250 let sc = SpacecraftUncertainty::builder()
251 .nominal(sc_seed)
252 .frame(LocalFrame::RIC)
253 .x_km(0.5)
254 .y_km(0.5)
255 .z_km(0.5)
256 .vx_km_s(5e-3)
257 .vy_km_s(5e-3)
258 .vz_km_s(5e-3)
259 .build();
260
261 // Build the filter initial estimate, which we will reuse in the filter.
262 let mut initial_estimate = sc.to_estimate()?;
263 initial_estimate.covar *= 3.0;
264
265 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
266
267 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
268 let process_noise = ProcessNoise3D::from_velocity_km_s(
269 &[1e-10, 1e-10, 1e-10],
270 1 * Unit::Hour,
271 10 * Unit::Minute,
272 None,
273 );
274
275 println!("{process_noise}");
276
277 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
278 let odp = SpacecraftKalmanOD::new(
279 setup,
280 KalmanVariant::ReferenceUpdate,
281 Some(ResidRejectCrit::default()),
282 proc_devices,
283 almanac.clone(),
284 )
285 .with_process_noise(process_noise);
286
287 let od_sol = odp.process_arc(initial_estimate, &arc)?;
288
289 let final_est = od_sol.estimates.last().unwrap();
290
291 println!("{final_est}");
292
293 let ric_err = traj_as_flown
294 .at(final_est.epoch())?
295 .orbit
296 .ric_difference(&final_est.orbital_state())?;
297 println!("== RIC at end ==");
298 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
299 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
300
301 println!(
302 "Num residuals rejected: #{}",
303 od_sol.rejected_residuals().len()
304 );
305 println!(
306 "Percentage within +/-3: {}",
307 od_sol.residual_ratio_within_threshold(3.0).unwrap()
308 );
309 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
310
311 od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
312
313 // In our case, we have the truth trajectory from NASA.
314 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
315 // Export the OD trajectory first.
316 let od_trajectory = od_sol.to_traj()?;
317 // Build the RIC difference.
318 od_trajectory.ric_diff_to_parquet(
319 &traj_as_flown,
320 "./04_lro_od_truth_error.parquet",
321 ExportCfg::default(),
322 )?;
323
324 Ok(())
325}
Source§impl SpacecraftUncertainty
impl SpacecraftUncertainty
Sourcepub fn to_estimate(&self) -> PhysicsResult<KfEstimate<Spacecraft>>
pub fn to_estimate(&self) -> PhysicsResult<KfEstimate<Spacecraft>>
Builds a Kalman filter estimate for a spacecraft state, ready to ingest into an OD Process.
Note: this function will rotate from the provided local frame into the inertial frame with the same central body.
Examples found in repository?
26fn main() -> Result<(), Box<dyn Error>> {
27 pel::init();
28 // Dynamics models require planetary constants and ephemerides to be defined.
29 // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
30 // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
31
32 // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
33 // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
34 let mut latest_jwst_ephem = MetaFile {
35 uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
36 crc32: None,
37 };
38 latest_jwst_ephem.process(true)?;
39
40 // Load this ephem in the general Almanac we're using for this analysis.
41 let almanac = Arc::new(
42 MetaAlmanac::latest()
43 .map_err(Box::new)?
44 .load_from_metafile(latest_jwst_ephem, true)?,
45 );
46
47 // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
48 // in the BSP. We need this ID in order to query the ephemeris.
49 const JWST_NAIF_ID: i32 = -170;
50 // Let's build a frame in the J2000 orientation centered on the JWST.
51 const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
52
53 // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
54 let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
55 println!("JWST defined from {earliest_epoch} to {latest_epoch}");
56 // Fetch the state, printing it in the Earth J2000 frame.
57 let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
58 println!("{jwst_orbit:x}");
59
60 // Build the spacecraft
61 // SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
62 // SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
63 let jwst = Spacecraft::builder()
64 .orbit(jwst_orbit)
65 .srp(SRPData {
66 area_m2: 21.197 * 14.162,
67 coeff_reflectivity: 1.56,
68 })
69 .mass(Mass::from_dry_mass(6200.0))
70 .build();
71
72 // Build up the spacecraft uncertainty builder.
73 // We can use the spacecraft uncertainty structure to build this up.
74 // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
75 // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
76 // Nyx can also estimate the deviation of the spacecraft parameters.
77 let jwst_uncertainty = SpacecraftUncertainty::builder()
78 .nominal(jwst)
79 .frame(LocalFrame::RIC)
80 .x_km(0.5)
81 .y_km(0.3)
82 .z_km(1.5)
83 .vx_km_s(1e-4)
84 .vy_km_s(0.6e-3)
85 .vz_km_s(3e-3)
86 .build();
87
88 println!("{jwst_uncertainty}");
89
90 // Build the Kalman filter estimate.
91 // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
92 // but this approach requires quite a bit more boilerplate code.
93 let jwst_estimate = jwst_uncertainty.to_estimate()?;
94
95 // Set up the spacecraft dynamics.
96 // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
97 // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
98
99 let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
100 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
101
102 // Finalize setting up the dynamics.
103 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
104
105 // Build the propagator set up to use for the whole analysis.
106 let setup = Propagator::default(dynamics);
107
108 // All of the analysis will use this duration.
109 let prediction_duration = 6.5 * Unit::Day;
110
111 // === Covariance mapping ===
112 // For the covariance mapping / prediction, we'll use the common orbit determination approach.
113 // This is done by setting up a spacecraft Kalman filter OD process, and predicting for the analysis duration.
114
115 // Build the propagation instance for the OD process.
116 let odp = SpacecraftKalmanOD::new(
117 setup.clone(),
118 KalmanVariant::DeviationTracking,
119 None,
120 BTreeMap::new(),
121 almanac.clone(),
122 );
123
124 // The prediction step is 1 minute by default, configured in the OD process, i.e. how often we want to know the covariance.
125 assert_eq!(odp.max_step, 1_i64.minutes());
126 // Finally, predict, and export the trajectory with covariance to a parquet file.
127 let od_sol = odp.predict_for(jwst_estimate, prediction_duration)?;
128 od_sol.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;
129
130 // === Monte Carlo framework ===
131 // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
132
133 let my_mc = MonteCarlo::new(
134 jwst, // Nominal state
135 jwst_estimate.to_random_variable()?,
136 "02_jwst".to_string(), // Scenario name
137 None, // No specific seed specified, so one will be drawn from the computer's entropy.
138 );
139
140 let num_runs = 5_000;
141 let rslts = my_mc.run_until_epoch(
142 setup,
143 almanac.clone(),
144 jwst.epoch() + prediction_duration,
145 num_runs,
146 );
147
148 assert_eq!(rslts.runs.len(), num_runs);
149 // Finally, export these results, computing the eclipse percentage for all of these results.
150
151 // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
152 let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
153 let umbra_event = eclipse_loc.to_umbra_event();
154 let penumbra_event = eclipse_loc.to_penumbra_event();
155
156 rslts.to_parquet(
157 "02_jwst_monte_carlo.parquet",
158 Some(vec![&umbra_event, &penumbra_event]),
159 ExportCfg::default(),
160 almanac,
161 )?;
162
163 Ok(())
164}
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 data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
45 .iter()
46 .collect();
47
48 let meta = data_folder.join("lro-dynamics.dhall");
49
50 // Load this ephem in the general Almanac we're using for this analysis.
51 let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
52 .map_err(Box::new)?
53 .process(true)
54 .map_err(Box::new)?;
55
56 let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
57 moon_pc.mu_km3_s2 = 4902.74987;
58 almanac.planetary_data.set_by_id(MOON, moon_pc)?;
59
60 let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
61 earth_pc.mu_km3_s2 = 398600.436;
62 almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
63
64 // Save this new kernel for reuse.
65 // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
66 almanac
67 .planetary_data
68 .save_as(&data_folder.join("lro-specific.pca"), true)?;
69
70 // Lock the almanac (an Arc is a read only structure).
71 let almanac = Arc::new(almanac);
72
73 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
74 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
75 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
76 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
77 let lro_frame = Frame::from_ephem_j2000(-85);
78
79 // To build the trajectory we need to provide a spacecraft template.
80 let sc_template = Spacecraft::builder()
81 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
82 .srp(SRPData {
83 // SRP configuration is arbitrary, but we will be estimating it anyway.
84 area_m2: 3.9 * 2.7,
85 coeff_reflectivity: 0.96,
86 })
87 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
88 .build();
89 // Now we can build the trajectory from the BSP file.
90 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
91 let traj_as_flown = Traj::from_bsp(
92 lro_frame,
93 MOON_J2000,
94 almanac.clone(),
95 sc_template,
96 5.seconds(),
97 Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
98 Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
99 Aberration::LT,
100 Some("LRO".to_string()),
101 )?;
102
103 println!("{traj_as_flown}");
104
105 // ====================== //
106 // === MODEL MATCHING === //
107 // ====================== //
108
109 // Set up the spacecraft dynamics.
110
111 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
112 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
113 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
114
115 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
116 // We're using the GRAIL JGGRX model.
117 let mut jggrx_meta = MetaFile {
118 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
119 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
120 };
121 // And let's download it if we don't have it yet.
122 jggrx_meta.process(true)?;
123
124 // Build the spherical harmonics.
125 // The harmonics must be computed in the body fixed frame.
126 // We're using the long term prediction of the Moon principal axes frame.
127 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
128 let sph_harmonics = Harmonics::from_stor(
129 almanac.frame_from_uid(moon_pa_frame)?,
130 HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
131 );
132
133 // Include the spherical harmonics into the orbital dynamics.
134 orbital_dyn.accel_models.push(sph_harmonics);
135
136 // We define the solar radiation pressure, using the default solar flux and accounting only
137 // for the eclipsing caused by the Earth and Moon.
138 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
139 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
140
141 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
142 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
143 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
144
145 println!("{dynamics}");
146
147 // Now we can build the propagator.
148 let setup = Propagator::default_dp78(dynamics.clone());
149
150 // For reference, let's build the trajectory with Nyx's models from that LRO state.
151 let (sim_final, traj_as_sim) = setup
152 .with(*traj_as_flown.first(), almanac.clone())
153 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
154
155 println!("SIM INIT: {:x}", traj_as_flown.first());
156 println!("SIM FINAL: {sim_final:x}");
157 // Compute RIC difference between SIM and LRO ephem
158 let sim_lro_delta = sim_final
159 .orbit
160 .ric_difference(&traj_as_flown.last().orbit)?;
161 println!("{traj_as_sim}");
162 println!(
163 "SIM v LRO - RIC Position (m): {:.3}",
164 sim_lro_delta.radius_km * 1e3
165 );
166 println!(
167 "SIM v LRO - RIC Velocity (m/s): {:.3}",
168 sim_lro_delta.velocity_km_s * 1e3
169 );
170
171 traj_as_sim.ric_diff_to_parquet(
172 &traj_as_flown,
173 "./04_lro_sim_truth_error.parquet",
174 ExportCfg::default(),
175 )?;
176
177 // ==================== //
178 // === OD SIMULATOR === //
179 // ==================== //
180
181 // 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
182 // and the truth LRO state.
183
184 // Therefore, we will actually run an estimation from a dispersed LRO state.
185 // The sc_seed is the true LRO state from the BSP.
186 let sc_seed = *traj_as_flown.first();
187
188 // Load the Deep Space Network ground stations.
189 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
190 let ground_station_file: PathBuf = [
191 env!("CARGO_MANIFEST_DIR"),
192 "examples",
193 "04_lro_od",
194 "dsn-network.yaml",
195 ]
196 .iter()
197 .collect();
198
199 let devices = GroundStation::load_named(ground_station_file)?;
200
201 let mut proc_devices = devices.clone();
202
203 // Increase the noise in the devices to accept more measurements.
204 for gs in proc_devices.values_mut() {
205 if let Some(noise) = &mut gs
206 .stochastic_noises
207 .as_mut()
208 .unwrap()
209 .get_mut(&MeasurementType::Range)
210 {
211 *noise.white_noise.as_mut().unwrap() *= 3.0;
212 }
213 }
214
215 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
216 // Nyx can build a tracking schedule for you based on the first station with access.
217 let trkconfg_yaml: PathBuf = [
218 env!("CARGO_MANIFEST_DIR"),
219 "examples",
220 "04_lro_od",
221 "tracking-cfg.yaml",
222 ]
223 .iter()
224 .collect();
225
226 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
227
228 // Build the tracking arc simulation to generate a "standard measurement".
229 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
230 devices.clone(),
231 traj_as_flown.clone(),
232 configs,
233 123, // Set a seed for reproducibility
234 )?;
235
236 trk.build_schedule(almanac.clone())?;
237 let arc = trk.generate_measurements(almanac.clone())?;
238 // Save the simulated tracking data
239 arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
240
241 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
242 println!("{arc}");
243
244 // Now that we have simulated measurements, we'll run the orbit determination.
245
246 // ===================== //
247 // === OD ESTIMATION === //
248 // ===================== //
249
250 let sc = SpacecraftUncertainty::builder()
251 .nominal(sc_seed)
252 .frame(LocalFrame::RIC)
253 .x_km(0.5)
254 .y_km(0.5)
255 .z_km(0.5)
256 .vx_km_s(5e-3)
257 .vy_km_s(5e-3)
258 .vz_km_s(5e-3)
259 .build();
260
261 // Build the filter initial estimate, which we will reuse in the filter.
262 let mut initial_estimate = sc.to_estimate()?;
263 initial_estimate.covar *= 3.0;
264
265 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
266
267 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
268 let process_noise = ProcessNoise3D::from_velocity_km_s(
269 &[1e-10, 1e-10, 1e-10],
270 1 * Unit::Hour,
271 10 * Unit::Minute,
272 None,
273 );
274
275 println!("{process_noise}");
276
277 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
278 let odp = SpacecraftKalmanOD::new(
279 setup,
280 KalmanVariant::ReferenceUpdate,
281 Some(ResidRejectCrit::default()),
282 proc_devices,
283 almanac.clone(),
284 )
285 .with_process_noise(process_noise);
286
287 let od_sol = odp.process_arc(initial_estimate, &arc)?;
288
289 let final_est = od_sol.estimates.last().unwrap();
290
291 println!("{final_est}");
292
293 let ric_err = traj_as_flown
294 .at(final_est.epoch())?
295 .orbit
296 .ric_difference(&final_est.orbital_state())?;
297 println!("== RIC at end ==");
298 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
299 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
300
301 println!(
302 "Num residuals rejected: #{}",
303 od_sol.rejected_residuals().len()
304 );
305 println!(
306 "Percentage within +/-3: {}",
307 od_sol.residual_ratio_within_threshold(3.0).unwrap()
308 );
309 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
310
311 od_sol.to_parquet("./04_lro_od_results.parquet", ExportCfg::default())?;
312
313 // In our case, we have the truth trajectory from NASA.
314 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
315 // Export the OD trajectory first.
316 let od_trajectory = od_sol.to_traj()?;
317 // Build the RIC difference.
318 od_trajectory.ric_diff_to_parquet(
319 &traj_as_flown,
320 "./04_lro_od_truth_error.parquet",
321 ExportCfg::default(),
322 )?;
323
324 Ok(())
325}
Sourcepub fn to_estimate_randomized(
&self,
seed: Option<u128>,
) -> PhysicsResult<KfEstimate<Spacecraft>>
pub fn to_estimate_randomized( &self, seed: Option<u128>, ) -> PhysicsResult<KfEstimate<Spacecraft>>
Returns an estimation whose nominal state is dispersed. Known bug in MultiVariate dispersion: https://github.com/nyx-space/nyx/issues/339
Examples found in repository?
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 let manifest_dir =
42 PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
43
44 let out = manifest_dir.join("data/04_output/");
45
46 let almanac = Arc::new(
47 Almanac::new(
48 &manifest_dir
49 .join("data/01_planetary/pck08.pca")
50 .to_string_lossy(),
51 )
52 .unwrap()
53 .load(
54 &manifest_dir
55 .join("data/01_planetary/de440s.bsp")
56 .to_string_lossy(),
57 )
58 .unwrap(),
59 );
60
61 let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
62 let moon_iau = almanac.frame_from_uid(IAU_MOON_FRAME).unwrap();
63
64 let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
65 let nrho = Orbit::cartesian(
66 166_473.631_302_239_7,
67 -274_715.487_253_382_7,
68 -211_233.210_176_686_7,
69 0.933_451_604_520_018_4,
70 0.436_775_046_841_900_9,
71 -0.082_211_021_250_348_95,
72 epoch,
73 eme2k,
74 );
75
76 let tx_nrho_sc = Spacecraft::from(nrho);
77
78 let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
79 println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
80
81 let bodies = vec![EARTH, SUN];
82 let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
83
84 let setup = Propagator::rk89(
85 dynamics,
86 IntegratorOptions::builder().max_step(0.5.minutes()).build(),
87 );
88
89 /* == Propagate the NRHO vehicle == */
90 let prop_time = 1.1 * state_luna.period().unwrap();
91
92 let (nrho_final, mut tx_traj) = setup
93 .with(tx_nrho_sc, almanac.clone())
94 .for_duration_with_traj(prop_time)
95 .unwrap();
96
97 tx_traj.name = Some("NRHO Tx SC".to_string());
98
99 println!("{tx_traj}");
100
101 /* == Propagate an LLO vehicle == */
102 let llo_orbit =
103 Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
104
105 let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
106
107 let (_, llo_traj) = setup
108 .with(llo_sc, almanac.clone())
109 .until_epoch_with_traj(nrho_final.epoch())
110 .unwrap();
111
112 // Export the subset of the first two hours.
113 llo_traj
114 .clone()
115 .filter_by_offset(..2.hours())
116 .to_parquet_simple(out.join("05_caps_llo_truth.pq"), almanac.clone())?;
117
118 /* == Setup the interlink == */
119
120 let mut measurement_types = IndexSet::new();
121 measurement_types.insert(MeasurementType::Range);
122 measurement_types.insert(MeasurementType::Doppler);
123
124 let mut stochastics = IndexMap::new();
125
126 let sa45_csac_allan_dev = 1e-11;
127
128 stochastics.insert(
129 MeasurementType::Range,
130 StochasticNoise::from_hardware_range_km(
131 sa45_csac_allan_dev,
132 10.0.seconds(),
133 link_specific::ChipRate::StandardT4B,
134 link_specific::SN0::Average,
135 ),
136 );
137
138 stochastics.insert(
139 MeasurementType::Doppler,
140 StochasticNoise::from_hardware_doppler_km_s(
141 sa45_csac_allan_dev,
142 10.0.seconds(),
143 link_specific::CarrierFreq::SBand,
144 link_specific::CN0::Average,
145 ),
146 );
147
148 let interlink = InterlinkTxSpacecraft {
149 traj: tx_traj,
150 measurement_types,
151 integration_time: None,
152 timestamp_noise_s: None,
153 ab_corr: Aberration::LT,
154 stochastic_noises: Some(stochastics),
155 };
156
157 // Devices are the transmitter, which is our NRHO vehicle.
158 let mut devices = BTreeMap::new();
159 devices.insert("NRHO Tx SC".to_string(), interlink);
160
161 let mut configs = BTreeMap::new();
162 configs.insert(
163 "NRHO Tx SC".to_string(),
164 TrkConfig::builder()
165 .strands(vec![Strand {
166 start: epoch,
167 end: nrho_final.epoch(),
168 }])
169 .build(),
170 );
171
172 let mut trk_sim =
173 TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
174 println!("{trk_sim}");
175
176 let trk_data = trk_sim.generate_measurements(almanac.clone()).unwrap();
177 println!("{trk_data}");
178
179 trk_data
180 .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
181 .unwrap();
182
183 // Run a truth OD where we estimate the LLO position
184 let llo_uncertainty = SpacecraftUncertainty::builder()
185 .nominal(llo_sc)
186 .x_km(1.0)
187 .y_km(1.0)
188 .z_km(1.0)
189 .vx_km_s(1e-3)
190 .vy_km_s(1e-3)
191 .vz_km_s(1e-3)
192 .build();
193
194 let mut proc_devices = devices.clone();
195
196 // Define the initial estimate, randomized, seed for reproducibility
197 let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
198 // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
199 initial_estimate.covar *= 2.5;
200
201 // Increase the noise in the devices to accept more measurements.
202
203 for link in proc_devices.values_mut() {
204 for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
205 *noise.white_noise.as_mut().unwrap() *= 3.0;
206 }
207 }
208
209 let init_err = initial_estimate
210 .orbital_state()
211 .ric_difference(&llo_orbit)
212 .unwrap();
213
214 println!("initial estimate:\n{initial_estimate}");
215 println!("RIC errors = {init_err}",);
216
217 let odp = InterlinkKalmanOD::new(
218 setup.clone(),
219 KalmanVariant::ReferenceUpdate,
220 Some(ResidRejectCrit::default()),
221 proc_devices,
222 almanac.clone(),
223 );
224
225 // Shrink the data to process.
226 let arc = trk_data.filter_by_offset(..2.hours());
227
228 let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
229
230 println!("{od_sol}");
231
232 od_sol
233 .to_parquet(
234 out.join(format!("05_caps_interlink_od_sol.pq")),
235 ExportCfg::default(),
236 )
237 .unwrap();
238
239 let od_traj = od_sol.to_traj().unwrap();
240
241 od_traj
242 .ric_diff_to_parquet(
243 &llo_traj,
244 out.join(format!("05_caps_interlink_llo_est_error.pq")),
245 ExportCfg::default(),
246 )
247 .unwrap();
248
249 let final_est = od_sol.estimates.last().unwrap();
250 assert!(final_est.within_3sigma(), "should be within 3 sigma");
251
252 println!("ESTIMATE\n{final_est:x}\n");
253 let truth = llo_traj.at(final_est.epoch()).unwrap();
254 println!("TRUTH\n{truth:x}");
255
256 let final_err = truth
257 .orbit
258 .ric_difference(&final_est.orbital_state())
259 .unwrap();
260 println!("ERROR {final_err}");
261
262 // Build the residuals versus reference plot.
263 let rvr_sol = odp
264 .process_arc(initial_estimate, &arc.resid_vs_ref_check())
265 .unwrap();
266
267 rvr_sol
268 .to_parquet(
269 out.join(format!("05_caps_interlink_resid_v_ref.pq")),
270 ExportCfg::default(),
271 )
272 .unwrap();
273
274 let final_rvr = rvr_sol.estimates.last().unwrap();
275
276 println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
277 println!(
278 "Pure prop error {:.3} m",
279 final_rvr
280 .orbital_state()
281 .ric_difference(&final_est.orbital_state())
282 .unwrap()
283 .rmag_km()
284 * 1e3
285 );
286
287 Ok(())
288}
Trait Implementations§
Source§impl Clone for SpacecraftUncertainty
impl Clone for SpacecraftUncertainty
Source§fn clone(&self) -> SpacecraftUncertainty
fn clone(&self) -> SpacecraftUncertainty
1.0.0 · Source§fn clone_from(&mut self, source: &Self)
fn clone_from(&mut self, source: &Self)
source
. Read moreSource§impl Debug for SpacecraftUncertainty
impl Debug for SpacecraftUncertainty
Source§impl Display for SpacecraftUncertainty
impl Display for SpacecraftUncertainty
impl Copy for SpacecraftUncertainty
Auto Trait Implementations§
impl Freeze for SpacecraftUncertainty
impl RefUnwindSafe for SpacecraftUncertainty
impl Send for SpacecraftUncertainty
impl Sync for SpacecraftUncertainty
impl Unpin for SpacecraftUncertainty
impl UnwindSafe for SpacecraftUncertainty
Blanket Implementations§
Source§impl<T> BorrowMut<T> for Twhere
T: ?Sized,
impl<T> BorrowMut<T> for Twhere
T: ?Sized,
Source§fn borrow_mut(&mut self) -> &mut T
fn borrow_mut(&mut self) -> &mut T
Source§impl<T> CloneToUninit for Twhere
T: Clone,
impl<T> CloneToUninit for Twhere
T: Clone,
§impl<T> Instrument for T
impl<T> Instrument for T
§fn instrument(self, span: Span) -> Instrumented<Self>
fn instrument(self, span: Span) -> Instrumented<Self>
§fn in_current_span(self) -> Instrumented<Self>
fn in_current_span(self) -> Instrumented<Self>
Source§impl<T> IntoEither for T
impl<T> IntoEither for T
Source§fn into_either(self, into_left: bool) -> Either<Self, Self>
fn into_either(self, into_left: bool) -> Either<Self, Self>
self
into a Left
variant of Either<Self, Self>
if into_left
is true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read moreSource§fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
self
into a Left
variant of Either<Self, Self>
if into_left(&self)
returns true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read more§impl<T> Pointable for T
impl<T> Pointable for T
§impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
§fn to_subset(&self) -> Option<SS>
fn to_subset(&self) -> Option<SS>
self
from the equivalent element of its
superset. Read more§fn is_in_subset(&self) -> bool
fn is_in_subset(&self) -> bool
self
is actually part of its subset T
(and can be converted to it).§fn to_subset_unchecked(&self) -> SS
fn to_subset_unchecked(&self) -> SS
self.to_subset
but without any property checks. Always succeeds.§fn from_subset(element: &SS) -> SP
fn from_subset(element: &SS) -> SP
self
to the equivalent element of its superset.