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    Spacecraft, State,
16    cosmic::{GuidanceMode, Mass, MetaAlmanac, Orbit, SRPData},
17    dynamics::{
18        GravityField, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
19        guidance::{Ruggiero, Thruster},
20    },
21    io::{ExportCfg, gravity::GravityFieldData},
22    mc::{MonteCarlo, MvnSpacecraft, StateDispersion},
23    md::prelude::{Objective, OrbitalElement, StateParameter},
24    propagators::{ErrorControl, IntegratorOptions, Propagator},
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 = GravityField::new(GravityFieldData::from_cof(
86        &jgm3_meta.uri,
87        8,
88        8,
89        almanac.frame_info(IAU_EARTH_FRAME)?,
90    )?);
91    orbital_dyn.accel_models.push(harmonics);
92
93    let srp_dyn = SolarPressure::default_flux(EARTH_J2000, &almanac)?;
94    let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn)
95        .with_guidance_law(ruggiero_ctrl.clone());
96
97    println!("{sc_dynamics}");
98
99    // Finally, let's use the Monte Carlo framework built into Nyx to propagate spacecraft.
100
101    // Let's start by defining the dispersion.
102    // 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.
103    // Note that additional validation on the MVN is in progress -- https://github.com/nyx-space/nyx/issues/339.
104    let mc_rv = MvnSpacecraft::new(
105        sc,
106        vec![StateDispersion::zero_mean(
107            StateParameter::Element(OrbitalElement::SemiMajorAxis),
108            3.0,
109        )],
110    )?;
111
112    let my_mc = MonteCarlo::new(
113        sc, // Nominal state
114        mc_rv,
115        "03_geo_sk".to_string(), // Scenario name
116        None, // No specific seed specified, so one will be drawn from the computer's entropy.
117    );
118
119    // Build the propagator setup.
120    let setup = Propagator::rk89(
121        sc_dynamics.clone(),
122        IntegratorOptions::builder()
123            .min_step(10.0_f64.seconds())
124            .error_ctrl(ErrorControl::RSSCartesianStep)
125            .build(),
126    );
127
128    let num_runs = 25;
129    let rslts = my_mc.run_until_epoch(setup, almanac.clone(), sc.epoch() + prop_time, num_runs);
130
131    assert_eq!(rslts.runs.len(), num_runs);
132
133    rslts.to_parquet("03_geo_sk.parquet", ExportCfg::default())?;
134
135    Ok(())
136}
03_geo_sk/stationkeeping.rs

03_geo_sk/
stationkeeping.rs
