Struct SolarPressure 

pub struct SolarPressure {
    pub phi: f64,
    pub shadow_model: ShadowModel,
    pub estimate: bool,
}

Computation of solar radiation pressure is based on STK: http://help.agi.com/stk/index.htm#gator/eq-solar.htm .

Fields

phi: f64

solar flux at 1 AU, in W/m^2

shadow_model: ShadowModel

estimate: bool

Set to true to estimate the coefficient of reflectivity

Implementations

impl SolarPressure

pub fn default_flux_raw( shadow_bodies: Vec<Frame>, almanac: &Almanac, ) -> Result<Self, DynamicsError>

Will set the solar flux at 1 AU to: Phi = 1367.0

pub fn default_flux( shadow_body: Frame, almanac: &Almanac, ) -> Result<Arc<Self>, DynamicsError>

Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0

Examples found in repository

nyx-core/examples/03_geo_analysis/stationkeeping.rs (line 94)

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();
    // Set up the dynamics like in the orbit raise.
    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    // Define the GEO orbit, and we're just going to maintain it very tightly.
    let earth_j2000 = almanac.frame_info(EARTH_J2000)?;
    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
    println!("{orbit:x}");

    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();

    // Set up the spacecraft dynamics like in the orbit raise example.

    let prop_time = 30.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,
        ),
    ];

    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
    println!("{ruggiero_ctrl}");

    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);

    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.
    };
    jgm3_meta.process(true)?;

    let harmonics = GravityField::new(GravityFieldData::from_cof(
        &jgm3_meta.uri,
        8,
        8,
        true,
        almanac.frame_info(IAU_EARTH_FRAME)?,
    )?);
    orbital_dyn.accel_models.push(harmonics);

    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;
    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
        .with_guidance_law(ruggiero_ctrl.clone());

    println!("{sc_dynamics}");

    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.

    // Let's start by defining the dispersion.
    // The MultivariateNormal structure allows us to define the dispersions in any of the orbital parameters, but these are applied directly in the Cartesian state space.
    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
    let mc_rv = MvnSpacecraft::new(
        sc,
        vec![StateDispersion::zero_mean(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            3.0,
        )],
    )?;

    let my_mc = MonteCarlo::new(
        sc, // Nominal state
        mc_rv,
        "03_geo_sk".to_string(), // Scenario name
        None, // No specific seed specified, so one will be drawn from the computer's entropy.
    );

    // Build the propagator setup.
    let setup = Propagator::rk89(
        sc_dynamics.clone(),
        IntegratorOptions::builder()
            .min_step(10.0_f64.seconds())
            .error_ctrl(ErrorControl::RSSCartesianStep)
            .build(),
    );

    let num_runs = 25;
    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);

    assert_eq!(rslts.runs.len(), num_runs);

    rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;

    Ok(())
}

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

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/01_orbit_prop/main.rs (line 139)

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(())
}

pub fn default_no_estimation( shadow_bodies: Vec<Frame>, almanac: &Almanac, ) -> Result<Arc<Self>, DynamicsError>

Accounts for the shadowing of only one body and will set the solar flux at 1 AU to: Phi = 1367.0

pub fn with_flux( flux_w_m2: f64, shadow_bodies: Vec<Frame>, almanac: &Almanac, ) -> Result<Arc<Self>, DynamicsError>

Must provide the flux in W/m^2

pub fn new( shadow_bodies: Vec<Frame>, almanac: &Almanac, ) -> Result<Arc<Self>, DynamicsError>

Solar radiation pressure force model accounting for the provided shadow bodies.

Examples found in repository

