Module nyx_space::dynamics

Provides several dynamics used for orbital mechanics and attitude dynamics, which can be elegantly combined.

Combining dynamics in a full spacecraft model.

extern crate nalgebra as na;
extern crate nyx_space as nyx;
// Warning: this is arguably a bad example: attitude dynamics very significantly
// faster than orbital mechanics. Hence you really should use different propagators
// for the attitude and orbital position and velocity.
use self::nyx::dynamics::Dynamics;
use self::nyx::dynamics::celestial::TwoBody;
use self::nyx::dynamics::momentum::AngularMom;
use self::nyx::celestia::EARTH;
use self::nyx::propagators::{CashKarp45, Options, Propagator};
use self::na::{Matrix3, U9, Vector3, Vector6, VectorN};

// In the following struct, we only store the dynamics because this is only a proof
// of concept. An engineer could add more useful information to this struct, such
// as a short cut to the position or an attitude.
#[derive(Copy, Clone)]
pub struct PosVelAttMom {
    pub twobody: TwoBody,
    pub momentum: AngularMom,
}

impl Dynamics for PosVelAttMom {
    type StateSize = U9;
    fn time(&self) -> f64 {
        // Both dynamical models have the same time because they share the propagator.
        self.twobody.time()
    }

    fn state(&self) -> VectorN<f64, Self::StateSize> {
        let twobody_state = self.twobody.state();
        let momentum_state = self.momentum.state();
        // We're channing both states to create a combined state.
        // The most important part here is make sure that the `state` and `set_state` handle the state in the same order.
        <VectorN<f64, U9>>::from_iterator(
            twobody_state.iter().chain(momentum_state.iter()).cloned(),
        )
    }

    fn set_state(&mut self, new_t: f64, new_state: &VectorN<f64, Self::StateSize>) {
        // HACK: Reconstructing the Vector6 from scratch because for some reason it isn't the correct type when using `fixed_slice`.
        // No doubt there's a more clever way to handle this, I just haven't figured it out yet.
        let mut pos_vel_vals = [0.0; 6];
        let mut mom_vals = [0.0; 3];
        for (i, val) in new_state.iter().enumerate() {
            if i < 6 {
                pos_vel_vals[i] = *val;
            } else {
                mom_vals[i - 6] = *val;
            }
        }
        self.twobody
            .set_state(new_t, &Vector6::from_row_slice(&pos_vel_vals));
        self.momentum
            .set_state(new_t, &Vector3::from_row_slice(&mom_vals));
    }

    fn eom(
        &self,
        _t: f64,
        state: &VectorN<f64, Self::StateSize>,
    ) -> VectorN<f64, Self::StateSize> {
        // Same issue as in `set_state`.
        let mut pos_vel_vals = [0.0; 6];
        let mut mom_vals = [0.0; 3];
        for (i, val) in state.iter().enumerate() {
            if i < 6 {
                pos_vel_vals[i] = *val;
            } else {
                mom_vals[i - 6] = *val;
            }
        }
        let dpos_vel_dt = self.twobody
            .eom(_t, &Vector6::from_row_slice(&pos_vel_vals));
        let domega_dt = self.momentum.eom(_t, &Vector3::from_row_slice(&mom_vals));
        <VectorN<f64, U9>>::from_iterator(dpos_vel_dt.iter().chain(domega_dt.iter()).cloned())
    }
}

// Let's initialize our combined dynamics.

fn main(){
    let dyn_twobody = TwoBody::around(
        &Vector6::from_row_slice(&[-2436.45, -2436.45, 6891.037, 5.088611, -5.088611, 0.0]),
        &EARTH,
    );

    let omega = Vector3::new(0.1, 0.4, -0.2);
    let tensor = Matrix3::new(10.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 2.0);
    let dyn_mom = AngularMom::from_tensor_matrix(&tensor, &omega);

    let mut full_model = PosVelAttMom {
        twobody: dyn_twobody,
        momentum: dyn_mom,
    };

    let init_momentum = full_model.momentum.momentum().norm();
    let mom_tolerance = 1e-8;

    // And now let's define the propagator and propagate for a short amount of time.
    let mut prop = Propagator::new::<CashKarp45>(&Options::with_adaptive_step(0.01, 30.0, 1e-12));

    // And propagate
    loop {
        let (t, state) = prop.derive(
            full_model.time(),
            &full_model.state(),
            |t_: f64, state_: &VectorN<f64, U9>| full_model.eom(t_, state_),
        );
        full_model.set_state(t, &state);
        if full_model.time() >= 3600.0 {
            println!("{:?}", prop.latest_details());
            println!("{}", full_model.state());
            let delta_mom =
                ((full_model.momentum.momentum().norm() - init_momentum) / init_momentum).abs();
            if delta_mom > mom_tolerance {
                panic!(
                    "angular momentum prop failed: momentum changed by {:e} (> {:e})",
                    delta_mom, mom_tolerance
                );
            }
            break;
        }
    }
}

Modules

celestial	The celestial module handles all Cartesian based dynamics.
momentum	The angular momentum module handles all angular momentum dynamics.

Traits

Dynamics

The Dynamics trait handles and stores any equation of motion and the state is integrated.