Update tests for new integrator method enum

Author: ChristopherRabotin

src/mc/montecarlo.rs

@@ -89,9 +89,9 @@ where
 
     /// Generate states and propagate each independently until a specific event is found `trigger` times.
     #[allow(clippy::needless_lifetimes)]
-    pub fn run_until_nth_event<'a, D, E, F>(
+    pub fn run_until_nth_event<D, E, F>(
         self,
-        prop: Propagator<'a, D, E>,
+        prop: Propagator<D, E>,
         almanac: Arc<Almanac>,
         max_duration: Duration,
         event: &F,
@@ -113,9 +113,9 @@ where
     /// Generate states and propagate each independently until a specific event is found `trigger` times.
     #[must_use = "Monte Carlo result must be used"]
     #[allow(clippy::needless_lifetimes)]
-    pub fn resume_run_until_nth_event<'a, D, E, F>(
+    pub fn resume_run_until_nth_event<D, E, F>(
         &self,
-        prop: Propagator<'a, D, E>,
+        prop: Propagator<D, E>,
         almanac: Arc<Almanac>,
         skip: usize,
         max_duration: Duration,
@@ -188,9 +188,9 @@ where
     /// Generate states and propagate each independently until a specific event is found `trigger` times.
     #[must_use = "Monte Carlo result must be used"]
     #[allow(clippy::needless_lifetimes)]
-    pub fn run_until_epoch<'a, D, E>(
+    pub fn run_until_epoch<D, E>(
         self,
-        prop: Propagator<'a, D, E>,
+        prop: Propagator<D, E>,
         almanac: Arc<Almanac>,
         end_epoch: Epoch,
         num_runs: usize,
@@ -209,9 +209,9 @@ where
     /// Resumes a Monte Carlo run by skipping the first `skip` items, generating states only after that, and propagate each independently until the specified epoch.
     #[must_use = "Monte Carlo result must be used"]
     #[allow(clippy::needless_lifetimes)]
-    pub fn resume_run_until_epoch<'a, D, E>(
+    pub fn resume_run_until_epoch<D, E>(
         &self,
-        prop: Propagator<'a, D, E>,
+        prop: Propagator<D, E>,
         almanac: Arc<Almanac>,
         skip: usize,
         end_epoch: Epoch,

src/md/opti/multipleshooting/altitude_heuristic.rs

@@ -42,7 +42,7 @@ impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> {
         node_count: usize,
         angular_velocity_deg_s: f64,
         body_frame: Frame,
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         almanac: Arc<Almanac>,
     ) -> Result<Self, MultipleShootingError> {
         if node_count < 3 {

src/md/opti/multipleshooting/equidistant_heuristic.rs

@@ -33,7 +33,7 @@ impl<'a, E: ErrorCtrl> MultipleShooting<'a, E, Node, 3, 3> {
         x0: Spacecraft,
         xf: Orbit,
         node_count: usize,
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
     ) -> Result<Self, TargetingError> {
         if node_count < 3 {
             error!("At least three nodes are needed for a multiple shooting optimization");

src/md/opti/multipleshooting/multishoot.rs

@@ -47,7 +47,7 @@ pub struct MultipleShooting<
     const OT: usize,
 > {
     /// The propagator setup (kind, stages, etc.)
-    pub prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+    pub prop: &'a Propagator<SpacecraftDynamics, E>,
     /// List of nodes of the optimal trajectory
     pub targets: Vec<T>,
     /// Starting point, must be a spacecraft equipped with a thruster
@@ -357,9 +357,9 @@ impl<T: MultishootNode<O>, const O: usize> fmt::Display for MultipleShootingSolu
 impl<T: MultishootNode<O>, const O: usize> MultipleShootingSolution<T, O> {
     /// Allows building the trajectories between different nodes
     /// This will rebuild the targeters and apply the solutions sequentially
-    pub fn build_trajectories<'a, E: ErrorCtrl>(
+    pub fn build_trajectories<E: ErrorCtrl>(
         &self,
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &Propagator<SpacecraftDynamics, E>,
         almanac: Arc<Almanac>,
     ) -> Result<Vec<ScTraj>, MultipleShootingError> {
         let mut trajz = Vec::with_capacity(self.nodes.len());

src/md/opti/optimizer.rs

@@ -35,7 +35,7 @@ use super::solution::TargeterSolution;
 #[derive(Clone)]
 pub struct Optimizer<'a, E: ErrorCtrl, const V: usize, const O: usize> {
     /// The propagator setup (kind, stages, etc.)
-    pub prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+    pub prop: &'a Propagator<SpacecraftDynamics, E>,
     /// The list of objectives of this targeter
     pub objectives: [Objective; O],
     /// An optional frame (and Cosm) to compute the objectives in.
@@ -68,7 +68,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> fmt::Display for Optimize
 impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> {
     /// Create a new Targeter which will apply an impulsive delta-v correction.
     pub fn delta_v(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         objectives: [Objective; O],
     ) -> Self {
         Self {
@@ -87,7 +87,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> {
 
     /// Create a new Targeter which will MOVE the position of the spacecraft at the correction epoch
     pub fn delta_r(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         objectives: [Objective; O],
     ) -> Self {
         Self {
@@ -105,10 +105,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> {
     }
 
     /// Create a new Targeter which will apply an impulsive delta-v correction on all components of the VNC frame. By default, max step is 0.5 km/s.
-    pub fn vnc(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
-        objectives: [Objective; O],
-    ) -> Self {
+    pub fn vnc(prop: &'a Propagator<SpacecraftDynamics, E>, objectives: [Objective; O]) -> Self {
         Self {
             prop,
             objectives,
@@ -127,7 +124,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 3, O> {
 impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 4, O> {
     /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
     pub fn thrust_dir(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         objectives: [Objective; O],
     ) -> Self {
         Self {
@@ -149,7 +146,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 4, O> {
 impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 7, O> {
     /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
     pub fn thrust_dir_rate(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         objectives: [Objective; O],
     ) -> Self {
         Self {
@@ -174,7 +171,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 7, O> {
 impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 10, O> {
     /// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
     pub fn thrust_profile(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         objectives: [Objective; O],
     ) -> Self {
         Self {
@@ -202,7 +199,7 @@ impl<'a, E: ErrorCtrl, const O: usize> Optimizer<'a, E, 10, O> {
 impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> {
     /// Create a new Targeter which will apply an impulsive delta-v correction.
     pub fn new(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         variables: [Variable; V],
         objectives: [Objective; O],
     ) -> Self {
@@ -218,7 +215,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> {
 
     /// Create a new Targeter which will apply an impulsive delta-v correction.
     pub fn in_frame(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         variables: [Variable; V],
         objectives: [Objective; O],
         objective_frame: Frame,
@@ -235,7 +232,7 @@ impl<'a, E: ErrorCtrl, const V: usize, const O: usize> Optimizer<'a, E, V, O> {
 
     /// Create a new Targeter which will apply an impulsive delta-v correction on the specified components of the VNC frame.
     pub fn vnc_with_components(
-        prop: &'a Propagator<'a, SpacecraftDynamics, E>,
+        prop: &'a Propagator<SpacecraftDynamics, E>,
         variables: [Variable; V],
         objectives: [Objective; O],
     ) -> Self {

src/propagators/instance.rs

@@ -48,7 +48,7 @@ where
     /// The state of this propagator instance
     pub state: D::StateType,
     /// The propagator setup (kind, stages, etc.)
-    pub prop: &'a Propagator<'a, D, E>,
+    pub prop: &'a Propagator<D, E>,
     /// Stores the details of the previous integration step
     pub details: IntegrationDetails,
     /// Should progress reports be logged
@@ -342,14 +342,14 @@ where
                 .context(DynamicsSnafu)?;
             self.k[0] = ki;
             let mut a_idx: usize = 0;
-            for i in 0..(self.prop.stages - 1) {
+            for i in 0..(self.prop.method.stages() - 1) {
                 // Let's compute the c_i by summing the relevant items from the list of coefficients.
                 // \sum_{j=1}^{i-1} a_ij  ∀ i ∈ [2, s]
                 let mut ci: f64 = 0.0;
                 // The wi stores the a_{s1} * k_1 + a_{s2} * k_2 + ... + a_{s, s-1} * k_{s-1} +
                 let mut wi = OVector::<f64, <D::StateType as State>::VecLength>::from_element(0.0);
                 for kj in &self.k[0..i + 1] {
-                    let a_ij = self.prop.a_coeffs[a_idx];
+                    let a_ij = self.prop.method.a_coeffs()[a_idx];
                     ci += a_ij;
                     wi += a_ij * kj;
                     a_idx += 1;
@@ -374,9 +374,9 @@ where
             let mut error_est =
                 OVector::<f64, <D::StateType as State>::VecLength>::from_element(0.0);
             for (i, ki) in self.k.iter().enumerate() {
-                let b_i = self.prop.b_coeffs[i];
+                let b_i = self.prop.method.b_coeffs()[i];
                 if !self.fixed_step {
-                    let b_i_star = self.prop.b_coeffs[i + self.prop.stages];
+                    let b_i_star = self.prop.method.b_coeffs()[i + self.prop.method.stages()];
                     error_est += step_size * (b_i - b_i_star) * ki;
                 }
                 next_state += step_size * b_i * ki;
@@ -407,7 +407,7 @@ where
                         let proposed_step = 0.9
                             * step_size
                             * (self.prop.opts.tolerance / self.details.error)
-                                .powf(1.0 / f64::from(self.prop.order));
+                                .powf(1.0 / f64::from(self.prop.method.order()));
                         step_size = if proposed_step > self.prop.opts.max_step.to_seconds() {
                             self.prop.opts.max_step.to_seconds()
                         } else {
@@ -424,7 +424,7 @@ where
                     let proposed_step = 0.9
                         * step_size
                         * (self.prop.opts.tolerance / self.details.error)
-                            .powf(1.0 / f64::from(self.prop.order - 1));
+                            .powf(1.0 / f64::from(self.prop.method.order() - 1));
                     step_size = if proposed_step < self.prop.opts.min_step.to_seconds() {
                         self.prop.opts.min_step.to_seconds()
                     } else {

src/propagators/propagator.rs

@@ -21,7 +21,7 @@ use std::sync::Arc;
 use anise::almanac::Almanac;
 
 use super::error_ctrl::{ErrorCtrl, RSSCartesianStep};
-use super::{Dormand78, IntegrationDetails, PropInstance, PropOpts, RK, RK89};
+use super::{IntegrationDetails, IntegratorMethod, PropInstance, PropOpts};
 use crate::dynamics::Dynamics;
 use crate::linalg::allocator::Allocator;
 use crate::linalg::{DefaultAllocator, OVector};
@@ -32,7 +32,7 @@ use crate::State;
 /// It is an EventTracker, without any event tracking. It includes the options, the integrator
 /// details of the previous step, and the set of coefficients used for the monomorphic instance.
 #[derive(Clone, Debug)]
-pub struct Propagator<'a, D: Dynamics, E: ErrorCtrl>
+pub struct Propagator<D: Dynamics, E: ErrorCtrl>
 where
     DefaultAllocator: Allocator<<D::StateType as State>::Size>
         + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size>
@@ -41,29 +41,23 @@ where
 {
     pub dynamics: D, // Stores the dynamics used. *Must* use this to get the latest values
     pub opts: PropOpts<E>, // Stores the integration options (tolerance, min/max step, init step, etc.)
-    pub(crate) order: u8,  // Order of the integrator
-    pub(crate) stages: usize, // Number of stages, i.e. how many times the derivatives will be called
-    pub(crate) a_coeffs: &'a [f64],
-    pub(crate) b_coeffs: &'a [f64],
+    pub method: IntegratorMethod,
 }
 
 /// The `Propagator` trait defines the functions of a propagator and of an event tracker.
-impl<'a, D: Dynamics, E: ErrorCtrl> Propagator<'a, D, E>
+impl<D: Dynamics, E: ErrorCtrl> Propagator<D, E>
 where
     DefaultAllocator: Allocator<<D::StateType as State>::Size>
         + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size>
         + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size>
         + Allocator<<D::StateType as State>::VecLength>,
 {
     /// Each propagator must be initialized with `new` which stores propagator information.
-    pub fn new<T: RK>(dynamics: D, opts: PropOpts<E>) -> Self {
+    pub fn new(dynamics: D, method: IntegratorMethod, opts: PropOpts<E>) -> Self {
         Self {
             dynamics,
             opts,
-            stages: T::STAGES,
-            order: T::ORDER,
-            a_coeffs: T::A_COEFFS,
-            b_coeffs: T::B_COEFFS,
+            method,
         }
     }
 
@@ -83,19 +77,19 @@ where
 
     /// An RK89 propagator (the default) with custom propagator options.
     pub fn rk89(dynamics: D, opts: PropOpts<E>) -> Self {
-        Self::new::<RK89>(dynamics, opts)
+        Self::new(dynamics, IntegratorMethod::RungeKutta89, opts)
     }
 
     /// A Dormand Prince 7-8 propagator with custom propagator options: it's about 20% faster than an RK98, and more stable in two body dynamics.
     /// WARNINGS: Dormand Prince may have issues with generating proper trajectories, leading to glitches in event finding.
     pub fn dp78(dynamics: D, opts: PropOpts<E>) -> Self {
-        Self::new::<Dormand78>(dynamics, opts)
+        Self::new(dynamics, IntegratorMethod::DormandPrince78, opts)
     }
 
-    pub fn with(&'a self, state: D::StateType, almanac: Arc<Almanac>) -> PropInstance<'a, D, E> {
+    pub fn with(&self, state: D::StateType, almanac: Arc<Almanac>) -> PropInstance<D, E> {
         // Pre-allocate the k used in the propagator
-        let mut k = Vec::with_capacity(self.stages + 1);
-        for _ in 0..self.stages {
+        let mut k = Vec::with_capacity(self.method.stages() + 1);
+        for _ in 0..self.method.stages() {
             k.push(OVector::<f64, <D::StateType as State>::VecLength>::zeros());
         }
         PropInstance {
@@ -115,7 +109,7 @@ where
     }
 }
 
-impl<'a, D: Dynamics> Propagator<'a, D, RSSCartesianStep>
+impl<D: Dynamics> Propagator<D, RSSCartesianStep>
 where
     DefaultAllocator: Allocator<<D::StateType as State>::Size>
         + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size>
@@ -124,13 +118,13 @@ where
 {
     /// Default propagator is an RK89 with the default PropOpts.
     pub fn default(dynamics: D) -> Self {
-        Self::new::<RK89>(dynamics, PropOpts::default())
+        Self::rk89(dynamics, PropOpts::default())
     }
 
     /// A default Dormand Prince 78 propagator with the default PropOpts.
     /// Faster and more stable than an RK89 (`default`) but seems to cause issues for event finding.
     /// WARNINGS: Dormand Prince may have issues with generating proper trajectories, leading to glitches in event finding.
     pub fn default_dp78(dynamics: D) -> Self {
-        Self::new::<Dormand78>(dynamics, PropOpts::default())
+        Self::dp78(dynamics, PropOpts::default())
     }
 }

src/propagators/rk_methods/mod.rs

@@ -84,7 +84,7 @@ impl IntegratorMethod {
         }
     }
 
-    /// Returns the stages of this integrator (as usize because it's used as indexing)
+    /// Returns the stages of this integrator, i.e. how many times the derivatives will be called
     pub const fn stages(self) -> usize {
         match self {
             Self::RungeKutta89 => RK89::STAGES,

tests/mission_design/orbitaldyn.rs

@@ -53,8 +53,9 @@ fn energy_conservation(almanac: Arc<Almanac>) {
         eme2k,
     );
 
-    let rk89_final = Propagator::new::<RK89>(
+    let rk89_final = Propagator::new(
         SpacecraftDynamics::new(OrbitalDynamics::two_body()),
+        IntegratorMethod::RungeKutta89,
         PropOpts::default(),
     )
     .with(Spacecraft::from(start_state), almanac.clone())
@@ -74,8 +75,9 @@ fn energy_conservation(almanac: Arc<Almanac>) {
     }
     println!();
 
-    let dp78_final = Propagator::new::<Dormand78>(
+    let dp78_final = Propagator::new(
         SpacecraftDynamics::new(OrbitalDynamics::two_body()),
+        IntegratorMethod::DormandPrince78,
         PropOpts::default(),
     )
     .with(start_state.into(), almanac)
@@ -1196,7 +1198,11 @@ fn val_cislunar_dynamics(almanac_gmat: Arc<Almanac>) {
     );
 
     let dynamics = SpacecraftDynamics::new(OrbitalDynamics::point_masses(vec![EARTH, SUN, MOON]));
-    let setup = Propagator::new::<RK4Fixed>(dynamics, PropOpts::with_fixed_step_s(0.5));
+    let setup = Propagator::new(
+        dynamics,
+        IntegratorMethod::RungeKutta4,
+        PropOpts::with_fixed_step_s(0.5),
+    );
     let mut prop = setup.with(state.into(), almanac);
     prop.for_duration(prop_time).unwrap();