Skip to main content

TrackingDataArc

Struct TrackingDataArc 

Source
pub struct TrackingDataArc {
    pub measurements: Vec<Measurement>,
    pub source: Option<String>,
    pub moduli: Option<IndexMap<MeasurementType, f64>>,
    pub force_reject: bool,
}
Expand description

Tracking data storing all of measurements as a B-Tree. It inherently does NOT support multiple concurrent measurements from several trackers.

§Measurement Moduli, e.g. range modulus

In the case of ranging, and possibly other data types, a code is used to measure the range to the spacecraft. The length of this code determines the ambiguity resolution, as per equation 9 in section 2.2.2.2 of the JPL DESCANSO, document 214, Pseudo-Noise and Regenerative Ranging. For example, using the JPL Range Code and a frequency range clock of 1 MHz, the range ambiguity is 75,660 km. In other words, as soon as the spacecraft is at a range of 75,660 + 1 km the JPL Range Code will report the vehicle to be at a range of 1 km. This is simply because the range code overlaps with itself, effectively loosing track of its own reference: it’s due to the phase shift of the signal “lapping” the original signal length.

            (Spacecraft)
            ^
            |    Actual Distance = 75,661 km
            |
0 km                                         75,660 km (Wrap-Around)
|-----------------------------------------------|
  When the "code length" is exceeded,
  measurements wrap back to 0.

So effectively:
    Observed code range = Actual range (mod 75,660 km)
    75,661 km → 1 km

Nyx can only resolve the range ambiguity if the tracking data specifies a modulus for this specific measurement type. For example, in the case of the JPL Range Code and a 1 MHz range clock, the ambiguity interval is 75,660 km.

The measurement used in the Orbit Determination Process then becomes the following, where // represents the Euclidian division.

k = computed_obs // ambiguity_interval
real_obs = measured_obs + k * modulus

Reference: JPL DESCANSO, document 214, Pseudo-Noise and Regenerative Ranging.

:type measurements: listMeasurement

Fields§

§measurements: Vec<Measurement>

All measurements in this data arc

§source: Option<String>

Source file if loaded from a file or saved to a file.

§moduli: Option<IndexMap<MeasurementType, f64>>

Optionally provide a map of modulos (e.g. the RANGE_MODULO of CCSDS TDM).

§force_reject: bool

Reject all of the measurements, useful for debugging passes.

Implementations§

Source§

impl TrackingDataArc

Source

pub fn from_tdm<P: AsRef<Path>>( path: P, aliases: Option<HashMap<String, String>>, ) -> Result<Self, InputOutputError>

Loads a tracking arc from its serialization in CCSDS TDM.

§Support level
  • Only the KVN format is supported.
  • Support is limited to orbit determination in “xGEO”, i.e. cislunar and deep space missions.
  • Only one metadata and data section per file is tested.
§Data types

Fully supported: - RANGE - DOPPLER_INSTANTANEOUS, DOPPLER_INTEGRATED - ANGLE_1 / ANGLE_2, as azimuth/elevation only

Partially supported: - TRANSMIT_FREQ / RECEIVE_FREQ : these will be converted to Doppler measurements using the TURNAROUND_NUMERATOR and TURNAROUND_DENOMINATOR in the TDM. The freq rate is not supported.

§Metadata support
§Mode

Only the MODE = SEQUENTIAL is supported.

§Time systems / time scales

All timescales supported by hifitime are supported here. This includes: UTC, TAI, GPS, TT, TDB, TAI, GST, QZSST, TL, TCL.

§Path

Only one way or two way data is supported, i.e. path must be either PATH n,m,n or PATH n,m.

Note that the actual indexes of the path are ignored.

§Participants

PARTICIPANT_1 must be the ground station / tracker. The second participant is ignored: the user must ensure that the Orbit Determination Process is properly configured and the proper arc is given.

§Turnaround ratio

The turnaround ratio is only accounted for when the data contains RECEIVE_FREQ and TRANSMIT_FREQ data.

§Range and modulus

Only kilometers are supported in range units. Range modulus is accounted for to compute range ambiguity.

Source

pub fn to_tdm_file<P: AsRef<Path>>( self, spacecraft_name: String, aliases: Option<HashMap<String, String>>, path: P, cfg: ExportCfg, ) -> Result<PathBuf, InputOutputError>

Store this tracking arc to a CCSDS TDM file, with optional metadata and a timestamp appended to the filename.

Source§

impl TrackingDataArc

Source

pub fn from_parquet<P: AsRef<Path>>(path: P) -> Result<Self, InputOutputError>

Loads a tracking arc from its serialization in parquet.

Warning: no metadata is read from the parquet file, even that written to it by Nyx.

Source

