Struct TrackingDataArc

pub struct TrackingDataArc {
    pub measurements: BTreeMap<Epoch, Measurement>,
    pub source: Option<String>,
    pub moduli: Option<IndexMap<MeasurementType, f64>>,
}

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.

Fields

measurements: BTreeMap<Epoch, 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).

Implementations

impl TrackingDataArc

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.

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.

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.

impl TrackingDataArc

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.

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

Store this tracking arc to a parquet file.

Examples found in repository

examples/04_lro_od/main.rs (line 223)

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();

    // ====================== //
    // === ALMANAC SET UP === //
    // ====================== //

    // Dynamics models require planetary constants and ephemerides to be defined.
    // Let's start by grabbing those by using ANISE's MetaAlmanac.

    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
        .iter()
        .collect();

    let meta = data_folder.join("lro-dynamics.dhall");

    // Load this ephem in the general Almanac we're using for this analysis.
    let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
        .map_err(Box::new)?
        .process(true)
        .map_err(Box::new)?;

    let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
    moon_pc.mu_km3_s2 = 4902.74987;
    almanac.planetary_data.set_by_id(MOON, moon_pc)?;

    let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
    earth_pc.mu_km3_s2 = 398600.436;
    almanac.planetary_data.set_by_id(EARTH, earth_pc)?;

    // Save this new kernel for reuse.
    // In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
    almanac
        .planetary_data
        .save_as(&data_folder.join("lro-specific.pca"), true)?;

    // Lock the almanac (an Arc is a read only structure).
    let almanac = Arc::new(almanac);

    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
    let lro_frame = Frame::from_ephem_j2000(-85);

    // To build the trajectory we need to provide a spacecraft template.
    let sc_template = Spacecraft::builder()
        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0)) // Launch masses
        .srp(SRPData {
            // SRP configuration is arbitrary, but we will be estimating it anyway.
            area_m2: 3.9 * 2.7,
            coeff_reflectivity: 0.96,
        })
        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
        .build();
    // Now we can build the trajectory from the BSP file.
    // We'll arbitrarily set the tracking arc to 24 hours with a five second time step.
    let traj_as_flown = Traj::from_bsp(
        lro_frame,
        MOON_J2000,
        almanac.clone(),
        sc_template,
        5.seconds(),
        Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
        Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
        Aberration::LT,
        Some("LRO".to_string()),
    )?;

    println!("{traj_as_flown}");

    // ====================== //
    // === MODEL MATCHING === //
    // ====================== //

    // Set up the spacecraft dynamics.

    // Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
    // The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);

    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
    // We're using the GRAIL JGGRX model.
    let mut jggrx_meta = MetaFile {
        uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
        crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
    };
    // And let's download it if we don't have it yet.
    jggrx_meta.process(true)?;

    // Build the spherical harmonics.
    // The harmonics must be computed in the body fixed frame.
    // We're using the long term prediction of the Moon principal axes frame.
    let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
    let sph_harmonics = Harmonics::from_stor(
        almanac.frame_from_uid(moon_pa_frame)?,
        HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
    );

    // Include the spherical harmonics into the orbital dynamics.
    orbital_dyn.accel_models.push(sph_harmonics);

    // We define the solar radiation pressure, using the default solar flux and accounting only
    // for the eclipsing caused by the Earth and Moon.
    // Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;

    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);

    println!("{dynamics}");

    // Now we can build the propagator.
    let setup = Propagator::default_dp78(dynamics.clone());

    // For reference, let's build the trajectory with Nyx's models from that LRO state.
    let (sim_final, traj_as_sim) = setup
        .with(*traj_as_flown.first(), almanac.clone())
        .until_epoch_with_traj(traj_as_flown.last().epoch())?;

    println!("SIM INIT:  {:x}", traj_as_flown.first());
    println!("SIM FINAL: {sim_final:x}");
    // Compute RIC difference between SIM and LRO ephem
    let sim_lro_delta = sim_final
        .orbit
        .ric_difference(&traj_as_flown.last().orbit)?;
    println!("{traj_as_sim}");
    println!(
        "SIM v LRO - RIC Position (m): {:.3}",
        sim_lro_delta.radius_km * 1e3
    );
    println!(
        "SIM v LRO - RIC Velocity (m/s): {:.3}",
        sim_lro_delta.velocity_km_s * 1e3
    );

    traj_as_sim.ric_diff_to_parquet(
        &traj_as_flown,
        "./04_lro_sim_truth_error.parquet",
        ExportCfg::default(),
    )?;

    // ==================== //
    // === OD SIMULATOR === //
    // ==================== //

    // 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
    // and the truth LRO state.

    // Therefore, we will actually run an estimation from a dispersed LRO state.
    // The sc_seed is the true LRO state from the BSP.
    let sc_seed = *traj_as_flown.first();

    // Load the Deep Space Network ground stations.
    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
    let ground_station_file: PathBuf = [
        env!("CARGO_MANIFEST_DIR"),
        "examples",
        "04_lro_od",
        "dsn-network.yaml",
    ]
    .iter()
    .collect();

    let devices = GroundStation::load_named(ground_station_file)?;

    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
    // Nyx can build a tracking schedule for you based on the first station with access.
    let trkconfg_yaml: PathBuf = [
        env!("CARGO_MANIFEST_DIR"),
        "examples",
        "04_lro_od",
        "tracking-cfg.yaml",
    ]
    .iter()
    .collect();

    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml)?;

    // Build the tracking arc simulation to generate a "standard measurement".
    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
        devices.clone(),
        traj_as_flown.clone(),
        configs,
    )?;

    trk.build_schedule(almanac.clone())?;
    let arc = trk.generate_measurements(almanac.clone())?;
    // Save the simulated tracking data
    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;

    // We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
    println!("{arc}");

    // Now that we have simulated measurements, we'll run the orbit determination.

    // ===================== //
    // === OD ESTIMATION === //
    // ===================== //

    let sc = SpacecraftUncertainty::builder()
        .nominal(sc_seed)
        .frame(LocalFrame::RIC)
        .x_km(0.5)
        .y_km(0.5)
        .z_km(0.5)
        .vx_km_s(5e-3)
        .vy_km_s(5e-3)
        .vz_km_s(5e-3)
        .build();

    // Build the filter initial estimate, which we will reuse in the filter.
    let initial_estimate = sc.to_estimate()?;

    println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");

    let kf = KF::new(
        // Increase the initial covariance to account for larger deviation.
        initial_estimate,
        // Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame.
        SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]),
    );

    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
    let mut odp = SpacecraftODProcess::ckf(
        setup.with(initial_estimate.state().with_stm(), almanac.clone()),
        kf,
        devices,
        Some(ResidRejectCrit::default()),
        almanac.clone(),
    );

    odp.process_arc(&arc)?;

    let ric_err = traj_as_flown
        .at(odp.estimates.last().unwrap().epoch())?
        .orbit
        .ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
    println!("== RIC at end ==");
    println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
    println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);

    odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;

    // In our case, we have the truth trajectory from NASA.
    // So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
    // Export the OD trajectory first.
    let od_trajectory = odp.to_traj()?;
    // Build the RIC difference.
    od_trajectory.ric_diff_to_parquet(
        &traj_as_flown,
        "./04_lro_od_truth_error.parquet",
        ExportCfg::default(),
    )?;

    Ok(())
}

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

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

impl TrackingDataArc

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

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

pub fn apply_moduli(&mut self)

Applies the moduli to each measurement, if defined.

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

Returns the unique list of aliases in this tracking data arc

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

Returns the unique measurement types in this tracking data arc

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

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

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

Returns the start epoch of this tracking arc

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

Returns the end epoch of this tracking arc

pub fn len(&self) -> usize

Returns the number of measurements in this data arc

pub fn is_empty(&self) -> bool

Returns whether this arc has no measurements.

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

Returns the minimum duration between two subsequent measurements. This is important to correctly set up the propagator and not miss any measurement.

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.

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

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

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

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.

Trait Implementations

impl Clone for TrackingDataArc

fn clone(&self) -> TrackingDataArc

Returns a copy of the value.

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

Performs copy-assignment from source.

impl Debug for TrackingDataArc

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

Formats the value using the given formatter.

impl Default for TrackingDataArc

fn default() -> TrackingDataArc

Returns the “default value” for a type.

impl Display for TrackingDataArc

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

Formats the value using the given formatter.

impl PartialEq for TrackingDataArc

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

Tests for self and other values to be equal, and is used by ==.

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

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.