nyx-core/examples/02_jwst_covar_monte_carlo/main.rs (line 100)

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.
    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.

    // Download the regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
    // Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
    let mut latest_jwst_ephem = MetaFile {
        uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
        crc32: None,
    };
    latest_jwst_ephem.process(true)?;

    // Load this ephem in the general Almanac we're using for this analysis.
    let almanac = Arc::new(
        MetaAlmanac::latest()
            .map_err(Box::new)?
            .load_from_metafile(latest_jwst_ephem, true)?,
    );

    // By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
    // in the BSP. We need this ID in order to query the ephemeris.
    const JWST_NAIF_ID: i32 = -170;
    // Let's build a frame in the J2000 orientation centered on the JWST.
    const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);

    // Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
    let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
    println!("JWST defined from {earliest_epoch} to {latest_epoch}");
    // Fetch the state, printing it in the Earth J2000 frame.
    let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
    println!("{jwst_orbit:x}");

    // Build the spacecraft
    // SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
    // SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
    let jwst = Spacecraft::builder()
        .orbit(jwst_orbit)
        .srp(SRPData {
            area_m2: 21.197 * 14.162,
            coeff_reflectivity: 1.56,
        })
        .mass(Mass::from_dry_mass(6200.0))
        .build();

    // Build up the spacecraft uncertainty builder.
    // We can use the spacecraft uncertainty structure to build this up.
    // We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
    // in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
    // Nyx can also estimate the deviation of the spacecraft parameters.
    let jwst_uncertainty = SpacecraftUncertainty::builder()
        .nominal(jwst)
        .frame(LocalFrame::RIC)
        .x_km(0.5)
        .y_km(0.3)
        .z_km(1.5)
        .vx_km_s(1e-4)
        .vy_km_s(0.6e-3)
        .vz_km_s(3e-3)
        .build();

    println!("{jwst_uncertainty}");

    // Build the Kalman filter estimate.
    // Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
    // but this approach requires quite a bit more boilerplate code.
    let jwst_estimate = jwst_uncertainty.to_estimate()?;

    // Set up the spacecraft dynamics.
    // We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
    // We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.

    let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
    let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], &almanac)?;

    // Finalize setting up the dynamics.
    let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);

    // Build the propagator set up to use for the whole analysis.
    let setup = Propagator::default(dynamics);

    // All of the analysis will use this duration.
    let prediction_duration = 6.5 * Unit::Day;

    // === Covariance mapping ===
    // For the covariance mapping / prediction, we'll use the common orbit determination approach.
    // This is done by setting up a spacecraft Kalman filter OD process, and predicting for the analysis duration.

    // Build the propagation instance for the OD process.
    let odp = SpacecraftKalmanOD::new(
        setup.clone(),
        KalmanVariant::DeviationTracking,
        None,
        BTreeMap::new(),
        almanac.clone(),
    );

    // The prediction step is 1 minute by default, configured in the OD process, i.e. how often we want to know the covariance.
    assert_eq!(odp.max_step, 1_i64.minutes());
    // Finally, predict, and export the trajectory with covariance to a parquet file.
    let od_sol = odp.predict_for(jwst_estimate, prediction_duration)?;
    od_sol.to_parquet("./02_jwst_covar_map.parquet", ExportCfg::default())?;

    // === Monte Carlo framework ===
    // Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.

    let my_mc = MonteCarlo::new(
        jwst, // Nominal state
        jwst_estimate.to_random_variable()?,
        "02_jwst".to_string(), // Scenario name
        None, // No specific seed specified, so one will be drawn from the computer's entropy.
    );

    let num_runs = 5_000;
    let rslts = my_mc.run_until_epoch(
        setup,
        almanac.clone(),
        jwst.epoch() + prediction_duration,
        num_runs,
    );

    assert_eq!(rslts.runs.len(), num_runs);
    // Finally, export these results, computing the eclipse percentage for all of these results.

    rslts.to_parquet("02_jwst_monte_carlo.parquet", ExportCfg::default())?;

    Ok(())
}

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

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)?;
    let arc = trk.generate_measurements(&almanac)?;
    // 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(SigmaRejection::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(())
}

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

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/04_lro_od/main.rs (line 150)

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,
        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)?;
    let arc = trk.generate_measurements(&almanac)?;
    // 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(SigmaRejection::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(())
}

