Skip to content

Commit ef56369

Browse files
ChristopherRabotinChristopherRabotin
committed
Add integration frame and serialization to integrator options
1 parent 450e821 commit ef56369

25 files changed

+236
-110
lines changed

‎Cargo.toml

+2-1
Original file line numberDiff line numberDiff line change
@@ -77,12 +77,13 @@ typed-builder = "0.20.0"
7777
pythonize = { version = "0.21", optional = true }
7878
snafu = { version = "0.8.3", features = ["backtrace"] }
7979
serde_dhall = "0.12"
80-
toml = "0.8.14"
80+
8181

8282
[dev-dependencies]
8383
polars = { version = "0.42.0", features = ["parquet"] }
8484
rstest = "0.22.0"
8585
pretty_env_logger = "0.5"
86+
toml = "0.8.14"
8687

8788
[build-dependencies]
8889
shadow-rs = "0.35.0"

‎src/cosmic/mod.rs

+6
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,12 @@ where
110110
fn set_value(&mut self, param: StateParameter, _val: f64) -> Result<(), StateError> {
111111
Err(StateError::Unavailable { param })
112112
}
113+
114+
/// Returns a copy of the orbit
115+
fn orbit(&self) -> Orbit;
116+
117+
/// Modifies this state's orbit
118+
fn set_orbit(&mut self, _orbit: Orbit) {}
113119
}
114120

