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