SpacecraftUncertainty

Struct SpacecraftUncertainty 

Source
pub struct SpacecraftUncertainty {
    pub nominal: Spacecraft,
    pub frame: Option<LocalFrame>,
    pub x_km: f64,
    pub y_km: f64,
    pub z_km: f64,
    pub vx_km_s: f64,
    pub vy_km_s: f64,
    pub vz_km_s: f64,
    pub coeff_reflectivity: f64,
    pub coeff_drag: f64,
    pub mass_kg: f64,
}
Expand description

Builds a spacecraft uncertainty in different local frames, dispersing any of the parameters of the spacecraft state.

§Usage

Use the TypeBuilder trait, e.g SpacecraftUncertainty::builder().nominal(spacecraft).frame(LocalFrame::RIC).x_km(0.5).y_km(0.5).z_km(0.5).build() to build an uncertainty on position in the RIC frame of 500 meters on R, I, and C, and zero on all other parameters (velocity components, Cr, Cd, mass).

Fields§

§nominal: Spacecraft§frame: Option<LocalFrame>§x_km: f64§y_km: f64§z_km: f64§vx_km_s: f64§vy_km_s: f64§vz_km_s: f64§coeff_reflectivity: f64§coeff_drag: f64§mass_kg: f64

Implementations§

Source§

impl SpacecraftUncertainty

Source

pub fn builder() -> SpacecraftUncertaintyBuilder<((), (), (), (), (), (), (), (), (), (), ())>

Create a builder for building SpacecraftUncertainty. On the builder, call .nominal(...), .frame(...)(optional), .x_km(...)(optional), .y_km(...)(optional), .z_km(...)(optional), .vx_km_s(...)(optional), .vy_km_s(...)(optional), .vz_km_s(...)(optional), .coeff_reflectivity(...)(optional), .coeff_drag(...)(optional), .mass_kg(...)(optional) to set the values of the fields. Finally, call .build() to create the instance of SpacecraftUncertainty.

