Add GEO drift, orbit raise, and station keeping example

.github/workflows/rust.yaml

@@ -85,6 +85,9 @@ jobs:
         run: |
           cargo run --example 01_orbit_prop --release
           cargo run --example 02_jwst --release
+          cargo run --example 03_geo_drift --release
+          cargo run --example 03_geo_raise --release
+          cargo run --example 03_geo_sk --release
 
   lints:
     name: Lints

Cargo.toml

@@ -105,3 +105,15 @@ path = "examples/01_orbit_prop/main.rs"
 [[example]]
 name = "02_jwst"
 path = "examples/02_jwst_covar_monte_carlo/main.rs"
+
+[[example]]
+name = "03_geo_drift"
+path = "examples/03_geo_analysis/drift.rs"
+
+[[example]]
+name = "03_geo_raise"
+path = "examples/03_geo_analysis/raise.rs"
+
+[[example]]
+name = "03_geo_sk"
+path = "examples/03_geo_analysis/stationkeeping.rs"

examples/03_geo_analysis/README.md

@@ -0,0 +1,126 @@
+# Geostationary spacecraft analysis
+
+GEO birds have specific station keeping requirements. In this preliminary study, we'll look at the drift rate inside the GEO box, the orbit raising from a GTO to a GEO orbit, and run a Monte Carlo of the station keeping requirements. All of this analysis is very preliminary, and any real-world mission would need much more thorough analysis.
+
+## Drift
+
+### Objective
+
+As documented in the [drift analysis code](./drift.rs), we're trying to gain insight into how quickly a vehicle would drift outside of an arbitrary box. The spacecraft is placed arbitrarily to be in a GEO orbit centered on Colorado, USA. We then propagate the original spacecraft position in high fidelity with the following dynamics:
+
+- Earth gravity field of 21x21, JGM3 model -- we don't need a super high field here because the spacecraft is ~42 Mm above the Earth
+- Solar radiation pressure, accounting for the shadow and penumbra cast by the Earth and the Moon at each time step
+- Point mass gravity of the Moon and Sun (in addition to the Earth, the central body)
+
+To run the [drift analysis](./drift.rs) example, just execute:
+```sh
+RUST_LOG=info cargo run --example 03_geo_drift --release
+```
+
+The vehicle is not under any active control in this example. We're really only trying to gain insight into the periodicity of the orbit.
+
+**Performance wise**, Nyx can propagate these high fidelity dynamics at about 560 days per minute on my 8 year old computer, or about 1.5 years per minute, even for much longer propagation times.
+
+```sh
+ INFO  nyx_space::propagators::instance   > Propagating for 1095 days 18 h until 2027-03-01T06:13:14 UTC
+ INFO  nyx_space::propagators::instance   > 	... current epoch 2025-09-22T07:51:27.645581934 UTC, remaing 524 days 22 h 21 min 46 s
+ INFO  nyx_space::propagators::instance   > 	... done in 1 min 55 s 661 ms 536 μs 822 ns
+(...)
+Longitude changed by -2.352 deg -- Box is 0.1 deg E-W
+Latitude changed by -2.637 deg -- Box is 0.05 deg N-S
+Altitude changed by -10.407 km -- Box is 30 km
+```
+
+### Analysis and verification
+
+#### GEO box
+
+The ITU allocates specific geostationary boxes to a number of states, who then grant their use to specific operators on the condition that they remain in that box.
+
+In this analysis, we're assessing how quickly an arbitrary GEO bird would drift out of the box. After running the `03_geo_drift` example, run the `plot_drift.py` and `plot_orbital_elements.py` (with the proper parquet file as an input) scripts.
+
+The latitude and longitude drift quite rapidly out of this arbitrary box of +/- 0.1 degrees East/West and +/- 0.05 degrees North/South.
+
+_Note:_ the drift analysis below was done with a much lighter spacecraft than in the two other examples.
+
+![Lat Long drift](./plots/drift-lat-long.png)
+
+![Orbital elements drift](./plots/drift-orbital-elements-3-years.png)
+
+Zooming into the drift of three weeks allows us to see a bit more of the short term trends. It also becomes more obvious why circular equatorial orbits usually use the true longitude instead of the true anomaly as a phasing angle: the former varies nearly-linearly whereas the latter is ill-defined.
+
+![Orbital elements short term drift](./plots/drift-orbital-elements-3-weeks.png)
+
+#### Eclipsing
+
+As expected for GEO birds, there are only a few eclipsing periods throughout a three year span. These eclipses are quite short, but we want to capture them. To do so, we export the penumbra event (look for `trajectory.to_parquet` in the code) at a 5 minute interval, ensuring that we interpolate the propagated trajectory every five minutes and evaluate the event at that same sample rate.
+
+![3D Traj and illumination](./plots/drift-3d-illumination.png)
+
+We can also use polars and plotly in Python directly to plot when these eclipses happen.
+
+```py
+>>> import polars as pl
+>>> import plotly.express as px
+>>> df = pl.read_parquet("03_geo_hf_prop.parquet")
+>>> px.line(df, x="Epoch (UTC)", y="penumbra event light-source: Sun J2000, shadows casted by: Earth J2000, Moon J2000")
+```
+
+We have daily eclipses during the equinoxes, as [verified by other sources](https://www.sws.bom.gov.au/Educational/5/4/3).
+
+> Between 28 February and 11 April, and between 2 September and 14 October, roughly 21 days either side equinoxes, satellites in geostationary orbits will pass through the shadow of the earth once every day.
+
+![Eclipses over 3 years](./plots/drift-eclipse.png)
+
+![Eclipse depth](./plots/drift-eclipse-depth.png)
+
+Refer to the [Eclipse section of the MathSpec](https://nyxspace.com/nyxspace/MathSpec/celestial/eclipse/) for computation details.
+
+## Orbit raise
+
+Most GEO spacecraft are deployed on a geostationary transfer orbit (GTO) and must rely on their own propulsion to reach their allocated GEO slot. Over the past fifteen years, we've seen a shift in the industry from high thrust birds to low thrust propulsion as they're many times more efficient. However, the optimization of the trajectory for low thrust propulsion is significantly more complicated. **Nyx allows for the development of custom guidance laws using the [`GuidanceLaw`](https://rustdoc.nyxspace.com/nyx_space/dynamics/guidance/trait.GuidanceLaw.html) trait.**
+
+As of version 2.0.0-rc, Nyx ships with the closed loop Ruggerio guidance law. This is a locally optimal law for low thrust propulsion. It allows targeting of all orbital elements _apart_ from the phasing (e.g. the true anomaly). As such, when applied to orbit raising, it is useful for reaching the vicinity of a desired orbital slot and maintaining control of the shape and inclination of an orbit, but more complicated methods are recommended to reach to exact ITU-allocated slot. **A more optimal guidance law for this problem is the Q-Law**, whose implementation in Nyx has been on the backburner for years. The Q-law guidance law is as optimal as a global trajectory optimizer where the phasing is left free, but runs extremely fast.
+
+To run the [orbit raise analysis](./raise.rs) example, just execute:
+```sh
+RUST_LOG=info cargo run --example 03_geo_raise --release
+```
+
+To build the following plots, use the `plot_3d_traj.py` script and the `plot_orbital_elements.py` scripts.
+
+![Orbital elements during orbit raise](./plots/raise-keplerian-oe.png)
+
+In the two follow plots, the colors correspond to the remaining fuel mass, thereby showing the fuel depletion over the orbit raise.
+
+![3D traj raise](./plots/raise-traj-3d.png)
+
+Looking edge-on shows how the inclination is changed over time.
+
+![3D traj raise edge-on](./plots/raise-traj-3d-edge-on.png)
+
+During the orbit raise, a low thrust vehicle will most likely not thrust when it's in the shadow. In this analysis, we specify that the vehicle should not thrust when it is in over 80% of penumbra, i.e. it can only generate _at most_ 20% of the power it would generate in full sunshine (depending on pointing of course).
+
+![3D traj illumination](./plots/raise-3d-illumination.png)
+
+## Station keeping
+
+As previously mentioned, station keeping is required for a GEO slot. In the third program, we look at the use of the Monte Carlo framework in Nyx which uses the multivariate normal distribution structure. We distribute the SMA of 25 spacecraft, propagate them for a two week period in high fidelity, and ensure that they Ruggiero guidance law is enabled and tight around the GEO box. **The Monte Carlo capabilities of Nyx are better demonstrated in the [02 JWST](../02_jwst_covar_monte_carlo/README.md) example.**
+
+To run the [station keeping Monte Carlo](./stationkeeping.rs) example, just execute:
+```sh
+RUST_LOG=info cargo run --example 03_geo_sk --release
+```
+
+Over a two week period, this two-ton spacecraft would need roughly 0.8 kg of fuel (if using the _NEXT-STEP_ engine, cf. the comments in the drift analysis code) +/- 0.1 kg for station keeping.
+
+![Fuel mass](./plots/sk-fuel-mass.png)
+
+The inclination plot shows when the guidance law turns on, and shows that we maintain the constraint within the specified objectives.
+
+![Inclination](./plots/sk-inc.png)
+![Eccentricity](./plots/sk-ecc.png)
+
+### Further analysis
+
+Additional analysis would run this Monte Carlo for longer and with many more spacecraft (upward of 100), and crucially ensure that the Ruggiero guidance law bounds correspond to the GEO box. Subsequently, one should implement the Q-Law guidance law for more fuel economy. Finally, the analysis should also include a variation on the tightness of the box, especially if the vehicle is equipped with a variable thrust engine as one may wish to drift less out of the SK box and keep the engine at a lower thrust level, or vice versa.

examples/03_geo_analysis/drift.rs

@@ -0,0 +1,217 @@
+extern crate log;
+extern crate nyx_space as nyx;
+extern crate pretty_env_logger as pel;
+
+use anise::{
+    almanac::metaload::MetaFile,
+    constants::{
+        celestial_objects::{MOON, SUN},
+        frames::{EARTH_J2000, IAU_EARTH_FRAME, MOON_J2000},
+    },
+};
+use hifitime::{Epoch, Unit};
+use nyx::{
+    cosmic::{eclipse::EclipseLocator, MetaAlmanac, Orbit, SrpConfig},
+    dynamics::{Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics},
+    io::{gravity::HarmonicsMem, ExportCfg},
+    propagators::Propagator,
+    Spacecraft, State,
+};
+use polars::{df, prelude::ParquetWriter};
+
+use std::fs::File;
+use std::{error::Error, sync::Arc};
+
+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_from_uid(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)
+        .dry_mass_kg(9.60)
+        .srp(SrpConfig {
+            area_m2: 10e-4,
+            cr: 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()?;
+
+    // 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 = Harmonics::from_stor(
+        almanac.frame_from_uid(IAU_EARTH_FRAME)?,
+        HarmonicsMem::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",
+        Some(vec![
+            &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event()
+        ]),
+        ExportCfg::builder().step(analysis_step).build(),
+        almanac.clone(),
+    )?;
+
+    // 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(())
+}