Trait Implementations

impl Clone for SolarPressure

fn clone(&self) -> SolarPressure

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl Debug for SolarPressure

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

Formats the value using the given formatter.

impl Default for SolarPressure

fn default() -> Self

Important: the default will FAIL at runtime if the shadow model is not manually defined with loaded frames.

impl DerefToPyAny for SolarPressure

impl<'de> Deserialize<'de> for SolarPressure

fn deserialize<__D>(__deserializer: __D) -> Result<Self, __D::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl Display for SolarPressure

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

Formats the value using the given formatter.

impl ForceModel for SolarPressure

fn estimation_index(&self) -> Option<usize>

If a parameter of this force model is stored in the spacecraft state, then this function should return the index where this parameter is being affected

fn eom( &self, ctx: &Spacecraft, almanac: &Almanac, ) -> Result<Vector3<f64>, DynamicsError>

Defines the equations of motion for this force model from the provided osculating state.

fn gradient( &self, ctx: &Spacecraft, almanac: &Almanac, ) -> Result<(Vector3<f64>, Matrix4x3<f64>), DynamicsError>

Force models must implement their partials, although those will only be called if the propagation requires the computation of the STM. The osc_ctx is the osculating context, i.e. it changes for each sub-step of the integrator. The last row corresponds to the partials of the parameter of this force model wrt the position, i.e. this only applies to conservative forces.

impl<'a, 'py> FromPyObject<'a, 'py> for SolarPressure

where Self: Clone,

type Error = PyClassGuardError<'a, 'py>

The type returned in the event of a conversion error.

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

Extracts Self from the bound smart pointer obj.

impl<'py> IntoPyObject<'py> for SolarPressure

type Target = SolarPressure

The Python output type

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

The smart pointer type to use.

type Error = PyErr

The type returned in the event of a conversion error.

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

Performs the conversion.

impl PyClass for SolarPressure

const NAME: &str = "SolarPressure"

Name of the class.

type Frozen = False

Whether the pyclass is frozen.

impl PyClassImpl for SolarPressure

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

Module which the class will be associated with.

const IS_BASETYPE: bool = false

#[pyclass(subclass)]

const IS_SUBCLASS: bool = false

#[pyclass(extends=…)]

const IS_MAPPING: bool = false

#[pyclass(mapping)]

const IS_SEQUENCE: bool = false

#[pyclass(sequence)]

const IS_IMMUTABLE_TYPE: bool = false

#[pyclass(immutable_type)]

const RAW_DOC: &'static CStr = /// Computation of solar radiation pressure is based on STK: <http://help.agi.com/stk/index.htm#gator/eq-solar.htm> .

Docstring for the class provided on the struct or enum.

const DOC: &'static CStr

Fully rendered class doc, including the text_signature if a constructor is defined.

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

Description of how this class is laid out in memory

type BaseType = PyAny

Base class

type ThreadChecker = NoopThreadChecker

This handles following two situations:

type Inventory = Pyo3MethodsInventoryForSolarPressure

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

Immutable or mutable

type Dict = PyClassDummySlot

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

type WeakRef = PyClassDummySlot

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

type BaseNativeType = PyAny

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

fn items_iter() -> PyClassItemsIter

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

impl PyClassNewTextSignature for SolarPressure

const TEXT_SIGNATURE: &'static str = "(shadow_bodies, almanac, flux_w_m2=..., estimate=True)"

impl PyTypeInfo for SolarPressure

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.

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.

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.

impl Serialize for SolarPressure

fn serialize<__S>(&self, __serializer: __S) -> Result<__S::Ok, __S::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.

impl StaticType for SolarPressure

where f64: StaticType, ShadowModel: StaticType, bool: StaticType,

fn static_type() -> SimpleType

Return the Dhall type that represents this type.