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::{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, OrbitalElement, 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_info(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(
58            StateParameter::Element(OrbitalElement::SemiMajorAxis),
59            42_165.0,
60            20.0,
61        ),
62        Objective::within_tolerance(
63            StateParameter::Element(OrbitalElement::Eccentricity),
64            0.001,
65            5e-5,
66        ),
67        Objective::within_tolerance(
68            StateParameter::Element(OrbitalElement::Inclination),
69            0.05,
70            1e-2,
71        ),
72    ];
73
74    let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, 0.2)?;
75    println!("{ruggiero_ctrl}");
76
77    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]);
78
79    let mut jgm3_meta = MetaFile {
80        uri: "http://public-data.nyxspace.com/nyx/models/JGM3.cof.gz".to_string(),
81        crc32: Some(0xF446F027), // Specifying the CRC32 avoids redownloading it if it's cached.
82    };
83    jgm3_meta.process(true)?;
84
85    let harmonics = Harmonics::from_stor(
86        almanac.frame_info(IAU_EARTH_FRAME)?,
87        HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?,
88    );
89    orbital_dyn.accel_models.push(harmonics);
90
91    let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?;
92    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
93        .with_guidance_law(ruggiero_ctrl.clone());
94
95    println!("{sc_dynamics}");
96
97    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
98
99    // Let's start by defining the dispersion.
100    // 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.
101    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
102    let mc_rv = MvnSpacecraft::new(
103        sc,
104        vec![StateDispersion::zero_mean(
105            StateParameter::Element(OrbitalElement::SemiMajorAxis),
106            3.0,
107        )],
108    )?;
109
110    let my_mc = MonteCarlo::new(
111        sc, // Nominal state
112        mc_rv,
113        "03_geo_sk".to_string(), // Scenario name
114        None, // No specific seed specified, so one will be drawn from the computer's entropy.
115    );
116
117    // Build the propagator setup.
118    let setup = Propagator::rk89(
119        sc_dynamics.clone(),
120        IntegratorOptions::builder()
121            .min_step(10.0_f64.seconds())
122            .error_ctrl(ErrorControl::RSSCartesianStep)
123            .build(),
124    );
125
126    let num_runs = 25;
127    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
128
129    assert_eq!(rslts.runs.len(), num_runs);
130
131    rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;
132
133    Ok(())
134}
03_geo_sk/stationkeeping.rs

03_geo_sk/
stationkeeping.rs