pub fn to_parquet_simple<P: AsRef<Path>>( &self, path: P, ) -> Result<PathBuf, InputOutputError>

Store this tracking arc to a parquet file.

Examples found in repository?
nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 179)
34fn main() -> Result<(), Box<dyn Error>> {
35    pel::init();
36
37    // ====================== //
38    // === ALMANAC SET UP === //
39    // ====================== //
40
41    let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
42
43    let out = manifest_dir.join("data/04_output/");
44
45    let almanac = Arc::new(
46        Almanac::new(
47            &manifest_dir
48                .join("data/01_planetary/pck08.pca")
49                .to_string_lossy(),
50        )
51        .unwrap()
52        .load(
53            &manifest_dir
54                .join("data/01_planetary/de440s.bsp")
55                .to_string_lossy(),
56        )
57        .unwrap(),
58    );
59
60    let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
61    let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();
62
63    let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
64    let nrho = Orbit::cartesian(
65        166_473.631_302_239_7,
66        -274_715.487_253_382_7,
67        -211_233.210_176_686_7,
68        0.933_451_604_520_018_4,
69        0.436_775_046_841_900_9,
70        -0.082_211_021_250_348_95,
71        epoch,
72        eme2k,
73    );
74
75    let tx_nrho_sc = Spacecraft::from(nrho);
76
77    let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
78    println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");
79
80    let bodies = vec![EARTH, SUN];
81    let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));
82
83    let setup = Propagator::rk89(
84        dynamics,
85        IntegratorOptions::builder().max_step(0.5.minutes()).build(),
86    );
87
88    /* == Propagate the NRHO vehicle == */
89    let prop_time = 1.1 * state_luna.period().unwrap();
90
91    let (nrho_final, mut tx_traj) = setup
92        .with(tx_nrho_sc, almanac.clone())
93        .for_duration_with_traj(prop_time)
94        .unwrap();
95
96    tx_traj.name = Some("NRHO Tx SC".to_string());
97
98    println!("{tx_traj}");
99
100    /* == Propagate an LLO vehicle == */
101    let llo_orbit =
102        Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();
103
104    let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();
105
106    let (_, llo_traj) = setup
107        .with(llo_sc, almanac.clone())
108        .until_epoch_with_traj(nrho_final.epoch())
109        .unwrap();
110
111    // Export the subset of the first two hours.
112    llo_traj
113        .clone()
114        .filter_by_offset(..2.hours())
115        .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;
116
117    /* == Setup the interlink == */
118
119    let mut measurement_types = IndexSet::new();
120    measurement_types.insert(MeasurementType::Range);
121    measurement_types.insert(MeasurementType::Doppler);
122
123    let mut stochastics = IndexMap::new();
124
125    let sa45_csac_allan_dev = 1e-11;
126
127    stochastics.insert(
128        MeasurementType::Range,
129        StochasticNoise::from_hardware_range_km(
130            sa45_csac_allan_dev,
131            10.0.seconds(),
132            link_specific::ChipRate::StandardT4B(),
133            link_specific::SN0::Average(),
134        ),
135    );
136
137    stochastics.insert(
138        MeasurementType::Doppler,
139        StochasticNoise::from_hardware_doppler_km_s(
140            sa45_csac_allan_dev,
141            10.0.seconds(),
142            link_specific::CarrierFreq::SBand(),
143            link_specific::CN0::Average(),
144        ),
145    );
146
147    let interlink = InterlinkTxSpacecraft {
148        traj: tx_traj,
149        measurement_types,
150        integration_time: None,
151        timestamp_noise_s: None,
152        ab_corr: Aberration::LT,
153        stochastic_noises: Some(stochastics),
154    };
155
156    // Devices are the transmitter, which is our NRHO vehicle.
157    let mut devices = BTreeMap::new();
158    devices.insert("NRHO Tx SC".to_string(), interlink);
159
160    let mut configs = BTreeMap::new();
161    configs.insert(
162        "NRHO Tx SC".to_string(),
163        TrkConfig::builder()
164            .strands(vec![Strand {
165                start: epoch,
166                end: nrho_final.epoch(),
167            }])
168            .build(),
169    );
170
171    let mut trk_sim =
172        TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
173    println!("{trk_sim}");
174
175    let trk_data = trk_sim.generate_measurements(&almanac).unwrap();
176    println!("{trk_data}");
177
178    trk_data
179        .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
180        .unwrap();
181
182    // Run a truth OD where we estimate the LLO position
183    let llo_uncertainty = SpacecraftUncertainty::builder()
184        .nominal(llo_sc)
185        .x_km(1.0)
186        .y_km(1.0)
187        .z_km(1.0)
188        .vx_km_s(1e-3)
189        .vy_km_s(1e-3)
190        .vz_km_s(1e-3)
191        .build();
192
193    let mut proc_devices = devices.clone();
194
195    // Define the initial estimate, randomized, seed for reproducibility
196    let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
197    // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
198    initial_estimate.covar *= 2.5;
199
200    // Increase the noise in the devices to accept more measurements.
201
202    for link in proc_devices.values_mut() {
203        for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
204            *noise.white_noise.as_mut().unwrap() *= 3.0;
205        }
206    }
207
208    let init_err = initial_estimate
209        .orbital_state()
210        .ric_difference(&llo_orbit)
211        .unwrap();
212
213    println!("initial estimate:\n{initial_estimate}");
214    println!("RIC errors = {init_err}",);
215
216    let odp = InterlinkKalmanOD::new(
217        setup.clone(),
218        KalmanVariant::ReferenceUpdate,
219        Some(SigmaRejection::default()),
220        proc_devices,
221        almanac.clone(),
222    );
223
224    // Shrink the data to process.
225    let arc = trk_data.filter_by_offset(..2.hours());
226
227    let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();
228
229    println!("{od_sol}");
230
231    od_sol
232        .to_parquet(
233            out.join("05_caps_interlink_od_sol.pq"),
234            ExportCfg::default(),
235        )
236        .unwrap();
237
238    let od_traj = od_sol.to_traj().unwrap();
239
240    od_traj
241        .ric_diff_to_parquet(
242            &llo_traj,
243            out.join("05_caps_interlink_llo_est_error.pq"),
244            ExportCfg::default(),
245        )
246        .unwrap();
247
248    let final_est = od_sol.estimates.last().unwrap();
249    assert!(final_est.within_3sigma(), "should be within 3 sigma");
250
251    println!("ESTIMATE\n{final_est:x}\n");
252    let truth = llo_traj.at(final_est.epoch()).unwrap();
253    println!("TRUTH\n{truth:x}");
254
255    let final_err = truth
256        .orbit
257        .ric_difference(&final_est.orbital_state())
258        .unwrap();
259    println!("ERROR {final_err}");
260
261    // Build the residuals versus reference plot.
262    let rvr_sol = odp
263        .process_arc(initial_estimate, &arc.resid_vs_ref_check())
264        .unwrap();
265
266    rvr_sol
267        .to_parquet(
268            out.join("05_caps_interlink_resid_v_ref.pq"),
269            ExportCfg::default(),
270        )
271        .unwrap();
272
273    let final_rvr = rvr_sol.estimates.last().unwrap();
274
275    println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
276    println!(
277        "Pure prop error {:.3} m",
278        final_rvr
279            .orbital_state()
280            .ric_difference(&final_est.orbital_state())
281            .unwrap()
282            .rmag_km()
283            * 1e3
284    );
285
286    Ok(())
287}
More examples
Hide additional examples
nyx-core/examples/06_lunar_orbit_determination/main.rs (line 160)
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}
nyx-core/examples/04_lro_od/main.rs (line 249)
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 output_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "../data", "04_output"]
46        .iter()
47        .collect();
48
49    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
50        .iter()
51        .collect();
52
53    let meta = data_folder.join("lro-dynamics.dhall");
54
55    // Load this ephem in the general Almanac we're using for this analysis.
56    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().as_ref())
57        .map_err(Box::new)?
58        .process(true)
59        .map_err(Box::new)?;
60
61    let mut moon_pc = almanac.get_planetary_data_from_id(MOON).unwrap();
62    moon_pc.mu_km3_s2 = 4902.74987;
63    almanac.set_planetary_data_from_id(MOON, moon_pc).unwrap();
64
65    let mut earth = almanac.get_planetary_data_from_id(EARTH).unwrap();
66    earth.mu_km3_s2 = 398600.436;
67    almanac.set_planetary_data_from_id(EARTH, earth).unwrap();
68
69    // Save this new kernel for reuse.
70    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
71    almanac
72        .planetary_data
73        .values()
74        .next()
75        .unwrap()
76        .save_as(&data_folder.join("lro-specific.pca"), true)?;
77
78    // Lock the almanac (an Arc is a read only structure).
79    let almanac = Arc::new(almanac);
80
81    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
82    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
83    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
84    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
85    let lro_frame = Frame::from_ephem_j2000(-85);
86
87    // To build the trajectory we need to provide a spacecraft template.
88    let sc_template = Spacecraft::builder()
89        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
90        .srp(SRPData {
91            // SRP configuration is arbitrary, but we will be estimating it anyway.
92            area_m2: 3.9 * 2.7,
93            coeff_reflectivity: 0.96,
94        })
95        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
96        .build();
97    // Now we can build the trajectory from the BSP file.
98    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
99    let traj_as_flown = Traj::from_bsp(
100        lro_frame,
101        MOON_J2000,
102        &almanac,
103        sc_template,
104        5.seconds(),
105        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
106        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
107        Aberration::LT,
108        Some("LRO".to_string()),
109    )?;
110
111    println!("{traj_as_flown}");
112
113    // ====================== //
114    // === MODEL MATCHING === //
115    // ====================== //
116
117    // Set up the spacecraft dynamics.
118
119    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
120    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
121    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
122
123    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
124    // We're using the GRAIL JGGRX model.
125    let mut jggrx_meta = MetaFile {
126        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
127        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
128    };
129    // And let's download it if we don't have it yet.
130    jggrx_meta.process(true)?;
131
132    // Build the spherical harmonics.
133    // The harmonics must be computed in the body fixed frame.
134    // We're using the long term prediction of the Moon principal axes frame.
135    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
136    let sph_harmonics = GravityField::new(GravityFieldData::from_shadr(
137        &jggrx_meta.uri,
138        80,
139        80,
140        almanac.frame_info(moon_pa_frame)?,
141    )?);
142
143    // Include the spherical harmonics into the orbital dynamics.
144    orbital_dyn.accel_models.push(sph_harmonics);
145
146    // We define the solar radiation pressure, using the default solar flux and accounting only
147    // for the eclipsing caused by the Earth and Moon.
148    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
149    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;
150
151    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
152    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
153    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
154
155    println!("{dynamics}");
156
157    // Now we can build the propagator.
158    let setup = Propagator::default_dp78(dynamics.clone());
159
160    // For reference, let's build the trajectory with Nyx's models from that LRO state.
161    let (sim_final, traj_as_sim) = setup
162        .with(*traj_as_flown.first(), almanac.clone())
163        .until_epoch_with_traj(traj_as_flown.last().epoch())?;
164
165    println!("SIM INIT:  {:x}", traj_as_flown.first());
166    println!("SIM FINAL: {sim_final:x}");
167    // Compute RIC difference between SIM and LRO ephem
168    let sim_lro_delta = sim_final
169        .orbit
170        .ric_difference(&traj_as_flown.last().orbit)?;
171    println!("{traj_as_sim}");
172    println!(
173        "SIM v LRO - RIC Position (m): {:.3}",
174        sim_lro_delta.radius_km * 1e3
175    );
176    println!(
177        "SIM v LRO - RIC Velocity (m/s): {:.3}",
178        sim_lro_delta.velocity_km_s * 1e3
179    );
180
181    traj_as_sim.ric_diff_to_parquet(
182        &traj_as_flown,
183        output_folder.join("./04_lro_sim_truth_error.parquet"),
184        ExportCfg::default(),
185    )?;
186
187    // ==================== //
188    // === OD SIMULATOR === //
189    // ==================== //
190
191    // 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
192    // and the truth LRO state.
193
194    // Therefore, we will actually run an estimation from a dispersed LRO state.
195    // The sc_seed is the true LRO state from the BSP.
196    let sc_seed = *traj_as_flown.first();
197
198    // Load the Deep Space Network ground stations.
199    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
200    let ground_station_file: PathBuf = [
201        env!("CARGO_MANIFEST_DIR"),
202        "examples",
203        "04_lro_od",
204        "dsn-network.yaml",
205    ]
206    .iter()
207    .collect();
208
209    let devices = GroundStation::load_named(ground_station_file)?;
210
211    let mut proc_devices = devices.clone();
212
213    // Increase the noise in the devices to accept more measurements.
214    for gs in proc_devices.values_mut() {
215        if let Some(noise) = &mut gs
216            .stochastic_noises
217            .as_mut()
218            .unwrap()
219            .get_mut(&MeasurementType::Range)
220        {
221            *noise.white_noise.as_mut().unwrap() *= 3.0;
222        }
223    }
224
225    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
226    // Nyx can build a tracking schedule for you based on the first station with access.
227    let trkconfg_yaml: PathBuf = [
228        env!("CARGO_MANIFEST_DIR"),
229        "examples",
230        "04_lro_od",
231        "tracking-cfg.yaml",
232    ]
233    .iter()
234    .collect();
235
236    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;
237
238    // Build the tracking arc simulation to generate a "standard measurement".
239    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
240        devices.clone(),
241        traj_as_flown.clone(),
242        configs,
243        123, // Set a seed for reproducibility
244    )?;
245
246    trk.build_schedule(&almanac)?;
247    let arc = trk.generate_measurements(&almanac)?;
248    // Save the simulated tracking data
249    arc.to_parquet_simple(output_folder.join("04_lro_simulated_tracking.parquet"))?;
250
251    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
252    println!("{arc}");
253
254    // Now that we have simulated measurements, we'll run the orbit determination.
255
256    // ===================== //
257    // === OD ESTIMATION === //
258    // ===================== //
259
260    let sc = SpacecraftUncertainty::builder()
261        .nominal(sc_seed)
262        .frame(LocalFrame::RIC)
263        .x_km(0.5)
264        .y_km(0.5)
265        .z_km(0.5)
266        .vx_km_s(5e-3)
267        .vy_km_s(5e-3)
268        .vz_km_s(5e-3)
269        .build();
270
271    // Build the filter initial estimate, which we will reuse in the filter.
272    let mut initial_estimate = sc.to_estimate()?;
273    initial_estimate.covar *= 3.0;
274
275    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
276
277    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
278    let process_noise = ProcessNoise3D::from_velocity_km_s(
279        &[1e-12, 1e-12, 1e-12],
280        1 * Unit::Hour,
281        10 * Unit::Minute,
282        None,
283    );
284
285    println!("{process_noise}");
286
287    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
288    let odp = SpacecraftKalmanOD::new(
289        setup,
290        KalmanVariant::ReferenceUpdate,
291        Some(SigmaRejection::default()),
292        proc_devices,
293        almanac.clone(),
294    )
295    .with_process_noise(process_noise);
296
297    let od_sol = odp.process_arc(initial_estimate, &arc)?;
298
299    let final_est = od_sol.estimates.last().unwrap();
300
301    println!("{final_est}");
302
303    let ric_err = traj_as_flown
304        .at(final_est.epoch())?
305        .orbit
306        .ric_difference(&final_est.orbital_state())?;
307    println!("== RIC at end ==");
308    println!("RIC Position (m): {:.3}", ric_err.radius_km * 1e3);
309    println!("RIC Velocity (m/s): {:.3}", ric_err.velocity_km_s * 1e3);
310
311    println!(
312        "Num residuals rejected: #{}",
313        od_sol.rejected_residuals().len()
314    );
315    println!(
316        "Percentage within +/-3: {}",
317        od_sol.residual_ratio_within_threshold(3.0).unwrap()
318    );
319    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());
320
321    od_sol.to_parquet(
322        output_folder.join("04_lro_od_results.parquet"),
323        ExportCfg::default(),
324    )?;
325
326    // Create the ephemeris
327    let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
328    let ephem_start = ephem.start_epoch().unwrap();
329    let ephem_end = ephem.end_epoch().unwrap();
330    // Check that the covariance is PSD throughout the ephemeris by interpolating it.
331    for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
332        ephem
333            .covar_at(
334                epoch,
335                anise::ephemerides::ephemeris::LocalFrame::RIC,
336                &almanac,
337            )
338            .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
339    }
340    // Export as BSP!
341    ephem
342        .write_spice_bsp(
343            -85,
344            output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
345            None,
346        )
347        .expect("could not built BSP");
348    let new_almanac = Almanac::default()
349        .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
350        .unwrap();
351    new_almanac.describe(None, None, None, None, None, None, None, None);
352    let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();
353
354    assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
355    assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);
356
357    // In our case, we have the truth trajectory from NASA.
358    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
359    // Export the OD trajectory first.
360    let od_trajectory = od_sol.to_traj()?;
361    // Build the RIC difference.
362    od_trajectory.ric_diff_to_parquet(
363        &traj_as_flown,
364        output_folder.join("04_lro_od_truth_error.parquet"),
365        ExportCfg::default(),
366    )?;
367
368    Ok(())
369}
Source

