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