Examples found in repository?
examples/02_jwst_covar_monte_carlo/main.rs (line 77)
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
Hide additional examples
examples/05_cislunar_spacecraft_link_od/main.rs (line 184)
34fn main() -> Result<(), Box<dyn Error>> {
35    pel::init();
36
37    // ====================== //
38    // === ALMANAC SET UP === //
39    // ====================== //
40
41    let manifest_dir =
42        PathBuf::from(std::env::var("CARGO_MANIFEST_DIR").unwrap_or(".".to_string()));
43
44    let out = manifest_dir.join("data/04_output/");
45
46    let almanac = Arc::new(
47        Almanac::new(
48            &manifest_dir
49                .join("data/01_planetary/pck08.pca")
50                .to_string_lossy(),
51        )
52        .unwrap()
53        .load(
54            &manifest_dir
55                .join("data/01_planetary/de440s.bsp")
56                .to_string_lossy(),
57        )
58        .unwrap(),
59    );
60
61    let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap();
62    let moon_iau = almanac.frame_from_uid(IAU_MOON_FRAME).unwrap();
63
64    let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
65    let nrho = Orbit::cartesian(
66        166_473.631_302_239_7,
67        -274_715.487_253_382_7,
68        -211_233.210_176_686_7,
69        0.933_451_604_520_018_4,
70        0.436_775_046_841_900_9,
71        -0.082_211_021_250_348_95,
72        epoch,
73        eme2k,
74    );
75
76    let tx_nrho_sc = Spacecraft::from(nrho);
77
78    let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
79    println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
80
81    let bodies = vec![EARTH, SUN];
82    let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
83
84    let setup = Propagator::rk89(
85        dynamics,
86        IntegratorOptions::builder().max_step(0.5.minutes()).build(),
87    );
88
89    /* == Propagate the NRHO vehicle == */
90    let prop_time = 1.1 * state_luna.period().unwrap();
91
92    let (nrho_final, mut tx_traj) = setup
93        .with(tx_nrho_sc, almanac.clone())
94        .for_duration_with_traj(prop_time)
95        .unwrap();
96
97    tx_traj.name = Some("NRHO Tx SC".to_string());
98
99    println!("{tx_traj}");
100
101    /* == Propagate an LLO vehicle == */
102    let llo_orbit =
103        Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
104
105    let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
106
107    let (_, llo_traj) = setup
108        .with(llo_sc, almanac.clone())
109        .until_epoch_with_traj(nrho_final.epoch())
110        .unwrap();
111
112    // Export the subset of the first two hours.
113    llo_traj
114        .clone()
115        .filter_by_offset(..2.hours())
116        .to_parquet_simple(out.join("05_caps_llo_truth.pq"), almanac.clone())?;
117
118    /* == Setup the interlink == */
119
120    let mut measurement_types = IndexSet::new();
121    measurement_types.insert(MeasurementType::Range);
122    measurement_types.insert(MeasurementType::Doppler);
123
124    let mut stochastics = IndexMap::new();
125
126    let sa45_csac_allan_dev = 1e-11;
127
128    stochastics.insert(
129        MeasurementType::Range,
130        StochasticNoise::from_hardware_range_km(
131            sa45_csac_allan_dev,
132            10.0.seconds(),
133            link_specific::ChipRate::StandardT4B,
134            link_specific::SN0::Average,
135        ),
136    );
137
138    stochastics.insert(
139        MeasurementType::Doppler,
140        StochasticNoise::from_hardware_doppler_km_s(
141            sa45_csac_allan_dev,
142            10.0.seconds(),
143            link_specific::CarrierFreq::SBand,
144            link_specific::CN0::Average,
145        ),
146    );
147
148    let interlink = InterlinkTxSpacecraft {
149        traj: tx_traj,
150        measurement_types,
151        integration_time: None,
152        timestamp_noise_s: None,
153        ab_corr: Aberration::LT,
154        stochastic_noises: Some(stochastics),
155    };
156
157    // Devices are the transmitter, which is our NRHO vehicle.
158    let mut devices = BTreeMap::new();
159    devices.insert("NRHO Tx SC".to_string(), interlink);
160
161    let mut configs = BTreeMap::new();
162    configs.insert(
163        "NRHO Tx SC".to_string(),
164        TrkConfig::builder()
165            .strands(vec![Strand {
166                start: epoch,
167                end: nrho_final.epoch(),
168            }])
169            .build(),
170    );
171
172    let mut trk_sim =
173        TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
174    println!("{trk_sim}");
175
176    let trk_data = trk_sim.generate_measurements(almanac.clone()).unwrap();
177    println!("{trk_data}");
178
179    trk_data
180        .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
181        .unwrap();
182
183    // Run a truth OD where we estimate the LLO position
184    let llo_uncertainty = SpacecraftUncertainty::builder()
185        .nominal(llo_sc)
186        .x_km(1.0)
187        .y_km(1.0)
188        .z_km(1.0)
189        .vx_km_s(1e-3)
190        .vy_km_s(1e-3)
191        .vz_km_s(1e-3)
192        .build();
193
194    let mut proc_devices = devices.clone();
195
196    // Define the initial estimate, randomized, seed for reproducibility
197    let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
198    // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
199    initial_estimate.covar *= 2.5;
200
201    // Increase the noise in the devices to accept more measurements.
202
203    for link in proc_devices.values_mut() {
204        for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
205            *noise.white_noise.as_mut().unwrap() *= 3.0;
206        }
207    }
208
209    let init_err = initial_estimate
210        .orbital_state()
211        .ric_difference(&llo_orbit)
212        .unwrap();
213
214    println!("initial estimate:\n{initial_estimate}");
215    println!("RIC errors = {init_err}",);
216
217    let odp = InterlinkKalmanOD::new(
218        setup.clone(),
219        KalmanVariant::ReferenceUpdate,
220        Some(ResidRejectCrit::default()),
221        proc_devices,
222        almanac.clone(),
223    );
224
225    // Shrink the data to process.
226    let arc = trk_data.filter_by_offset(..2.hours());
227
228    let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
229
230    println!("{od_sol}");
231
232    od_sol
233        .to_parquet(
234            out.join(format!("05_caps_interlink_od_sol.pq")),
235            ExportCfg::default(),
236        )
237        .unwrap();
238
239    let od_traj = od_sol.to_traj().unwrap();
240
241    od_traj
242        .ric_diff_to_parquet(
243            &llo_traj,
244            out.join(format!("05_caps_interlink_llo_est_error.pq")),
245            ExportCfg::default(),
246        )
247        .unwrap();
248
249    let final_est = od_sol.estimates.last().unwrap();
250    assert!(final_est.within_3sigma(), "should be within 3 sigma");
251
252    println!("ESTIMATE\n{final_est:x}\n");
253    let truth = llo_traj.at(final_est.epoch()).unwrap();
254    println!("TRUTH\n{truth:x}");
255
256    let final_err = truth
257        .orbit
258        .ric_difference(&final_est.orbital_state())
259        .unwrap();
260    println!("ERROR {final_err}");
261
262    // Build the residuals versus reference plot.
263    let rvr_sol = odp
264        .process_arc(initial_estimate, &arc.resid_vs_ref_check())
265        .unwrap();
266
267    rvr_sol
268        .to_parquet(
269            out.join(format!("05_caps_interlink_resid_v_ref.pq")),
270            ExportCfg::default(),
271        )
272        .unwrap();
273
274    let final_rvr = rvr_sol.estimates.last().unwrap();
275
276    println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
277    println!(
278        "Pure prop error {:.3} m",
279        final_rvr
280            .orbital_state()
281            .ric_difference(&final_est.orbital_state())
282            .unwrap()
283            .rmag_km()
284            * 1e3
285    );
286
287    Ok(())
288}
examples/04_lro_od/main.rs (line 250)
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}
Source§