pub fn to_parquet<P: AsRef<Path>>( &self, path: P, cfg: ExportCfg, ) -> Result<PathBuf, InputOutputError>

Store this tracking arc to a parquet file, with optional metadata and a timestamp appended to the filename.

Source§

impl TrackingDataArc

Source

pub fn sort(&mut self)

Sort these measurements by epoch :rtype: None

Source

pub fn start_epoch(&self) -> Option<Epoch>

Returns the start epoch of this tracking arc :rtype: Epoch | None

Source

pub fn end_epoch(&self) -> Option<Epoch>

Returns the end epoch of this tracking arc :rtype: Epoch | None

Source

pub fn duration(&self) -> Option<Duration>

Returns the duration this tracking arc :rtype: Duration | None

Source

pub fn len(&self) -> usize

Returns the number of measurements in this data arc :rtype: int

Source

pub fn is_empty(&self) -> bool

Returns whether this arc has no measurements. :rtype: bool

Source

pub fn min_duration_sep(&self) -> Option<Duration>

Returns the minimum duration between two subsequent measurements. :rtype: Duration | None

Source

pub fn set_moduli(&mut self, msr_type: MeasurementType, modulus: f64)

Set (or overwrites) the modulus of the provided measurement type.

:type msr_type: MeasurementType :type modulus: float :rtype: None

