03_geo_sk/
stationkeeping.rs

1#![doc = include_str!("./README.md")]
2extern crate log;
3extern crate nyx_space as nyx;
4extern crate pretty_env_logger as pel;
5
6use anise::{
7    almanac::metaload::MetaFile,
8    constants::{
9        celestial_objects::{MOON, SUN},
10        frames::{EARTH_J2000, IAU_EARTH_FRAME},
11    },
12};
13use hifitime::{Epoch, TimeUnits, Unit};
14use nyx::{
15    cosmic::{eclipse::EclipseLocator, GuidanceMode, Mass, MetaAlmanac, Orbit, SRPData},
16    dynamics::{
17        guidance::{Ruggiero, Thruster},
18        Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
19    },
20    io::{gravity::HarmonicsMem, ExportCfg},
21    mc::{MonteCarlo, MvnSpacecraft, StateDispersion},
22    md::{prelude::Objective, StateParameter},
23    propagators::{ErrorControl, IntegratorOptions, Propagator},
24    Spacecraft, State,
25};
26use std::{error::Error, sync::Arc};
27
28fn main() -> Result<(), Box<dyn Error>> {
29    pel::init();
30    // Set up the dynamics like in the orbit raise.
31    let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?);
32    let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14);
33
34    // Define the GEO orbit, and we're just going to maintain it very tightly.
35    let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?;
36    let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?;
37    println!("{orbit:x}");
38
39    let sc = Spacecraft::builder()
40        .orbit(orbit)
41        .mass(Mass::from_dry_and_prop_masses(1000.0, 1000.0)) // 1000 kg of dry mass and prop, totalling 2.0 tons
42        .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
43        .thruster(Thruster {
44            // "NEXT-STEP" row in Table 2
45            isp_s: 4435.0,
46            thrust_N: 0.472,
47        })
48        .mode(GuidanceMode::Thrust) // Start thrusting immediately.
49        .build();
50
51    // Set up the spacecraft dynamics like in the orbit raise example.
52
53    let prop_time = 30.0 * Unit::Day;
54
55    // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089.
56    let objectives = &[
57        Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km
58        Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5),
59        Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2),
60    ];
61
62    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
63    println!("{ruggiero_ctrl}");
64
65    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
66
67    let mut jgm3_meta = MetaFile {
68        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
69        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
70    };
71    jgm3_meta.process(true)?;
72
73    let harmonics = Harmonics::from_stor(
74        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
75        HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
76    );
77    orbital_dyn.accel_models.push(harmonics);
78
79    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
80    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
81        .with_guidance_law(ruggiero_ctrl.clone());
82
83    println!("{sc_dynamics}");
84
85    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
86
87    // Let's start by defining the dispersion.
88    // 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.
89    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
90    let mc_rv = MvnSpacecraft::new(
91        sc,
92        vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
93    )?;
94
95    let my_mc = MonteCarlo::new(
96        sc, // Nominal state
97        mc_rv,
98        "03_geo_sk".to_string(), // Scenario name
99        None, // No specific seed specified, so one will be drawn from the computer's entropy.
100    );
101
102    // Build the propagator setup.
103    let setup = Propagator::rk89(
104        sc_dynamics.clone(),
105        IntegratorOptions::builder()
106            .min_step(10.0_f64.seconds())
107            .error_ctrl(ErrorControl::RSSCartesianStep)
108            .build(),
109    );
110
111    let num_runs = 25;
112    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
113
114    assert_eq!(rslts.runs.len(), num_runs);
115
116    // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
117
118    rslts.to_parquet(
119        "03_geo_sk.parquet",
120        Some(vec![
121            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
122        ]),
123        ExportCfg::default(),
124        almanac,
125    )?;
126
127    Ok(())
128}
03_geo_sk/stationkeeping.rs

03_geo_sk/
stationkeeping.rs
