Struct PropInstance 

pub struct PropInstance<'a, D: Dynamics>
where
    DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,
{
    pub state: D::StateType,
    pub prop: &'a Propagator<D>,
    pub details: IntegrationDetails,
    pub log_progress: bool,
    /* private fields */
}

A Propagator allows propagating a set of dynamics forward or backward in time. It is an EventTracker, without any event tracking. It includes the options, the integrator details of the previous step, and the set of coefficients used for the monomorphic instance.

Fields

state: D::StateType

The state of this propagator instance

prop: &'a Propagator<D>

The propagator setup (kind, stages, etc.)

details: IntegrationDetails

Stores the details of the previous integration step

log_progress: bool

Should progress reports be logged

Implementations

impl<D: Dynamics> PropInstance<'_, D>

where DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,

pub fn until_event( &mut self, max_duration: Duration, event: &Event, event_frame: Option<Frame>, ) -> Result<(D::StateType, Traj<D::StateType>), PropagationError>

where <DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send, D::StateType: Interpolatable,

Propagates the dynamics until the specified event has occurred one, or until max_duration is reached. Refer to [until_nth_event] for details.

pub fn until_nth_event( &mut self, max_duration: Duration, event: &Event, event_frame: Option<Frame>, trigger: usize, ) -> Result<(D::StateType, Traj<D::StateType>), PropagationError>

where <DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send, D::StateType: Interpolatable,

Propagates the dynamics until the specified event has occurred trigger times, or until max_duration is reached.

This method monitors the provided event during propagation. Once the event condition is met trigger number of times (e.g., set trigger to 1 for the first occurrence), the propagation stops at the end of that integration step.

A root-finding algorithm (Brent’s method) is then used to locate the exact time of the event within the final integration step. The returned state corresponds to this precise event time, interpolated from the trajectory.

Arguments

• max_duration - The maximum duration to propagate if the event is not triggered the requested number of times.
• event - The event definition (scalar expression and condition) to monitor.
• trigger - The 1-based index of the event occurrence to stop at (e.g. 1 for the first crossing, 2 for the second).

Returns

A tuple containing:

1. The interpolated state exactly at the moment the $n$-th event occurred.
2. The full trajectory recorded up to the end of the propagation step where the event occurred.

Errors

• PropagationError::NthEventError: Returned if max_duration is reached before the event was triggered trigger times.
• PropagationError::TrajectoryEvent: Returned if the interpolation of the event state fails.
• PropagationError::Analysis: Returned if the event evaluation fails during the search.

impl<D: Dynamics> PropInstance<'_, D>

where DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength>,

pub fn quiet(self) -> Self

Sets this instance to not log progress

pub fn verbose(self) -> Self

Sets this instance to log progress

pub fn set_step(&mut self, step_size: Duration, fixed: bool)

Allows setting the step size of the propagator

pub fn for_duration( &mut self, duration: Duration, ) -> Result<D::StateType, PropagationError>

This method propagates the provided Dynamics for the provided duration.

pub fn for_duration_with_channel( &mut self, duration: Duration, tx_chan: Sender<D::StateType>, ) -> Result<D::StateType, PropagationError>

This method propagates the provided Dynamics for the provided duration and publishes each state on the channel.

pub fn until_epoch( &mut self, end_time: Epoch, ) -> Result<D::StateType, PropagationError>

Propagates the provided Dynamics until the provided epoch. Returns the end state.

pub fn until_epoch_with_channel( &mut self, end_time: Epoch, tx_chan: Sender<D::StateType>, ) -> Result<D::StateType, PropagationError>

Propagates the provided Dynamics until the provided epoch and publishes states on the provided channel. Returns the end state.

pub fn for_duration_with_traj( &mut self, duration: Duration, ) -> Result<(D::StateType, Traj<D::StateType>), PropagationError>

where <DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send, D::StateType: Interpolatable,

Propagates the provided Dynamics for the provided duration and generate the trajectory of these dynamics on its own thread. Returns the end state and the trajectory.

Examples found in repository

nyx-core/examples/03_geo_analysis/raise.rs (line 143)

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

    // Dynamics models require planetary constants and ephemerides to be defined.
    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
    // This will automatically download the DE440s planetary ephemeris,
    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
    // planetary constants kernels.
    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
    // references to many functions.
    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
    // Fetch the EME2000 frame from the Almabac
    let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
    // Define the orbit epoch
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    // Build the spacecraft itself.
    // Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf
    // for the "next gen" SEP characteristics.

    // GTO start
    let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k);

    let sc = Spacecraft::builder()
        .orbit(orbit)
        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
        .srp(SRPData::from_area(3.0 * 6.0)) // Assuming 1 kW/m^2 or 18 kW, giving a margin of 4.35 kW for on-propulsion consumption
        .thruster(Thruster {
            // "NEXT-STEP" row in Table 2
            isp_s: 4435.0,
            thrust_N: 0.472,
        })
        .mode(GuidanceMode::Thrust) // Start thrusting immediately.
        .build();

    let prop_time = 180.0 * Unit::Day;

    // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
    let objectives = &[
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            42_165.0,
            20.0,
        ),
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::Eccentricity),
            0.001,
            5e-5,
        ),
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::Inclination),
            0.05,
            1e-2,
        ),
    ];

    // Ensure that we only thrust if we have more than 20% illumination.
    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2).unwrap();
    println!("{ruggiero_ctrl}");

    // Define the high fidelity dynamics

    // Set up the spacecraft dynamics.

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

    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
    // We're using the JGM3 model here, which is the default in GMAT.
    let mut jgm3_meta = MetaFile {
        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
    };
    // And let's download it if we don't have it yet.
    jgm3_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 Earth centered Earth fixed frame, IAU Earth.
    let harmonics = GravityField::new(
        GravityFieldData::from_cof(
            &jgm3_meta.uri,
            8,
            8,
            true,
            almanac.frame_info(IAU_EARTH_FRAME)?,
        )
        .unwrap(),
    );

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

    // We define the solar radiation pressure, using the default solar flux and accounting only
    // for the eclipsing caused by the Earth.
    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;

    // 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 sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
        .with_guidance_law(ruggiero_ctrl.clone());

    println!("{orbit:x}");

    // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low.
    let (final_state, traj) = Propagator::rk89(
        sc_dynamics.clone(),
        IntegratorOptions::builder()
            .min_step(10.0_f64.seconds())
            .error_ctrl(ErrorControl::RSSCartesianStep)
            .build(),
    )
    .with(sc, almanac.clone())
    .for_duration_with_traj(prop_time)?;

    let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;
    println!("{:x}", final_state.orbit);
    println!("prop usage: {prop_usage:.3} kg");

    // Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise.
    traj.to_parquet("./03_geo_raise.parquet", ExportCfg::default())?;

    for status_line in ruggiero_ctrl.status(&final_state) {
        println!("{status_line}");
    }

    ruggiero_ctrl
        .achieved(&final_state)
        .expect("objective not achieved");

    Ok(())
}

nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 93)

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

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

    let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));

    let out = manifest_dir.join("data/04_output/");

    let almanac = Arc::new(
        Almanac::new(
            &manifest_dir
                .join("data/01_planetary/pck08.pca")
                .to_string_lossy(),
        )
        .unwrap()
        .load(
            &manifest_dir
                .join("data/01_planetary/de440s.bsp")
                .to_string_lossy(),
        )
        .unwrap(),
    );

    let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
    let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();

    let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
    let nrho = Orbit::cartesian(
        166_473.631_302_239_7,
        -274_715.487_253_382_7,
        -211_233.210_176_686_7,
        0.933_451_604_520_018_4,
        0.436_775_046_841_900_9,
        -0.082_211_021_250_348_95,
        epoch,
        eme2k,
    );

    let tx_nrho_sc = Spacecraft::from(nrho);

    let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
    println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");

    let bodies = vec![EARTH, SUN];
    let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));

    let setup = Propagator::rk89(
        dynamics,
        IntegratorOptions::builder().max_step(0.5.minutes()).build(),
    );

    /* == Propagate the NRHO vehicle == */
    let prop_time = 1.1 * state_luna.period().unwrap();

    let (nrho_final, mut tx_traj) = setup
        .with(tx_nrho_sc, almanac.clone())
        .for_duration_with_traj(prop_time)
        .unwrap();

    tx_traj.name = Some("NRHO Tx SC".to_string());

    println!("{tx_traj}");

    /* == Propagate an LLO vehicle == */
    let llo_orbit =
        Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();

    let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();

    let (_, llo_traj) = setup
        .with(llo_sc, almanac.clone())
        .until_epoch_with_traj(nrho_final.epoch())
        .unwrap();

    // Export the subset of the first two hours.
    llo_traj
        .clone()
        .filter_by_offset(..2.hours())
        .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;

    /* == Setup the interlink == */

    let mut measurement_types = IndexSet::new();
    measurement_types.insert(MeasurementType::Range);
    measurement_types.insert(MeasurementType::Doppler);

    let mut stochastics = IndexMap::new();

    let sa45_csac_allan_dev = 1e-11;

    stochastics.insert(
        MeasurementType::Range,
        StochasticNoise::from_hardware_range_km(
            sa45_csac_allan_dev,
            10.0.seconds(),
            link_specific::ChipRate::StandardT4B(),
            link_specific::SN0::Average(),
        ),
    );

    stochastics.insert(
        MeasurementType::Doppler,
        StochasticNoise::from_hardware_doppler_km_s(
            sa45_csac_allan_dev,
            10.0.seconds(),
            link_specific::CarrierFreq::SBand(),
            link_specific::CN0::Average(),
        ),
    );

    let interlink = InterlinkTxSpacecraft {
        traj: tx_traj,
        measurement_types,
        integration_time: None,
        timestamp_noise_s: None,
        ab_corr: Aberration::LT,
        stochastic_noises: Some(stochastics),
    };

    // Devices are the transmitter, which is our NRHO vehicle.
    let mut devices = BTreeMap::new();
    devices.insert("NRHO Tx SC".to_string(), interlink);

    let mut configs = BTreeMap::new();
    configs.insert(
        "NRHO Tx SC".to_string(),
        TrkConfig::builder()
            .strands(vec![Strand {
                start: epoch,
                end: nrho_final.epoch(),
            }])
            .build(),
    );

    let mut trk_sim =
        TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
    println!("{trk_sim}");

    let trk_data = trk_sim.generate_measurements(almanac.clone()).unwrap();
    println!("{trk_data}");

    trk_data
        .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
        .unwrap();

    // Run a truth OD where we estimate the LLO position
    let llo_uncertainty = SpacecraftUncertainty::builder()
        .nominal(llo_sc)
        .x_km(1.0)
        .y_km(1.0)
        .z_km(1.0)
        .vx_km_s(1e-3)
        .vy_km_s(1e-3)
        .vz_km_s(1e-3)
        .build();

    let mut proc_devices = devices.clone();

    // Define the initial estimate, randomized, seed for reproducibility
    let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
    // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
    initial_estimate.covar *= 2.5;

    // Increase the noise in the devices to accept more measurements.

    for link in proc_devices.values_mut() {
        for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
            *noise.white_noise.as_mut().unwrap() *= 3.0;
        }
    }

    let init_err = initial_estimate
        .orbital_state()
        .ric_difference(&llo_orbit)
        .unwrap();

    println!("initial estimate:\n{initial_estimate}");
    println!("RIC errors = {init_err}",);

    let odp = InterlinkKalmanOD::new(
        setup.clone(),
        KalmanVariant::ReferenceUpdate,
        Some(ResidRejectCrit::default()),
        proc_devices,
        almanac.clone(),
    );

    // Shrink the data to process.
    let arc = trk_data.filter_by_offset(..2.hours());

    let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();

    println!("{od_sol}");

    od_sol
        .to_parquet(
            out.join("05_caps_interlink_od_sol.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let od_traj = od_sol.to_traj().unwrap();

    od_traj
        .ric_diff_to_parquet(
            &llo_traj,
            out.join("05_caps_interlink_llo_est_error.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let final_est = od_sol.estimates.last().unwrap();
    assert!(final_est.within_3sigma(), "should be within 3 sigma");

    println!("ESTIMATE\n{final_est:x}\n");
    let truth = llo_traj.at(final_est.epoch()).unwrap();
    println!("TRUTH\n{truth:x}");

    let final_err = truth
        .orbit
        .ric_difference(&final_est.orbital_state())
        .unwrap();
    println!("ERROR {final_err}");

    // Build the residuals versus reference plot.
    let rvr_sol = odp
        .process_arc(initial_estimate, &arc.resid_vs_ref_check())
        .unwrap();

    rvr_sol
        .to_parquet(
            out.join("05_caps_interlink_resid_v_ref.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let final_rvr = rvr_sol.estimates.last().unwrap();

    println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
    println!(
        "Pure prop error {:.3} m",
        final_rvr
            .orbital_state()
            .ric_difference(&final_est.orbital_state())
            .unwrap()
            .rmag_km()
            * 1e3
    );

    Ok(())
}

nyx-core/examples/06_lunar_orbit_determination/main.rs (line 131)

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",
        "06_lunar_orbit_determination",
    ]
    .iter()
    .collect();

    let meta = data_folder.join("metaalmanac.dhall");

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

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

    // Build a nominal trajectory
    // TODO: Switch this to a sequence once the OD over a spacecraft sequence is implemented.

    let epoch = Epoch::from_gregorian_utc_at_noon(2024, 2, 29);
    let moon_j2000 = almanac.frame_info(MOON_J2000)?;

    // To build the trajectory we need to provide a spacecraft template.
    let orbiter = Spacecraft::builder()
        .mass(Mass::from_dry_and_prop_masses(1018.0, 900.0))
        .srp(SRPData {
            area_m2: 3.9 * 2.7,
            coeff_reflectivity: 0.96,
        })
        .orbit(Orbit::try_keplerian_altitude(
            150.0, 0.00212, 33.6, 45.0, 45.0, 0.0, epoch, moon_j2000,
        )?) // Setting a zero orbit here because it's just a template
        .build();

    // ========================== //
    // === BUILD NOMINAL TRAJ === //
    // ========================== //

    // 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 = GravityField::new(GravityFieldData::from_shadr(
        &jggrx_meta.uri,
        80,
        80,
        true,
        almanac.frame_info(moon_pa_frame)?,
    )?);

    // 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![MOON_J2000], &almanac)?;

    // 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}");

    let setup = Propagator::rk89(dynamics.clone(), IntegratorOptions::default());

    let truth_traj = setup
        .with(orbiter, almanac.clone())
        .for_duration_with_traj(Unit::Day * 2)?
        .1;

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

    // 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 = data_folder.join("dsn-network.yaml");
    let devices = GroundStation::load_named(ground_station_file)?;

    let proc_devices = devices.clone();

    // 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 configs: BTreeMap<String, TrkConfig> =
        TrkConfig::load_named(data_folder.join("tracking-cfg.yaml"))?;

    // Build the tracking arc simulation to generate a "standard measurement".
    let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::with_seed(
        devices.clone(),
        truth_traj.clone(),
        configs,
        123, // Set a seed for reproducibility
    )?;

    trk.build_schedule(almanac.clone())?;
    let arc = trk.generate_measurements(almanac.clone())?;
    // Save the simulated tracking data
    arc.to_parquet_simple("./data/04_output/06_lunar_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(orbiter)
        .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{orbiter:x}\n{initial_estimate}");

    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
    let process_noise = ProcessNoise3D::from_velocity_km_s(
        &[1e-14, 1e-14, 1e-14],
        1 * Unit::Hour,
        10 * Unit::Minute,
        None,
    );

    println!("{process_noise}");

    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
    let odp = SpacecraftKalmanScalarOD::new(
        setup,
        KalmanVariant::ReferenceUpdate,
        Some(ResidRejectCrit::default()),
        proc_devices,
        almanac.clone(),
    )
    .with_process_noise(process_noise);

    let od_sol = odp.process_arc(initial_estimate, &arc)?;

    let final_est = od_sol.estimates.last().unwrap();

    println!("{final_est}");

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

    println!(
        "Num residuals rejected: #{}",
        od_sol.rejected_residuals().len()
    );
    println!(
        "Percentage within +/-3: {}",
        od_sol.residual_ratio_within_threshold(3.0).unwrap()
    );
    println!("Whitened residuals normal? {}", od_sol.is_normal(None)?);
    println!("NIS test success? {}", od_sol.is_nis_consistent(None)?);

    od_sol.to_parquet(
        "./data/04_output/06_lunar_od_results.parquet",
        ExportCfg::default(),
    )?;

    let od_trajectory = od_sol.to_traj()?;
    // Build the RIC difference.
    od_trajectory.ric_diff_to_parquet(
        &truth_traj,
        "./data/04_output/06_lunar_od_truth_error.parquet",
        ExportCfg::default(),
    )?;

    Ok(())
}

pub fn until_epoch_with_traj( &mut self, end_time: Epoch, ) -> Result<(D::StateType, Traj<D::StateType>), PropagationError>

where <DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send, D::StateType: Interpolatable,

Propagates the provided Dynamics until the provided epoch and generate the trajectory of these dynamics on its own thread. Returns the end state and the trajectory.

Examples found in repository

nyx-core/examples/05_cislunar_spacecraft_link_od/main.rs (line 108)

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

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

    let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));

    let out = manifest_dir.join("data/04_output/");

    let almanac = Arc::new(
        Almanac::new(
            &manifest_dir
                .join("data/01_planetary/pck08.pca")
                .to_string_lossy(),
        )
        .unwrap()
        .load(
            &manifest_dir
                .join("data/01_planetary/de440s.bsp")
                .to_string_lossy(),
        )
        .unwrap(),
    );

    let eme2k = almanac.frame_info(EARTH_J2000).unwrap();
    let moon_iau = almanac.frame_info(IAU_MOON_FRAME).unwrap();

    let epoch = Epoch::from_gregorian_tai(2021, 5, 29, 19, 51, 16, 852_000);
    let nrho = Orbit::cartesian(
        166_473.631_302_239_7,
        -274_715.487_253_382_7,
        -211_233.210_176_686_7,
        0.933_451_604_520_018_4,
        0.436_775_046_841_900_9,
        -0.082_211_021_250_348_95,
        epoch,
        eme2k,
    );

    let tx_nrho_sc = Spacecraft::from(nrho);

    let state_luna = almanac.transform_to(nrho, MOON_J2000, None).unwrap();
    println!("Start state (dynamics: Earth, Moon, Sun gravity):\n{state_luna}");

    let bodies = vec![EARTH, SUN];
    let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(bodies));

    let setup = Propagator::rk89(
        dynamics,
        IntegratorOptions::builder().max_step(0.5.minutes()).build(),
    );

    /* == Propagate the NRHO vehicle == */
    let prop_time = 1.1 * state_luna.period().unwrap();

    let (nrho_final, mut tx_traj) = setup
        .with(tx_nrho_sc, almanac.clone())
        .for_duration_with_traj(prop_time)
        .unwrap();

    tx_traj.name = Some("NRHO Tx SC".to_string());

    println!("{tx_traj}");

    /* == Propagate an LLO vehicle == */
    let llo_orbit =
        Orbit::try_keplerian_altitude(110.0, 1e-4, 90.0, 0.0, 0.0, 0.0, epoch, moon_iau).unwrap();

    let llo_sc = Spacecraft::builder().orbit(llo_orbit).build();

    let (_, llo_traj) = setup
        .with(llo_sc, almanac.clone())
        .until_epoch_with_traj(nrho_final.epoch())
        .unwrap();

    // Export the subset of the first two hours.
    llo_traj
        .clone()
        .filter_by_offset(..2.hours())
        .to_parquet_simple(out.join("05_caps_llo_truth.pq"))?;

    /* == Setup the interlink == */

    let mut measurement_types = IndexSet::new();
    measurement_types.insert(MeasurementType::Range);
    measurement_types.insert(MeasurementType::Doppler);

    let mut stochastics = IndexMap::new();

    let sa45_csac_allan_dev = 1e-11;

    stochastics.insert(
        MeasurementType::Range,
        StochasticNoise::from_hardware_range_km(
            sa45_csac_allan_dev,
            10.0.seconds(),
            link_specific::ChipRate::StandardT4B(),
            link_specific::SN0::Average(),
        ),
    );

    stochastics.insert(
        MeasurementType::Doppler,
        StochasticNoise::from_hardware_doppler_km_s(
            sa45_csac_allan_dev,
            10.0.seconds(),
            link_specific::CarrierFreq::SBand(),
            link_specific::CN0::Average(),
        ),
    );

    let interlink = InterlinkTxSpacecraft {
        traj: tx_traj,
        measurement_types,
        integration_time: None,
        timestamp_noise_s: None,
        ab_corr: Aberration::LT,
        stochastic_noises: Some(stochastics),
    };

    // Devices are the transmitter, which is our NRHO vehicle.
    let mut devices = BTreeMap::new();
    devices.insert("NRHO Tx SC".to_string(), interlink);

    let mut configs = BTreeMap::new();
    configs.insert(
        "NRHO Tx SC".to_string(),
        TrkConfig::builder()
            .strands(vec![Strand {
                start: epoch,
                end: nrho_final.epoch(),
            }])
            .build(),
    );

    let mut trk_sim =
        TrackingArcSim::with_seed(devices.clone(), llo_traj.clone(), configs, 0).unwrap();
    println!("{trk_sim}");

    let trk_data = trk_sim.generate_measurements(almanac.clone()).unwrap();
    println!("{trk_data}");

    trk_data
        .to_parquet_simple(out.clone().join("nrho_interlink_msr.pq"))
        .unwrap();

    // Run a truth OD where we estimate the LLO position
    let llo_uncertainty = SpacecraftUncertainty::builder()
        .nominal(llo_sc)
        .x_km(1.0)
        .y_km(1.0)
        .z_km(1.0)
        .vx_km_s(1e-3)
        .vy_km_s(1e-3)
        .vz_km_s(1e-3)
        .build();

    let mut proc_devices = devices.clone();

    // Define the initial estimate, randomized, seed for reproducibility
    let mut initial_estimate = llo_uncertainty.to_estimate_randomized(Some(0)).unwrap();
    // Inflate the covariance -- https://github.com/nyx-space/nyx/issues/339
    initial_estimate.covar *= 2.5;

    // Increase the noise in the devices to accept more measurements.

    for link in proc_devices.values_mut() {
        for noise in &mut link.stochastic_noises.as_mut().unwrap().values_mut() {
            *noise.white_noise.as_mut().unwrap() *= 3.0;
        }
    }

    let init_err = initial_estimate
        .orbital_state()
        .ric_difference(&llo_orbit)
        .unwrap();

    println!("initial estimate:\n{initial_estimate}");
    println!("RIC errors = {init_err}",);

    let odp = InterlinkKalmanOD::new(
        setup.clone(),
        KalmanVariant::ReferenceUpdate,
        Some(ResidRejectCrit::default()),
        proc_devices,
        almanac.clone(),
    );

    // Shrink the data to process.
    let arc = trk_data.filter_by_offset(..2.hours());

    let od_sol = odp.process_arc(initial_estimate, &arc).unwrap();

    println!("{od_sol}");

    od_sol
        .to_parquet(
            out.join("05_caps_interlink_od_sol.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let od_traj = od_sol.to_traj().unwrap();

    od_traj
        .ric_diff_to_parquet(
            &llo_traj,
            out.join("05_caps_interlink_llo_est_error.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let final_est = od_sol.estimates.last().unwrap();
    assert!(final_est.within_3sigma(), "should be within 3 sigma");

    println!("ESTIMATE\n{final_est:x}\n");
    let truth = llo_traj.at(final_est.epoch()).unwrap();
    println!("TRUTH\n{truth:x}");

    let final_err = truth
        .orbit
        .ric_difference(&final_est.orbital_state())
        .unwrap();
    println!("ERROR {final_err}");

    // Build the residuals versus reference plot.
    let rvr_sol = odp
        .process_arc(initial_estimate, &arc.resid_vs_ref_check())
        .unwrap();

    rvr_sol
        .to_parquet(
            out.join("05_caps_interlink_resid_v_ref.pq"),
            ExportCfg::default(),
        )
        .unwrap();

    let final_rvr = rvr_sol.estimates.last().unwrap();

    println!("RMAG error {:.3} m", final_err.rmag_km() * 1e3);
    println!(
        "Pure prop error {:.3} m",
        final_rvr
            .orbital_state()
            .ric_difference(&final_est.orbital_state())
            .unwrap()
            .rmag_km()
            * 1e3
    );

    Ok(())
}

nyx-core/examples/03_geo_analysis/drift.rs (line 119)

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();
    // Dynamics models require planetary constants and ephemerides to be defined.
    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
    // This will automatically download the DE440s planetary ephemeris,
    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
    // planetary constants kernels.
    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
    // references to many functions.
    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
    // Define the orbit epoch
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    // Define the orbit.
    // First we need to fetch the Earth J2000 from information from the Almanac.
    // This allows the frame to include the gravitational parameters and the shape of the Earth,
    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
    // by loading a different set of planetary constants.
    let earth_j2000 = almanac.frame_info(EARTH_J2000)?;

    // Placing this GEO bird just above Colorado.
    // In theory, the eccentricity is zero, but in practice, it's about 1e-5 to 1e-6 at best.
    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
    // Print in in Keplerian form.
    println!("{orbit:x}");

    let state_bf = almanac.transform_to(orbit, IAU_EARTH_FRAME, None)?;
    let (orig_lat_deg, orig_long_deg, orig_alt_km) = state_bf.latlongalt()?;

    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
    // models such as solar radiation pressure.

    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
    let sc = Spacecraft::builder()
        .orbit(orbit)
        .mass(Mass::from_dry_mass(9.60))
        .srp(SRPData {
            area_m2: 10e-4,
            coeff_reflectivity: 1.1,
        })
        .build();
    println!("{sc:x}");

    // Set up the spacecraft dynamics.

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

    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
    // We're using the JGM3 model here, which is the default in GMAT.
    let mut jgm3_meta = MetaFile {
        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
    };
    // And let's download it if we don't have it yet.
    jgm3_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 Earth centered Earth fixed frame, IAU Earth.
    let harmonics_21x21 = GravityField::new(
        GravityFieldData::from_cof(
            &jgm3_meta.uri,
            21,
            21,
            true,
            almanac.frame_info(IAU_EARTH_FRAME)?,
        )
        .unwrap(),
    );

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

    // We define the solar radiation pressure, using the default solar flux and accounting only
    // for the eclipsing caused by the Earth and Moon.
    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;

    // 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}");

    // Finally, let's propagate this orbit to the same epoch as above.
    // The first returned value is the spacecraft state at the final epoch.
    // The second value is the full trajectory where the step size is variable step used by the propagator.
    let (future_sc, trajectory) = Propagator::default(dynamics)
        .with(sc, almanac.clone())
        .until_epoch_with_traj(epoch + Unit::Century * 0.03)?;

    println!("=== High fidelity propagation ===");
    println!(
        "SMA changed by {:.3} km",
        orbit.sma_km()? - future_sc.orbit.sma_km()?
    );
    println!(
        "ECC changed by {:.6}",
        orbit.ecc()? - future_sc.orbit.ecc()?
    );
    println!(
        "INC changed by {:.3e} deg",
        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
    );
    println!(
        "RAAN changed by {:.3} deg",
        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
    );
    println!(
        "AOP changed by {:.3} deg",
        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
    );
    println!(
        "TA changed by {:.3} deg",
        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
    );

    // We also have access to the full trajectory throughout the propagation.
    println!("{trajectory}");

    println!("Spacecraft params after 3 years without active control:\n{future_sc:x}");

    // With the trajectory, let's build a few data products.

    // 1. Export the trajectory as a parquet file, which includes the Keplerian orbital elements.

    let analysis_step = Unit::Minute * 5;

    trajectory.to_parquet(
        "./03_geo_hf_prop.parquet",
        ExportCfg::builder().step(analysis_step).build(),
    )?;

    // 2. Compute the latitude, longitude, and altitude throughout the trajectory by rotating the spacecraft position into the Earth body fixed frame.

    // We iterate over the trajectory, grabbing a state every two minutes.
    let mut offset_s = vec![];
    let mut epoch_str = vec![];
    let mut longitude_deg = vec![];
    let mut latitude_deg = vec![];
    let mut altitude_km = vec![];

    for state in trajectory.every(analysis_step) {
        // Convert the GEO bird state into the body fixed frame, and keep track of its latitude, longitude, and altitude.
        // These define the GEO stationkeeping box.

        let this_epoch = state.epoch();

        offset_s.push((this_epoch - orbit.epoch).to_seconds());
        epoch_str.push(this_epoch.to_isoformat());

        let state_bf = almanac.transform_to(state.orbit, IAU_EARTH_FRAME, None)?;
        let (lat_deg, long_deg, alt_km) = state_bf.latlongalt()?;
        longitude_deg.push(long_deg);
        latitude_deg.push(lat_deg);
        altitude_km.push(alt_km);
    }

    println!(
        "Longitude changed by {:.3} deg -- Box is 0.1 deg E-W",
        orig_long_deg - longitude_deg.last().unwrap()
    );

    println!(
        "Latitude changed by {:.3} deg -- Box is 0.05 deg N-S",
        orig_lat_deg - latitude_deg.last().unwrap()
    );

    println!(
        "Altitude changed by {:.3} km -- Box is 30 km",
        orig_alt_km - altitude_km.last().unwrap()
    );

    // Build the station keeping data frame.
    let mut sk_df = df!(
        "Offset (s)" => offset_s.clone(),
        "Epoch (UTC)" => epoch_str.clone(),
        "Longitude E-W (deg)" => longitude_deg,
        "Latitude N-S (deg)" => latitude_deg,
        "Altitude (km)" => altitude_km,

    )?;

    // Create a file to write the Parquet to
    let file = File::create("./03_geo_lla.parquet").expect("Could not create file");

    // Create a ParquetWriter and write the DataFrame to the file
    ParquetWriter::new(file).finish(&mut sk_df)?;

    Ok(())
}

nyx-core/examples/01_orbit_prop/main.rs (line 152)

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();
    // Dynamics models require planetary constants and ephemerides to be defined.
    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
    // This will automatically download the DE440s planetary ephemeris,
    // the daily-updated Earth Orientation Parameters, the high fidelity Moon orientation
    // parameters (for the Moon Mean Earth and Moon Principal Axes frames), and the PCK11
    // planetary constants kernels.
    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
    // Note that we place the Almanac into an Arc so we can clone it cheaply and provide read-only
    // references to many functions.
    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
    // Define the orbit epoch
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    // Define the orbit.
    // First we need to fetch the Earth J2000 from information from the Almanac.
    // This allows the frame to include the gravitational parameters and the shape of the Earth,
    // defined as a tri-axial ellipoid. Note that this shape can be changed manually or in the Almanac
    // by loading a different set of planetary constants.
    let earth_j2000 = almanac.frame_info(EARTH_J2000)?;

    let orbit =
        Orbit::try_keplerian_altitude(300.0, 0.015, 68.5, 65.2, 75.0, 0.0, epoch, earth_j2000)?;
    // Print in in Keplerian form.
    println!("{orbit:x}");

    // There are two ways to propagate an orbit. We can make a quick approximation assuming only two-body
    // motion. This is a useful first order approximation but it isn't used in real-world applications.

    // This approach is a feature of ANISE.
    let future_orbit_tb = orbit.at_epoch(epoch + Unit::Day * 3)?;
    println!("{future_orbit_tb:x}");

    // Two body propagation relies solely on Kepler's laws, so only the true anomaly will change.
    println!(
        "SMA changed by {:.3e} km",
        orbit.sma_km()? - future_orbit_tb.sma_km()?
    );
    println!(
        "ECC changed by {:.3e}",
        orbit.ecc()? - future_orbit_tb.ecc()?
    );
    println!(
        "INC changed by {:.3e} deg",
        orbit.inc_deg()? - future_orbit_tb.inc_deg()?
    );
    println!(
        "RAAN changed by {:.3e} deg",
        orbit.raan_deg()? - future_orbit_tb.raan_deg()?
    );
    println!(
        "AOP changed by {:.3e} deg",
        orbit.aop_deg()? - future_orbit_tb.aop_deg()?
    );
    println!(
        "TA changed by {:.3} deg",
        orbit.ta_deg()? - future_orbit_tb.ta_deg()?
    );

    // Nyx is used for high fidelity propagation, not Keplerian propagation as above.
    // Nyx only propagates Spacecraft at the moment, which allows it to account for acceleration
    // models such as solar radiation pressure.

    // Let's build a cubesat sized spacecraft, with an SRP area of 10 cm^2 and a mass of 9.6 kg.
    let sc = Spacecraft::builder()
        .orbit(orbit)
        .mass(Mass::from_dry_mass(9.60))
        .srp(SRPData {
            area_m2: 10e-4,
            coeff_reflectivity: 1.1,
        })
        .build();
    println!("{sc:x}");

    // Set up the spacecraft dynamics.

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

    // We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
    // We're using the JGM3 model here, which is the default in GMAT.
    let mut jgm3_meta = MetaFile {
        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
    };
    // And let's download it if we don't have it yet.
    jgm3_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 Earth centered Earth fixed frame, IAU Earth.
    let harmonics_21x21 = GravityField::new(
        GravityFieldData::from_cof(
            &jgm3_meta.uri,
            21,
            21,
            true,
            almanac.frame_info(IAU_EARTH_FRAME)?,
        )
        .unwrap(),
    );

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

    // We define the solar radiation pressure, using the default solar flux and accounting only
    // for the eclipsing caused by the Earth.
    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;

    // 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}");

    // Finally, let's propagate this orbit to the same epoch as above.
    // The first returned value is the spacecraft state at the final epoch.
    // The second value is the full trajectory where the step size is variable step used by the propagator.
    let (future_sc, trajectory) = Propagator::default(dynamics)
        .with(sc, almanac.clone())
        .until_epoch_with_traj(future_orbit_tb.epoch)?;

    println!("=== High fidelity propagation ===");
    println!(
        "SMA changed by {:.3} km",
        orbit.sma_km()? - future_sc.orbit.sma_km()?
    );
    println!(
        "ECC changed by {:.6}",
        orbit.ecc()? - future_sc.orbit.ecc()?
    );
    println!(
        "INC changed by {:.3e} deg",
        orbit.inc_deg()? - future_sc.orbit.inc_deg()?
    );
    println!(
        "RAAN changed by {:.3} deg",
        orbit.raan_deg()? - future_sc.orbit.raan_deg()?
    );
    println!(
        "AOP changed by {:.3} deg",
        orbit.aop_deg()? - future_sc.orbit.aop_deg()?
    );
    println!(
        "TA changed by {:.3} deg",
        orbit.ta_deg()? - future_sc.orbit.ta_deg()?
    );

    // We also have access to the full trajectory throughout the propagation.
    println!("{trajectory}");

    // With the trajectory, let's build a few data products.

    // 1. Export the trajectory as a CCSDS OEM version 2.0 file and as a parquet file, which includes the Keplerian orbital elements.

    trajectory.to_oem_file(
        "./01_cubesat_hf_prop.oem",
        "CUBESAT-ID".to_string(),
        Some("Nyx Space".to_string()),
        Some("CUBESAT".to_string()),
        ExportCfg::builder().step(Unit::Minute * 2).build(),
    )?;

    trajectory.to_parquet_with_cfg(
        "./01_cubesat_hf_prop.parquet",
        ExportCfg::builder().step(Unit::Minute * 2).build(),
    )?;

    // 2. Compare the difference in the radial-intrack-crosstrack frame between the high fidelity
    // and Keplerian propagation. The RIC frame is commonly used to compute the difference in position
    // and velocity of different spacecraft.
    // 3. Compute the azimuth, elevation, range, and range-rate data of that spacecraft as seen from Boulder, CO, USA.

    let boulder_station = GroundStation::from_point(
        "Boulder, CO, USA".to_string(),
        40.014984,   // latitude in degrees
        -105.270546, // longitude in degrees
        1.6550,      // altitude in kilometers
        almanac.frame_info(IAU_EARTH_FRAME)?,
    );

    // We iterate over the trajectory, grabbing a state every two minutes.
    let mut offset_s = vec![];
    let mut epoch_str = vec![];
    let mut ric_x_km = vec![];
    let mut ric_y_km = vec![];
    let mut ric_z_km = vec![];
    let mut ric_vx_km_s = vec![];
    let mut ric_vy_km_s = vec![];
    let mut ric_vz_km_s = vec![];

    let mut azimuth_deg = vec![];
    let mut elevation_deg = vec![];
    let mut range_km = vec![];
    let mut range_rate_km_s = vec![];
    for state in trajectory.every(Unit::Minute * 2) {
        // Try to compute the Keplerian/two body state just in time.
        // This method occasionally fails to converge on an appropriate true anomaly
        // from the mean anomaly. If that happens, we just skip this state.
        // The high fidelity and Keplerian states diverge continuously, and we're curious
        // about the divergence in this quick analysis.
        let this_epoch = state.epoch();
        match orbit.at_epoch(this_epoch) {
            Ok(tb_then) => {
                offset_s.push((this_epoch - orbit.epoch).to_seconds());
                epoch_str.push(format!("{this_epoch}"));
                // Compute the two body state just in time.
                let ric = state.orbit.ric_difference(&tb_then)?;
                ric_x_km.push(ric.radius_km.x);
                ric_y_km.push(ric.radius_km.y);
                ric_z_km.push(ric.radius_km.z);
                ric_vx_km_s.push(ric.velocity_km_s.x);
                ric_vy_km_s.push(ric.velocity_km_s.y);
                ric_vz_km_s.push(ric.velocity_km_s.z);

                // Compute the AER data for each state.
                let aer = almanac.azimuth_elevation_range_sez(
                    state.orbit,
                    boulder_station.to_orbit(this_epoch, &almanac)?,
                    None,
                    None,
                )?;
                azimuth_deg.push(aer.azimuth_deg);
                elevation_deg.push(aer.elevation_deg);
                range_km.push(aer.range_km);
                range_rate_km_s.push(aer.range_rate_km_s);
            }
            Err(e) => warn!("{} {e}", state.epoch()),
        };
    }

    // Build the data frames.
    let ric_df = df!(
        "Offset (s)" => offset_s.clone(),
        "Epoch" => epoch_str.clone(),
        "RIC X (km)" => ric_x_km,
        "RIC Y (km)" => ric_y_km,
        "RIC Z (km)" => ric_z_km,
        "RIC VX (km/s)" => ric_vx_km_s,
        "RIC VY (km/s)" => ric_vy_km_s,
        "RIC VZ (km/s)" => ric_vz_km_s,
    )?;

    println!("RIC difference at start\n{}", ric_df.head(Some(10)));
    println!("RIC difference at end\n{}", ric_df.tail(Some(10)));

    let aer_df = df!(
        "Offset (s)" => offset_s.clone(),
        "Epoch" => epoch_str.clone(),
        "azimuth (deg)" => azimuth_deg,
        "elevation (deg)" => elevation_deg,
        "range (km)" => range_km,
        "range rate (km/s)" => range_rate_km_s,
    )?;

    // Finally, let's see when the spacecraft is visible, assuming 15 degrees minimum elevation.
    let mask = aer_df
        .column("elevation (deg)")?
        .gt(&Column::Scalar(ScalarColumn::new(
            "elevation mask (deg)".into(),
            Scalar::new(DataType::Float64, AnyValue::Float64(15.0)),
            offset_s.len(),
        )))?;
    let cubesat_visible = aer_df.filter(&mask)?;

    println!("{cubesat_visible}");

    Ok(())
}

nyx-core/examples/04_lro_od/main.rs (line 164)

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 output_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "../data", "04_output"]
        .iter()
        .collect();

    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().as_ref())
        .map_err(Box::new)?
        .process(true)
        .map_err(Box::new)?;

    let mut moon_pc = almanac.get_planetary_data_from_id(MOON).unwrap();
    moon_pc.mu_km3_s2 = 4902.74987;
    almanac.set_planetary_data_from_id(MOON, moon_pc).unwrap();

    let mut earth = almanac.get_planetary_data_from_id(EARTH).unwrap();
    earth.mu_km3_s2 = 398600.436;
    almanac.set_planetary_data_from_id(EARTH, earth).unwrap();

    // 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
        .values()
        .next()
        .unwrap()
        .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 = GravityField::new(GravityFieldData::from_shadr(
        &jggrx_meta.uri,
        80,
        80,
        true,
        almanac.frame_info(moon_pa_frame)?,
    )?);

    // 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)?;

    // 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,
        output_folder.join("./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)?;

    let mut proc_devices = devices.clone();

    // Increase the noise in the devices to accept more measurements.
    for gs in proc_devices.values_mut() {
        if let Some(noise) = &mut gs
            .stochastic_noises
            .as_mut()
            .unwrap()
            .get_mut(&MeasurementType::Range)
        {
            *noise.white_noise.as_mut().unwrap() *= 3.0;
        }
    }

    // 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>::with_seed(
        devices.clone(),
        traj_as_flown.clone(),
        configs,
        123, // Set a seed for reproducibility
    )?;

    trk.build_schedule(almanac.clone())?;
    let arc = trk.generate_measurements(almanac.clone())?;
    // Save the simulated tracking data
    arc.to_parquet_simple(output_folder.join("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 mut initial_estimate = sc.to_estimate()?;
    initial_estimate.covar *= 3.0;

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

    // Build the SNC in the Moon J2000 frame, specified as a velocity noise over time.
    let process_noise = ProcessNoise3D::from_velocity_km_s(
        &[1e-12, 1e-12, 1e-12],
        1 * Unit::Hour,
        10 * Unit::Minute,
        None,
    );

    println!("{process_noise}");

    // We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
    let odp = SpacecraftKalmanOD::new(
        setup,
        KalmanVariant::ReferenceUpdate,
        Some(ResidRejectCrit::default()),
        proc_devices,
        almanac.clone(),
    )
    .with_process_noise(process_noise);

    let od_sol = odp.process_arc(initial_estimate, &arc)?;

    let final_est = od_sol.estimates.last().unwrap();

    println!("{final_est}");

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

    println!(
        "Num residuals rejected: #{}",
        od_sol.rejected_residuals().len()
    );
    println!(
        "Percentage within +/-3: {}",
        od_sol.residual_ratio_within_threshold(3.0).unwrap()
    );
    println!("Ratios normal? {}", od_sol.is_normal(None).unwrap());

    od_sol.to_parquet(
        output_folder.join("04_lro_od_results.parquet"),
        ExportCfg::default(),
    )?;

    // Create the ephemeris
    let ephem = od_sol.to_ephemeris("LRO rebuilt".to_string());
    let ephem_start = ephem.start_epoch().unwrap();
    let ephem_end = ephem.end_epoch().unwrap();
    // Check that the covariance is PSD throughout the ephemeris by interpolating it.
    for epoch in TimeSeries::inclusive(ephem_start, ephem_end, Unit::Minute * 5) {
        ephem
            .covar_at(
                epoch,
                anise::ephemerides::ephemeris::LocalFrame::RIC,
                &almanac,
            )
            .unwrap_or_else(|e| panic!("covar not PSD at {epoch}: {e}"));
    }
    // Export as BSP!
    ephem
        .write_spice_bsp(
            -85,
            output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap(),
            None,
        )
        .expect("could not built BSP");
    let new_almanac = Almanac::default()
        .load(output_folder.join("04_lro_rebuilt.bsp").to_str().unwrap())
        .unwrap();
    new_almanac.describe(None, None, None, None, None, None, None, None);
    let (spk_start, spk_end) = new_almanac.spk_domain(-85).unwrap();

    assert!((ephem_start - spk_start).abs() < Unit::Microsecond * 1);
    assert!((ephem_end - spk_end).abs() < Unit::Microsecond * 1);

    // 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 = od_sol.to_traj()?;
    // Build the RIC difference.
    od_trajectory.ric_diff_to_parquet(
        &traj_as_flown,
        output_folder.join("04_lro_od_truth_error.parquet"),
        ExportCfg::default(),
    )?;

    Ok(())
}

pub fn single_step(&mut self) -> Result<(), PropagationError>

Take a single propagator step and emit the result on the TX channel (if enabled)

pub fn latest_details(&self) -> IntegrationDetails

Copy the details of the latest integration step.