Source

pub fn apply_moduli(&mut self)

Applies the moduli to each measurement, if defined. :rtype: None

Source

pub fn downsample(&self, target_step: Duration) -> Self

Downsamples the tracking data to a lower frequency using a simple moving average low-pass filter followed by decimation, returning new TrackingDataArc with downsampled measurements.

It provides a computationally efficient approach to reduce the sampling rate while mitigating aliasing effects.

§Algorithm
  1. A simple moving average filter is applied as a low-pass filter.
  2. Decimation is performed by selecting every Nth sample after filtering.
§Advantages
  • Computationally efficient, suitable for large datasets common in spaceflight applications.
  • Provides basic anti-aliasing, crucial for preserving signal integrity in orbit determination and tracking.
  • Maintains phase information, important for accurate timing in spacecraft state estimation.
§Limitations
  • The frequency response is not as sharp as more sophisticated filters (e.g., FIR, IIR).
  • May not provide optimal stopband attenuation for high-precision applications.
§Considerations for Spaceflight Applications
  • Suitable for initial data reduction in ground station tracking pipelines.
  • Adequate for many orbit determination and tracking tasks where computational speed is prioritized.
  • For high-precision applications (e.g., interplanetary navigation), consider using more advanced filtering techniques.

:type target_step: Duration :rtype: TrackingDataArc

