Struct SpacecraftDynamics 

pub struct SpacecraftDynamics {
    pub orbital_dyn: OrbitalDynamics,
    pub force_models: Vec<Arc<dyn ForceModel>>,
    pub guid_law: Option<Arc<dyn GuidanceLaw>>,
    pub decrement_mass: bool,
}

A generic spacecraft dynamics with associated force models, guidance law, and flag specifying whether to decrement the prop mass or not. Note: when developing new guidance laws, it is recommended to not enable prop decrement until the guidance law seems to work without proper physics. Note: if the spacecraft runs out of prop, the propagation segment will return an error.

Fields

orbital_dyn: OrbitalDynamics

force_models: Vec<Arc<dyn ForceModel>>

guid_law: Option<Arc<dyn GuidanceLaw>>

decrement_mass: bool

Implementations

impl SpacecraftDynamics

pub fn from_guidance_law( orbital_dyn: OrbitalDynamics, guid_law: Arc<dyn GuidanceLaw>, ) -> Self

Initialize a Spacecraft with a set of orbital dynamics and a propulsion subsystem. By default, the mass of the vehicle will be decremented as propellant is consumed.

pub fn from_guidance_law_no_decr( orbital_dyn: OrbitalDynamics, guid_law: Arc<dyn GuidanceLaw>, ) -> Self

Initialize a Spacecraft with a set of orbital dynamics and a propulsion subsystem. Will not decrement the prop mass as propellant is consumed.

pub fn new(orbital_dyn: OrbitalDynamics) -> Self

Initialize a Spacecraft with a set of orbital dynamics and with SRP enabled.

Examples found in repository

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

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

pub fn from_model( orbital_dyn: OrbitalDynamics, force_model: Arc<dyn ForceModel>, ) -> Self

Initialize new spacecraft dynamics with the provided orbital mechanics and with the provided force model.

Examples found in repository

nyx-core/examples/03_geo_analysis/raise_optim.rs (line 178)

fn evaluate_weights(
    weights: &[f32],
    prop_time_days: f64,
    state: Arc<SharedState>,
) -> Result<(f64, f64), Box<dyn Error>> {
    let ηthresholds: Vec<f64> = weights.iter().map(|w| *w as f64).collect();

    let eme2k = state.almanac.frame_info(EARTH_J2000).unwrap();
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    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))
        .srp(SRPData::from_area(3.0 * 6.0))
        .thruster(Thruster {
            isp_s: 4435.0,
            thrust_N: 0.472,
        })
        .mode(GuidanceMode::Thrust)
        .build();

    let prop_time = prop_time_days * Unit::Day;

    let objectives = &[
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            30_000.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 kluever_ctrl = Kluever::from_max_eclipse(objectives, &weights_f64, 0.2);
    let ctrl = Ruggiero::from_ηthresholds(objectives, &ηthresholds, sc)?;

    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
    orbital_dyn.accel_models.push(state.harmonics.clone());

    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, state.srp_dyn.clone())
        .with_guidance_law(ctrl.clone());

    let (final_state, _traj) = Propagator::rk89(
        sc_dynamics.clone(),
        IntegratorOptions::builder()
            .min_step(10.0_f64.seconds())
            .tolerance(1e-8)
            .error_ctrl(ErrorControl::RSSCartesianStep)
            .build(),
    )
    .with(sc, state.almanac.clone())
    .for_duration_with_traj(prop_time)?;

    let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;

    let mut penalty = 0.0;
    for obj in objectives {
        let (achieved, error) = obj.assess(&final_state)?;
        if !achieved {
            penalty += error.abs();
        }
        info!("{obj} error: {error:.3}, achieved? {achieved}");
    }

    info!("{ηthresholds:?} -> {prop_usage:.3} kg\tpenalty = {penalty:.3}");

    Ok((prop_usage, penalty * 1000.0))
}

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

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::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true)?,
    );
    orbital_dyn.accel_models.push(harmonics);

    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, almanac.clone())?;
    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 123)

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::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true).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.clone())?;

    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
    let 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/02_jwst_covar_monte_carlo/main.rs (line 103)

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.clone())?;

    // 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 120)

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

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

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

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

    println!("{dynamics}");

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

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

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::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 21, 21, true).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.clone())?;

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

    println!("{dynamics}");

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

Additional examples can be found in:

• nyx-core/examples/01_orbit_prop/main.rs
• nyx-core/examples/04_lro_od/main.rs

pub fn from_models( orbital_dyn: OrbitalDynamics, force_models: Vec<Arc<dyn ForceModel>>, ) -> Self

Initialize new spacecraft dynamics with a vector of force models.