impl SpacecraftUncertainty

Source

pub fn to_estimate(&self) -> PhysicsResult<KfEstimate<Spacecraft>>

Builds a Kalman filter estimate for a spacecraft state, ready to ingest into an OD Process.

Note: this function will rotate from the provided local frame into the inertial frame with the same central body.

Examples found in repository?
examples/02_jwst_covar_monte_carlo/main.rs (line 93)
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
Hide additional examples
examples/04_lro_od/main.rs (line 262)
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}
Source

pub fn to_estimate_randomized( &self, seed: Option<u128>, ) -> PhysicsResult<KfEstimate<Spacecraft>>

Returns an estimation whose nominal state is dispersed. Known bug in MultiVariate dispersion: https://github.com/nyx-space/nyx/issues/339

Examples found in repository?
examples/05_cislunar_spacecraft_link_od/main.rs (line 197)
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}

Trait Implementations§

Source§

impl Clone for SpacecraftUncertainty

Source§

fn clone(&self) -> SpacecraftUncertainty

Returns a duplicate of the value. Read more
1.0.0 · Source§

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source. Read more
Source§

impl Debug for SpacecraftUncertainty

Source§

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter. Read more
Source§

impl Display for SpacecraftUncertainty

Source§

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter. Read more
Source§

impl Copy for SpacecraftUncertainty

Auto Trait Implementations§

Blanket Implementations§

Source§

impl<T> Any for T
where T: 'static + ?Sized,

Source§

fn type_id(&self) -> TypeId

Gets the TypeId of self. Read more
Source§

impl<T> Borrow<T> for T
where T: ?Sized,

Source§

fn borrow(&self) -> &T

Immutably borrows from an owned value. Read more
Source§

impl<T> BorrowMut<T> for T
where T: ?Sized,

Source§

fn borrow_mut(&mut self) -> &mut T

Mutably borrows from an owned value. Read more
Source§

impl<T> CloneToUninit for T
where T: Clone,

Source§

unsafe fn clone_to_uninit(&self, dest: *mut u8)

