diff --git a/.github/workflows/rust.yaml b/.github/workflows/rust.yaml index b5f49de6..ae945ff4 100644 --- a/.github/workflows/rust.yaml +++ b/.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 diff --git a/Cargo.toml b/Cargo.toml index 572119e4..199f0cd7 100644 --- a/Cargo.toml +++ b/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" diff --git a/examples/03_geo_analysis/README.md b/examples/03_geo_analysis/README.md new file mode 100644 index 00000000..6bef6fca --- /dev/null +++ b/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. \ No newline at end of file diff --git a/examples/03_geo_analysis/drift.rs b/examples/03_geo_analysis/drift.rs new file mode 100644 index 00000000..41d6503a --- /dev/null +++ b/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> { + 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(()) +} diff --git a/examples/03_geo_analysis/plot_3d_traj.py b/examples/03_geo_analysis/plot_3d_traj.py new file mode 100644 index 00000000..7dab2ffd --- /dev/null +++ b/examples/03_geo_analysis/plot_3d_traj.py @@ -0,0 +1,98 @@ +import argparse +from datetime import datetime +import polars as pl +import plotly.graph_objs as go +from plotly.subplots import make_subplots +import numpy as np + + +def build_sphere(size, num_points=500, opacity=1.0): + """ + Source: https://python.plainenglish.io/how-to-create-a-3d-model-of-the-solar-system-with-plotly-in-python-2a74e672b771 + """ + + color = (86, 180, 233) # #56b4e9 + + theta = np.linspace(0, 2 * np.pi, num_points) + phi = np.linspace(0, np.pi, num_points) + + # Set up coordinates for points on the sphere + x0 = size * np.outer(np.cos(theta), np.sin(phi)) + y0 = size * np.outer(np.sin(theta), np.sin(phi)) + z0 = size * np.outer(np.ones(num_points), np.cos(phi)) + + # Set up trace + trace = go.Surface( + x=x0, + y=y0, + z=z0, + colorscale=[ + [0, f"rgba({int(color[0])}, {int(color[1])}, {int(color[2])}, {opacity})"], + [1, f"rgba({int(color[0])}, {int(color[1])}, {int(color[2])}, {opacity})"], + ], + hoverinfo="none", + ) + trace.update(showscale=False) + + return trace + + +def plot_traj( + df: pl.DataFrame, + colored_by="fuel_mass (kg)", + color_descr="fuel mass (kg)", + scale=1.0 +): + """ + Plot a trajectory in 3D + + Args: + dfs (pandas.DataFrame): The data frame containing the trajectory (or a list thereof) + title (str): The title of the plot + html_out (str, optional): The path to save the HTML to. Defaults to None. + copyright (str, optional): The copyright to display on the plot. Defaults to None. + fig (plotly.graph_objects.Figure, optional): The figure to add the trajectory to. Defaults to None. + center (str, optional): The name of the center object, e.g. `Luna` (Moon). Defaults to "Earth". + show (bool, optional): Whether to show the plot. Defaults to True. If set to false, the figure will be returned. + """ + + df = df.with_columns(pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f")).sort( + "Epoch (UTC)", descending=False + ) + + fig = go.Figure() + + earth_radius = 6378.136300 + + traces = [build_sphere(earth_radius, opacity=0.8)] + + traces += [ + go.Scatter3d( + x=df["x (km)"], + y=df["y (km)"], + z=df["z (km)"], + mode="lines", + line=dict( + color=df[colored_by]*scale, + colorscale="Viridis", + colorbar=dict(title=color_descr), + ), + name="GEO", + ) + ] + + # Now that we have the data, let's plot it + fig.add_traces(traces) + fig.update_layout(scene_aspectmode="data") + + fig.show() + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="3D traj plotter") + parser.add_argument("pq", type=str, help="Path to the parquet file") + args = parser.parse_args() + + df = pl.read_parquet(args.pq) + plot_traj(df) + plot_traj(df, colored_by="penumbra event light-source: Sun J2000, shadows casted by: Earth J2000, Moon J2000", color_descr="Illumination %", scale=100.0) diff --git a/examples/03_geo_analysis/plot_drift.py b/examples/03_geo_analysis/plot_drift.py new file mode 100644 index 00000000..d040b2e9 --- /dev/null +++ b/examples/03_geo_analysis/plot_drift.py @@ -0,0 +1,101 @@ +import polars as pl +import plotly.graph_objs as go +from plotly.subplots import make_subplots + +if __name__ == "__main__": + df_lla = pl.read_parquet("./03_geo_lla.parquet") + df_lla = df_lla.with_columns( + pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f") + ).sort("Epoch (UTC)", descending=False) + print(df_lla.describe()) + + # Create subplots + fig = make_subplots( + rows=2, + cols=1, + shared_xaxes=True, + vertical_spacing=0.1, + subplot_titles=("Latitude vs Epoch", "Longitude vs Epoch"), + ) + + # Scattergl plot for Latitude vs Epoch + fig.add_trace( + go.Scattergl( + x=df_lla["Epoch (UTC)"], + y=df_lla["Latitude N-S (deg)"], + mode="markers+lines", + name="Latitude", + marker=dict(color="blue"), + ), + row=1, + col=1, + ) + + # Scattergl plot for Longitude vs Epoch + fig.add_trace( + go.Scattergl( + x=df_lla["Epoch (UTC)"], + y=df_lla["Longitude E-W (deg)"], + mode="markers+lines", + name="Longitude", + marker=dict(color="blue"), + ), + row=2, + col=1, + ) + + # Add vertical line to the Latitude plot + first_latitude = df_lla["Latitude N-S (deg)"][0] + lower_lat_deg = first_latitude - 0.05 + higher_lat_deg = first_latitude + 0.05 + + fig.add_hline( + lower_lat_deg, + row=1, + col=1, + annotation_text=f"{lower_lat_deg:.3} - 0.05 deg", + annotation_position="top right", + line=dict(color="Red", width=1, dash="dash"), + ) + fig.add_hline( + higher_lat_deg, + row=1, + col=1, + annotation_text=f"{lower_lat_deg:.3} + 0.05 deg", + annotation_position="top right", + line=dict(color="Red", width=1, dash="dash"), + ) + + # Add vertical line to the longitude plot + first_longitude = df_lla["Longitude E-W (deg)"][0] + lower_long_deg = first_longitude - 0.1 + higher_long_deg = first_longitude + 0.1 + + fig.add_hline( + lower_long_deg, + row=2, + col=1, + annotation_text=f"{lower_long_deg:.3} - 0.1 deg", + annotation_position="top right", + line=dict(color="Red", width=1, dash="dash"), + ) + fig.add_hline( + higher_long_deg, + row=2, + col=1, + annotation_text=f"{lower_lat_deg:.3} + 0.1 deg", + annotation_position="top right", + line=dict(color="Red", width=1, dash="dash"), + ) + + # Update layout + fig.update_layout( + title_text="Latitude and Longitude vs Epoch", + xaxis_title="Epoch", + yaxis_title="Latitude N-S (deg)", + xaxis2_title="Epoch", + yaxis2_title="Longitude E-W (deg)", + ) + + # Show the plot + fig.show() diff --git a/examples/03_geo_analysis/plot_orbital_elements.py b/examples/03_geo_analysis/plot_orbital_elements.py new file mode 100644 index 00000000..8475bc5d --- /dev/null +++ b/examples/03_geo_analysis/plot_orbital_elements.py @@ -0,0 +1,62 @@ +import argparse +import polars as pl +import plotly.graph_objs as go +from plotly.subplots import make_subplots + + +def plot_orbit_elements( + df: pl.DataFrame, +): + """ + Plots the orbital elements: SMA, ECC, INC, RAAN, AOP, TA, True Longitude, AOL + """ + + 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)", + ] + + fig = make_subplots( + rows=4, + cols=2, + subplot_titles=columns, + shared_xaxes=True, + vertical_spacing=0.1, + ) + + row_i = 0 + col_i = 0 + + for col in columns: + name = col + + fig.add_trace( + go.Scattergl(x=df["Epoch (UTC)"], y=df[col], name=name), + row=row_i + 1, + col=col_i + 1, + ) + col_i = (col_i + 1) % 2 + if col_i == 0: + row_i = (row_i + 1) % 4 + + fig.show() + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="3D traj plotter") + parser.add_argument("pq", type=str, help="Path to the parquet file") + args = parser.parse_args() + + df = pl.read_parquet(args.pq) + + plot_orbit_elements(df) diff --git a/examples/03_geo_analysis/plot_sk_mc.py b/examples/03_geo_analysis/plot_sk_mc.py new file mode 100644 index 00000000..29420129 --- /dev/null +++ b/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() diff --git a/examples/03_geo_analysis/plots/drift-3d-illumination.png b/examples/03_geo_analysis/plots/drift-3d-illumination.png new file mode 100644 index 00000000..786fa49f Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-3d-illumination.png differ diff --git a/examples/03_geo_analysis/plots/drift-eclipse-depth.png b/examples/03_geo_analysis/plots/drift-eclipse-depth.png new file mode 100644 index 00000000..20114c05 Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-eclipse-depth.png differ diff --git a/examples/03_geo_analysis/plots/drift-eclipse.png b/examples/03_geo_analysis/plots/drift-eclipse.png new file mode 100644 index 00000000..596553f7 Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-eclipse.png differ diff --git a/examples/03_geo_analysis/plots/drift-lat-long.png b/examples/03_geo_analysis/plots/drift-lat-long.png new file mode 100644 index 00000000..47050178 Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-lat-long.png differ diff --git a/examples/03_geo_analysis/plots/drift-orbital-elements-3-weeks.png b/examples/03_geo_analysis/plots/drift-orbital-elements-3-weeks.png new file mode 100644 index 00000000..df7ef793 Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-orbital-elements-3-weeks.png differ diff --git a/examples/03_geo_analysis/plots/drift-orbital-elements-3-years.png b/examples/03_geo_analysis/plots/drift-orbital-elements-3-years.png new file mode 100644 index 00000000..f95ed8bb Binary files /dev/null and b/examples/03_geo_analysis/plots/drift-orbital-elements-3-years.png differ diff --git a/examples/03_geo_analysis/plots/raise-3d-illumination-2.png b/examples/03_geo_analysis/plots/raise-3d-illumination-2.png new file mode 100644 index 00000000..efc6eefe Binary files /dev/null and b/examples/03_geo_analysis/plots/raise-3d-illumination-2.png differ diff --git a/examples/03_geo_analysis/plots/raise-3d-illumination.png b/examples/03_geo_analysis/plots/raise-3d-illumination.png new file mode 100644 index 00000000..e7dc2937 Binary files /dev/null and b/examples/03_geo_analysis/plots/raise-3d-illumination.png differ diff --git a/examples/03_geo_analysis/plots/raise-keplerian-oe.png b/examples/03_geo_analysis/plots/raise-keplerian-oe.png new file mode 100644 index 00000000..9cb6c837 Binary files /dev/null and b/examples/03_geo_analysis/plots/raise-keplerian-oe.png differ diff --git a/examples/03_geo_analysis/plots/raise-traj-3d-edge-on.png b/examples/03_geo_analysis/plots/raise-traj-3d-edge-on.png new file mode 100644 index 00000000..fd1cdfed Binary files /dev/null and b/examples/03_geo_analysis/plots/raise-traj-3d-edge-on.png differ diff --git a/examples/03_geo_analysis/plots/raise-traj-3d.png b/examples/03_geo_analysis/plots/raise-traj-3d.png new file mode 100644 index 00000000..a362e618 Binary files /dev/null and b/examples/03_geo_analysis/plots/raise-traj-3d.png differ diff --git a/examples/03_geo_analysis/plots/sk-ecc.png b/examples/03_geo_analysis/plots/sk-ecc.png new file mode 100644 index 00000000..1c4089e0 Binary files /dev/null and b/examples/03_geo_analysis/plots/sk-ecc.png differ diff --git a/examples/03_geo_analysis/plots/sk-fuel-mass.png b/examples/03_geo_analysis/plots/sk-fuel-mass.png new file mode 100644 index 00000000..8221f6e6 Binary files /dev/null and b/examples/03_geo_analysis/plots/sk-fuel-mass.png differ diff --git a/examples/03_geo_analysis/plots/sk-inc.png b/examples/03_geo_analysis/plots/sk-inc.png new file mode 100644 index 00000000..e067b2db Binary files /dev/null and b/examples/03_geo_analysis/plots/sk-inc.png differ diff --git a/examples/03_geo_analysis/raise.rs b/examples/03_geo_analysis/raise.rs new file mode 100644 index 00000000..063db4a7 --- /dev/null +++ b/examples/03_geo_analysis/raise.rs @@ -0,0 +1,154 @@ +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}, + }, +}; +use hifitime::{Epoch, TimeUnits, Unit}; +use nyx::{ + cosmic::{ + eclipse::{EclipseLocator, EclipseState}, + GuidanceMode, MetaAlmanac, Orbit, SrpConfig, + }, + dynamics::{ + guidance::{GuidanceLaw, Ruggiero, Thruster}, + Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics, + }, + io::{gravity::HarmonicsMem, ExportCfg}, + md::{prelude::Objective, StateParameter}, + propagators::{PropOpts, Propagator, RSSCartesianStep}, + Spacecraft, +}; +use std::{error::Error, sync::Arc}; + +fn main() -> Result<(), Box> { + 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)?); + // Fetch the EME2000 frame from the Almabac + let eme2k = almanac.frame_from_uid(EARTH_J2000).unwrap(); + // Define the orbit epoch + let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14); + + // Build the spacecraft itself. + // Using slide 6 of https://aerospace.org/sites/default/files/2018-11/Davis-Mayberry_HPSEP_11212018.pdf + // for the "next gen" SEP characteristics. + + // GTO start + let orbit = Orbit::keplerian(24505.9, 0.725, 7.05, 0.0, 0.0, 0.0, epoch, eme2k); + + let sc = Spacecraft::builder() + .orbit(orbit) + .dry_mass_kg(1000.0) // 1000 kg of dry mass + .fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons + .srp(SrpConfig::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 + .thruster(Thruster { + // "NEXT-STEP" row in Table 2 + isp_s: 4435.0, + thrust_N: 0.472, + }) + .mode(GuidanceMode::Thrust) // Start thrusting immediately. + .build(); + + let prop_time = 180.0 * Unit::Day; + + // Define the guidance law -- we're just using a Ruggiero controller as demonstrated in AAS-2004-5089. + let objectives = &[ + Objective::within_tolerance(StateParameter::SMA, 42_165.0, 20.0), + Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5), + Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2), + ]; + + // Ensure that we only thrust if we have more than 20% illumination. + let ruggiero_ctrl = + Ruggiero::from_max_eclipse(objectives, sc, EclipseState::Penumbra(0.2)).unwrap(); + println!("{ruggiero_ctrl}"); + + // Define the high fidelity dynamics + + // 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 = Harmonics::from_stor( + almanac.frame_from_uid(IAU_EARTH_FRAME)?, + HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true).unwrap(), + ); + + // Include the spherical harmonics into the orbital dynamics. + orbital_dyn.accel_models.push(harmonics); + + // We define the solar radiation pressure, using the default solar flux and accounting only + // for the eclipsing caused by the Earth. + let srp_dyn = SolarPressure::default(EARTH_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 sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn) + .with_guidance_law(ruggiero_ctrl.clone()); + + println!("{:x}", orbit); + + // We specify a minimum step in the propagator because the Ruggiero control would otherwise drive this step very low. + let (final_state, traj) = Propagator::rk89( + sc_dynamics.clone(), + PropOpts::builder() + .min_step(10.0_f64.seconds()) + .error_ctrl(RSSCartesianStep {}) + .build(), + ) + .with(sc, almanac.clone()) + .for_duration_with_traj(prop_time)?; + + let fuel_usage = sc.fuel_mass_kg - final_state.fuel_mass_kg; + println!("{:x}", final_state.orbit); + println!("fuel usage: {:.3} kg", fuel_usage); + + // Finally, export the results for analysis, including the penumbra percentage throughout the orbit raise. + traj.to_parquet( + "./03_geo_raise.parquet", + Some(vec![ + &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event() + ]), + ExportCfg::default(), + almanac, + )?; + + for status_line in ruggiero_ctrl.status(&final_state) { + println!("{status_line}"); + } + + ruggiero_ctrl + .achieved(&final_state) + .expect("objective not achieved"); + + Ok(()) +} diff --git a/examples/03_geo_analysis/requirements.txt b/examples/03_geo_analysis/requirements.txt new file mode 100644 index 00000000..fcd2d3ec --- /dev/null +++ b/examples/03_geo_analysis/requirements.txt @@ -0,0 +1,37 @@ +asttokens==2.4.1 +attrs==23.2.0 +decorator==5.1.1 +executing==2.0.1 +iniconfig==2.0.0 +ipython==8.25.0 +jedi==0.19.1 +matplotlib-inline==0.1.7 +maturin==1.6.0 +numpy==1.26.4 +packaging==24.1 +pandas==2.1.4 +parso==0.8.4 +pexpect==4.9.0 +pip==24.0 +plotly==5.16.1 +pluggy==1.5.0 +polars==0.20.31 +prompt-toolkit==3.0.47 +ptyprocess==0.7.0 +pure-eval==0.2.2 +pyarrow==13.0.0 +pygments==2.18.0 +pytest==7.2.2 +python-dateutil==2.9.0.post0 +python-slugify==8.0.4 +pytz==2024.1 +ruff==0.5.0 +scipy==1.11.4 +six==1.16.0 +stack-data==0.6.3 +tenacity==8.3.0 +text-unidecode==1.3 +traitlets==5.14.3 +typing-extensions==4.12.2 +tzdata==2024.1 +wcwidth==0.2.13 diff --git a/examples/03_geo_analysis/stationkeeping.rs b/examples/03_geo_analysis/stationkeeping.rs new file mode 100644 index 00000000..442b5667 --- /dev/null +++ b/examples/03_geo_analysis/stationkeeping.rs @@ -0,0 +1,131 @@ +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}, + }, +}; +use hifitime::{Epoch, TimeUnits, Unit}; +use nyx::{ + cosmic::{ + eclipse::{EclipseLocator, EclipseState}, + GuidanceMode, MetaAlmanac, Orbit, SrpConfig, + }, + dynamics::{ + guidance::{Ruggiero, Thruster}, + Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics, + }, + io::{gravity::HarmonicsMem, ExportCfg}, + mc::{MonteCarlo, MultivariateNormal, StateDispersion}, + md::{prelude::Objective, StateParameter}, + propagators::{PropOpts, Propagator, RSSCartesianStep}, + Spacecraft, State, +}; +use std::{error::Error, sync::Arc}; + +fn main() -> Result<(), Box> { + pel::init(); + // Set up the dynamics like in the orbit raise. + let almanac = Arc::new(MetaAlmanac::latest().map_err(Box::new)?); + let epoch = Epoch::from_gregorian_utc_hms(2024, 2, 29, 12, 13, 14); + + // Define the GEO orbit, and we're just going to maintain it very tightly. + let earth_j2000 = almanac.frame_from_uid(EARTH_J2000)?; + let orbit = Orbit::try_keplerian(42164.0, 1e-5, 0., 163.0, 75.0, 0.0, epoch, earth_j2000)?; + println!("{orbit:x}"); + + let sc = Spacecraft::builder() + .orbit(orbit) + .dry_mass_kg(1000.0) // 1000 kg of dry mass + .fuel_mass_kg(1000.0) // 1000 kg of fuel, totalling 2.0 tons + .srp(SrpConfig::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 + .thruster(Thruster { + // "NEXT-STEP" row in Table 2 + isp_s: 4435.0, + thrust_N: 0.472, + }) + .mode(GuidanceMode::Thrust) // Start thrusting immediately. + .build(); + + // Set up the spacecraft dynamics like in the orbit raise example. + + 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 = &[ + Objective::within_tolerance(StateParameter::SMA, 42_164.0, 5.0), // 5 km + Objective::within_tolerance(StateParameter::Eccentricity, 0.001, 5e-5), + Objective::within_tolerance(StateParameter::Inclination, 0.05, 1e-2), + ]; + + let ruggiero_ctrl = Ruggiero::from_max_eclipse(objectives, sc, EclipseState::Penumbra(0.2))?; + println!("{ruggiero_ctrl}"); + + let mut orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN]); + + 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. + }; + jgm3_meta.process()?; + + let harmonics = Harmonics::from_stor( + almanac.frame_from_uid(IAU_EARTH_FRAME)?, + HarmonicsMem::from_cof(&jgm3_meta.uri, 8, 8, true)?, + ); + orbital_dyn.accel_models.push(harmonics); + + let srp_dyn = SolarPressure::default(EARTH_J2000, almanac.clone())?; + let sc_dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn) + .with_guidance_law(ruggiero_ctrl.clone()); + + println!("{sc_dynamics}"); + + // 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)], + )?; + + let my_mc = MonteCarlo::new( + sc, // Nominal state + mc_rv, + "03_geo_sk".to_string(), // Scenario name + None, // No specific seed specified, so one will be drawn from the computer's entropy. + ); + + // Build the propagator setup. + let setup = Propagator::rk89( + sc_dynamics.clone(), + PropOpts::builder() + .min_step(10.0_f64.seconds()) + .error_ctrl(RSSCartesianStep {}) + .build(), + ); + + 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); + + // For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra. + + rslts.to_parquet( + "03_geo_sk.parquet", + Some(vec![ + &EclipseLocator::cislunar(almanac.clone()).to_penumbra_event() + ]), + ExportCfg::default(), + almanac, + )?; + + Ok(()) +} diff --git a/src/cosmic/eclipse.rs b/src/cosmic/eclipse.rs index 93c4b8c3..1f5c1976 100644 --- a/src/cosmic/eclipse.rs +++ b/src/cosmic/eclipse.rs @@ -244,7 +244,7 @@ impl EventEvaluator for PenumbraEvent { { EclipseState::Umbra => Ok(0.0), EclipseState::Visibilis => Ok(1.0), - EclipseState::Penumbra(val) => Ok(val - 1.0), + EclipseState::Penumbra(val) => Ok(val), } } @@ -252,7 +252,7 @@ impl EventEvaluator for PenumbraEvent { fn epoch_precision(&self) -> Duration { 0.1 * Unit::Second } - /// Finds the slightest penumbra within 2%(i.e. 98% in visibility) + /// Finds the slightest penumbra within 2% (i.e. 98% in visibility) fn value_precision(&self) -> f64 { 0.02 } @@ -368,13 +368,13 @@ pub fn eclipse_state( let shadow_area = circ_seg_area(r_eb_prime, d1) + circ_seg_area(r_ls_prime, d2); if shadow_area.is_nan() { - warn!( + error!( "Shadow area is NaN! Please file a bug with initial states, eclipsing bodies, etc." ); return Ok(EclipseState::Umbra); } // Compute the nominal area of the light source - let nominal_area = std::f64::consts::PI * r_ls_prime.powi(2); + let nominal_area = core::f64::consts::PI * r_ls_prime.powi(2); // And return the percentage (between 0 and 1) of the eclipse. Ok(EclipseState::Penumbra(1.0 - shadow_area / nominal_area)) } else { diff --git a/src/dynamics/guidance/finiteburns.rs b/src/dynamics/guidance/finiteburns.rs index 88398daa..d2336f96 100644 --- a/src/dynamics/guidance/finiteburns.rs +++ b/src/dynamics/guidance/finiteburns.rs @@ -1,3 +1,4 @@ +use anise::prelude::Almanac; /* Nyx, blazing fast astrodynamics Copyright (C) 2018-onwards Christopher Rabotin @@ -93,7 +94,7 @@ impl GuidanceLaw for FiniteBurns { } } - fn next(&self, sc: &mut Spacecraft) { + fn next(&self, sc: &mut Spacecraft, _almanac: Arc) { // Grab the last maneuver if let Some(last_mnvr) = self.mnvrs.last() { // If the last maneuver ends before the current epoch, switch back into coast diff --git a/src/dynamics/guidance/mnvr.rs b/src/dynamics/guidance/mnvr.rs index 1a211611..61336295 100644 --- a/src/dynamics/guidance/mnvr.rs +++ b/src/dynamics/guidance/mnvr.rs @@ -25,9 +25,11 @@ use crate::linalg::Vector3; use crate::polyfit::CommonPolynomial; use crate::time::{Epoch, Unit}; use crate::State; +use anise::prelude::Almanac; use hifitime::{Duration, TimeUnits}; use snafu::ResultExt; use std::fmt; +use std::sync::Arc; /// Mnvr defined a single maneuver. Direction MUST be in the VNC frame (Velocity / Normal / Cross). /// It may be used with a maneuver scheduler. @@ -254,7 +256,7 @@ impl GuidanceLaw for Mnvr { } } - fn next(&self, sc: &mut Spacecraft) { + fn next(&self, sc: &mut Spacecraft, _almanac: Arc) { let next_mode = if sc.epoch() >= self.start && sc.epoch() < self.end { GuidanceMode::Thrust } else { diff --git a/src/dynamics/guidance/mod.rs b/src/dynamics/guidance/mod.rs index 56c9b764..b28ab856 100644 --- a/src/dynamics/guidance/mod.rs +++ b/src/dynamics/guidance/mod.rs @@ -22,6 +22,7 @@ use crate::linalg::Vector3; use anise::astro::PhysicsResult; use anise::errors::PhysicsError; use anise::math::rotation::DCM; +use anise::prelude::Almanac; use serde::{Deserialize, Serialize}; mod finiteburns; @@ -35,6 +36,7 @@ pub use ruggiero::{Objective, Ruggiero, StateParameter}; use snafu::Snafu; use std::fmt; +use std::sync::Arc; #[cfg(feature = "python")] use pyo3::prelude::*; @@ -78,7 +80,7 @@ pub trait GuidanceLaw: fmt::Display + Send + Sync { fn throttle(&self, osc_state: &Spacecraft) -> Result; /// Updates the state of the BaseSpacecraft for the next maneuver, e.g. prepares the controller for the next maneuver - fn next(&self, next_state: &mut Spacecraft); + fn next(&self, next_state: &mut Spacecraft, almanac: Arc); /// Returns whether this thrust control has been achieved, if it has an objective fn achieved(&self, _osc_state: &Spacecraft) -> Result { diff --git a/src/dynamics/guidance/ruggiero.rs b/src/dynamics/guidance/ruggiero.rs index 68f8d76f..d953b0c6 100644 --- a/src/dynamics/guidance/ruggiero.rs +++ b/src/dynamics/guidance/ruggiero.rs @@ -16,12 +16,14 @@ along with this program. If not, see . */ +use anise::prelude::Almanac; use snafu::ResultExt; use super::{ unit_vector_from_plane_angles, GuidStateSnafu, GuidanceError, GuidanceLaw, GuidanceMode, GuidancePhysicsSnafu, NyxError, Orbit, Spacecraft, Vector3, }; +use crate::cosmic::eclipse::{EclipseLocator, EclipseState}; pub use crate::md::objective::Objective; pub use crate::md::StateParameter; use crate::State; @@ -36,6 +38,8 @@ pub struct Ruggiero { pub objectives: [Option; 5], /// Stores the minimum efficiency to correct a given orbital element, defaults to zero (i.e. always correct) pub ηthresholds: [f64; 5], + /// If define, coast until vehicle is out of the provided eclipse state. + pub max_eclipse: Option, init_state: Spacecraft, } @@ -45,13 +49,14 @@ pub struct Ruggiero { impl Ruggiero { /// Creates a new Ruggiero locally optimal control as an Arc /// Note: this returns an Arc so it can be plugged into the Spacecraft dynamics directly. - pub fn new(objectives: &[Objective], initial: Spacecraft) -> Result, NyxError> { - Self::with_ηthresholds(objectives, &[0.0; 5], initial) + pub fn simple(objectives: &[Objective], initial: Spacecraft) -> Result, NyxError> { + Self::from_ηthresholds(objectives, &[0.0; 5], initial) } - /// Creates a new Ruggiero locally optimal control as an Arc + /// Creates a new Ruggiero locally optimal control with the provided efficiency threshold. + /// If the efficiency to correct the mapped orbital element is greater than the threshold, then the control law will be applied to this orbital element. /// Note: this returns an Arc so it can be plugged into the Spacecraft dynamics directly. - pub fn with_ηthresholds( + pub fn from_ηthresholds( objectives: &[Objective], ηthresholds: &[f64], initial: Spacecraft, @@ -100,9 +105,61 @@ impl Ruggiero { objectives: objs, init_state: initial, ηthresholds: eff, + max_eclipse: None, })) } + /// Creates a new Ruggiero locally optimal control as an Arc + /// Note: this returns an Arc so it can be plugged into the Spacecraft dynamics directly. + pub fn from_max_eclipse( + objectives: &[Objective], + initial: Spacecraft, + max_eclipse: EclipseState, + ) -> Result, NyxError> { + let mut objs: [Option; 5] = [None, None, None, None, None]; + let eff: [f64; 5] = [0.0; 5]; + if objectives.len() > 5 || objectives.is_empty() { + return Err(NyxError::GuidanceConfigError { + msg: format!( + "Must provide between 1 and 5 objectives (included), provided {}", + objectives.len() + ), + }); + } + + for (i, obj) in objectives.iter().enumerate() { + if [ + StateParameter::SMA, + StateParameter::Eccentricity, + StateParameter::Inclination, + StateParameter::RAAN, + StateParameter::AoP, + ] + .contains(&obj.parameter) + { + objs[i] = Some(*obj); + } else { + return Err(NyxError::GuidanceConfigError { + msg: format!("Objective {} not supported in Ruggerio", obj.parameter), + }); + } + } + for i in 0..objectives.len() { + objs[i] = Some(objectives[i]); + } + Ok(Arc::new(Self { + objectives: objs, + init_state: initial, + ηthresholds: eff, + max_eclipse: Some(max_eclipse), + })) + } + + /// Sets the maximum eclipse during which we can thrust. + pub fn set_max_eclipse(&mut self, max_eclipse: EclipseState) { + self.max_eclipse = Some(max_eclipse); + } + /// Returns the efficiency η ∈ [0; 1] of correcting a specific orbital element at the provided osculating orbit pub fn efficency(parameter: &StateParameter, osc_orbit: &Orbit) -> Result { let e = osc_orbit.ecc().context(GuidancePhysicsSnafu { @@ -183,11 +240,40 @@ impl Ruggiero { .abs() } } + + /// Returns whether the guidance law has achieved all goals + pub fn status(&self, state: &Spacecraft) -> Vec { + self.objectives + .iter() + .flatten() + .map(|obj| { + let (ok, err) = obj.assess(state).unwrap(); + format!( + "{} achieved: {}\t error = {:.5} {}", + obj, + ok, + err, + obj.parameter.unit() + ) + }) + .collect::>() + } } impl fmt::Display for Ruggiero { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "Ruggiero with {} objectives", self.objectives.len()) + let obj_msg = self + .objectives + .iter() + .flatten() + .map(|obj| format!("{obj}")) + .collect::>(); + write!( + f, + "Ruggiero Controller (max eclipse: {:?}): \n {}", + self.max_eclipse, + obj_msg.join("\n") + ) } } @@ -196,7 +282,7 @@ impl GuidanceLaw for Ruggiero { fn achieved(&self, state: &Spacecraft) -> Result { for obj in self.objectives.iter().flatten() { if !obj - .assess_raw(state.value(obj.parameter).context(GuidStateSnafu)?) + .assess_value(state.value(obj.parameter).context(GuidStateSnafu)?) .0 { return Ok(false); @@ -321,29 +407,40 @@ impl GuidanceLaw for Ruggiero { // Either thrust full power or not at all fn throttle(&self, sc: &Spacecraft) -> Result { if sc.mode() == GuidanceMode::Thrust { - for (i, obj) in self.objectives.iter().flatten().enumerate() { - let weight = self.weighting(obj, sc, self.ηthresholds[i]); - if weight.abs() > 0.0 { - return Ok(1.0); - } + if self.direction(sc)?.norm() > 0.0 { + Ok(1.0) + } else { + Ok(0.0) } - Ok(0.0) } else { Ok(0.0) } } /// Update the state for the next iteration - fn next(&self, sc: &mut Spacecraft) { + fn next(&self, sc: &mut Spacecraft, almanac: Arc) { if sc.mode() != GuidanceMode::Inhibit { if !self.achieved(sc).unwrap() { - if sc.mode() == GuidanceMode::Coast { - info!("enabling steering: {:x}", sc.orbit); + // Check eclipse state if applicable. + if let Some(max_eclipse) = self.max_eclipse { + let locator = EclipseLocator::cislunar(almanac.clone()); + if locator + .compute(sc.orbit, almanac) + .expect("cannot compute eclipse") + > max_eclipse + { + // Coast in eclipse + sc.mode = GuidanceMode::Coast; + } else { + sc.mode = GuidanceMode::Thrust; + } + } else if sc.mode() == GuidanceMode::Coast { + debug!("enabling steering: {:x}", sc.orbit); } sc.mut_mode(GuidanceMode::Thrust); } else { if sc.mode() == GuidanceMode::Thrust { - info!("disabling steering: {:x}", sc.orbit); + debug!("disabling steering: {:x}", sc.orbit); } sc.mut_mode(GuidanceMode::Coast); } @@ -367,7 +464,7 @@ fn ruggiero_weight() { Objective::within_tolerance(StateParameter::Eccentricity, 0.01, 5e-5), ]; - let ruggiero = Ruggiero::new(objectives, sc).unwrap(); + let ruggiero = Ruggiero::simple(objectives, sc).unwrap(); // 7301.597157 201.699933 0.176016 -0.202974 7.421233 0.006476 298.999726 let osc = Orbit::new( 7_303.253_461_441_64f64, diff --git a/src/dynamics/spacecraft.rs b/src/dynamics/spacecraft.rs index 18a401da..9e20f021 100644 --- a/src/dynamics/spacecraft.rs +++ b/src/dynamics/spacecraft.rs @@ -237,7 +237,7 @@ impl Dynamics for SpacecraftDynamics { fn finally( &self, next_state: Self::StateType, - _almanac: Arc, + almanac: Arc, ) -> Result { if next_state.fuel_mass_kg < 0.0 { error!("negative fuel mass at {}", next_state.epoch()); @@ -249,7 +249,7 @@ impl Dynamics for SpacecraftDynamics { if let Some(guid_law) = &self.guid_law { let mut state = next_state; // Update the control mode - guid_law.next(&mut state); + guid_law.next(&mut state, almanac.clone()); Ok(state) } else { Ok(next_state) diff --git a/src/mc/montecarlo.rs b/src/mc/montecarlo.rs index efb2b1a6..2d63efe8 100644 --- a/src/mc/montecarlo.rs +++ b/src/mc/montecarlo.rs @@ -251,6 +251,7 @@ where |(arc_prop, tx), (index, dispersed_state)| { let result = arc_prop .with(dispersed_state.state, almanac.clone()) + .quiet() .until_epoch_with_traj(end_epoch); // Build a single run result diff --git a/src/md/objective.rs b/src/md/objective.rs index 2c721ec0..e1b96b9e 100644 --- a/src/md/objective.rs +++ b/src/md/objective.rs @@ -17,7 +17,7 @@ */ use super::StateParameter; -use crate::cosmic::OrbitPartial; +use crate::{errors::StateError, Spacecraft, State}; use std::fmt; /// Defines a state parameter event finder @@ -59,13 +59,13 @@ impl Objective { } /// Returns whether this objective has been achieved, and the associated parameter error. - pub fn assess(&self, achieved: OrbitPartial) -> (bool, f64) { - self.assess_raw(achieved.real()) + pub fn assess(&self, achieved: &Spacecraft) -> Result<(bool, f64), StateError> { + Ok(self.assess_value(achieved.value(self.parameter)?)) } /// Returns whether this objective has been achieved, and the associated parameter error. /// Warning: the parameter `achieved` must be in the same unit as the objective. - pub fn assess_raw(&self, achieved: f64) -> (bool, f64) { + pub fn assess_value(&self, achieved: f64) -> (bool, f64) { let param_err = self.multiplicative_factor * (self.desired_value - achieved) + self.additive_factor; @@ -85,14 +85,15 @@ impl fmt::LowerHex for Objective { write!( f, - "{:?} → {:.prec$} ", + "{:?} → {:.prec$} {}", self.parameter, self.desired_value, - prec = max_obj_tol + self.parameter.unit(), + prec = max_obj_tol, )?; if self.tolerance.abs() < 1e-1 { - write!(f, "(± {:.1e})", self.tolerance) + write!(f, " (± {:.1e})", self.tolerance) } else { write!(f, " (± {:.2})", self.tolerance) } diff --git a/src/md/opti/raphson_finite_diff.rs b/src/md/opti/raphson_finite_diff.rs index 4f1570a7..88b9112a 100644 --- a/src/md/opti/raphson_finite_diff.rs +++ b/src/md/opti/raphson_finite_diff.rs @@ -279,7 +279,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { let achieved = partial.real(); - let (ok, param_err) = obj.assess_raw(achieved); + let (ok, param_err) = obj.assess_value(achieved); if !ok { converged = false; } diff --git a/src/md/opti/raphson_hyperdual.rs b/src/md/opti/raphson_hyperdual.rs index 4a3b877f..35da9725 100644 --- a/src/md/opti/raphson_hyperdual.rs +++ b/src/md/opti/raphson_hyperdual.rs @@ -188,7 +188,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> { let achieved = xf_partial.real(); - let (ok, param_err) = obj.assess_raw(achieved); + let (ok, param_err) = obj.assess_value(achieved); if !ok { converged = false; } diff --git a/src/od/process/mod.rs b/src/od/process/mod.rs index 200e0b51..6b4de35b 100644 --- a/src/od/process/mod.rs +++ b/src/od/process/mod.rs @@ -134,7 +134,7 @@ where ) -> Self { let init_state = prop.state; Self { - prop, + prop: prop.quiet(), kf, estimates: Vec::with_capacity(10_000), residuals: Vec::with_capacity(10_000), diff --git a/src/propagators/instance.rs b/src/propagators/instance.rs index 4f235ff0..9ca52c95 100644 --- a/src/propagators/instance.rs +++ b/src/propagators/instance.rs @@ -52,6 +52,8 @@ where pub prop: &'a Propagator<'a, D, E>, /// Stores the details of the previous integration step pub details: IntegrationDetails, + /// Should progress reports be logged + pub log_progress: bool, pub(crate) almanac: Arc, pub(crate) step_size: Duration, // Stores the adapted step for the _next_ call pub(crate) fixed_step: bool, @@ -66,6 +68,18 @@ where + Allocator::Size, ::Size> + Allocator::VecLength>, { + /// Sets this instance to not log progress + pub fn quiet(mut self) -> Self { + self.log_progress = false; + self + } + + /// Sets this instance to log progress + pub fn verbose(mut self) -> Self { + self.log_progress = true; + self + } + /// Allows setting the step size of the propagator pub fn set_step(&mut self, step_size: Duration, fixed: bool) { self.step_size = step_size; @@ -85,9 +99,10 @@ where #[cfg(not(target_arch = "wasm32"))] let tick = Instant::now(); - let log_progress = duration.abs() >= 7 * Unit::Day; + #[cfg(not(target_arch = "wasm32"))] + let mut prev_tick = Instant::now(); - if log_progress { + if self.log_progress { // Prevent the print spam for orbit determination cases info!("Propagating for {} until {}", duration, stop_time); } @@ -111,7 +126,7 @@ where // No propagation necessary #[cfg(not(target_arch = "wasm32"))] { - if log_progress { + if self.log_progress { let tock: Duration = tick.elapsed().into(); debug!("Done in {}", tock); } @@ -141,14 +156,27 @@ where #[cfg(not(target_arch = "wasm32"))] { - if log_progress { + if self.log_progress { let tock: Duration = tick.elapsed().into(); - debug!("Done in {}", tock); + info!("\t... done in {}", tock); } } return Ok(self.state); } else { + #[cfg(not(target_arch = "wasm32"))] + { + if self.log_progress { + let tock: Duration = prev_tick.elapsed().into(); + if tock.to_unit(Unit::Minute) > 1.0 { + // Report status every minute + let cur_epoch = self.state.epoch(); + let dur_to_go = (stop_time - cur_epoch).floor(Unit::Second * 1); + info!("\t... current epoch {}, remaing {}", cur_epoch, dur_to_go); + prev_tick = Instant::now(); + } + } + } self.single_step()?; // Publish to channel if provided if let Some(ref chan) = maybe_tx_chan { diff --git a/src/propagators/options.rs b/src/propagators/options.rs index cc28519c..1c52e1ab 100644 --- a/src/propagators/options.rs +++ b/src/propagators/options.rs @@ -21,6 +21,7 @@ use std::fmt; use crate::time::{Duration, Unit}; use super::{ErrorCtrl, RSSCartesianStep}; +use typed_builder::TypedBuilder; /// PropOpts stores the integrator options, including the minimum and maximum step sizes, and the /// max error size. @@ -29,15 +30,22 @@ use super::{ErrorCtrl, RSSCartesianStep}; /// methods. To use a fixed step integrator, initialize the options using `with_fixed_step`, and /// use whichever adaptive step integrator is desired. For example, initializing an RK45 with /// fixed step options will lead to an RK4 being used instead of an RK45. -#[derive(Clone, Copy, Debug)] +#[derive(Clone, Copy, Debug, TypedBuilder)] +#[builder(doc)] pub struct PropOpts { + #[builder(default_code = "60.0 * Unit::Second")] pub init_step: Duration, + #[builder(default_code = "0.001 * Unit::Second")] pub min_step: Duration, + #[builder(default_code = "2700.0 * Unit::Second")] pub max_step: Duration, + #[builder(default = 1e-12)] pub tolerance: f64, + #[builder(default = 50)] pub attempts: u8, + #[builder(default = false)] pub fixed_step: bool, - pub _errctrl: E, + pub error_ctrl: E, } impl PropOpts { @@ -47,7 +55,7 @@ impl PropOpts { min_step: Duration, max_step: Duration, tolerance: f64, - errctrl: E, + error_ctrl: E, ) -> Self { PropOpts { init_step: max_step, @@ -56,16 +64,21 @@ impl PropOpts { tolerance, attempts: 50, fixed_step: false, - _errctrl: errctrl, + error_ctrl, } } - pub fn with_adaptive_step_s(min_step: f64, max_step: f64, tolerance: f64, errctrl: E) -> Self { + pub fn with_adaptive_step_s( + min_step: f64, + max_step: f64, + tolerance: f64, + error_ctrl: E, + ) -> Self { Self::with_adaptive_step( min_step * Unit::Second, max_step * Unit::Second, tolerance, - errctrl, + error_ctrl, ) } @@ -116,7 +129,7 @@ impl PropOpts { tolerance: 0.0, fixed_step: true, attempts: 0, - _errctrl: RSSCartesianStep {}, + error_ctrl: RSSCartesianStep {}, } } @@ -151,7 +164,7 @@ impl Default for PropOpts { tolerance: 1e-12, attempts: 50, fixed_step: false, - _errctrl: RSSCartesianStep {}, + error_ctrl: RSSCartesianStep {}, } } } diff --git a/src/propagators/propagator.rs b/src/propagators/propagator.rs index 363948cc..9e0c0039 100644 --- a/src/propagators/propagator.rs +++ b/src/propagators/propagator.rs @@ -106,6 +106,7 @@ where error: 0.0, attempts: 1, }, + log_progress: true, almanac, step_size: self.opts.init_step, fixed_step: self.opts.fixed_step, diff --git a/tests/propagation/trajectory.rs b/tests/propagation/trajectory.rs index 2f99e6cb..1e3f01de 100644 --- a/tests/propagation/trajectory.rs +++ b/tests/propagation/trajectory.rs @@ -270,7 +270,7 @@ fn traj_spacecraft(almanac: Arc) { 5e-3, )]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); // Build the spacecraft state let fuel_mass = 67.0; @@ -301,7 +301,7 @@ fn traj_spacecraft(almanac: Arc) { for mut sc_state in traj.every(1 * Unit::Day) { // We need to evaluate the mode of this state because the trajectory does not store discrete information - ruggiero_ctrl.next(&mut sc_state); + ruggiero_ctrl.next(&mut sc_state, almanac.clone()); if sc_state.mode() != prev_mode { println!( "Mode changed from {:?} to {:?} @ {}", @@ -316,7 +316,7 @@ fn traj_spacecraft(almanac: Arc) { for epoch in TimeSeries::inclusive(start_dt, start_dt + prop_time, 1 * Unit::Day) { // Note: the `evaluate` function will return a Result which prevents a panic if you request something out of the ephemeris let mut sc_state = traj.at(epoch).unwrap(); - ruggiero_ctrl.next(&mut sc_state); + ruggiero_ctrl.next(&mut sc_state, almanac.clone()); if sc_state.mode() != prev_mode { println!( "Mode changed from {:?} to {:?} @ {}", diff --git a/tests/propulsion/closedloop_multi_oe_ruggiero.rs b/tests/propulsion/closedloop_multi_oe_ruggiero.rs index d23703cf..feae05a8 100644 --- a/tests/propulsion/closedloop_multi_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_multi_oe_ruggiero.rs @@ -55,7 +55,7 @@ fn qlaw_as_ruggiero_case_a(almanac: Arc) { Event::within_tolerance(StateParameter::Eccentricity, 0.01, 5e-5), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let dry_mass = 1.0; let fuel_mass = 299.0; @@ -117,7 +117,7 @@ fn qlaw_as_ruggiero_case_b(almanac: Arc) { Objective::within_tolerance(StateParameter::Inclination, 0.05, 5e-3), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 1999.9; let dry_mass = 0.1; @@ -172,7 +172,7 @@ fn qlaw_as_ruggiero_case_c(almanac: Arc) { Objective::within_tolerance(StateParameter::Eccentricity, 0.7, 5e-5), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 299.9; let dry_mass = 0.1; @@ -230,7 +230,7 @@ fn qlaw_as_ruggiero_case_d(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 360.0 - 90.0, 5e-3), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -290,7 +290,7 @@ fn qlaw_as_ruggiero_case_e(almanac: Arc) { Objective::within_tolerance(StateParameter::AoP, 180.0, 5e-3), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 1999.9; let dry_mass = 0.1; @@ -346,7 +346,7 @@ fn qlaw_as_ruggiero_case_f(almanac: Arc) { let objectives = &[Objective::new(StateParameter::Eccentricity, 0.15)]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -406,7 +406,7 @@ fn ruggiero_iepc_2011_102(almanac: Arc) { Objective::within_tolerance(StateParameter::Eccentricity, 0.011, 5e-5), ]; - let ruggiero_ctrl = Ruggiero::new(objectives, orbit.into()).unwrap(); + let ruggiero_ctrl = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; diff --git a/tests/propulsion/closedloop_single_oe_ruggiero.rs b/tests/propulsion/closedloop_single_oe_ruggiero.rs index d74fb853..6116c134 100644 --- a/tests/propulsion/closedloop_single_oe_ruggiero.rs +++ b/tests/propulsion/closedloop_single_oe_ruggiero.rs @@ -43,7 +43,7 @@ fn rugg_sma(almanac: Arc) { 1.0, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -96,7 +96,7 @@ fn rugg_sma_regress_threshold(almanac: Arc) { let fuel_mass = 67.0; let dry_mass = 300.0; for (threshold, expected_fuel_usage) in &[(0.9, 16.9), (0.0, 21.3)] { - let guid_law = Ruggiero::with_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); + let guid_law = Ruggiero::from_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); let sc_state = Spacecraft::from_thruster(orbit, dry_mass, fuel_mass, lowt, GuidanceMode::Thrust); @@ -147,7 +147,7 @@ fn rugg_sma_decr(almanac: Arc) { 1.0, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -202,7 +202,7 @@ fn rugg_inc(almanac: Arc) { 5e-3, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -258,7 +258,7 @@ fn rugg_inc_threshold(almanac: Arc) { 5e-3, )]; - let guid_law = Ruggiero::with_ηthresholds(objectives, &[0.9], orbit.into()).unwrap(); + let guid_law = Ruggiero::from_ηthresholds(objectives, &[0.9], orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -313,7 +313,7 @@ fn rugg_inc_decr(almanac: Arc) { 5e-3, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -368,7 +368,7 @@ fn rugg_ecc(almanac: Arc) { 5e-5, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -424,7 +424,7 @@ fn rugg_ecc_regress_threshold(almanac: Arc) { )]; for (threshold, expected_fuel_usage) in &[(0.9, 8.2), (0.0, 10.37)] { - let guid_law = Ruggiero::with_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); + let guid_law = Ruggiero::from_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); let sc_state = Spacecraft::from_thruster(orbit, dry_mass, fuel_mass, lowt, GuidanceMode::Thrust); @@ -478,7 +478,7 @@ fn rugg_ecc_decr(almanac: Arc) { 5e-5, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -535,7 +535,7 @@ fn rugg_aop(almanac: Arc) { 5e-3, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -591,7 +591,7 @@ fn rugg_aop_decr(almanac: Arc) { 5e-3, )]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -644,7 +644,7 @@ fn rugg_raan(almanac: Arc) { // Define the objectives let objectives = &[Objective::within_tolerance(StateParameter::RAAN, 5.0, 5e-5)]; - let guid_law = Ruggiero::new(objectives, orbit.into()).unwrap(); + let guid_law = Ruggiero::simple(objectives, orbit.into()).unwrap(); let fuel_mass = 67.0; let dry_mass = 300.0; @@ -696,7 +696,7 @@ fn rugg_raan_regress_threshold(almanac: Arc) { let objectives = &[Objective::within_tolerance(StateParameter::RAAN, 5.0, 5e-5)]; for (threshold, expected_fuel_usage) in &[(0.9, 14.787), (0.0, 22.189)] { - let guid_law = Ruggiero::with_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); + let guid_law = Ruggiero::from_ηthresholds(objectives, &[*threshold], orbit.into()).unwrap(); let sc_state = Spacecraft::from_thruster(orbit, dry_mass, fuel_mass, lowt, GuidanceMode::Thrust);