Add GEO station keeping example

examples/03_geo_analysis/README.md

@@ -1,6 +1,6 @@
 # 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.
+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
 
@@ -101,4 +101,26 @@ Looking edge-on shows how the inclination is changed over time.
 
 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)
+![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/plot_sk_mc.py

@@ -0,0 +1,25 @@
+import polars as pl
+import plotly.graph_objs as go
+
+
+if __name__ == "__main__":
+    df = pl.read_parquet("03_geo_sk.parquet")
+
+    df = df.with_columns(pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f")).sort(
+        "Epoch (UTC)", descending=False
+    )
+
+    columns = [
+        "sma (km)",
+        "ecc",
+        "inc (deg)",
+        "raan (deg)",
+        "aop (deg)",
+        "ta (deg)",
+        "aol (deg)",
+        "tlong (deg)",
+        "fuel_mass (kg)",
+    ]
+
+    for col in columns:
+        go.Figure(data=[go.Scattergl(x=df["Epoch (UTC)"], y=df[col], name=col, opacity=0.05, showlegend=True,text=df["Monte Carlo Run Index"])]).show()

examples/03_geo_analysis/stationkeeping.rs

@@ -53,7 +53,7 @@ fn main() -> Result<(), Box<dyn Error>> {
 
     // Set up the spacecraft dynamics like in the orbit raise example.
 
-    let prop_time = 60.0 * Unit::Day;
+    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 = &[
@@ -88,14 +88,11 @@ fn main() -> Result<(), Box<dyn Error>> {
     // 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 = MultivariateNormal::new(
         sc,
-        vec![
-            StateDispersion::zero_mean(StateParameter::SMA, 3.0),
-            StateDispersion::zero_mean(StateParameter::Inclination, 0.05),
-        ],
+        vec![StateDispersion::zero_mean(StateParameter::SMA, 3.0)],
     )?;
 
     let my_mc = MonteCarlo::new(
@@ -114,7 +111,7 @@ fn main() -> Result<(), Box<dyn Error>> {
             .build(),
     );
 
-    let num_runs = 10;
+    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);