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?
nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 114)
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
42
43 let out = manifest_dir.join("data/04_output/");
44
45 let almanac = Arc::new(
46 Almanac::new(
47 &manifest_dir
48 .join("data/01_planetary/pck08.pca")
49 .to_string_lossy(),
50 )
51 .unwrap()
52 .load(
53 &manifest_dir
54 .join("data/01_planetary/de440s.bsp")
55 .to_string_lossy(),
56 )
57 .unwrap(),
58 );
59
60 let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
61 let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();
62
63 let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
64 let nrho = Orbit::cartesian(
65 166_473.631_302_239_7,
66 -274_715.487_253_382_7,
67 -211_233.210_176_686_7,
68 0.933_451_604_520_018_4,
69 0.436_775_046_841_900_9,
70 -0.082_211_021_250_348_95,
71 epoch,
72 eme2k,
73 );
74
75 let tx_nrho_sc = Spacecraft::from(nrho);
76
77 let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
78 println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
79
80 let bodies = vec![EARTH, SUN];
81 let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
82
83 let setup = Propagator::rk89(
84 dynamics,
85 IntegratorOptions::builder().max_step(0.5.minutes()).build(),
86 );
87
88 /* == Propagate the NRHO vehicle == */
89 let prop_time = 1.1 * state_luna.period().unwrap();
90
91 let (nrho_final, mut tx_traj) = setup
92 .with(tx_nrho_sc, almanac.clone())
93 .for_duration_with_traj(prop_time)
94 .unwrap();
95
96 tx_traj.name = Some("NRHO Tx SC".to_string());
97
98 println!("{tx_traj}");
99
100 /* == Propagate an LLO vehicle == */
101 let llo_orbit =
102 Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
103
104 let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
105
106 let (_, llo_traj) = setup
107 .with(llo_sc, almanac.clone())
108 .until_epoch_with_traj(nrho_final.epoch())
109 .unwrap();
110
111 // Export the subset of the first two hours.
112 llo_traj
113 .clone()
114 .filter_by_offset(..2.hours())
115 .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;
116
117 /* == Setup the interlink == */
118
119 let mut measurement_types = IndexSet::new();
120 measurement_types.insert(MeasurementType::Range);
121 measurement_types.insert(MeasurementType::Doppler);
122
123 let mut stochastics = IndexMap::new();
124
125 let sa45_csac_allan_dev = 1e-11;
126
127 stochastics.insert(
128 MeasurementType::Range,
129 StochasticNoise::from_hardware_range_km(
130 sa45_csac_allan_dev,
131 10.0.seconds(),
132 link_specific::ChipRate::StandardT4B(),
133 link_specific::SN0::Average(),
134 ),
135 );
136
137 stochastics.insert(
138 MeasurementType::Doppler,
139 StochasticNoise::from_hardware_doppler_km_s(
140 sa45_csac_allan_dev,
141 10.0.seconds(),
142 link_specific::CarrierFreq::SBand(),
143 link_specific::CN0::Average(),
144 ),
145 );
146
147 let interlink = InterlinkTxSpacecraft {
148 traj: tx_traj,
149 measurement_types,
150 integration_time: None,
151 timestamp_noise_s: None,
152 ab_corr: Aberration::LT,
153 stochastic_noises: Some(stochastics),
154 };
155
156 // Devices are the transmitter, which is our NRHO vehicle.
157 let mut devices = BTreeMap::new();
158 devices.insert("NRHO Tx SC".to_string(), interlink);
159
160 let mut configs = BTreeMap::new();
161 configs.insert(
162 "NRHO Tx SC".to_string(),
163 TrkConfig::builder()
164 .strands(vec![Strand {
165 start: epoch,
166 end: nrho_final.epoch(),
167 }])
168 .build(),
169 );
170
171 let mut trk_sim =
172 TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
173 println!("{trk_sim}");
174
175 let trk_data = trk_sim.generate_measurements(&almanac).unwrap();
176 println!("{trk_data}");
177
178 trk_data
179 .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
180 .unwrap();
181
182 // Run a truth OD where we estimate the LLO position
183 let llo_uncertainty = SpacecraftUncertainty::builder()
184 .nominal(llo_sc)
185 .x_km(1.0)
186 .y_km(1.0)
187 .z_km(1.0)
188 .vx_km_s(1e-3)
189 .vy_km_s(1e-3)
190 .vz_km_s(1e-3)
191 .build();
192
193 let mut proc_devices = devices.clone();
194
195 // Define the initial estimate, randomized, seed for reproducibility
196 let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
197 // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
198 initial_estimate.covar *= 2.5;
199
200 // Increase the noise in the devices to accept more measurements.
201
202 for link in proc_devices.values_mut() {
203 for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
204 *noise.white_noise.as_mut().unwrap() *= 3.0;
205 }
206 }
207
208 let init_err = initial_estimate
209 .orbital_state()
210 .ric_difference(&llo_orbit)
211 .unwrap();
212
213 println!("initial estimate:\n{initial_estimate}");
214 println!("RIC errors = {init_err}",);
215
216 let odp = InterlinkKalmanOD::new(
217 setup.clone(),
218 KalmanVariant::ReferenceUpdate,
219 Some(SigmaRejection::default()),
220 proc_devices,
221 almanac.clone(),
222 );
223
224 // Shrink the data to process.
225 let arc = trk_data.filter_by_offset(..2.hours());
226
227 let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
228
229 println!("{od_sol}");
230
231 od_sol
232 .to_parquet(
233 out.join("05_caps_interlink_od_sol.pq"),
234 ExportCfg::default(),
235 )
236 .unwrap();
237
238 let od_traj = od_sol.to_traj().unwrap();
239
240 od_traj
241 .ric_diff_to_parquet(
242 &llo_traj,
243 out.join("05_caps_interlink_llo_est_error.pq"),
244 ExportCfg::default(),
245 )
246 .unwrap();
247
248 let final_est = od_sol.estimates.last().unwrap();
249 assert!(final_est.within_3sigma(), "should be within 3 sigma");
250
251 println!("ESTIMATE\n{final_est:x}\n");
252 let truth = llo_traj.at(final_est.epoch()).unwrap();
253 println!("TRUTH\n{truth:x}");
254
255 let final_err = truth
256 .orbit
257 .ric_difference(&final_est.orbital_state())
258 .unwrap();
259 println!("ERROR {final_err}");
260
261 // Build the residuals versus reference plot.
262 let rvr_sol = odp
263 .process_arc(initial_estimate, &arc.resid_vs_ref_check())
264 .unwrap();
265
266 rvr_sol
267 .to_parquet(
268 out.join("05_caps_interlink_resid_v_ref.pq"),
269 ExportCfg::default(),
270 )
271 .unwrap();
272
273 let final_rvr = rvr_sol.estimates.last().unwrap();
274
275 println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
276 println!(
277 "Pure prop error {:.3} m",
278 final_rvr
279 .orbital_state()
280 .ric_difference(&final_est.orbital_state())
281 .unwrap()
282 .rmag_km()
283 * 1e3
284 );
285
286 Ok(())
287}fn minutes(self) -> Duration
fn minutes(self) -> Duration
Examples found in repository?
nyx-core/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)?;
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 rslts.to_parquet("02_jwst_monte_carlo.parquet", ExportCfg::default())?;
152
153 Ok(())
154}More examples
nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 85)
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
42
43 let out = manifest_dir.join("data/04_output/");
44
45 let almanac = Arc::new(
46 Almanac::new(
47 &manifest_dir
48 .join("data/01_planetary/pck08.pca")
49 .to_string_lossy(),
50 )
51 .unwrap()
52 .load(
53 &manifest_dir
54 .join("data/01_planetary/de440s.bsp")
55 .to_string_lossy(),
56 )
57 .unwrap(),
58 );
59
60 let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
61 let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();
62
63 let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
64 let nrho = Orbit::cartesian(
65 166_473.631_302_239_7,
66 -274_715.487_253_382_7,
67 -211_233.210_176_686_7,
68 0.933_451_604_520_018_4,
69 0.436_775_046_841_900_9,
70 -0.082_211_021_250_348_95,
71 epoch,
72 eme2k,
73 );
74
75 let tx_nrho_sc = Spacecraft::from(nrho);
76
77 let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
78 println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
79
80 let bodies = vec![EARTH, SUN];
81 let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
82
83 let setup = Propagator::rk89(
84 dynamics,
85 IntegratorOptions::builder().max_step(0.5.minutes()).build(),
86 );
87
88 /* == Propagate the NRHO vehicle == */
89 let prop_time = 1.1 * state_luna.period().unwrap();
90
91 let (nrho_final, mut tx_traj) = setup
92 .with(tx_nrho_sc, almanac.clone())
93 .for_duration_with_traj(prop_time)
94 .unwrap();
95
96 tx_traj.name = Some("NRHO Tx SC".to_string());
97
98 println!("{tx_traj}");
99
100 /* == Propagate an LLO vehicle == */
101 let llo_orbit =
102 Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
103
104 let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
105
106 let (_, llo_traj) = setup
107 .with(llo_sc, almanac.clone())
108 .until_epoch_with_traj(nrho_final.epoch())
109 .unwrap();
110
111 // Export the subset of the first two hours.
112 llo_traj
113 .clone()
114 .filter_by_offset(..2.hours())
115 .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;
116
117 /* == Setup the interlink == */
118
119 let mut measurement_types = IndexSet::new();
120 measurement_types.insert(MeasurementType::Range);
121 measurement_types.insert(MeasurementType::Doppler);
122
123 let mut stochastics = IndexMap::new();
124
125 let sa45_csac_allan_dev = 1e-11;
126
127 stochastics.insert(
128 MeasurementType::Range,
129 StochasticNoise::from_hardware_range_km(
130 sa45_csac_allan_dev,
131 10.0.seconds(),
132 link_specific::ChipRate::StandardT4B(),
133 link_specific::SN0::Average(),
134 ),
135 );
136
137 stochastics.insert(
138 MeasurementType::Doppler,
139 StochasticNoise::from_hardware_doppler_km_s(
140 sa45_csac_allan_dev,
141 10.0.seconds(),
142 link_specific::CarrierFreq::SBand(),
143 link_specific::CN0::Average(),
144 ),
145 );
146
147 let interlink = InterlinkTxSpacecraft {
148 traj: tx_traj,
149 measurement_types,
150 integration_time: None,
151 timestamp_noise_s: None,
152 ab_corr: Aberration::LT,
153 stochastic_noises: Some(stochastics),
154 };
155
156 // Devices are the transmitter, which is our NRHO vehicle.
157 let mut devices = BTreeMap::new();
158 devices.insert("NRHO Tx SC".to_string(), interlink);
159
160 let mut configs = BTreeMap::new();
161 configs.insert(
162 "NRHO Tx SC".to_string(),
163 TrkConfig::builder()
164 .strands(vec![Strand {
165 start: epoch,
166 end: nrho_final.epoch(),
167 }])
168 .build(),
169 );
170
171 let mut trk_sim =
172 TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
173 println!("{trk_sim}");
174
175 let trk_data = trk_sim.generate_measurements(&almanac).unwrap();
176 println!("{trk_data}");
177
178 trk_data
179 .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
180 .unwrap();
181
182 // Run a truth OD where we estimate the LLO position
183 let llo_uncertainty = SpacecraftUncertainty::builder()
184 .nominal(llo_sc)
185 .x_km(1.0)
186 .y_km(1.0)
187 .z_km(1.0)
188 .vx_km_s(1e-3)
189 .vy_km_s(1e-3)
190 .vz_km_s(1e-3)
191 .build();
192
193 let mut proc_devices = devices.clone();
194
195 // Define the initial estimate, randomized, seed for reproducibility
196 let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
197 // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
198 initial_estimate.covar *= 2.5;
199
200 // Increase the noise in the devices to accept more measurements.
201
202 for link in proc_devices.values_mut() {
203 for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
204 *noise.white_noise.as_mut().unwrap() *= 3.0;
205 }
206 }
207
208 let init_err = initial_estimate
209 .orbital_state()
210 .ric_difference(&llo_orbit)
211 .unwrap();
212
213 println!("initial estimate:\n{initial_estimate}");
214 println!("RIC errors = {init_err}",);
215
216 let odp = InterlinkKalmanOD::new(
217 setup.clone(),
218 KalmanVariant::ReferenceUpdate,
219 Some(SigmaRejection::default()),
220 proc_devices,
221 almanac.clone(),
222 );
223
224 // Shrink the data to process.
225 let arc = trk_data.filter_by_offset(..2.hours());
226
227 let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
228
229 println!("{od_sol}");
230
231 od_sol
232 .to_parquet(
233 out.join("05_caps_interlink_od_sol.pq"),
234 ExportCfg::default(),
235 )
236 .unwrap();
237
238 let od_traj = od_sol.to_traj().unwrap();
239
240 od_traj
241 .ric_diff_to_parquet(
242 &llo_traj,
243 out.join("05_caps_interlink_llo_est_error.pq"),
244 ExportCfg::default(),
245 )
246 .unwrap();
247
248 let final_est = od_sol.estimates.last().unwrap();
249 assert!(final_est.within_3sigma(), "should be within 3 sigma");
250
251 println!("ESTIMATE\n{final_est:x}\n");
252 let truth = llo_traj.at(final_est.epoch()).unwrap();
253 println!("TRUTH\n{truth:x}");
254
255 let final_err = truth
256 .orbit
257 .ric_difference(&final_est.orbital_state())
258 .unwrap();
259 println!("ERROR {final_err}");
260
261 // Build the residuals versus reference plot.
262 let rvr_sol = odp
263 .process_arc(initial_estimate, &arc.resid_vs_ref_check())
264 .unwrap();
265
266 rvr_sol
267 .to_parquet(
268 out.join("05_caps_interlink_resid_v_ref.pq"),
269 ExportCfg::default(),
270 )
271 .unwrap();
272
273 let final_rvr = rvr_sol.estimates.last().unwrap();
274
275 println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
276 println!(
277 "Pure prop error {:.3} m",
278 final_rvr
279 .orbital_state()
280 .ric_difference(&final_est.orbital_state())
281 .unwrap()
282 .rmag_km()
283 * 1e3
284 );
285
286 Ok(())
287}fn seconds(self) -> Duration
fn seconds(self) -> Duration
Examples found in repository?
nyx-core/examples/03_geo_analysis/stationkeeping.rs (line 123)
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_info(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(
58 StateParameter::Element(OrbitalElement::SemiMajorAxis),
59 42_165.0,
60 20.0,
61 ),
62 Objective::within_tolerance(
63 StateParameter::Element(OrbitalElement::Eccentricity),
64 0.001,
65 5e-5,
66 ),
67 Objective::within_tolerance(
68 StateParameter::Element(OrbitalElement::Inclination),
69 0.05,
70 1e-2,
71 ),
72 ];
73
74 let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
75 println!("{ruggiero_ctrl}");
76
77 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
78
79 let mut jgm3_meta = MetaFile {
80 uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
81 crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
82 };
83 jgm3_meta.process(true)?;
84
85 let harmonics = GravityField::new(GravityFieldData::from_cof(
86 &jgm3_meta.uri,
87 8,
88 8,
89 almanac.frame_info(IAU_EARTH_FRAME)?,
90 )?);
91 orbital_dyn.accel_models.push(harmonics);
92
93 let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;
94 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
95 .with_guidance_law(ruggiero_ctrl.clone());
96
97 println!("{sc_dynamics}");
98
99 // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
100
101 // Let's start by defining the dispersion.
102 // 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.
103 // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
104 let mc_rv = MvnSpacecraft::new(
105 sc,
106 vec![StateDispersion::zero_mean(
107 StateParameter::Element(OrbitalElement::SemiMajorAxis),
108 3.0,
109 )],
110 )?;
111
112 let my_mc = MonteCarlo::new(
113 sc, // Nominal state
114 mc_rv,
115 "03_geo_sk".to_string(), // Scenario name
116 None, // No specific seed specified, so one will be drawn from the computer's entropy.
117 );
118
119 // Build the propagator setup.
120 let setup = Propagator::rk89(
121 sc_dynamics.clone(),
122 IntegratorOptions::builder()
123 .min_step(10.0_f64.seconds())
124 .error_ctrl(ErrorControl::RSSCartesianStep)
125 .build(),
126 );
127
128 let num_runs = 25;
129 let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
130
131 assert_eq!(rslts.runs.len(), num_runs);
132
133 rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;
134
135 Ok(())
136}More examples
nyx-core/examples/03_geo_analysis/raise.rs (line 132)
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_info(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(
69 StateParameter::Element(OrbitalElement::SemiMajorAxis),
70 42_165.0,
71 20.0,
72 ),
73 Objective::within_tolerance(
74 StateParameter::Element(OrbitalElement::Eccentricity),
75 0.001,
76 5e-5,
77 ),
78 Objective::within_tolerance(
79 StateParameter::Element(OrbitalElement::Inclination),
80 0.05,
81 1e-2,
82 ),
83 ];
84
85 // Ensure that we only thrust if we have more than 20% illumination.
86 let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
87 println!("{ruggiero_ctrl}");
88
89 // Define the high fidelity dynamics
90
91 // Set up the spacecraft dynamics.
92
93 // Specify that the orbital dynamics must account for the graviational pull of the Moon and the Sun.
94 // The gravity of the Earth will also be accounted for since the spaceraft in an Earth orbit.
95 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
96
97 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
98 // We're using the JGM3 model here, which is the default in GMAT.
99 let mut jgm3_meta = MetaFile {
100 uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
101 crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
102 };
103 // And let's download it if we don't have it yet.
104 jgm3_meta.process(true)?;
105
106 // Build the spherical harmonics.
107 // The harmonics must be computed in the body fixed frame.
108 // We're using the long term prediction of the Earth centered Earth fixed frame, IAU Earth.
109 let harmonics = GravityField::new(
110 GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, almanac.frame_info(IAU_EARTH_FRAME)?)
111 .unwrap(),
112 );
113
114 // Include the spherical harmonics into the orbital dynamics.
115 orbital_dyn.accel_models.push(harmonics);
116
117 // We define the solar radiation pressure, using the default solar flux and accounting only
118 // for the eclipsing caused by the Earth.
119 let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;
120
121 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
122 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
123 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
124 .with_guidance_law(ruggiero_ctrl.clone());
125
126 println!("{orbit:x}");
127
128 // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
129 let (final_state, traj) = Propagator::rk89(
130 sc_dynamics.clone(),
131 IntegratorOptions::builder()
132 .min_step(10.0_f64.seconds())
133 .error_ctrl(ErrorControl::RSSCartesianStep)
134 .build(),
135 )
136 .with(sc, almanac.clone())
137 .for_duration_with_traj(prop_time)?;
138
139 let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
140 println!("{:x}", final_state.orbit);
141 println!("prop usage: {prop_usage:.3} kg");
142
143 // Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
144 traj.to_parquet("./03_geo_raise.parquet", ExportCfg::default())?;
145
146 for status_line in ruggiero_ctrl.status(&final_state) {
147 println!("{status_line}");
148 }
149
150 ruggiero_ctrl
151 .achieved(&final_state)
152 .expect("objective not achieved");
153
154 Ok(())
155}nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 131)
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
42
43 let out = manifest_dir.join("data/04_output/");
44
45 let almanac = Arc::new(
46 Almanac::new(
47 &manifest_dir
48 .join("data/01_planetary/pck08.pca")
49 .to_string_lossy(),
50 )
51 .unwrap()
52 .load(
53 &manifest_dir
54 .join("data/01_planetary/de440s.bsp")
55 .to_string_lossy(),
56 )
57 .unwrap(),
58 );
59
60 let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
61 let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();
62
63 let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
64 let nrho = Orbit::cartesian(
65 166_473.631_302_239_7,
66 -274_715.487_253_382_7,
67 -211_233.210_176_686_7,
68 0.933_451_604_520_018_4,
69 0.436_775_046_841_900_9,
70 -0.082_211_021_250_348_95,
71 epoch,
72 eme2k,
73 );
74
75 let tx_nrho_sc = Spacecraft::from(nrho);
76
77 let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
78 println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
79
80 let bodies = vec![EARTH, SUN];
81 let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
82
83 let setup = Propagator::rk89(
84 dynamics,
85 IntegratorOptions::builder().max_step(0.5.minutes()).build(),
86 );
87
88 /* == Propagate the NRHO vehicle == */
89 let prop_time = 1.1 * state_luna.period().unwrap();
90
91 let (nrho_final, mut tx_traj) = setup
92 .with(tx_nrho_sc, almanac.clone())
93 .for_duration_with_traj(prop_time)
94 .unwrap();
95
96 tx_traj.name = Some("NRHO Tx SC".to_string());
97
98 println!("{tx_traj}");
99
100 /* == Propagate an LLO vehicle == */
101 let llo_orbit =
102 Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
103
104 let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
105
106 let (_, llo_traj) = setup
107 .with(llo_sc, almanac.clone())
108 .until_epoch_with_traj(nrho_final.epoch())
109 .unwrap();
110
111 // Export the subset of the first two hours.
112 llo_traj
113 .clone()
114 .filter_by_offset(..2.hours())
115 .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;
116
117 /* == Setup the interlink == */
118
119 let mut measurement_types = IndexSet::new();
120 measurement_types.insert(MeasurementType::Range);
121 measurement_types.insert(MeasurementType::Doppler);
122
123 let mut stochastics = IndexMap::new();
124
125 let sa45_csac_allan_dev = 1e-11;
126
127 stochastics.insert(
128 MeasurementType::Range,
129 StochasticNoise::from_hardware_range_km(
130 sa45_csac_allan_dev,
131 10.0.seconds(),
132 link_specific::ChipRate::StandardT4B(),
133 link_specific::SN0::Average(),
134 ),
135 );
136
137 stochastics.insert(
138 MeasurementType::Doppler,
139 StochasticNoise::from_hardware_doppler_km_s(
140 sa45_csac_allan_dev,
141 10.0.seconds(),
142 link_specific::CarrierFreq::SBand(),
143 link_specific::CN0::Average(),
144 ),
145 );
146
147 let interlink = InterlinkTxSpacecraft {
148 traj: tx_traj,
149 measurement_types,
150 integration_time: None,
151 timestamp_noise_s: None,
152 ab_corr: Aberration::LT,
153 stochastic_noises: Some(stochastics),
154 };
155
156 // Devices are the transmitter, which is our NRHO vehicle.
157 let mut devices = BTreeMap::new();
158 devices.insert("NRHO Tx SC".to_string(), interlink);
159
160 let mut configs = BTreeMap::new();
161 configs.insert(
162 "NRHO Tx SC".to_string(),
163 TrkConfig::builder()
164 .strands(vec![Strand {
165 start: epoch,
166 end: nrho_final.epoch(),
167 }])
168 .build(),
169 );
170
171 let mut trk_sim =
172 TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
173 println!("{trk_sim}");
174
175 let trk_data = trk_sim.generate_measurements(&almanac).unwrap();
176 println!("{trk_data}");
177
178 trk_data
179 .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
180 .unwrap();
181
182 // Run a truth OD where we estimate the LLO position
183 let llo_uncertainty = SpacecraftUncertainty::builder()
184 .nominal(llo_sc)
185 .x_km(1.0)
186 .y_km(1.0)
187 .z_km(1.0)
188 .vx_km_s(1e-3)
189 .vy_km_s(1e-3)
190 .vz_km_s(1e-3)
191 .build();
192
193 let mut proc_devices = devices.clone();
194
195 // Define the initial estimate, randomized, seed for reproducibility
196 let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
197 // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
198 initial_estimate.covar *= 2.5;
199
200 // Increase the noise in the devices to accept more measurements.
201
202 for link in proc_devices.values_mut() {
203 for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
204 *noise.white_noise.as_mut().unwrap() *= 3.0;
205 }
206 }
207
208 let init_err = initial_estimate
209 .orbital_state()
210 .ric_difference(&llo_orbit)
211 .unwrap();
212
213 println!("initial estimate:\n{initial_estimate}");
214 println!("RIC errors = {init_err}",);
215
216 let odp = InterlinkKalmanOD::new(
217 setup.clone(),
218 KalmanVariant::ReferenceUpdate,
219 Some(SigmaRejection::default()),
220 proc_devices,
221 almanac.clone(),
222 );
223
224 // Shrink the data to process.
225 let arc = trk_data.filter_by_offset(..2.hours());
226
227 let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
228
229 println!("{od_sol}");
230
231 od_sol
232 .to_parquet(
233 out.join("05_caps_interlink_od_sol.pq"),
234 ExportCfg::default(),
235 )
236 .unwrap();
237
238 let od_traj = od_sol.to_traj().unwrap();
239
240 od_traj
241 .ric_diff_to_parquet(
242 &llo_traj,
243 out.join("05_caps_interlink_llo_est_error.pq"),
244 ExportCfg::default(),
245 )
246 .unwrap();
247
248 let final_est = od_sol.estimates.last().unwrap();
249 assert!(final_est.within_3sigma(), "should be within 3 sigma");
250
251 println!("ESTIMATE\n{final_est:x}\n");
252 let truth = llo_traj.at(final_est.epoch()).unwrap();
253 println!("TRUTH\n{truth:x}");
254
255 let final_err = truth
256 .orbit
257 .ric_difference(&final_est.orbital_state())
258 .unwrap();
259 println!("ERROR {final_err}");
260
261 // Build the residuals versus reference plot.
262 let rvr_sol = odp
263 .process_arc(initial_estimate, &arc.resid_vs_ref_check())
264 .unwrap();
265
266 rvr_sol
267 .to_parquet(
268 out.join("05_caps_interlink_resid_v_ref.pq"),
269 ExportCfg::default(),
270 )
271 .unwrap();
272
273 let final_rvr = rvr_sol.estimates.last().unwrap();
274
275 println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
276 println!(
277 "Pure prop error {:.3} m",
278 final_rvr
279 .orbital_state()
280 .ric_difference(&final_est.orbital_state())
281 .unwrap()
282 .rmag_km()
283 * 1e3
284 );
285
286 Ok(())
287}nyx-core/examples/04_lro_od/main.rs (line 85)
34fn main() -> Result<(), Box<dyn Error>> {
35 pel::init();
36
37 // ====================== //
38 // === ALMANAC SET UP === //
39 // ====================== //
40
41 // Dynamics models require planetary constants and ephemerides to be defined.
42 // Let's start by grabbing those by using ANISE's MetaAlmanac.
43
44 let output_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "../data", "04_output"]
45 .iter()
46 .collect();
47
48 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
49 .iter()
50 .collect();
51
52 let meta = data_folder.join("lro-dynamics.dhall");
53
54 // Load this ephem in the general Almanac we're using for this analysis.
55 let almanac = Arc::new(
56 MetaAlmanac::new(meta.to_string_lossy().as_ref())
57 .map_err(Box::new)?
58 .process(true)
59 .map_err(Box::new)?,
60 );
61
62 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
63 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
64 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
65 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
66 let lro_frame = Frame::from_ephem_j2000(-85);
67
68 // To build the trajectory we need to provide a spacecraft template.
69 let sc_template = Spacecraft::builder()
70 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
71 .srp(SRPData {
72 // SRP configuration is arbitrary, but we will be estimating it anyway.
73 area_m2: 3.9 * 2.7,
74 coeff_reflectivity: 0.96,
75 })
76 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
77 .build();
78 // Now we can build the trajectory from the BSP file.
79 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
80 let traj_as_flown = Traj::from_bsp(
81 lro_frame,
82 MOON_J2000,
83 &almanac,
84 sc_template,
85 5.seconds(),
86 Some(Epoch::from_str("2024-01-01 01:00:00 UTC")?),
87 Some(Epoch::from_str("2024-01-02 01:00:00 UTC")?),
88 None,
89 Some("LRO".to_string()),
90 )?;
91
92 println!("{traj_as_flown}");
93
94 // ====================== //
95 // === MODEL MATCHING === //
96 // ====================== //
97
98 // Set up the spacecraft dynamics.
99
100 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
101 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
102 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
103
104 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
105 // We're using the GRAIL JGGRX model.
106 let mut jggrx_meta = MetaFile {
107 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
108 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
109 };
110 // And let's download it if we don't have it yet.
111 jggrx_meta.process(true)?;
112
113 // Build the spherical harmonics.
114 // The harmonics must be computed in the body fixed frame.
115 // We're using the long term prediction of the Moon principal axes frame.
116 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
117 let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
118 &jggrx_meta.uri,
119 81,
120 81,
121 almanac.frame_info(moon_pa_frame)?,
122 )?);
123
124 // Include the spherical harmonics into the orbital dynamics.
125 orbital_dyn.accel_models.push(sph_harmonics);
126
127 // We define the solar radiation pressure, using the default solar flux and accounting only
128 // for the eclipsing caused by the Earth and Moon.
129 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
130 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;
131
132 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
133 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
134 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
135
136 println!("{dynamics}");
137
138 // Now we can build the propagator.
139 let setup = Propagator::default_dp78(dynamics.clone());
140
141 // For reference, let's build the trajectory with Nyx's models from that LRO state.
142 let (sim_final, traj_as_sim) = setup
143 .with(*traj_as_flown.first(), almanac.clone())
144 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
145
146 println!("SIM INIT: {:x}", traj_as_flown.first());
147 println!("SIM FINAL: {sim_final:x}");
148 // Compute RIC difference between SIM and LRO ephem
149 let sim_lro_delta = sim_final
150 .orbit
151 .ric_difference(&traj_as_flown.last().orbit)?;
152 println!("{traj_as_sim}");
153 println!(
154 "SIM v LRO - RIC Position (m): {:.3}",
155 sim_lro_delta.radius_km * 1e3
156 );
157 println!(
158 "SIM v LRO - RIC Velocity (m/s): {:.3}",
159 sim_lro_delta.velocity_km_s * 1e3
160 );
161
162 traj_as_sim.ric_diff_to_parquet(
163 &traj_as_flown,
164 output_folder.join("./04_lro_sim_truth_error.parquet"),
165 ExportCfg::default(),
166 )?;
167
168 // ==================== //
169 // === OD SIMULATOR === //
170 // ==================== //
171
172 // After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
173 // and the truth LRO state.
174
175 // Therefore, we will actually run an estimation from a dispersed LRO state.
176 // The sc_seed is the true LRO state from the BSP.
177 let sc_seed = *traj_as_flown.first();
178
179 // Load the Deep Space Network ground stations.
180 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
181 let ground_station_file: PathBuf = [
182 env!("CARGO_MANIFEST_DIR"),
183 "examples",
184 "04_lro_od",
185 "dsn-network.yaml",
186 ]
187 .iter()
188 .collect();
189
190 let devices = GroundStation::load_named(ground_station_file)?;
191
192 let proc_devices = devices.clone();
193
194 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
195 // Nyx can build a tracking schedule for you based on the first station with access.
196 let trkconfg_yaml: PathBuf = [
197 env!("CARGO_MANIFEST_DIR"),
198 "examples",
199 "04_lro_od",
200 "tracking-cfg.yaml",
201 ]
202 .iter()
203 .collect();
204
205 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
206
207 // Build the tracking arc simulation to generate a "standard measurement".
208 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
209 devices.clone(),
210 traj_as_flown.clone(),
211 configs,
212 123, // Set a seed for reproducibility
213 )?;
214
215 trk.build_schedule(&almanac)?;
216 let arc = trk.generate_measurements(&almanac)?;
217 // Save the simulated tracking data
218 arc.to_parquet_simple(output_folder.join("04_lro_simulated_tracking.parquet"))?;
219
220 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
221 println!("{arc}");
222
223 // Now that we have simulated measurements, we'll run the orbit determination.
224
225 // ===================== //
226 // === OD ESTIMATION === //
227 // ===================== //
228
229 let sc = SpacecraftUncertainty::builder()
230 .nominal(sc_seed)
231 .frame(LocalFrame::RIC)
232 .x_km(0.5)
233 .y_km(0.5)
234 .z_km(0.5)
235 .vx_km_s(5e-3)
236 .vy_km_s(5e-3)
237 .vz_km_s(5e-3)
238 .build();
239
240 // Build the filter initial estimate, which we will reuse in the filter.
241 let mut initial_estimate = sc.to_estimate()?;
242 initial_estimate.covar *= 1.5;
243
244 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
245
246 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
247 let process_noise = ProcessNoise3D::from_velocity_km_s(
248 &[5e-13, 5e-13, 5e-13],
249 1 * Unit::Hour,
250 10 * Unit::Minute,
251 None,
252 );
253
254 println!("{process_noise}");
255
256 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
257 let odp = SpacecraftKalmanOD::new(
258 setup,
259 KalmanVariant::ReferenceUpdate,
260 Some(SigmaRejection::default()),
261 proc_devices,
262 almanac.clone(),
263 )
264 .with_process_noise(process_noise);
265
266 let od_sol = odp.process_arc(initial_estimate, &arc)?;
267
268 let final_est = od_sol.estimates.last().unwrap();
269
270 println!("{final_est}");
271
272 let ric_err = traj_as_flown
273 .at(final_est.epoch())?
274 .orbit
275 .ric_difference(&final_est.orbital_state())?;
276 println!("== RIC at end ==");
277 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
278 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
279
280 println!(
281 "Num residuals rejected: #{}",
282 od_sol.rejected_residuals().len()
283 );
284 println!(
285 "Percentage within +/-3: {}",
286 od_sol.residual_ratio_within_threshold(3.0).unwrap()
287 );
288 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
289 od_sol.nis_consistency(None)?.log();
290
291 od_sol.to_parquet(
292 output_folder.join("04_lro_od_results.parquet"),
293 ExportCfg::default(),
294 )?;
295
296 // Create the ephemeris
297 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
298 let ephem_start = ephem.start_epoch().unwrap();
299 let ephem_end = ephem.end_epoch().unwrap();
300 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
301 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
302 ephem
303 .covar_at(
304 epoch,
305 anise::ephemerides::ephemeris::LocalFrame::RIC,
306 &almanac,
307 )
308 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
309 }
310 // Export as BSP!
311 ephem
312 .write_spice_bsp(
313 -85,
314 output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
315 None,
316 )
317 .expect("could not built BSP");
318 let new_almanac = Almanac::default()
319 .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
320 .unwrap();
321 new_almanac.describe(None, None, None, None, None, None, None, None);
322 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
323
324 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
325 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
326
327 // In our case, we have the truth trajectory from NASA.
328 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
329 // Export the OD trajectory first.
330 let od_trajectory = od_sol.to_traj()?;
331 // Build the RIC difference.
332 od_trajectory.ric_diff_to_parquet(
333 &traj_as_flown,
334 output_folder.join("04_lro_od_truth_error.parquet"),
335 ExportCfg::default(),
336 )?;
337
338 Ok(())
339}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".