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

06_seq_od/
main.rs