Source

pub fn chunk(&self, max_duration: Duration) -> Vec<TrackingDataArc>

Splits a long tracking data arc into smaller chunks, each up to max_duration long.

:type max_duration: Duration :rtype: listTrackingDataArc

Source§

impl TrackingDataArc

Source

pub fn unique_aliases(&self) -> IndexSet<String>

Returns the unique list of aliases in this tracking data arc

Source

pub fn unique_types(&self) -> IndexSet<MeasurementType>

Returns the unique measurement types in this tracking data arc

Source

pub fn unique(&self) -> (IndexSet<String>, IndexSet<MeasurementType>)

Returns the unique trackers and unique measurement types in this data arc

Source

pub fn filter_by_epoch<R: RangeBounds<Epoch>>(self, bound: R) -> Self

Returns a new tracking arc that only contains measurements that fall within the given epoch range.

Executes in O(N) time strictly due to memory shifting, requiring zero new allocations.

Source

pub fn filter_by_offset<R: RangeBounds<Duration>>(self, bound: R) -> Self

Returns a new tracking arc that only contains measurements that fall within the given offset from the first epoch. For example, a bound of 30.minutes()..90.minutes() will only read measurements from the start of the arc + 30 minutes until start + 90 minutes.

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

pub fn filter_by_tracker(self, tracker: String) -> Self