pub fn guidance_achieved( &self, state: &Spacecraft, ) -> Result<bool, GuidanceError>

A shortcut to spacecraft.guid_law if a guidance law is defined for these dynamics

pub fn with_guidance_law(&self, guid_law: Arc<dyn GuidanceLaw>) -> Self

Clone these spacecraft dynamics and update the control to the one provided.

Examples found in repository

nyx-core/examples/03_geo_analysis/raise_optim.rs (line 179)

fn evaluate_weights(
    weights: &[f32],
    prop_time_days: f64,
    state: Arc<SharedState>,
) -> Result<(f64, f64), Box<dyn Error>> {
    let ηthresholds: Vec<f64> = weights.iter().map(|w| *w as f64).collect();

    let eme2k = state.almanac.frame_info(EARTH_J2000).unwrap();
    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);

    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))
        .srp(SRPData::from_area(3.0 * 6.0))
        .thruster(Thruster {
            isp_s: 4435.0,
            thrust_N: 0.472,
        })
        .mode(GuidanceMode::Thrust)
        .build();

    let prop_time = prop_time_days * Unit::Day;

    let objectives = &[
        Objective::within_tolerance(
            StateParameter::Element(OrbitalElement::SemiMajorAxis),
            30_000.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 kluever_ctrl = Kluever::from_max_eclipse(objectives, &weights_f64, 0.2);
    let ctrl = Ruggiero::from_ηthresholds(objectives, &ηthresholds, sc)?;

    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
    orbital_dyn.accel_models.push(state.harmonics.clone());

    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, state.srp_dyn.clone())
        .with_guidance_law(ctrl.clone());

    let (final_state, _traj) = Propagator::rk89(
        sc_dynamics.clone(),
        IntegratorOptions::builder()
            .min_step(10.0_f64.seconds())
            .tolerance(1e-8)
            .error_ctrl(ErrorControl::RSSCartesianStep)
            .build(),
    )
    .with(sc, state.almanac.clone())
    .for_duration_with_traj(prop_time)?;

    let prop_usage = sc.mass.prop_mass_kg - final_state.mass.prop_mass_kg;

    let mut penalty = 0.0;
    for obj in objectives {
        let (achieved, error) = obj.assess(&final_state)?;
        if !achieved {
            penalty += error.abs();
        }
        info!("{obj} error: {error:.3}, achieved? {achieved}");
    }

    info!("{ηthresholds:?} -> {prop_usage:.3} kg\tpenalty = {penalty:.3}");

    Ok((prop_usage, penalty * 1000.0))
}

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

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::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true)?,
    );
    orbital_dyn.accel_models.push(harmonics);

    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, almanac.clone())?;
    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 124)

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::from_stor(
        almanac.frame_info(IAU_EARTH_FRAME)?,
        GravityFieldData::from_cof(&jgm3_meta.uri, 8, 8, true).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.clone())?;

    // Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
    // acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
    let 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(())
}

Trait Implementations

impl Clone for SpacecraftDynamics

fn clone(&self) -> SpacecraftDynamics

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl Debug for SpacecraftDynamics

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

Formats the value using the given formatter.

impl Display for SpacecraftDynamics

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

Formats the value using the given formatter.

impl Dynamics for SpacecraftDynamics

type HyperdualSize = Const<9>

The state of the associated hyperdual state, almost always StateType + U1

type StateType = Spacecraft

fn finally( &self, next_state: Self::StateType, almanac: Arc<Almanac>, ) -> Result<Self::StateType, DynamicsError>

Optionally performs some final changes after each successful integration of the equations of motion. For example, this can be used to update the Guidance mode. NOTE: This function is also called just prior to very first integration step in order to update the initial state if needed.

fn eom( &self, delta_t_s: f64, state: &OVector<f64, Const<90>>, ctx: &Self::StateType, almanac: Arc<Almanac>, ) -> Result<OVector<f64, Const<90>>, DynamicsError>

Defines the equations of motion for these dynamics, or a combination of provided dynamics. The time delta_t is in seconds PAST the context epoch. The state vector is the state which changes for every intermediate step of the integration. The state context is the state of what is being propagated, it should allow rebuilding a new state context from the provided state vector.

fn dual_eom( &self, delta_t_s: f64, ctx: &Self::StateType, almanac: Arc<Almanac>, ) -> Result<(OVector<f64, Const<9>>, OMatrix<f64, Const<9>, Const<9>>), DynamicsError>

Defines the equations of motion for Dual numbers for these dynamics. All dynamics need to allow for automatic differentiation. However, if differentiation is not supported, then the dynamics should prevent initialization with a context which has an STM defined.