06_seq_od/
main.rs

1// #![doc = include_str!("./README.md")]
2extern crate log;
3extern crate nyx_space as nyx;
4extern crate pretty_env_logger as pel;
5
6use anise::{
7    almanac::metaload::MetaFile,
8    constants::{
9        celestial_objects::{EARTH, JUPITER_BARYCENTER, SUN},
10        frames::{MOON_J2000, MOON_PA_FRAME},
11    },
12};
13use hifitime::{Epoch, Unit};
14use nyx::{
15    Orbit, Spacecraft,
16    cosmic::{Mass, MetaAlmanac, SRPData},
17    dynamics::{
18        GravityField, OrbitalDynamics, SolarPressure, SpacecraftDynamics, guidance::LocalFrame,
19    },
20    io::{ConfigRepr, ExportCfg},
21    md::prelude::GravityFieldData,
22    od::{
23        GroundStation, SpacecraftKalmanScalarOD,
24        prelude::{KalmanVariant, TrackingArcSim, TrkConfig},
25        process::{Estimate, NavSolution, SigmaRejection, SpacecraftUncertainty},
26        snc::ProcessNoise3D,
27    },
28    propagators::{IntegratorOptions, Propagator},
29};
30
31use std::{collections::BTreeMap, error::Error, path::PathBuf, sync::Arc};
32
33// TODO: Convert this to a Spacecraft Sequence
34
35fn main() -> Result<(), Box<dyn Error>> {
36    pel::init();
37
38    // ====================== //
39    // === ALMANAC SET UP === //
40    // ====================== //
41
42    // Dynamics models require planetary constants and ephemerides to be defined.
43    // Let's start by grabbing those by using ANISE's MetaAlmanac.
44
45    let data_folder: PathBuf = [
46        env!("CARGO_MANIFEST_DIR"),
47        "examples",
48        "06_lunar_orbit_determination",
49    ]
50    .iter()
51    .collect();
52
53    let meta = data_folder.join("metaalmanac.dhall");
54
55    // Load this ephem in the general Almanac we're using for this analysis.
56    let almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
57        .map_err(Box::new)?
58        .process(true)
59        .map_err(Box::new)?;
60
61    // Lock the almanac (an Arc is a read only structure).
62    let almanac = Arc::new(almanac);
63
64    // Build a nominal trajectory
65    // TODO: Switch this to a sequence once the OD over a spacecraft sequence is implemented.
66
67    let epoch = Epoch::from_gregorian_utc_at_noon(2024, 2, 29);
68    let moon_j2000 = almanac.frame_info(MOON_J2000)?;
69
70    // To build the trajectory we need to provide a spacecraft template.
71    let orbiter = Spacecraft::builder()
72        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0))
73        .srp(SRPData {
74            area_m2: 3.9 * 2.7,
75            coeff_reflectivity: 0.96,
76        })
77        .orbit(Orbit::try_keplerian_altitude(
78            150.0, 0.00212, 33.6, 45.0, 45.0, 0.0, epoch, moon_j2000,
79        )?) // Setting a zero orbit here because it's just a template
80        .build();
81
82    // ========================== //
83    // === BUILD NOMINAL TRAJ === //
84    // ========================== //
85
86    // Set up the spacecraft dynamics.
87
88    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
89    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
90    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
91
92    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
93    // We're using the GRAIL JGGRX model.
94    let mut jggrx_meta = MetaFile {
95        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
96        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
97    };
98    // And let's download it if we don't have it yet.
99    jggrx_meta.process(true)?;
100
101    // Build the spherical harmonics.
102    // The harmonics must be computed in the body fixed frame.
103    // We're using the long term prediction of the Moon principal axes frame.
104    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
105    let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
106        &jggrx_meta.uri,
107        80,
108        80,
109        true,
110        almanac.frame_info(moon_pa_frame)?,
111    )?);
112
113    // Include the spherical harmonics into the orbital dynamics.
114    orbital_dyn.accel_models.push(sph_harmonics);
115
116    // We define the solar radiation pressure, using the default solar flux and accounting only
117    // for the eclipsing caused by the Earth and Moon.
118    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
119    let srp_dyn = SolarPressure::new(vec![MOON_J2000], &almanac)?;
120
121    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
122    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
123    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
124
125    println!("{dynamics}");
126
127    let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::default());
128
129    let truth_traj = setup
130        .with(orbiter, almanac.clone())
131        .for_duration_with_traj(Unit::Day * 2)?
132        .1;
133
134    // ==================== //
135    // === OD SIMULATOR === //
136    // ==================== //
137
138    // Load the Deep Space Network ground stations.
139    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
140    let ground_station_file = data_folder.join("dsn-network.yaml");
141    let devices = GroundStation::load_named(ground_station_file)?;
142
143    let proc_devices = devices.clone();
144
145    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
146    // Nyx can build a tracking schedule for you based on the first station with access.
147    let configs: BTreeMap<String, TrkConfig> =
148        TrkConfig::load_named(data_folder.join("tracking-cfg.yaml"))?;
149
150    // Build the tracking arc simulation to generate a "standard measurement".
151    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
152        devices.clone(),
153        truth_traj.clone(),
154        configs,
155        123, // Set a seed for reproducibility
156    )?;
157
158    trk.build_schedule(&almanac)?;
159    let arc = trk.generate_measurements(&almanac)?;
160    // Save the simulated tracking data
161    arc.to_parquet_simple("./data/04_output/06_lunar_simulated_tracking.parquet")?;
162
163    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
164    println!("{arc}");
165
166    // Now that we have simulated measurements, we'll run the orbit determination.
167
168    // ===================== //
169    // === OD ESTIMATION === //
170    // ===================== //
171
172    let sc = SpacecraftUncertainty::builder()
173        .nominal(orbiter)
174        .frame(LocalFrame::RIC)
175        .x_km(0.5)
176        .y_km(0.5)
177        .z_km(0.5)
178        .vx_km_s(5e-3)
179        .vy_km_s(5e-3)
180        .vz_km_s(5e-3)
181        .build();
182
183    // Build the filter initial estimate, which we will reuse in the filter.
184    let initial_estimate = sc.to_estimate()?;
185
186    println!("== FILTER STATE ==\n{orbiter:x}\n{initial_estimate}");
187
188    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
189    let process_noise = ProcessNoise3D::from_velocity_km_s(
190        &[1e-14, 1e-14, 1e-14],
191        1 * Unit::Hour,
192        10 * Unit::Minute,
193        None,
194    );
195
196    println!("{process_noise}");
197
198    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
199    let odp = SpacecraftKalmanScalarOD::new(
200        setup,
201        KalmanVariant::ReferenceUpdate,
202        Some(SigmaRejection::default()),
203        proc_devices,
204        almanac.clone(),
205    )
206    .with_process_noise(process_noise);
207
208    let od_sol = odp.process_arc(initial_estimate, &arc)?;
209
210    let final_est = od_sol.estimates.last().unwrap();
211
212    println!("{final_est}");
213
214    let ric_err = truth_traj
215        .at(final_est.epoch())?
216        .orbit
217        .ric_difference(&final_est.orbital_state())?;
218    println!("== RIC at end ==");
219    println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
220    println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
221
222    println!(
223        "Num residuals rejected: #{}",
224        od_sol.rejected_residuals().len()
225    );
226    println!(
227        "Percentage within +/-3: {}",
228        od_sol.residual_ratio_within_threshold(3.0).unwrap()
229    );
230    println!("Whitened residuals normal? {}", od_sol.is_normal(None)?);
231    println!("NIS test success? {}", od_sol.is_nis_consistent(None)?);
232
233    od_sol.to_parquet(
234        "./data/04_output/06_lunar_od_results.parquet",
235        ExportCfg::default(),
236    )?;
237
238    let od_trajectory = od_sol.to_traj()?;
239    // Build the RIC difference.
240    od_trajectory.ric_diff_to_parquet(
241        &truth_traj,
242        "./data/04_output/06_lunar_od_truth_error.parquet",
243        ExportCfg::default(),
244    )?;
245
246    Ok(())
247}
06_seq_od/main.rs

06_seq_od/
main.rs
