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.clone()).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(ResidRejectCrit::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.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 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.clone()).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(ResidRejectCrit::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/raise_optim.rs (line 189)
135fn evaluate_weights(
136 weights: &[f32],
137 prop_time_days: f64,
138 state: Arc<SharedState>,
139) -> Result<(f64, f64), Box<dyn Error>> {
140 let ηthresholds: Vec<f64> = weights.iter().map(|w| *w as f64).collect();
141
142 let eme2k = state.almanac.frame_info(EARTH_J2000).unwrap();
143 let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
144
145 let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);
146
147 let sc = Spacecraft::builder()
148 .orbit(orbit)
149 .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0))
150 .srp(SRPData::from_area(3.0 * 6.0))
151 .thruster(Thruster {
152 isp_s: 4435.0,
153 thrust_N: 0.472,
154 })
155 .mode(GuidanceMode::Thrust)
156 .build();
157
158 let prop_time = prop_time_days * Unit::Day;
159
160 let objectives = &[
161 Objective::within_tolerance(
162 StateParameter::Element(OrbitalElement::SemiMajorAxis),
163 30_000.0,
164 20.0,
165 ),
166 Objective::within_tolerance(
167 StateParameter::Element(OrbitalElement::Eccentricity),
168 0.001,
169 5e-5,
170 ),
171 Objective::within_tolerance(
172 StateParameter::Element(OrbitalElement::Inclination),
173 0.05,
174 1e-2,
175 ),
176 ];
177
178 let ctrl = Ruggiero::from_ηthresholds(objectives, &ηthresholds, sc)?;
179
180 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
181 orbital_dyn.accel_models.push(state.harmonics.clone());
182
183 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, state.srp_dyn.clone())
184 .with_guidance_law(ctrl.clone());
185
186 let (final_state, _traj) = Propagator::rk89(
187 sc_dynamics.clone(),
188 IntegratorOptions::builder()
189 .min_step(10.0_f64.seconds())
190 .tolerance(1e-8)
191 .error_ctrl(ErrorControl::RSSCartesianStep)
192 .build(),
193 )
194 .with(sc, state.almanac.clone())
195 .for_duration_with_traj(prop_time)?;
196
197 let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
198
199 let mut penalty = 0.0;
200 for obj in objectives {
201 let (achieved, error) = obj.assess(&final_state)?;
202 if !achieved {
203 penalty += error.abs();
204 }
205 info!("{obj} error: {error:.3}, achieved? {achieved}");
206 }
207
208 info!("{ηthresholds:?} -> {prop_usage:.3} kg\tpenalty = {penalty:.3}");
209
210 Ok((prop_usage, penalty * 1000.0))
211}More examples
nyx-core/examples/03_geo_analysis/stationkeeping.rs (line 121)
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::from_stor(
86 almanac.frame_info(IAU_EARTH_FRAME)?,
87 GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true)?,
88 );
89 orbital_dyn.accel_models.push(harmonics);
90
91 let srp_dyn = SolarPressure::default_flux(EARTH_J2000, almanac.clone())?;
92 let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
93 .with_guidance_law(ruggiero_ctrl.clone());
94
95 println!("{sc_dynamics}");
96
97 // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
98
99 // Let's start by defining the dispersion.
100 // 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.
101 // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
102 let mc_rv = MvnSpacecraft::new(
103 sc,
104 vec![StateDispersion::zero_mean(
105 StateParameter::Element(OrbitalElement::SemiMajorAxis),
106 3.0,
107 )],
108 )?;
109
110 let my_mc = MonteCarlo::new(
111 sc, // Nominal state
112 mc_rv,
113 "03_geo_sk".to_string(), // Scenario name
114 None, // No specific seed specified, so one will be drawn from the computer's entropy.
115 );
116
117 // Build the propagator setup.
118 let setup = Propagator::rk89(
119 sc_dynamics.clone(),
120 IntegratorOptions::builder()
121 .min_step(10.0_f64.seconds())
122 .error_ctrl(ErrorControl::RSSCartesianStep)
123 .build(),
124 );
125
126 let num_runs = 25;
127 let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
128
129 assert_eq!(rslts.runs.len(), num_runs);
130
131 rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;
132
133 Ok(())
134}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::from_stor(
110 almanac.frame_info(IAU_EARTH_FRAME)?,
111 GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true).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.clone())?;
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.clone()).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(ResidRejectCrit::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 100)
35fn main() -> Result<(), Box<dyn Error>> {
36 pel::init();
37
38 // ====================== //
39 // === ALMANAC SET UP === //
40 // ====================== //
41
42 // Dynamics models require planetary constants and ephemerides to be defined.
43 // Let's start by grabbing those by using ANISE's MetaAlmanac.
44
45 let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
46 .iter()
47 .collect();
48
49 let meta = data_folder.join("lro-dynamics.dhall");
50
51 // Load this ephem in the general Almanac we're using for this analysis.
52 let mut almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
53 .map_err(Box::new)?
54 .process(true)
55 .map_err(Box::new)?;
56
57 let mut moon_pc = almanac.get_planetary_data_from_id(MOON).unwrap();
58 moon_pc.mu_km3_s2 = 4902.74987;
59 almanac.set_planetary_data_from_id(MOON, moon_pc).unwrap();
60
61 let mut earth = almanac.get_planetary_data_from_id(EARTH).unwrap();
62 earth.mu_km3_s2 = 398600.436;
63 almanac.set_planetary_data_from_id(EARTH, earth).unwrap();
64
65 // Save this new kernel for reuse.
66 // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
67 almanac
68 .planetary_data
69 .values()
70 .next()
71 .unwrap()
72 .save_as(&data_folder.join("lro-specific.pca"), true)?;
73
74 // Lock the almanac (an Arc is a read only structure).
75 let almanac = Arc::new(almanac);
76
77 // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
78 // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
79 // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
80 // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
81 let lro_frame = Frame::from_ephem_j2000(-85);
82
83 // To build the trajectory we need to provide a spacecraft template.
84 let sc_template = Spacecraft::builder()
85 .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
86 .srp(SRPData {
87 // SRP configuration is arbitrary, but we will be estimating it anyway.
88 area_m2: 3.9 * 2.7,
89 coeff_reflectivity: 0.96,
90 })
91 .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
92 .build();
93 // Now we can build the trajectory from the BSP file.
94 // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
95 let traj_as_flown = Traj::from_bsp(
96 lro_frame,
97 MOON_J2000,
98 almanac.clone(),
99 sc_template,
100 5.seconds(),
101 Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
102 Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
103 Aberration::LT,
104 Some("LRO".to_string()),
105 )?;
106
107 println!("{traj_as_flown}");
108
109 // ====================== //
110 // === MODEL MATCHING === //
111 // ====================== //
112
113 // Set up the spacecraft dynamics.
114
115 // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
116 // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
117 let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
118
119 // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
120 // We're using the GRAIL JGGRX model.
121 let mut jggrx_meta = MetaFile {
122 uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
123 crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
124 };
125 // And let's download it if we don't have it yet.
126 jggrx_meta.process(true)?;
127
128 // Build the spherical harmonics.
129 // The harmonics must be computed in the body fixed frame.
130 // We're using the long term prediction of the Moon principal axes frame.
131 let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
132 let sph_harmonics = GravityField::from_stor(
133 almanac.frame_info(moon_pa_frame)?,
134 GravityFieldData::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
135 );
136
137 // Include the spherical harmonics into the orbital dynamics.
138 orbital_dyn.accel_models.push(sph_harmonics);
139
140 // We define the solar radiation pressure, using the default solar flux and accounting only
141 // for the eclipsing caused by the Earth and Moon.
142 // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
143 let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
144
145 // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
146 // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
147 let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
148
149 println!("{dynamics}");
150
151 // Now we can build the propagator.
152 let setup = Propagator::default_dp78(dynamics.clone());
153
154 // For reference, let's build the trajectory with Nyx's models from that LRO state.
155 let (sim_final, traj_as_sim) = setup
156 .with(*traj_as_flown.first(), almanac.clone())
157 .until_epoch_with_traj(traj_as_flown.last().epoch())?;
158
159 println!("SIM INIT: {:x}", traj_as_flown.first());
160 println!("SIM FINAL: {sim_final:x}");
161 // Compute RIC difference between SIM and LRO ephem
162 let sim_lro_delta = sim_final
163 .orbit
164 .ric_difference(&traj_as_flown.last().orbit)?;
165 println!("{traj_as_sim}");
166 println!(
167 "SIM v LRO - RIC Position (m): {:.3}",
168 sim_lro_delta.radius_km * 1e3
169 );
170 println!(
171 "SIM v LRO - RIC Velocity (m/s): {:.3}",
172 sim_lro_delta.velocity_km_s * 1e3
173 );
174
175 traj_as_sim.ric_diff_to_parquet(
176 &traj_as_flown,
177 "./data/04_output/04_lro_sim_truth_error.parquet",
178 ExportCfg::default(),
179 )?;
180
181 // ==================== //
182 // === OD SIMULATOR === //
183 // ==================== //
184
185 // 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
186 // and the truth LRO state.
187
188 // Therefore, we will actually run an estimation from a dispersed LRO state.
189 // The sc_seed is the true LRO state from the BSP.
190 let sc_seed = *traj_as_flown.first();
191
192 // Load the Deep Space Network ground stations.
193 // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
194 let ground_station_file: PathBuf = [
195 env!("CARGO_MANIFEST_DIR"),
196 "examples",
197 "04_lro_od",
198 "dsn-network.yaml",
199 ]
200 .iter()
201 .collect();
202
203 let devices = GroundStation::load_named(ground_station_file)?;
204
205 let mut proc_devices = devices.clone();
206
207 // Increase the noise in the devices to accept more measurements.
208 for gs in proc_devices.values_mut() {
209 if let Some(noise) = &mut gs
210 .stochastic_noises
211 .as_mut()
212 .unwrap()
213 .get_mut(&MeasurementType::Range)
214 {
215 *noise.white_noise.as_mut().unwrap() *= 3.0;
216 }
217 }
218
219 // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
220 // Nyx can build a tracking schedule for you based on the first station with access.
221 let trkconfg_yaml: PathBuf = [
222 env!("CARGO_MANIFEST_DIR"),
223 "examples",
224 "04_lro_od",
225 "tracking-cfg.yaml",
226 ]
227 .iter()
228 .collect();
229
230 let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
231
232 // Build the tracking arc simulation to generate a "standard measurement".
233 let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
234 devices.clone(),
235 traj_as_flown.clone(),
236 configs,
237 123, // Set a seed for reproducibility
238 )?;
239
240 trk.build_schedule(almanac.clone())?;
241 let arc = trk.generate_measurements(almanac.clone())?;
242 // Save the simulated tracking data
243 arc.to_parquet_simple("./data/04_output/04_lro_simulated_tracking.parquet")?;
244
245 // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
246 println!("{arc}");
247
248 // Now that we have simulated measurements, we'll run the orbit determination.
249
250 // ===================== //
251 // === OD ESTIMATION === //
252 // ===================== //
253
254 let sc = SpacecraftUncertainty::builder()
255 .nominal(sc_seed)
256 .frame(LocalFrame::RIC)
257 .x_km(0.5)
258 .y_km(0.5)
259 .z_km(0.5)
260 .vx_km_s(5e-3)
261 .vy_km_s(5e-3)
262 .vz_km_s(5e-3)
263 .build();
264
265 // Build the filter initial estimate, which we will reuse in the filter.
266 let mut initial_estimate = sc.to_estimate()?;
267 initial_estimate.covar *= 3.0;
268
269 println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
270
271 // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
272 let process_noise = ProcessNoise3D::from_velocity_km_s(
273 &[1e-12, 1e-12, 1e-12],
274 1 * Unit::Hour,
275 10 * Unit::Minute,
276 None,
277 );
278
279 println!("{process_noise}");
280
281 // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
282 let odp = SpacecraftKalmanOD::new(
283 setup,
284 KalmanVariant::ReferenceUpdate,
285 Some(ResidRejectCrit::default()),
286 proc_devices,
287 almanac.clone(),
288 )
289 .with_process_noise(process_noise);
290
291 let od_sol = odp.process_arc(initial_estimate, &arc)?;
292
293 let final_est = od_sol.estimates.last().unwrap();
294
295 println!("{final_est}");
296
297 let ric_err = traj_as_flown
298 .at(final_est.epoch())?
299 .orbit
300 .ric_difference(&final_est.orbital_state())?;
301 println!("== RIC at end ==");
302 println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
303 println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
304
305 println!(
306 "Num residuals rejected: #{}",
307 od_sol.rejected_residuals().len()
308 );
309 println!(
310 "Percentage within +/-3: {}",
311 od_sol.residual_ratio_within_threshold(3.0).unwrap()
312 );
313 println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
314
315 od_sol.to_parquet(
316 "./data/04_output/04_lro_od_results.parquet",
317 ExportCfg::default(),
318 )?;
319
320 // Create the ephemeris
321 let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
322 let ephem_start = ephem.start_epoch().unwrap();
323 let ephem_end = ephem.end_epoch().unwrap();
324 // Check that the covariance is PSD throughout the ephemeris by interpolating it.
325 for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
326 ephem
327 .covar_at(
328 epoch,
329 anise::ephemerides::ephemeris::LocalFrame::RIC,
330 &almanac,
331 )
332 .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
333 }
334 // Export as BSP!
335 ephem
336 .write_spice_bsp(-85, "./data/04_output/04_lro_rebuilt.bsp", None)
337 .expect("could not built BSP");
338 let new_almanac = Almanac::default()
339 .load("./data/04_output/04_lro_rebuilt.bsp")
340 .unwrap();
341 new_almanac.describe(None, None, None, None, None, None, None, None);
342 let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
343
344 assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
345 assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
346
347 // In our case, we have the truth trajectory from NASA.
348 // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
349 // Export the OD trajectory first.
350 let od_trajectory = od_sol.to_traj()?;
351 // Build the RIC difference.
352 od_trajectory.ric_diff_to_parquet(
353 &traj_as_flown,
354 "./data/04_output/04_lro_od_truth_error.parquet",
355 ExportCfg::default(),
356 )?;
357
358 Ok(())
359}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".