🔬This is a nightly-only experimental API. (clone_to_uninit)
Performs copy-assignment from self to dest. Read more
Source§

impl<T> From<T> for T

Source§

fn from(t: T) -> T

Returns the argument unchanged.

§

impl<T> Instrument for T

§

fn instrument(self, span: Span) -> Instrumented<Self>

Instruments this type with the provided [Span], returning an Instrumented wrapper. Read more
§

fn in_current_span(self) -> Instrumented<Self>

Instruments this type with the current Span, returning an Instrumented wrapper. Read more
Source§

impl<T, U> Into<U> for T
where U: From<T>,

Source§

fn into(self) -> U

Calls U::from(self).

That is, this conversion is whatever the implementation of From<T> for U chooses to do.

Source§

impl<T> IntoEither for T

Source§

fn into_either(self, into_left: bool) -> Either<Self, Self>

Converts 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 more
Source§

fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
where F: FnOnce(&Self) -> bool,

Converts 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

§

const ALIGN: usize

The alignment of pointer.
§

type Init = T

The type for initializers.
§

unsafe fn init(init: <T as Pointable>::Init) -> usize

Initializes a with the given initializer. Read more
§

unsafe fn deref<'a>(ptr: usize) -> &'a T

Dereferences the given pointer. Read more
§

unsafe fn deref_mut<'a>(ptr: usize) -> &'a mut T

Mutably dereferences the given pointer. Read more
§

unsafe fn drop(ptr: usize)

Drops the object pointed to by the given pointer. Read more
Source§

impl<T> Same for T

Source§

type Output = T

Should always be Self
§

impl<SS, SP> SupersetOf<SS> for SP
where SS: SubsetOf<SP>,

§

fn to_subset(&self) -> Option<SS>

The inverse inclusion map: attempts to construct self from the equivalent element of its superset. Read more
§

fn is_in_subset(&self) -> bool

Checks if self is actually part of its subset T (and can be converted to it).
§

fn to_subset_unchecked(&self) -> SS

Use with care! Same as self.to_subset but without any property checks. Always succeeds.
§

fn from_subset(element: &SS) -> SP

The inclusion map: converts self to the equivalent element of its superset.
Source§

impl<T> ToOwned for T
where T: Clone,

Source§

type Owned = T

The resulting type after obtaining ownership.
Source§

fn to_owned(&self) -> T

Creates owned data from borrowed data, usually by cloning. Read more
Source§

fn clone_into(&self, target: &mut T)

Uses borrowed data to replace owned data, usually by cloning. Read more
Source§

impl<T> ToString for T
where T: Display + ?Sized,

Source§

fn to_string(&self) -> String

Converts the given value to a String. Read more
Source§

impl<T, U> TryFrom<U> for T
where U: Into<T>,

Source§

type Error = Infallible

The type returned in the event of a conversion error.
Source§

fn try_from(value: U) -> Result<T, <T as TryFrom<U>>::Error>

Performs the conversion.
Source§

impl<T, U> TryInto<U> for T
where U: TryFrom<T>,

Source§

type Error = <U as TryFrom<T>>::Error

The type returned in the event of a conversion error.
Source§

fn try_into(self) -> Result<U, <U as TryFrom<T>>::Error>

Performs the conversion.
§

impl<V, T> VZip<V> for T
where V: MultiLane<T>,

§

fn vzip(self) -> V

§

impl<T> WithSubscriber for T

§

fn with_subscriber<S>(self, subscriber: S) -> WithDispatch<Self>
where S: Into<Dispatch>,

Attaches the provided Subscriber to this type, returning a [WithDispatch] wrapper. Read more
§

fn with_current_subscriber(self) -> WithDispatch<Self>

Attaches the current default Subscriber to this type, returning a [WithDispatch] wrapper. Read more
§

impl<T> Allocation for T
where T: RefUnwindSafe + Send + Sync,

§

impl<T> ErasedDestructor for T
where T: 'static,