115121
pub fn assert_orbit_eq_or_abs(left: &Orbit, right: &Orbit, epsilon: f64, msg: &str) {

‎src/cosmic/spacecraft.rs

+8
Original file line numberDiff line numberDiff line change
@@ -761,6 +761,14 @@ impl State for Spacecraft {
761761
fn unset_stm(&mut self) {
762762
self.stm = None;
763763
}
764+
765+
fn orbit(&self) -> Orbit {
766+
self.orbit
767+
}
768+
769+
fn set_orbit(&mut self, orbit: Orbit) {
770+
self.orbit = orbit;
771+
}
764772
}
765773

766774
impl Add<OVector<f64, Const<6>>> for Spacecraft {

‎src/dynamics/mod.rs

+1
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,7 @@ pub trait AccelModel: Send + Sync + fmt::Display {
172172

173173
/// Stores dynamical model errors
174174
#[derive(Debug, PartialEq, Snafu)]
175+
#[snafu(visibility(pub(crate)))]
175176
pub enum DynamicsError {
176177
/// Fuel exhausted at the provided spacecraft state
177178
#[snafu(display("fuel exhausted at {sc}"))]

‎src/md/mod.rs

+2-2
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ use snafu::prelude::*;
2525

2626
pub mod prelude {
2727
pub use super::{
28-
optimizer::*,
28+
targeter::*,
2929
trajectory::{ExportCfg, Interpolatable, Traj},
3030
Event, ScTraj, StateParameter,
3131
};
@@ -52,7 +52,7 @@ pub use events::{Event, EventEvaluator};
5252

5353
pub mod objective;
5454
pub mod opti;
55-
pub use opti::optimizer;
55+
pub use opti::targeter;
5656
pub type ScTraj = trajectory::Traj<Spacecraft>;
5757
// pub type Ephemeris = trajectory::Traj<Orbit>;
5858

‎src/md/opti/mod.rs

+4-4
Original file line numberDiff line numberDiff line change
@@ -19,16 +19,16 @@
1919
// pub mod convert_impulsive;
2020
pub mod multipleshooting;
2121
pub use multipleshooting::{ctrlnodes, multishoot};
22-
/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem.
23-
// #[cfg(feature = "broken-donotuse")]
24-
// pub mod minimize_lm;
25-
pub mod optimizer;
2622
/// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via finite differencing.
2723
pub mod raphson_finite_diff;
2824
/// Uses a [Newton Raphson](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) method where the Jacobian is computed via hyperdual numbers.
2925
pub mod raphson_hyperdual;
3026
pub mod solution;
3127
pub mod target_variable;
28+
/// Uses a Levenberg Marquardt minimizer to solve the damped least squares problem.
29+
// #[cfg(feature = "broken-donotuse")]
30+
// pub mod minimize_lm;
31+
pub mod targeter;
3232

3333
#[derive(Copy, Clone, Debug, PartialEq, Eq)]
3434
pub enum DiffMethod {

‎src/md/opti/multipleshooting/multishoot.rs

+7-7
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ pub use super::CostFunction;
2222
use super::{MultipleShootingError, TargetingSnafu};
2323
use crate::linalg::{DMatrix, DVector, SVector};
2424
use crate::md::opti::solution::TargeterSolution;
25-
use crate::md::optimizer::Optimizer;
25+
use crate::md::targeter::Targeter;
2626
use crate::md::{prelude::*, TargetingError};
2727
use crate::pseudo_inverse;
2828
use crate::{Orbit, Spacecraft};
@@ -81,7 +81,7 @@ impl<'a, T: MultishootNode<OT>, const VT: usize, const OT: usize> MultipleShooti
8181
/* ***
8282
** 1. Solve the delta-v differential corrector between each node
8383
** *** */
84-
let tgt = Optimizer {
84+
let tgt = Targeter {
8585
prop: self.prop,
8686
objectives: self.targets[i].into(),
8787
variables: self.variables,
@@ -125,7 +125,7 @@ impl<'a, T: MultishootNode<OT>, const VT: usize, const OT: usize> MultipleShooti
125125
** Note that because the first initial_state is x0, the i-th "initial state"
126126
** is the initial state to reach the i-th node.
127127
** *** */
128-
let inner_tgt_a = Optimizer::delta_v(self.prop, next_node);
128+
let inner_tgt_a = Targeter::delta_v(self.prop, next_node);
129129
let inner_sol_a = inner_tgt_a
130130
.try_achieve_dual(
131131
initial_states[i],
@@ -151,7 +151,7 @@ impl<'a, T: MultishootNode<OT>, const VT: usize, const OT: usize> MultipleShooti
151151
/* ***
152152
** 2.C. Rerun the targeter from the new state at the perturbed node to the next unpertubed node
153153
** *** */
154-
let inner_tgt_b = Optimizer::delta_v(self.prop, self.targets[i + 1].into());
154+
let inner_tgt_b = Targeter::delta_v(self.prop, self.targets[i + 1].into());
155155
let inner_sol_b = inner_tgt_b
156156
.try_achieve_dual(
157157
inner_sol_a.achieved_state,
@@ -241,7 +241,7 @@ impl<'a, T: MultishootNode<OT>, const VT: usize, const OT: usize> MultipleShooti
241241

242242
for (i, node) in self.targets.iter().enumerate() {
243243
// Run the unpertubed targeter
244-
let tgt = Optimizer::delta_v(self.prop, (*node).into());
244+
let tgt = Targeter::delta_v(self.prop, (*node).into());
245245
let sol = tgt
246246
.try_achieve_dual(
247247
initial_states[i],
@@ -355,8 +355,8 @@ impl<T: MultishootNode<O>, const O: usize> MultipleShootingSolution<T, O> {
355355
) -> Result<Vec<ScTraj>, MultipleShootingError> {
356356
let mut trajz = Vec::with_capacity(self.nodes.len());
357357

358-
for (i, node) in self.nodes.iter().enumerate() {
359-
let (_, traj) = Optimizer::delta_v(prop, (*node).into())
358+
for (i, node) in self.nodes.iter().copied().enumerate() {
359+
let (_, traj) = Targeter::delta_v(prop, node.into())
360360
.apply_with_traj(&self.solutions[i], almanac.clone())
361361
.context(TargetingSnafu { segment: i })?;
362362
trajz.push(traj);

‎src/md/opti/raphson_finite_diff.rs

+2-2
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
along with this program. If not, see <https://www.gnu.org/licenses/>.
1717
*/
1818

19-
use super::optimizer::Optimizer;
19+
use super::targeter::Targeter;
2020
use super::solution::TargeterSolution;
2121
use crate::cosmic::{AstroAlmanacSnafu, AstroPhysicsSnafu};
2222
use crate::dynamics::guidance::{GuidanceError, LocalFrame, Mnvr};
@@ -33,7 +33,7 @@ use snafu::{ensure, ResultExt};
3333
#[cfg(not(target_arch = "wasm32"))]
3434
use std::time::Instant;
3535

36-
impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> {
36+
impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> {
3737
/// Differential correction using finite differencing
3838
#[allow(clippy::comparison_chain)]
3939
pub fn try_achieve_fd(

‎src/md/opti/raphson_hyperdual.rs

+1-1
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,7 @@ use crate::utils::are_eigenvalues_stable;
3030
#[cfg(not(target_arch = "wasm32"))]
3131
use std::time::Instant;
3232

33-
impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> {
33+
impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> {
3434
/// Differential correction using hyperdual numbers for the objectives
3535
#[allow(clippy::comparison_chain)]
3636
pub fn try_achieve_dual(

‎src/md/opti/optimizer.rs ‎src/md/opti/targeter.rs

+7-9
Original file line numberDiff line numberDiff line change
@@ -30,11 +30,9 @@ use std::fmt;
3030

3131
use super::solution::TargeterSolution;
3232

33-
// TODO(now): rename to Differential Controller
34-
3533
/// An optimizer structure with V control variables and O objectives.
3634
#[derive(Clone)]
37-
pub struct Optimizer<'a, const V: usize, const O: usize> {
35+
pub struct Targeter<'a, const V: usize, const O: usize> {
3836
/// The propagator setup (kind, stages, etc.)
3937
pub prop: &'a Propagator<SpacecraftDynamics>,
4038
/// The list of objectives of this targeter
@@ -50,7 +48,7 @@ pub struct Optimizer<'a, const V: usize, const O: usize> {
5048
pub iterations: usize,
5149
}
5250

53-
impl<'a, const V: usize, const O: usize> fmt::Display for Optimizer<'a, V, O> {
51+
impl<'a, const V: usize, const O: usize> fmt::Display for Targeter<'a, V, O> {
5452
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
5553
let mut objmsg = String::from("");
5654
for obj in &self.objectives {
@@ -66,7 +64,7 @@ impl<'a, const V: usize, const O: usize> fmt::Display for Optimizer<'a, V, O> {
6664
}
6765
}
6866

69-
impl<'a, const O: usize> Optimizer<'a, 3, O> {
67+
impl<'a, const O: usize> Targeter<'a, 3, O> {
7068
/// Create a new Targeter which will apply an impulsive delta-v correction.
7169
pub fn delta_v(prop: &'a Propagator<SpacecraftDynamics>, objectives: [Objective; O]) -> Self {
7270
Self {
@@ -116,7 +114,7 @@ impl<'a, const O: usize> Optimizer<'a, 3, O> {
116114
}
117115
}
118116

119-
impl<'a, const O: usize> Optimizer<'a, 4, O> {
117+
impl<'a, const O: usize> Targeter<'a, 4, O> {
120118
/// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
121119
pub fn thrust_dir(
122120
prop: &'a Propagator<SpacecraftDynamics>,
@@ -138,7 +136,7 @@ impl<'a, const O: usize> Optimizer<'a, 4, O> {
138136
}
139137
}
140138

141-
impl<'a, const O: usize> Optimizer<'a, 7, O> {
139+
impl<'a, const O: usize> Targeter<'a, 7, O> {
142140
/// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
143141
pub fn thrust_dir_rate(
144142
prop: &'a Propagator<SpacecraftDynamics>,
@@ -163,7 +161,7 @@ impl<'a, const O: usize> Optimizer<'a, 7, O> {
163161
}
164162
}
165163

166-
impl<'a, const O: usize> Optimizer<'a, 10, O> {
164+
impl<'a, const O: usize> Targeter<'a, 10, O> {
167165
/// Create a new Targeter which will apply a continuous thrust for the whole duration of the segment
168166
pub fn thrust_profile(
169167
prop: &'a Propagator<SpacecraftDynamics>,
@@ -191,7 +189,7 @@ impl<'a, const O: usize> Optimizer<'a, 10, O> {
191189
}
192190
}
193191

194-
impl<'a, const V: usize, const O: usize> Optimizer<'a, V, O> {
192+
impl<'a, const V: usize, const O: usize> Targeter<'a, V, O> {
195193
/// Create a new Targeter which will apply an impulsive delta-v correction.
196194
pub fn new(
197195
prop: &'a Propagator<SpacecraftDynamics>,

‎src/md/trajectory/interpolatable.rs

-7
Original file line numberDiff line numberDiff line change
@@ -47,9 +47,6 @@ where
4747

4848
/// List of state parameters that will be exported to a trajectory file in addition to the epoch (provided in this different formats).
4949
fn export_params() -> Vec<StateParameter>;
50-
51-
/// Returns the orbit
52-
fn orbit(&self) -> &Orbit;
5350
}
5451

5552
impl Interpolatable for Spacecraft {
@@ -156,8 +153,4 @@ impl Interpolatable for Spacecraft {
156153
]
157154
.concat()
158155
}
159-
160-
fn orbit(&self) -> &Orbit {
161-
&self.orbit
162-
}
163156
}

‎src/md/trajectory/traj.rs

+2-2
Original file line numberDiff line numberDiff line change
@@ -445,8 +445,8 @@ where
445445
// Build an array of all the RIC differences
446446
let mut ric_diff = Vec::with_capacity(other_states.len());
447447
for (ii, other_state) in other_states.iter().enumerate() {
448-
let self_orbit = *self_states[ii].orbit();
449-
let other_orbit = *other_state.orbit();
448+
let self_orbit = self_states[ii].orbit();
449+
let other_orbit = other_state.orbit();
450450

451451
let this_ric_diff = self_orbit.ric_difference(&other_orbit).map_err(Box::new)?;
452452

‎src/propagators/error_ctrl.rs

+6
Original file line numberDiff line numberDiff line change
@@ -186,6 +186,12 @@ impl ErrorControl {
186186
}
187187
}
188188

189+
impl Default for ErrorControl {
190+
fn default() -> Self {
191+
Self::RSSCartesianStep
192+
}
193+
}
194+
189195
/// An RSS step error control which effectively computes the L2 norm of the provided Vector of size 3
190196
///
191197
/// Note that this error controller should be preferably be used only with slices of a state with the same units.

‎src/propagators/instance.rs

+59-6
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@
1717
*/
1818

1919
use super::{DynamicsSnafu, IntegrationDetails, PropagationError, Propagator};
20-
use crate::dynamics::Dynamics;
20+
use crate::dynamics::{Dynamics, DynamicsAlmanacSnafu};
2121
use crate::linalg::allocator::Allocator;
2222
use crate::linalg::{DefaultAllocator, OVector};
2323
use crate::md::trajectory::{Interpolatable, Traj};
@@ -94,11 +94,6 @@ where
9494
}
9595
let stop_time = self.state.epoch() + duration;
9696

97-
#[cfg(not(target_arch = "wasm32"))]
98-
let tick = Instant::now();
99-
#[cfg(not(target_arch = "wasm32"))]
100-
let mut prev_tick = Instant::now();
101-
10297
if self.log_progress {
10398
// Prevent the print spam for orbit determination cases
10499
info!("Propagating for {} until {}", duration, stop_time);
@@ -114,6 +109,39 @@ where
114109
if backprop {
115110
self.step_size = -self.step_size; // Invert the step size
116111
}
112+
113+
// Transform the state if needed
114+
let mut original_frame = None;
115+
if let Some(integration_frame) = self.prop.opts.integration_frame {
116+
if integration_frame != self.state.orbit().frame {
117+
original_frame = Some(self.state.orbit().frame);
118+
let mut new_orbit = self
119+
.almanac
120+
.transform_to(self.state.orbit(), integration_frame, None)
121+
.context(DynamicsAlmanacSnafu {
122+
action: "transforming state into desired integration frame",
123+
})
124+
.context(DynamicsSnafu)?;
125+
// If the integration frame has parameters, we set them here.
126+
if let Some(mu_km3_s2) = integration_frame.mu_km3_s2 {
127+
new_orbit.frame.mu_km3_s2 = Some(mu_km3_s2);
128+
}
129+
// If the integration frame has parameters, we set them here.
130+
if let Some(shape) = integration_frame.shape {
131+
new_orbit.frame.shape = Some(shape);
132+
}
133+
if self.log_progress {
134+
info!("State transformed to the integration frame {integration_frame}");
135+
}
136+
self.state.set_orbit(new_orbit);
137+
}
138+
}
139+
140+
#[cfg(not(target_arch = "wasm32"))]
141+
let tick = Instant::now();
142+
#[cfg(not(target_arch = "wasm32"))]
143+
let mut prev_tick = Instant::now();
144+
117145
loop {
118146
let epoch = self.state.epoch();
119147
if (!backprop && epoch + self.step_size > stop_time)
@@ -128,6 +156,19 @@ where
128156
debug!("Done in {}", tock);
129157
}
130158
}
159+
160+
// Rotate back if needed
161+
if let Some(original_frame) = original_frame {
162+
let new_orbit = self
163+
.almanac
164+
.transform_to(self.state.orbit(), original_frame, None)
165+
.context(DynamicsAlmanacSnafu {
166+
action: "transforming state from desired integration frame",
167+
})
168+
.context(DynamicsSnafu)?;
169+
self.state.set_orbit(new_orbit);
170+
}
171+
131172
return Ok(self.state);
132173
}
133174
// Take one final step of exactly the needed duration until the stop time
@@ -159,6 +200,18 @@ where
159200
}
160201
}
161202

203+
// Rotate back if needed
204+
if let Some(original_frame) = original_frame {
205+
let new_orbit = self
206+
.almanac
207+
.transform_to(self.state.orbit(), original_frame, None)
208+
.context(DynamicsAlmanacSnafu {
209+
action: "transforming state from desired integration frame",
210+
})
211+
.context(DynamicsSnafu)?;
212+
self.state.set_orbit(new_orbit);
213+
}
214+
162215
return Ok(self.state);
163216
} else {
164217
#[cfg(not(target_arch = "wasm32"))]

0 commit comments

Comments
 (0)