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        almanac.frame_info(moon_pa_frame)?,
110    )?);
111
112    // Include the spherical harmonics into the orbital dynamics.
113    orbital_dyn.accel_models.push(sph_harmonics);
114
115    // We define the solar radiation pressure, using the default solar flux and accounting only
116    // for the eclipsing caused by the Earth and Moon.
117    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
118    let srp_dyn = SolarPressure::new(vec![MOON_J2000], &almanac)?;
119
120    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
121    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
122    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
123
124    println!("{dynamics}");
125
126    let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::default());
127
128    let truth_traj = setup
129        .with(orbiter, almanac.clone())
130        .for_duration_with_traj(Unit::Day * 2)?
131        .1;
132
133    // ==================== //
134    // === OD SIMULATOR === //
135    // ==================== //
136
137    // Load the Deep Space Network ground stations.
138    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
139    let ground_station_file = data_folder.join("dsn-network.yaml");
140    let devices = GroundStation::load_named(ground_station_file)?;
141
142    let proc_devices = devices.clone();
143
144    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
145    // Nyx can build a tracking schedule for you based on the first station with access.
146    let configs: BTreeMap<String, TrkConfig> =
147        TrkConfig::load_named(data_folder.join("tracking-cfg.yaml"))?;
148
149    // Build the tracking arc simulation to generate a "standard measurement".
150    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
151        devices.clone(),
152        truth_traj.clone(),
153        configs,
154        123, // Set a seed for reproducibility
155    )?;
156
157    trk.build_schedule(&almanac)?;
158    let arc = trk.generate_measurements(&almanac)?;
159    // Save the simulated tracking data
160    arc.to_parquet_simple("./data/04_output/06_lunar_simulated_tracking.parquet")?;
161
162    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
163    println!("{arc}");
164
165    // Now that we have simulated measurements, we'll run the orbit determination.
166
167    // ===================== //
168    // === OD ESTIMATION === //
169    // ===================== //
170
171    let sc = SpacecraftUncertainty::builder()
172        .nominal(orbiter)
173        .frame(LocalFrame::RIC)
174        .x_km(0.5)
175        .y_km(0.5)
176        .z_km(0.5)
177        .vx_km_s(5e-3)
178        .vy_km_s(5e-3)
179        .vz_km_s(5e-3)
180        .build();
181
182    // Build the filter initial estimate, which we will reuse in the filter.
183    let initial_estimate = sc.to_estimate()?;
184
185    println!("== FILTER STATE ==\n{orbiter:x}\n{initial_estimate}");
186
187    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
188    let process_noise = ProcessNoise3D::from_velocity_km_s(
189        &[1e-14, 1e-14, 1e-14],
190        1 * Unit::Hour,
191        10 * Unit::Minute,
192        None,
193    );
194
195    println!("{process_noise}");
196
197    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
198    let odp = SpacecraftKalmanScalarOD::new(
199        setup,
200        KalmanVariant::ReferenceUpdate,
201        Some(SigmaRejection::default()),
202        proc_devices,
203        almanac.clone(),
204    )
205    .with_process_noise(process_noise);
206
207    let od_sol = odp.process_arc(initial_estimate, &arc)?;
208
209    let final_est = od_sol.estimates.last().unwrap();
210
211    println!("{final_est}");
212
213    let ric_err = truth_traj
214        .at(final_est.epoch())?
215        .orbit
216        .ric_difference(&final_est.orbital_state())?;
217    println!("== RIC at end ==");
218    println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
219    println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
220
221    println!(
222        "Num residuals rejected: #{}",
223        od_sol.rejected_residuals().len()
224    );
225    println!(
226        "Percentage within +/-3: {}",
227        od_sol.residual_ratio_within_threshold(3.0).unwrap()
228    );
229    println!("Whitened residuals normal? {}", od_sol.is_normal(None)?);
230    println!("NIS consistency: {}", od_sol.nis_consistency(None)?);
231
232    od_sol.to_parquet(
233        "./data/04_output/06_lunar_od_results.parquet",
234        ExportCfg::default(),
235    )?;
236
237    let od_trajectory = od_sol.to_traj()?;
238    // Build the RIC difference.
239    od_trajectory.ric_diff_to_parquet(
240        &truth_traj,
241        "./data/04_output/06_lunar_od_truth_error.parquet",
242        ExportCfg::default(),
243    )?;
244
245    Ok(())
246}
06_seq_od/main.rs

06_seq_od/
main.rs
