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, MOON, SUN},
10        frames::{EARTH_J2000, MOON_J2000, MOON_PA_FRAME},
11    },
12    prelude::Almanac,
13};
14use hifitime::{Epoch, TimeSeries, TimeUnits, Unit};
15use nyx::{
16    cosmic::{Aberration, Frame, Mass, MetaAlmanac, SRPData},
17    dynamics::{
18        guidance::LocalFrame, GravityField, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
19    },
20    io::{ConfigRepr, ExportCfg},
21    md::prelude::{GravityFieldData, Traj},
22    od::{
23        msr::MeasurementType,
24        prelude::{KalmanVariant, TrackingArcSim, TrkConfig},
25        process::{Estimate, NavSolution, ResidRejectCrit, SpacecraftUncertainty},
26        snc::ProcessNoise3D,
27        GroundStation, SpacecraftKalmanOD, SpacecraftKalmanScalarOD,
28    },
29    propagators::{IntegratorOptions, Propagator},
30    Orbit, Spacecraft, State,
31};
32
33use std::{collections::BTreeMap, error::Error, path::PathBuf, str::FromStr, sync::Arc};
34
35// TODO: Convert this to a Spacecraft Sequence
36
37fn main() -> Result<(), Box<dyn Error>> {
38    pel::init();
39
40    // ====================== //
41    // === ALMANAC SET UP === //
42    // ====================== //
43
44    // Dynamics models require planetary constants and ephemerides to be defined.
45    // Let's start by grabbing those by using ANISE's MetaAlmanac.
46
47    let data_folder: PathBuf = [
48        env!("CARGO_MANIFEST_DIR"),
49        "examples",
50        "06_lunar_orbit_determination",
51    ]
52    .iter()
53    .collect();
54
55    let meta = data_folder.join("metaalmanac.dhall");
56
57    // Load this ephem in the general Almanac we're using for this analysis.
58    let almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
59        .map_err(Box::new)?
60        .process(true)
61        .map_err(Box::new)?;
62
63    // Lock the almanac (an Arc is a read only structure).
64    let almanac = Arc::new(almanac);
65
66    // Build a nominal trajectory
67    // TODO: Switch this to a sequence once the OD over a spacecraft sequence is implemented.
68
69    let epoch = Epoch::from_gregorian_utc_at_noon(2024, 2, 29);
70    let moon_j2000 = almanac.frame_info(MOON_J2000)?;
71
72    // To build the trajectory we need to provide a spacecraft template.
73    let orbiter = Spacecraft::builder()
74        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0))
75        .srp(SRPData {
76            area_m2: 3.9 * 2.7,
77            coeff_reflectivity: 0.96,
78        })
79        .orbit(Orbit::try_keplerian_altitude(
80            150.0, 0.00212, 33.6, 45.0, 45.0, 0.0, epoch, moon_j2000,
81        )?) // Setting a zero orbit here because it's just a template
82        .build();
83
84    // ========================== //
85    // === BUILD NOMINAL TRAJ === //
86    // ========================== //
87
88    // Set up the spacecraft dynamics.
89
90    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
91    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
92    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
93
94    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
95    // We're using the GRAIL JGGRX model.
96    let mut jggrx_meta = MetaFile {
97        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
98        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
99    };
100    // And let's download it if we don't have it yet.
101    jggrx_meta.process(true)?;
102
103    // Build the spherical harmonics.
104    // The harmonics must be computed in the body fixed frame.
105    // We're using the long term prediction of the Moon principal axes frame.
106    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
107    let sph_harmonics = GravityField::from_stor(
108        almanac.frame_info(moon_pa_frame)?,
109        GravityFieldData::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
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.clone())?;
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.clone())?;
158    let arc = trk.generate_measurements(almanac.clone())?;
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-13, 1e-13, 1e-13],
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(ResidRejectCrit::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!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
230
231    od_sol.to_parquet(
232        "./data/04_output/06_lunar_od_results.parquet",
233        ExportCfg::default(),
234    )?;
235
236    let od_trajectory = od_sol.to_traj()?;
237    // Build the RIC difference.
238    od_trajectory.ric_diff_to_parquet(
239        &truth_traj,
240        "./data/04_output/06_lunar_od_truth_error.parquet",
241        ExportCfg::default(),
242    )?;
243
244    Ok(())
245}
06_seq_od/main.rs

06_seq_od/
main.rs