Returns a new tracking arc that only contains measurements from the desired tracker.

Source

pub fn filter_by_measurement_type(self, included_type: MeasurementType) -> Self

Returns a new tracking arc that only contains measurements of the provided type.

Source

pub fn exclude_tracker(self, excluded_tracker: String) -> Self

Returns a new tracking arc that contains measurements from all trackers except the one provided

Source

pub fn exclude_by_epoch<R: RangeBounds<Epoch>>(self, bound: R) -> Self

Returns a new tracking arc that excludes measurements within the given epoch range.

Executes an in-place O(N) memory shift with zero heap allocations.

Source

pub fn exclude_measurement_type(self, excluded_type: MeasurementType) -> Self

Returns a new tracking arc that contains measurements from all trackers except the one provided

Source

pub fn reject_by_epoch<R: RangeBounds<Epoch>>(self, bound: R) -> Self

Marks measurements within the given epoch range as rejected.

Operates in O(log N) for bound resolution and O(K) for iteration, where K is the slice length.

Source

pub fn reject_by_tracker(self, tracker: &str) -> Self

Marks measurements from the provided tracker as rejected. Requires an O(N) scan. The parameter is downgraded to &str to prevent heap allocations.

Source

pub fn resid_vs_ref_check(self) -> Self

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

Trait Implementations§

Source§

impl Add for TrackingDataArc

Source§

type Output = TrackingDataArc

The resulting type after applying the + operator.
Source§

fn add(self, rhs: Self) -> Self::Output

Performs the + operation. Read more
Source§

impl AddAssign for TrackingDataArc

Source§

fn add_assign(&mut self, rhs: Self)

Performs the += operation. Read more
Source§

impl Clone for TrackingDataArc

Source§

fn clone(&self) -> TrackingDataArc

Returns a duplicate of the value. Read more
1.0.0 (const: unstable) · Source§

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

Performs copy-assignment from source. Read more
Source§

impl Debug for TrackingDataArc

Source§

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

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

