From ef56369a3e392fef10c0f66d080f4bed284bd727 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 29 Sep 2024 23:10:51 -0600 Subject: [PATCH] Add integration frame and serialization to integrator options --- Cargo.toml | 3 +- src/cosmic/mod.rs | 6 ++ src/cosmic/spacecraft.rs | 8 ++ src/dynamics/mod.rs | 1 + src/md/mod.rs | 4 +- src/md/opti/mod.rs | 8 +- src/md/opti/multipleshooting/multishoot.rs | 14 +-- src/md/opti/raphson_finite_diff.rs | 4 +- src/md/opti/raphson_hyperdual.rs | 2 +- src/md/opti/{optimizer.rs => targeter.rs} | 16 ++-- src/md/trajectory/interpolatable.rs | 7 -- src/md/trajectory/traj.rs | 4 +- src/propagators/error_ctrl.rs | 6 ++ src/propagators/instance.rs | 65 +++++++++++-- src/propagators/options.rs | 91 ++++++++++++------- src/propagators/rk_methods/dormand.rs | 4 +- src/propagators/rk_methods/mod.rs | 6 +- src/propagators/rk_methods/rk.rs | 6 +- src/propagators/rk_methods/verner.rs | 2 +- tests/mission_design/force_models.rs | 37 +++++++- tests/mission_design/targeter/b_plane.rs | 8 +- tests/mission_design/targeter/finite_burns.rs | 10 +- tests/mission_design/targeter/multi_oe.rs | 8 +- tests/mission_design/targeter/multi_oe_vnc.rs | 6 +- tests/mission_design/targeter/single_oe.rs | 20 ++-- 25 files changed, 236 insertions(+), 110 deletions(-) rename src/md/opti/{optimizer.rs => targeter.rs} (96%) diff --git a/Cargo.toml b/Cargo.toml index 909b56b7..fa895fa1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -77,12 +77,13 @@ typed-builder = "0.20.0" pythonize = { version = "0.21", optional = true } snafu = { version = "0.8.3", features = ["backtrace"] } serde_dhall = "0.12" -toml = "0.8.14" + [dev-dependencies] polars = { version = "0.42.0", features = ["parquet"] } rstest = "0.22.0" pretty_env_logger = "0.5" +toml = "0.8.14" [build-dependencies] shadow-rs = "0.35.0" diff --git a/src/cosmic/mod.rs b/src/cosmic/mod.rs index da23739f..4352e9b7 100644 --- a/src/cosmic/mod.rs +++ b/src/cosmic/mod.rs @@ -110,6 +110,12 @@ where fn set_value(&mut self, param: StateParameter, _val: f64) -> Result<(), StateError> { Err(StateError::Unavailable { param }) } + + /// Returns a copy of the orbit + fn orbit(&self) -> Orbit; + + /// Modifies this state's orbit + fn set_orbit(&mut self, _orbit: Orbit) {} } pub fn assert_orbit_eq_or_abs(left: &Orbit, right: &Orbit, epsilon: f64, msg: &str) { diff --git a/src/cosmic/spacecraft.rs b/src/cosmic/spacecraft.rs index eb7d0204..fcc4d16e 100644 --- a/src/cosmic/spacecraft.rs +++ b/src/cosmic/spacecraft.rs @@ -761,6 +761,14 @@ impl State for Spacecraft { fn unset_stm(&mut self) { self.stm = None; } + + fn orbit(&self) -> Orbit { + self.orbit + } + + fn set_orbit(&mut self, orbit: Orbit) { + self.orbit = orbit; + } } impl Add>> for Spacecraft { diff --git a/src/dynamics/mod.rs b/src/dynamics/mod.rs index 73f44e8a..a356a5ee 100644 --- a/src/dynamics/mod.rs +++ b/src/dynamics/mod.rs @@ -172,6 +172,7 @@ pub trait AccelModel: Send + Sync + fmt::Display { /// Stores dynamical model errors #[derive(Debug, PartialEq, Snafu)] +#[snafu(visibility(pub(crate)))] pub enum DynamicsError { /// Fuel exhausted at the provided spacecraft state #[snafu(display("fuel exhausted at {sc}"))] diff --git a/src/md/mod.rs b/src/md/mod.rs index ec7a6377..7045d5b7 100644 --- a/src/md/mod.rs +++ b/src/md/mod.rs @@ -25,7 +25,7 @@ use snafu::prelude::*; pub mod prelude { pub use super::{ - optimizer::*, + targeter::*, trajectory::{ExportCfg, Interpolatable, Traj}, Event, ScTraj, StateParameter, }; @@ -52,7 +52,7 @@ pub use events::{Event, EventEvaluator}; pub mod objective; pub mod opti; -pub use opti::optimizer; +pub use opti::targeter; pub type ScTraj = trajectory::Traj; // pub type Ephemeris = trajectory::Traj; diff --git a/src/md/opti/mod.rs b/src/md/opti/mod.rs index 7f09afa0..1524a795 100644 --- a/src/md/opti/mod.rs +++ b/src/md/opti/mod.rs @@ -19,16 +19,16 @@ // pub mod convert_impulsive; pub mod multipleshooting; pub use multipleshooting::{ctrlnodes, multishoot}; -/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem. -// #[cfg(feature = "broken-donotuse")] -// pub mod minimize_lm; -pub mod optimizer; /// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via finite differencing. pub mod raphson_finite_diff; /// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via hyperdual numbers. pub mod raphson_hyperdual; pub mod solution; pub mod target_variable; +/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem. +// #[cfg(feature = "broken-donotuse")] +// pub mod minimize_lm; +pub mod targeter; #[derive(Copy, Clone, Debug, PartialEq, Eq)] pub enum DiffMethod { diff --git a/src/md/opti/multipleshooting/multishoot.rs b/src/md/opti/multipleshooting/multishoot.rs index 579bf3e7..c31ac95b 100644 --- a/src/md/opti/multipleshooting/multishoot.rs +++ b/src/md/opti/multipleshooting/multishoot.rs @@ -22,7 +22,7 @@ pub use super::CostFunction; use super::{MultipleShootingError, TargetingSnafu}; use crate::linalg::{DMatrix, DVector, SVector}; use crate::md::opti::solution::TargeterSolution; -use crate::md::optimizer::Optimizer; +use crate::md::targeter::Targeter; use crate::md::{prelude::*, TargetingError}; use crate::pseudo_inverse; use crate::{Orbit, Spacecraft}; @@ -81,7 +81,7 @@ impl<'a, T: MultishootNode, const VT: usize, const OT: usize> MultipleShooti /* *** ** 1. Solve the delta-v differential corrector between each node ** *** */ - let tgt = Optimizer { + let tgt = Targeter { prop: self.prop, objectives: self.targets[i].into(), variables: self.variables, @@ -125,7 +125,7 @@ impl<'a, T: MultishootNode, const VT: usize, const OT: usize> MultipleShooti ** Note that because the first initial_state is x0, the i-th "initial state" ** is the initial state to reach the i-th node. ** *** */ - let inner_tgt_a = Optimizer::delta_v(self.prop, next_node); + let inner_tgt_a = Targeter::delta_v(self.prop, next_node); let inner_sol_a = inner_tgt_a .try_achieve_dual( initial_states[i], @@ -151,7 +151,7 @@ impl<'a, T: MultishootNode, const VT: usize, const OT: usize> MultipleShooti /* *** ** 2.C. Rerun the targeter from the new state at the perturbed node to the next unpertubed node ** *** */ - let inner_tgt_b = Optimizer::delta_v(self.prop, self.targets[i + 1].into()); + let inner_tgt_b = Targeter::delta_v(self.prop, self.targets[i + 1].into()); let inner_sol_b = inner_tgt_b .try_achieve_dual( inner_sol_a.achieved_state, @@ -241,7 +241,7 @@ impl<'a, T: MultishootNode, const VT: usize, const OT: usize> MultipleShooti for (i, node) in self.targets.iter().enumerate() { // Run the unpertubed targeter - let tgt = Optimizer::delta_v(self.prop, (*node).into()); + let tgt = Targeter::delta_v(self.prop, (*node).into()); let sol = tgt .try_achieve_dual( initial_states[i], @@ -355,8 +355,8 @@ impl, const O: usize> MultipleShootingSolution { ) -> Result, MultipleShootingError> { let mut trajz = Vec::with_capacity(self.nodes.len()); - for (i, node) in self.nodes.iter().enumerate() { - let (_, traj) = Optimizer::delta_v(prop, (*node).into()) + for (i, node) in self.nodes.iter().copied().enumerate() { + let (_, traj) = Targeter::delta_v(prop, node.into()) .apply_with_traj(&self.solutions[i], almanac.clone()) .context(TargetingSnafu { segment: i })?; trajz.push(traj); diff --git a/src/md/opti/raphson_finite_diff.rs b/src/md/opti/raphson_finite_diff.rs index 5d3f3104..15aa5558 100644 --- a/src/md/opti/raphson_finite_diff.rs +++ b/src/md/opti/raphson_finite_diff.rs @@ -16,7 +16,7 @@ along with this program. If not, see . */ -use super::optimizer::Optimizer; +use super::targeter::Targeter; use super::solution::TargeterSolution; use crate::cosmic::{AstroAlmanacSnafu, AstroPhysicsSnafu}; use crate::dynamics::guidance::{GuidanceError, LocalFrame, Mnvr}; @@ -33,7 +33,7 @@ use snafu::{ensure, ResultExt}; #[cfg(not(target_arch = "wasm32"))] use std::time::Instant; -impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Differential correction using finite differencing #[allow(clippy::comparison_chain)] pub fn try_achieve_fd( diff --git a/src/md/opti/raphson_hyperdual.rs b/src/md/opti/raphson_hyperdual.rs index 7d4ffa30..eda4ca32 100644 --- a/src/md/opti/raphson_hyperdual.rs +++ b/src/md/opti/raphson_hyperdual.rs @@ -30,7 +30,7 @@ use crate::utils::are_eigenvalues_stable; #[cfg(not(target_arch = "wasm32"))] use std::time::Instant; -impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Differential correction using hyperdual numbers for the objectives #[allow(clippy::comparison_chain)] pub fn try_achieve_dual( diff --git a/src/md/opti/optimizer.rs b/src/md/opti/targeter.rs similarity index 96% rename from src/md/opti/optimizer.rs rename to src/md/opti/targeter.rs index dda79528..389e4edf 100644 --- a/src/md/opti/optimizer.rs +++ b/src/md/opti/targeter.rs @@ -30,11 +30,9 @@ use std::fmt; use super::solution::TargeterSolution; -// TODO(now): rename to Differential Controller - /// An optimizer structure with V control variables and O objectives. #[derive(Clone)] -pub struct Optimizer<'a, const V: usize, const O: usize> { +pub struct Targeter<'a, const V: usize, const O: usize> { /// The propagator setup (kind, stages, etc.) pub prop: &'a Propagator, /// The list of objectives of this targeter @@ -50,7 +48,7 @@ pub struct Optimizer<'a, const V: usize, const O: usize> { pub iterations: usize, } -impl<'a, const V: usize, const O: usize> fmt::Display for Optimizer<'a, V, O> { +impl<'a, const V: usize, const O: usize> fmt::Display for Targeter<'a, V, O> { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { let mut objmsg = String::from(""); for obj in &self.objectives { @@ -66,7 +64,7 @@ impl<'a, const V: usize, const O: usize> fmt::Display for Optimizer<'a, V, O> { } } -impl<'a, const O: usize> Optimizer<'a, 3, O> { +impl<'a, const O: usize> Targeter<'a, 3, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. pub fn delta_v(prop: &'a Propagator, objectives: [Objective; O]) -> Self { Self { @@ -116,7 +114,7 @@ impl<'a, const O: usize> Optimizer<'a, 3, O> { } } -impl<'a, const O: usize> Optimizer<'a, 4, O> { +impl<'a, const O: usize> Targeter<'a, 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, @@ -138,7 +136,7 @@ impl<'a, const O: usize> Optimizer<'a, 4, O> { } } -impl<'a, const O: usize> Optimizer<'a, 7, O> { +impl<'a, const O: usize> Targeter<'a, 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, @@ -163,7 +161,7 @@ impl<'a, const O: usize> Optimizer<'a, 7, O> { } } -impl<'a, const O: usize> Optimizer<'a, 10, O> { +impl<'a, const O: usize> Targeter<'a, 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, @@ -191,7 +189,7 @@ impl<'a, const O: usize> Optimizer<'a, 10, O> { } } -impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> { +impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> { /// Create a new Targeter which will apply an impulsive delta-v correction. pub fn new( prop: &'a Propagator, diff --git a/src/md/trajectory/interpolatable.rs b/src/md/trajectory/interpolatable.rs index 32e54641..b845cd5b 100644 --- a/src/md/trajectory/interpolatable.rs +++ b/src/md/trajectory/interpolatable.rs @@ -47,9 +47,6 @@ where /// List of state parameters that will be exported to a trajectory file in addition to the epoch (provided in this different formats). fn export_params() -> Vec; - - /// Returns the orbit - fn orbit(&self) -> &Orbit; } impl Interpolatable for Spacecraft { @@ -156,8 +153,4 @@ impl Interpolatable for Spacecraft { ] .concat() } - - fn orbit(&self) -> &Orbit { - &self.orbit - } } diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs index 58ba2f45..01a2c302 100644 --- a/src/md/trajectory/traj.rs +++ b/src/md/trajectory/traj.rs @@ -445,8 +445,8 @@ where // Build an array of all the RIC differences let mut ric_diff = Vec::with_capacity(other_states.len()); for (ii, other_state) in other_states.iter().enumerate() { - let self_orbit = *self_states[ii].orbit(); - let other_orbit = *other_state.orbit(); + let self_orbit = self_states[ii].orbit(); + let other_orbit = other_state.orbit(); let this_ric_diff = self_orbit.ric_difference(&other_orbit).map_err(Box::new)?; diff --git a/src/propagators/error_ctrl.rs b/src/propagators/error_ctrl.rs index 0fc2efde..d3c0b232 100644 --- a/src/propagators/error_ctrl.rs +++ b/src/propagators/error_ctrl.rs @@ -186,6 +186,12 @@ impl ErrorControl { } } +impl Default for ErrorControl { + fn default() -> Self { + Self::RSSCartesianStep + } +} + /// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3 /// /// Note that this error controller should be preferably be used only with slices of a state with the same units. diff --git a/src/propagators/instance.rs b/src/propagators/instance.rs index 7214a69d..26a770bf 100644 --- a/src/propagators/instance.rs +++ b/src/propagators/instance.rs @@ -17,7 +17,7 @@ */ use super::{DynamicsSnafu, IntegrationDetails, PropagationError, Propagator}; -use crate::dynamics::Dynamics; +use crate::dynamics::{Dynamics, DynamicsAlmanacSnafu}; use crate::linalg::allocator::Allocator; use crate::linalg::{DefaultAllocator, OVector}; use crate::md::trajectory::{Interpolatable, Traj}; @@ -94,11 +94,6 @@ where } let stop_time = self.state.epoch() + duration; - #[cfg(not(target_arch = "wasm32"))] - let tick = Instant::now(); - #[cfg(not(target_arch = "wasm32"))] - let mut prev_tick = Instant::now(); - if self.log_progress { // Prevent the print spam for orbit determination cases info!("Propagating for {} until {}", duration, stop_time); @@ -114,6 +109,39 @@ where if backprop { self.step_size = -self.step_size; // Invert the step size } + + // Transform the state if needed + let mut original_frame = None; + if let Some(integration_frame) = self.prop.opts.integration_frame { + if integration_frame != self.state.orbit().frame { + original_frame = Some(self.state.orbit().frame); + let mut new_orbit = self + .almanac + .transform_to(self.state.orbit(), integration_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state into desired integration frame", + }) + .context(DynamicsSnafu)?; + // If the integration frame has parameters, we set them here. + if let Some(mu_km3_s2) = integration_frame.mu_km3_s2 { + new_orbit.frame.mu_km3_s2 = Some(mu_km3_s2); + } + // If the integration frame has parameters, we set them here. + if let Some(shape) = integration_frame.shape { + new_orbit.frame.shape = Some(shape); + } + if self.log_progress { + info!("State transformed to the integration frame {integration_frame}"); + } + self.state.set_orbit(new_orbit); + } + } + + #[cfg(not(target_arch = "wasm32"))] + let tick = Instant::now(); + #[cfg(not(target_arch = "wasm32"))] + let mut prev_tick = Instant::now(); + loop { let epoch = self.state.epoch(); if (!backprop && epoch + self.step_size > stop_time) @@ -128,6 +156,19 @@ where debug!("Done in {}", tock); } } + + // Rotate back if needed + if let Some(original_frame) = original_frame { + let new_orbit = self + .almanac + .transform_to(self.state.orbit(), original_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state from desired integration frame", + }) + .context(DynamicsSnafu)?; + self.state.set_orbit(new_orbit); + } + return Ok(self.state); } // Take one final step of exactly the needed duration until the stop time @@ -159,6 +200,18 @@ where } } + // Rotate back if needed + if let Some(original_frame) = original_frame { + let new_orbit = self + .almanac + .transform_to(self.state.orbit(), original_frame, None) + .context(DynamicsAlmanacSnafu { + action: "transforming state from desired integration frame", + }) + .context(DynamicsSnafu)?; + self.state.set_orbit(new_orbit); + } + return Ok(self.state); } else { #[cfg(not(target_arch = "wasm32"))] diff --git a/src/propagators/options.rs b/src/propagators/options.rs index 2d8e59c8..8ec73d49 100644 --- a/src/propagators/options.rs +++ b/src/propagators/options.rs @@ -21,16 +21,17 @@ use std::fmt; use crate::time::{Duration, Unit}; use super::ErrorControl; +use anise::frames::Frame; +use serde::{Deserialize, Serialize}; use typed_builder::TypedBuilder; -/// PropOpts stores the integrator options, including the minimum and maximum step sizes, and the -/// max error size. +/// Stores the integrator options, including the minimum and maximum step sizes, and the central body to perform the integration. /// /// Note that different step sizes and max errors are only used for adaptive /// 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, TypedBuilder)] +#[derive(Clone, Copy, Debug, TypedBuilder, Serialize, Deserialize, PartialEq)] #[builder(doc)] pub struct IntegratorOptions { #[builder(default_code = "60.0 * Unit::Second")] @@ -45,8 +46,12 @@ pub struct IntegratorOptions { pub attempts: u8, #[builder(default = false)] pub fixed_step: bool, + #[builder(default)] pub error_ctrl: ErrorControl, - // TODO(now): Add an optional central body that will swap central bodies on init + /// If a frame is specified and the propagator state is in a different frame, it it changed to this frame prior to integration. + /// Note, when setting this, it's recommended to call `strip` on the Frame. + #[builder(default, setter(strip_option))] + pub integration_frame: Option, } impl IntegratorOptions { @@ -66,6 +71,7 @@ impl IntegratorOptions { attempts: 50, fixed_step: false, error_ctrl, + integration_frame: None, } } @@ -94,6 +100,7 @@ impl IntegratorOptions { fixed_step: true, attempts: 0, error_ctrl: ErrorControl::RSSCartesianStep, + integration_frame: None, } } @@ -164,37 +171,55 @@ impl Default for IntegratorOptions { attempts: 50, fixed_step: false, error_ctrl: ErrorControl::RSSCartesianStep, + integration_frame: None, } } } -#[test] -fn test_options() { - let opts = IntegratorOptions::with_fixed_step_s(1e-1); - assert_eq!(opts.min_step, 1e-1 * Unit::Second); - assert_eq!(opts.max_step, 1e-1 * Unit::Second); - assert!(opts.tolerance.abs() < f64::EPSILON); - assert!(opts.fixed_step); - - let opts = IntegratorOptions::with_adaptive_step_s(1e-2, 10.0, 1e-12, ErrorControl::RSSStep); - assert_eq!(opts.min_step, 1e-2 * Unit::Second); - assert_eq!(opts.max_step, 10.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert!(!opts.fixed_step); - - let opts: IntegratorOptions = Default::default(); - assert_eq!(opts.init_step, 60.0 * Unit::Second); - assert_eq!(opts.min_step, 0.001 * Unit::Second); - assert_eq!(opts.max_step, 2700.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert_eq!(opts.attempts, 50); - assert!(!opts.fixed_step); - - let opts = IntegratorOptions::with_max_step(1.0 * Unit::Second); - assert_eq!(opts.init_step, 1.0 * Unit::Second); - assert_eq!(opts.min_step, 0.001 * Unit::Second); - assert_eq!(opts.max_step, 1.0 * Unit::Second); - assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); - assert_eq!(opts.attempts, 50); - assert!(!opts.fixed_step); +#[cfg(test)] +mod ut_integr_opts { + use hifitime::Unit; + + use crate::propagators::{ErrorControl, IntegratorOptions}; + + #[test] + fn test_options() { + let opts = IntegratorOptions::with_fixed_step_s(1e-1); + assert_eq!(opts.min_step, 1e-1 * Unit::Second); + assert_eq!(opts.max_step, 1e-1 * Unit::Second); + assert!(opts.tolerance.abs() < f64::EPSILON); + assert!(opts.fixed_step); + + let opts = + IntegratorOptions::with_adaptive_step_s(1e-2, 10.0, 1e-12, ErrorControl::RSSStep); + assert_eq!(opts.min_step, 1e-2 * Unit::Second); + assert_eq!(opts.max_step, 10.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert!(!opts.fixed_step); + + let opts: IntegratorOptions = Default::default(); + assert_eq!(opts.init_step, 60.0 * Unit::Second); + assert_eq!(opts.min_step, 0.001 * Unit::Second); + assert_eq!(opts.max_step, 2700.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert_eq!(opts.attempts, 50); + assert!(!opts.fixed_step); + + let opts = IntegratorOptions::with_max_step(1.0 * Unit::Second); + assert_eq!(opts.init_step, 1.0 * Unit::Second); + assert_eq!(opts.min_step, 0.001 * Unit::Second); + assert_eq!(opts.max_step, 1.0 * Unit::Second); + assert!((opts.tolerance - 1e-12).abs() < f64::EPSILON); + assert_eq!(opts.attempts, 50); + assert!(!opts.fixed_step); + } + + #[test] + fn test_serde() { + let opts = IntegratorOptions::default(); + let serialized = toml::to_string(&opts).unwrap(); + println!("{serialized}"); + let deserd: IntegratorOptions = toml::from_str(&serialized).unwrap(); + assert_eq!(deserd, opts); + } } diff --git a/src/propagators/rk_methods/dormand.rs b/src/propagators/rk_methods/dormand.rs index b124c96d..31948135 100644 --- a/src/propagators/rk_methods/dormand.rs +++ b/src/propagators/rk_methods/dormand.rs @@ -19,7 +19,7 @@ use super::RK; /// `Dormand45` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). -pub struct Dormand45 {} +pub(crate) struct Dormand45 {} impl RK for Dormand45 { const ORDER: u8 = 5; @@ -68,7 +68,7 @@ impl RK for Dormand45 { /// `Dormand78` is a [Dormand-Prince integrator](https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method). /// /// Coefficients taken from GMAT `src/base/propagator/PrinceDormand78.cpp`. -pub struct Dormand78 {} +pub(crate) struct Dormand78 {} impl RK for Dormand78 { const ORDER: u8 = 8; diff --git a/src/propagators/rk_methods/mod.rs b/src/propagators/rk_methods/mod.rs index dea49f02..f6500385 100644 --- a/src/propagators/rk_methods/mod.rs +++ b/src/propagators/rk_methods/mod.rs @@ -24,11 +24,11 @@ use serde::Serialize; use crate::io::ConfigError; -pub use self::rk::*; +use self::rk::*; mod dormand; -pub use self::dormand::*; +use self::dormand::*; mod verner; -pub use self::verner::*; +use self::verner::*; use super::PropagationError; diff --git a/src/propagators/rk_methods/rk.rs b/src/propagators/rk_methods/rk.rs index fe78c203..282b64cc 100644 --- a/src/propagators/rk_methods/rk.rs +++ b/src/propagators/rk_methods/rk.rs @@ -19,7 +19,7 @@ use super::RK; /// `CashKarp45` is a [Runge Kutta Cash Karp integrator](https://en.wikipedia.org/wiki/Cash%E2%80%93Karp_method). -pub struct CashKarp45 {} +pub(crate) struct CashKarp45 {} impl RK for CashKarp45 { const ORDER: u8 = 5; @@ -61,7 +61,7 @@ impl RK for CashKarp45 { /// /// If initialized with an `PropOpts.with_adaptive_step`, the variable step will **not** be taken into consideration. #[allow(clippy::upper_case_acronyms)] -pub struct RK4Fixed {} +pub(crate) struct RK4Fixed {} impl RK for RK4Fixed { const ORDER: u8 = 4; @@ -86,7 +86,7 @@ const SQRT6: f64 = 2.449_489_742_783_178; /// /// Coefficients taken from GMAT `src/base/propagator/RungeKutta89.cpp`. #[allow(clippy::upper_case_acronyms)] -pub struct RK89 {} +pub(crate) struct RK89 {} impl RK for RK89 { const ORDER: u8 = 9; diff --git a/src/propagators/rk_methods/verner.rs b/src/propagators/rk_methods/verner.rs index 9318d93a..0adb5dbe 100644 --- a/src/propagators/rk_methods/verner.rs +++ b/src/propagators/rk_methods/verner.rs @@ -21,7 +21,7 @@ use super::RK; /// `Verner56` is an RK Verner integrator of order 5-6. /// /// Coefficients taken from [here (PDF)](http://people.math.sfu.ca/~jverner/classify.1992.ps). -pub struct Verner56 {} +pub(crate) struct Verner56 {} impl RK for Verner56 { const ORDER: u8 = 6; diff --git a/tests/mission_design/force_models.rs b/tests/mission_design/force_models.rs index 41fcf65c..4139c793 100644 --- a/tests/mission_design/force_models.rs +++ b/tests/mission_design/force_models.rs @@ -1,5 +1,6 @@ extern crate nyx_space as nyx; +use anise::constants::frames::MOON_J2000; use nyx::cosmic::{Orbit, Spacecraft}; use nyx::dynamics::{Drag, OrbitalDynamics, SolarPressure, SpacecraftDynamics}; use nyx::linalg::Vector6; @@ -8,6 +9,7 @@ use nyx::time::{Epoch, Unit}; use nyx::utils::rss_orbit_vec_errors; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; +use nyx_space::md::prelude::Interpolatable; use rstest::*; use std::sync::Arc; @@ -150,7 +152,7 @@ fn srp_earth_meo_ecc_inc(almanac: Arc) { let sc = Spacecraft::from_srp_defaults(orbit, dry_mass, 16.0); - let setup = Propagator::default(sc_dyn); + let setup = Propagator::default(sc_dyn.clone()); let mut prop = setup.with(sc, almanac.clone()); let final_state = prop.for_duration(prop_time).unwrap(); @@ -176,6 +178,39 @@ fn srp_earth_meo_ecc_inc(almanac: Arc) { assert!(err_r < 2e-3, "position error too large for SRP"); assert!(err_v < 1e-6, "velocity error too large for SRP"); + // Test that we can specify an integration frame in the options and that the result is the same. + // Note that we also test here that we're setting the GM and shape of the integration frame as defined + // and not just grabbing that data from the almanac. + let orbit = almanac.transform_to(orbit, MOON_J2000, None).unwrap(); + println!("{:x}", orbit); + + let sc_moon = Spacecraft::from_srp_defaults(orbit, dry_mass, 16.0); + + let mut setup_moon = Propagator::default(sc_dyn); + setup_moon.opts.integration_frame = Some(eme2k); + let mut prop = setup_moon.with(sc_moon, almanac.clone()); + let final_moon_state = prop.for_duration(prop_time).unwrap(); + assert!( + final_moon_state.frame().ephem_origin_match(MOON_J2000), + "expected a result in the Moon frame" + ); + println!("{}", final_moon_state); + + let final_earth_orbit = almanac + .transform_to(final_moon_state.orbit, EARTH_J2000, None) + .unwrap(); + + let (fw_err_r, fw_err_v) = + rss_orbit_vec_errors(&final_earth_orbit.to_cartesian_pos_vel(), &rslt); + println!( + "Error accumulated in ecc+inc MEO (with penumbras) over {} after integration frame swap: {:.6} m \t{:.6} m/s", + prop_time, + fw_err_r * 1e3, + fw_err_v * 1e3 + ); + assert!(dbg!((fw_err_r - err_r).abs() / err_r) < 1e-9); + assert!(dbg!((fw_err_v - err_v).abs() / err_v) < 1e-12); + // Compare the case with the hyperdual EOMs (computation uses another part of the code) let mut prop = setup.with(sc.with_stm(), almanac); let final_state_dual = prop.for_duration(prop_time).unwrap(); diff --git a/tests/mission_design/targeter/b_plane.rs b/tests/mission_design/targeter/b_plane.rs index c1889c02..3f65ddf3 100644 --- a/tests/mission_design/targeter/b_plane.rs +++ b/tests/mission_design/targeter/b_plane.rs @@ -4,7 +4,7 @@ use anise::constants::celestial_objects::JUPITER_BARYCENTER; use anise::constants::celestial_objects::MOON; use anise::constants::celestial_objects::SUN; use anise::constants::frames::MOON_J2000; -use nyx::md::optimizer::*; +use nyx::md::targeter::*; use nyx::md::prelude::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -48,7 +48,7 @@ fn tgt_b_plane_earth_gravity_assist_no_propagation(almanac: Arc) { let b_plane_tgt = BPlaneTarget::from_bt_br(13135.7982982557, 5022.26511510685); - let tgt = Optimizer::delta_v(&prop, b_plane_tgt.to_objectives()); + let tgt = Targeter::delta_v(&prop, b_plane_tgt.to_objectives()); let sol = tgt .try_achieve_from(spacecraft, epoch, epoch, almanac.clone()) @@ -117,7 +117,7 @@ fn tgt_b_plane_lunar_transfer(almanac: Arc) { // GMAT truth with central differencing: 1.15740867962, -0.576350387399, 0.632247251449 // GMAT truth with forward differencing: 1.33490412071, -0.5447988683, 1.77697094604 (approach currently in Nyx) - let tgt = Optimizer::in_frame( + let tgt = Targeter::in_frame( &prop, [ Variable { @@ -209,7 +209,7 @@ fn tgt_b_plane_earth_gravity_assist_with_propagation(almanac: Arc) { let b_plane_tgt = BPlaneTarget::from_bt_br(13135.7982982557, 5022.26511510685); - let tgt = Optimizer::delta_v(&prop, b_plane_tgt.to_objectives()); + let tgt = Targeter::delta_v(&prop, b_plane_tgt.to_objectives()); let sol = tgt .try_achieve_from(prior_sc, prior_sc.epoch(), epoch, almanac.clone()) diff --git a/tests/mission_design/targeter/finite_burns.rs b/tests/mission_design/targeter/finite_burns.rs index 104099fa..e4b92d4d 100644 --- a/tests/mission_design/targeter/finite_burns.rs +++ b/tests/mission_design/targeter/finite_burns.rs @@ -4,7 +4,7 @@ use anise::constants::celestial_objects::{JUPITER_BARYCENTER, MOON, SUN}; use hifitime::TimeUnits; use nyx::dynamics::guidance::{LocalFrame, Mnvr, Thruster}; use nyx::linalg::Vector3; -use nyx::md::optimizer::*; +use nyx::md::targeter::*; use nyx::md::prelude::*; use crate::propagation::GMAT_EARTH_GM; @@ -52,7 +52,7 @@ fn thrust_dir_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_dir(&setup, objectives); + let tgt = Targeter::thrust_dir(&setup, objectives); println!("{}", tgt); @@ -99,7 +99,7 @@ fn thrust_dir_rate_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_dir_rate(&setup, objectives); + let tgt = Targeter::thrust_dir_rate(&setup, objectives); println!("{}", tgt); @@ -148,7 +148,7 @@ fn thrust_profile_tgt_sma_aop_raan(almanac: Arc) { Objective::within_tolerance(StateParameter::RAAN, 60.000182, 1e-3), ]; - let tgt = Optimizer::thrust_profile(&setup, objectives); + let tgt = Targeter::thrust_profile(&setup, objectives); println!("{}", tgt); @@ -222,7 +222,7 @@ fn val_tgt_finite_burn(almanac: Arc) { let sc_no_thrust = SpacecraftDynamics::new(orbital_dyn); let mut prop_no_thrust = Propagator::default(sc_no_thrust); prop_no_thrust.set_max_step(mnvr0.duration()); - let impulsive_tgt = Optimizer::delta_v( + let impulsive_tgt = Targeter::delta_v( &prop_no_thrust, [ Objective::within_tolerance(StateParameter::X, sc_xf_desired.orbit.radius_km.x, 1e-5), diff --git a/tests/mission_design/targeter/multi_oe.rs b/tests/mission_design/targeter/multi_oe.rs index 77caa13c..cb46838a 100644 --- a/tests/mission_design/targeter/multi_oe.rs +++ b/tests/mission_design/targeter/multi_oe.rs @@ -2,7 +2,7 @@ extern crate nyx_space as nyx; // use nyx::dynamics::guidance::Mnvr; use nyx::dynamics::guidance::Thruster; -use nyx::md::optimizer::*; +use nyx::md::targeter::*; use nyx::md::prelude::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -40,7 +40,7 @@ fn tgt_c3_decl(almanac: Arc) { Objective::within_tolerance(StateParameter::C3, -5.0, 0.5), ]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -106,7 +106,7 @@ fn conv_tgt_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -209,7 +209,7 @@ fn tgt_hd_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { diff --git a/tests/mission_design/targeter/multi_oe_vnc.rs b/tests/mission_design/targeter/multi_oe_vnc.rs index d265320a..75a44e65 100644 --- a/tests/mission_design/targeter/multi_oe_vnc.rs +++ b/tests/mission_design/targeter/multi_oe_vnc.rs @@ -1,6 +1,6 @@ extern crate nyx_space as nyx; -use nyx::md::optimizer::*; +use nyx::md::targeter::*; use nyx::md::prelude::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -38,7 +38,7 @@ fn tgt_vnc_c3_decl(almanac: Arc) { Objective::within_tolerance(StateParameter::C3, -5.0, 0.5), ]; - let tgt = Optimizer::vnc(&setup, objectives); + let tgt = Targeter::vnc(&setup, objectives); println!("{}", tgt); @@ -86,7 +86,7 @@ fn tgt_vnc_sma_ecc(almanac: Arc) { Objective::within_tolerance(StateParameter::SMA, 8100.0, 0.1), ]; - let tgt = Optimizer::vnc_with_components( + let tgt = Targeter::vnc_with_components( &setup, [ Variable { diff --git a/tests/mission_design/targeter/single_oe.rs b/tests/mission_design/targeter/single_oe.rs index b7736655..ac6cfebe 100644 --- a/tests/mission_design/targeter/single_oe.rs +++ b/tests/mission_design/targeter/single_oe.rs @@ -3,7 +3,7 @@ extern crate nyx_space as nyx; use anise::constants::celestial_objects::JUPITER_BARYCENTER; use anise::constants::celestial_objects::MOON; use anise::constants::celestial_objects::SUN; -use nyx::md::optimizer::*; +use nyx::md::targeter::*; use nyx::md::prelude::*; use anise::{constants::frames::EARTH_J2000, prelude::Almanac}; @@ -43,7 +43,7 @@ fn tgt_sma_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -103,7 +103,7 @@ fn tgt_sma_from_peri_fd(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -161,7 +161,7 @@ fn tgt_hd_sma_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::SMA, xf_desired_sma)]; - let mut tgt = Optimizer::delta_v(&setup, objectives); + let mut tgt = Targeter::delta_v(&setup, objectives); tgt.iterations = 5; println!("{}", tgt); @@ -268,7 +268,7 @@ fn tgt_ecc_from_apo(almanac: Arc) { let xf_desired_ecc = 0.4; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -342,7 +342,7 @@ fn tgt_ecc_from_peri(almanac: Arc) { let xf_desired_ecc = 0.4; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -416,7 +416,7 @@ fn tgt_raan_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::RAAN, xf_desired_raan)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -469,7 +469,7 @@ fn tgt_raan_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::RAAN, xf_desired_raan)]; - let tgt = Optimizer::new( + let tgt = Targeter::new( &setup, [ Variable { @@ -543,7 +543,7 @@ fn tgt_aop_from_apo(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::AoP, xf_desired_aop)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt); @@ -596,7 +596,7 @@ fn tgt_aop_from_peri(almanac: Arc) { // Define the objective let objectives = [Objective::new(StateParameter::AoP, xf_desired_aop)]; - let tgt = Optimizer::delta_v(&setup, objectives); + let tgt = Targeter::delta_v(&setup, objectives); println!("{}", tgt);