Trait TimeUnits
pub trait TimeUnits: Copy + Mul<Unit, Output = Duration> {
// Provided methods
fn centuries(self) -> Duration { ... }
fn weeks(self) -> Duration { ... }
fn days(self) -> Duration { ... }
fn hours(self) -> Duration { ... }
fn minutes(self) -> Duration { ... }
fn seconds(self) -> Duration { ... }
fn milliseconds(self) -> Duration { ... }
fn microseconds(self) -> Duration { ... }
fn nanoseconds(self) -> Duration { ... }
}Expand description
A trait to automatically convert some primitives to a duration
#[cfg(feature = "std")]
{
use hifitime::prelude::*;
use std::str::FromStr;
assert_eq!(Duration::from_str("1 d").unwrap(), 1.days());
assert_eq!(Duration::from_str("10.598 days").unwrap(), 10.598.days());
assert_eq!(Duration::from_str("10.598 min").unwrap(), 10.598.minutes());
assert_eq!(Duration::from_str("10.598 us").unwrap(), 10.598.microseconds());
assert_eq!(Duration::from_str("10.598 seconds").unwrap(), 10.598.seconds());
assert_eq!(Duration::from_str("10.598 nanosecond").unwrap(), 10.598.nanoseconds());
}Provided Methods§
fn centuries(self) -> Duration
fn weeks(self) -> Duration
fn days(self) -> Duration
fn hours(self) -> Duration
fn hours(self) -> Duration
Examples found in repository?
examples/05_cislunar_spacecraft_link_od/main.rs (line 115)
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}fn minutes(self) -> Duration
fn minutes(self) -> Duration
Examples found in repository?
examples/02_jwst_covar_monte_carlo/main.rs (line 125)
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
examples/05_cislunar_spacecraft_link_od/main.rs (line 86)
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}fn seconds(self) -> Duration
fn seconds(self) -> Duration
Examples found in repository?
examples/03_geo_analysis/stationkeeping.rs (line 106)
28fn main() -> Result<(), Box<dyn Error>> {
29 pel::init();
30 // Set up the dynamics like in the orbit raise.
31 let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
32 let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
33
34 // Define the GEO orbit, and we're just going to maintain it very tightly.
35 let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
36 let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
37 println!("{orbit:x}");
38
39 let sc = Spacecraft::builder()
40 .orbit(orbit)
41 .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
42 .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
43 .thruster(Thruster {
44 // "NEXT-STEP" row in Table 2
45 isp_s: 4435.0,
46 thrust_N: 0.472,
47 })
48 .mode(GuidanceMode::Thrust) // Start thrusting immediately.
49 .build();
50
51 // Set up the spacecraft dynamics like in the orbit raise example.
52
53 let prop_time = 30.0 * Unit::Day;
54
55 // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
56 let objectives = &[
57 Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km
58 Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
59 Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
60 ];
61
62 let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
63 println!("{ruggiero_ctrl}");
64
65 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
66
67 let mut jgm3_meta = MetaFile {
68 uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
69 crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
70 };
71 jgm3_meta.process(true)?;
72
73 let harmonics = Harmonics::from_stor(
74 almanac.frame_from_uid(IAU_EARTH_FRAME)?,
75 HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
76 );
77 orbital_dyn.accel_models.push(harmonics);
78
79 let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
80 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
81 .with_guidance_law(ruggiero_ctrl.clone());
82
83 println!("{sc_dynamics}");
84
85 // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
86
87 // Let's start by defining the dispersion.
88 // The MultivariateNormal structure allows us to define the dispersions in any of the orbital parameters, but these are applied directly in the Cartesian state space.
89 // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
90 let mc_rv = MvnSpacecraft::new(
91 sc,
92 vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
93 )?;
94
95 let my_mc = MonteCarlo::new(
96 sc, // Nominal state
97 mc_rv,
98 "03_geo_sk".to_string(), // Scenario name
99 None, // No specific seed specified, so one will be drawn from the computer's entropy.
100 );
101
102 // Build the propagator setup.
103 let setup = Propagator::rk89(
104 sc_dynamics.clone(),
105 IntegratorOptions::builder()
106 .min_step(10.0_f64.seconds())
107 .error_ctrl(ErrorControl::RSSCartesianStep)
108 .build(),
109 );
110
111 let num_runs = 25;
112 let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
113
114 assert_eq!(rslts.runs.len(), num_runs);
115
116 // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
117
118 rslts.to_parquet(
119 "03_geo_sk.parquet",
120 Some(vec![
121 &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
122 ]),
123 ExportCfg::default(),
124 almanac,
125 )?;
126
127 Ok(())
128}More examples
examples/03_geo_analysis/raise.rs (line 120)
27fn main() -> Result<(), Box<dyn Error>> {
28 pel::init();
29
30 // Dynamics models require planetary constants and ephemerides to be defined.
31 // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
32 // This will automatically download the DE440s planetary ephemeris,
33 // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
34 // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
35 // planetary constants kernels.
36 // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
37 // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
38 // references to many functions.
39 let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
40 // Fetch the EME2000 frame from the Almabac
41 let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
42 // Define the orbit epoch
43 let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
44
45 // Build the spacecraft itself.
46 // Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf
47 // for the "next gen" SEP characteristics.
48
49 // GTO start
50 let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
51
52 let sc = Spacecraft::builder()
53 .orbit(orbit)
54 .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
55 .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
56 .thruster(Thruster {
57 // "NEXT-STEP" row in Table 2
58 isp_s: 4435.0,
59 thrust_N: 0.472,
60 })
61 .mode(GuidanceMode::Thrust) // Start thrusting immediately.
62 .build();
63
64 let prop_time = 180.0 * Unit::Day;
65
66 // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
67 let objectives = &[
68 Objective::within_tolerance(StateParameter::SMA, 42_165.0, 20.0),
69 Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
70 Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
71 ];
72
73 // Ensure that we only thrust if we have more than 20% illumination.
74 let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
75 println!("{ruggiero_ctrl}");
76
77 // Define the high fidelity dynamics
78
79 // Set up the spacecraft dynamics.
80
81 // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
82 // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
83 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
84
85 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
86 // We're using the JGM3 model here, which is the default in GMAT.
87 let mut jgm3_meta = MetaFile {
88 uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
89 crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
90 };
91 // And let's download it if we don't have it yet.
92 jgm3_meta.process(true)?;
93
94 // Build the spherical harmonics.
95 // The harmonics must be computed in the body fixed frame.
96 // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
97 let harmonics = Harmonics::from_stor(
98 almanac.frame_from_uid(IAU_EARTH_FRAME)?,
99 HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true).unwrap(),
100 );
101
102 // Include the spherical harmonics into the orbital dynamics.
103 orbital_dyn.accel_models.push(harmonics);
104
105 // We define the solar radiation pressure, using the default solar flux and accounting only
106 // for the eclipsing caused by the Earth.
107 let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
108
109 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
110 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
111 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
112 .with_guidance_law(ruggiero_ctrl.clone());
113
114 println!("{orbit:x}");
115
116 // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
117 let (final_state, traj) = Propagator::rk89(
118 sc_dynamics.clone(),
119 IntegratorOptions::builder()
120 .min_step(10.0_f64.seconds())
121 .error_ctrl(ErrorControl::RSSCartesianStep)
122 .build(),
123 )
124 .with(sc, almanac.clone())
125 .for_duration_with_traj(prop_time)?;
126
127 let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
128 println!("{:x}", final_state.orbit);
129 println!("prop usage: {prop_usage:.3} kg");
130
131 // Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
132 traj.to_parquet(
133 "./03_geo_raise.parquet",
134 Some(vec![
135 &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
136 ]),
137 ExportCfg::default(),
138 almanac,
139 )?;
140
141 for status_line in ruggiero_ctrl.status(&final_state) {
142 println!("{status_line}");
143 }
144
145 ruggiero_ctrl
146 .achieved(&final_state)
147 .expect("objective not achieved");
148
149 Ok(())
150}examples/05_cislunar_spacecraft_link_od/main.rs (line 132)
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}examples/04_lro_od/main.rs (line 96)
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}fn milliseconds(self) -> Duration
fn microseconds(self) -> Duration
fn nanoseconds(self) -> Duration
Dyn Compatibility§
This trait is not dyn compatible.
In older versions of Rust, dyn compatibility was called "object safety", so this trait is not object safe.