impl Default for TrackingDataArc

Source§

fn default() -> TrackingDataArc

Returns the “default value” for a type. Read more
Source§

impl DerefToPyAny for TrackingDataArc

Source§

impl Display for TrackingDataArc

Source§

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

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

impl<'a, 'py> FromPyObject<'a, 'py> for TrackingDataArc
where Self: Clone,

Source§

type Error = PyClassGuardError<'a, 'py>

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

fn extract( obj: Borrowed<'a, 'py, PyAny>, ) -> Result<Self, <Self as FromPyObject<'a, 'py>>::Error>

Extracts Self from the bound smart pointer obj. Read more
Source§

impl<'py> IntoPyObject<'py> for TrackingDataArc

Source§

type Target = TrackingDataArc

The Python output type
Source§

type Output = Bound<'py, <TrackingDataArc as IntoPyObject<'py>>::Target>

The smart pointer type to use. Read more
Source§

type Error = PyErr

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

fn into_pyobject( self, py: Python<'py>, ) -> Result<<Self as IntoPyObject<'_>>::Output, <Self as IntoPyObject<'_>>::Error>

Performs the conversion.
Source§

impl PartialEq for TrackingDataArc

Source§

fn eq(&self, other: &Self) -> bool

Equality operator ==. Read more
1.0.0 (const: unstable) · Source§

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=. Read more
Source§

impl PyClass for TrackingDataArc

Source§

const NAME: &str = "TrackingDataArc"

Name of the class. Read more
Source§

type Frozen = False

Whether the pyclass is frozen. Read more
Source§

impl PyClassImpl for TrackingDataArc

Source§

const MODULE: Option<&str> = ::core::option::Option::None

Module which the class will be associated with. Read more
Source§

const IS_BASETYPE: bool = false

#[pyclass(subclass)]
Source§

const IS_SUBCLASS: bool = false

#[pyclass(extends=…)]
Source§

const IS_MAPPING: bool = false

#[pyclass(mapping)]
Source§

const IS_SEQUENCE: bool = false

#[pyclass(sequence)]
Source§

const IS_IMMUTABLE_TYPE: bool = false

#[pyclass(immutable_type)]
Source§

const RAW_DOC: &'static CStr = /// Tracking data storing all of measurements as a B-Tree. /// It inherently does NOT support multiple concurrent measurements from several trackers. /// /// # Measurement Moduli, e.g. range modulus /// /// In the case of ranging, and possibly other data types, a code is used to measure the range to the spacecraft. The length of this code /// determines the ambiguity resolution, as per equation 9 in section 2.2.2.2 of the JPL DESCANSO, document 214, _Pseudo-Noise and Regenerative Ranging_. /// For example, using the JPL Range Code and a frequency range clock of 1 MHz, the range ambiguity is 75,660 km. In other words, /// as soon as the spacecraft is at a range of 75,660 + 1 km the JPL Range Code will report the vehicle to be at a range of 1 km. /// This is simply because the range code overlaps with itself, effectively loosing track of its own reference: /// it's due to the phase shift of the signal "lapping" the original signal length. /// /// ```text /// (Spacecraft) /// ^ /// | Actual Distance = 75,661 km /// | /// 0 km 75,660 km (Wrap-Around) /// |-----------------------------------------------| /// When the "code length" is exceeded, /// measurements wrap back to 0. /// /// So effectively: /// Observed code range = Actual range (mod 75,660 km) /// 75,661 km → 1 km /// /// ``` /// /// Nyx can only resolve the range ambiguity if the tracking data specifies a modulus for this specific measurement type. /// For example, in the case of the JPL Range Code and a 1 MHz range clock, the ambiguity interval is 75,660 km. /// /// The measurement used in the Orbit Determination Process then becomes the following, where `//` represents the [Euclidian division](https://doc.rust-lang.org/std/primitive.f64.html#method.div_euclid). /// /// ```text /// k = computed_obs // ambiguity_interval /// real_obs = measured_obs + k * modulus /// ``` /// /// Reference: JPL DESCANSO, document 214, _Pseudo-Noise and Regenerative Ranging_. /// /// :type measurements: list[Measurement]

Docstring for the class provided on the struct or enum. Read more
Source§

const DOC: &'static CStr

Fully rendered class doc, including the text_signature if a constructor is defined. Read more
Source§

type Layout = <<TrackingDataArc as PyClassImpl>::BaseNativeType as PyClassBaseType>::Layout<TrackingDataArc>

Description of how this class is laid out in memory
Source§

type BaseType = PyAny

Base class
Source§

type ThreadChecker = NoopThreadChecker

This handles following two situations: Read more
Source§

type Inventory = Pyo3MethodsInventoryForTrackingDataArc

Source§

type PyClassMutability = <<PyAny as PyClassBaseType>::PyClassMutability as PyClassMutability>::MutableChild

Immutable or mutable
Source§

type Dict = PyClassDummySlot

Specify this class has #[pyclass(dict)] or not.
Source§

type WeakRef = PyClassDummySlot

Specify this class has #[pyclass(weakref)] or not.
Source§

type BaseNativeType = PyAny

The closest native ancestor. This is PyAny by default, and when you declare #[pyclass(extends=PyDict)], it’s PyDict.
Source§

fn items_iter() -> PyClassItemsIter

Source§

fn lazy_type_object() -> &'static LazyTypeObject<Self>

§

fn dict_offset() -> Option<PyObjectOffset>

Used to provide the dictoffset slot (equivalent to tp_dictoffset)
§

fn weaklist_offset() -> Option<PyObjectOffset>

Used to provide the weaklistoffset slot (equivalent to tp_weaklistoffset
Source§

impl PyClassNewTextSignature for TrackingDataArc

Source§

const TEXT_SIGNATURE: &'static str = "(measurements)"

Source§

impl PyTypeInfo for TrackingDataArc

Source§

const NAME: &str = <Self as ::pyo3::PyClass>::NAME

👎Deprecated since 0.28.0:

prefer using ::type_object(py).name() to get the correct runtime value

Class name.
Source§

const MODULE: Option<&str> = <Self as ::pyo3::impl_::pyclass::PyClassImpl>::MODULE

👎Deprecated since 0.28.0:

prefer using ::type_object(py).module() to get the correct runtime value

Module name, if any.
Source§

fn type_object_raw(py: Python<'_>) -> *mut PyTypeObject

Returns the PyTypeObject instance for this type.
§

fn type_object(py: Python<'_>) -> Bound<'_, PyType>

Returns the safe abstraction over the type object.
§

fn is_type_of(object: &Bound<'_, PyAny>) -> bool

Checks if object is an instance of this type or a subclass of this type.
§

fn is_exact_type_of(object: &Bound<'_, PyAny>) -> bool

Checks if object is an instance of this type.

Auto Trait Implementations§

Blanket Implementations§

§

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

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
§

impl<ST, DT> CastableFrom<ST, Initialized, Initialized> for DT
where ST: ?Sized, DT: ?Sized,

§

impl<ST, DT> CastableFrom<ST, Uninit, Uninit> for DT
where ST: ?Sized, DT: ?Sized,

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
§

impl<T, Right> ClosedAdd<Right> for T
where T: Add<Right, Output = T> + AddAssign<Right>,

§

impl<T, Right> ClosedAddAssign<Right> for T
where T: ClosedAdd<Right> + AddAssign<Right>,

Source§

impl<T> From<T> for T

Source§

fn from(t: T) -> T

Returns the argument unchanged.

§

impl<'py, T> FromPyObjectOwned<'py> for T
where T: for<'a> FromPyObject<'a, 'py>,

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<'py, T> IntoPyObjectExt<'py> for T
where T: IntoPyObject<'py>,

§

fn into_bound_py_any(self, py: Python<'py>) -> Result<Bound<'py, PyAny>, PyErr>

Converts self into an owned Python object, dropping type information.
§

fn into_py_any(self, py: Python<'py>) -> Result<Py<PyAny>, PyErr>

Converts self into an owned Python object, dropping type information and unbinding it from the 'py lifetime.
§

fn into_pyobject_or_pyerr(self, py: Python<'py>) -> Result<Self::Output, PyErr>

Converts self into a Python object. 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
§

impl<T> PyErrArguments for T
where T: for<'py> IntoPyObject<'py> + Send + Sync,

§

fn arguments(self, py: Python<'_>) -> Py<PyAny>

Arguments for exception
§

impl<T> PyTypeCheck for T
where T: PyTypeInfo,

§

const NAME: &'static str = T::NAME

👎Deprecated since 0.27.0:

Use ::classinfo_object() instead and format the type name at runtime. Note that using built-in cast features is often better than manual PyTypeCheck usage.

Name of self. This is used in error messages, for example.
§

fn type_check(object: &Bound<'_, PyAny>) -> bool

Checks if object is an instance of Self, which may include a subtype. Read more
§

fn classinfo_object(py: Python<'_>) -> Bound<'_, PyAny>

Returns the expected type as a possible argument for the isinstance and issubclass function. Read more
§

impl<T> Read<Exclusive, BecauseExclusive> for T
where T: ?Sized,

Source§

impl<T> Same for T

Source§

type Output = T

Should always be Self
Source§

impl<T> Scalar for T
where T: 'static + Clone + PartialEq + Debug,

§

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<T> Ungil for T
where T: Send,