Permalink
Browse files

Initial work on RKF54. Panics are runtime because I can't get an elem…

…ent as an f64.

Maybe is there a more memory efficient way to do this?
  • Loading branch information...
ChristopherRabotin committed Feb 13, 2018
1 parent 331015b commit 1dbaa888a08d0a3bf843b0442aa2c5340067600e
Showing with 155 additions and 4 deletions.
  1. +71 −0 src/propagators/classic_rk.rs
  2. +14 −4 src/propagators/mod.rs
  3. +2 −0 tests/lib.rs
  4. +68 −0 tests/propagators.rs
@@ -0,0 +1,71 @@
extern crate nalgebra as na;

use self::na::DMatrix;
pub use super::{Options, RK};

pub struct RKF54 {
opts: Options,
}

/// RKF54 is a [Runge Kutta Ferhlberg integrator](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method).
impl RK for RKF54 {
fn a_coeffs() -> DMatrix<f64> {
DMatrix::from_row_slice(
5,
5,
&[
1.0 / 4.0,
0.0,
0.0,
0.0,
0.0, // end row 1
3.0 / 32.0,
9.0 / 32.0,
0.0,
0.0,
0.0, // end row 2
1932.0 / 2197.0,
-7200.0 / 2197.0,
7296.0 / 2197.0,
0.0,
0.0, // end row 3
439.0 / 216.0,
-8.0,
3680.0 / 513.0,
-845.0 / 4104.0,
0.0, // end row 4
-8.0 / 27.0,
2.0,
-3544.0 / 2565.0,
1859.0 / 4104.0,
-11.0 / 40.0, // end row 5
],
)
}
fn b_coeffs() -> DMatrix<f64> {
DMatrix::from_row_slice(
2,
6,
&[
16.0 / 135.0,
0.0,
6656.0 / 12825.0,
28561.0 / 56430.0,
-9.0 / 50.0,
2.0 / 55.0,
25.0 / 216.0,
0.0,
1408.0 / 2565.0,
2197.0 / 4104.0,
-1.0 / 5.0,
0.0,
],
)
}
fn from_options(opts: Options) -> RKF54 {
RKF54 { opts: opts }
}
fn options(self) -> Options {
self.opts
}
}
@@ -2,6 +2,10 @@ extern crate nalgebra as na;

use self::na::{DMatrix, DVector};

// Re-Export
mod classic_rk;
pub use self::classic_rk::*;

/// The `RK` trait defines a Runge Kutta
pub trait RK
where
@@ -13,7 +17,7 @@ where
/// Returns a Matrix corresponding to the b_i and b^*_i coefficients of the Butcher table for that RK
fn b_coeffs() -> DMatrix<f64>;
fn from_options(opts: Options) -> Self;
fn options(&self) -> Options;
fn options(self) -> Options;
}

#[derive(Clone)]
@@ -53,7 +57,7 @@ impl Propagator {
let ci: f64 = self.a_coeffs.row(i).iter().sum();
let mut wi = DVector::from_element(self.state.shape().0 as usize, 0.0);
for j in 0..k.len() {
let a_ij = self.a_coeffs.slice((i, j), (1, 1));
let a_ij = self.a_coeffs.slice((i, j), (1, 1)).0;
wi += a_ij * &k[j];
}
let ki = d_xdt(
@@ -66,8 +70,8 @@ impl Propagator {
let mut next_state = self.state.clone();
let mut next_state_star = self.state.clone();
for i in 0..k.len() {
let b_ij = self.a_coeffs.slice((i, 0), (1, 1));
let b_ij_star = self.a_coeffs.slice((i, 1), (1, 1));
let b_ij = self.b_coeffs.slice((i, 0), (1, 1));
let b_ij_star = self.b_coeffs.slice((i, 1), (1, 1));
next_state += self.opts.min_step * b_ij * &k[i];
next_state_star += self.opts.min_step * b_ij_star * &k[i];
}
@@ -147,3 +151,9 @@ fn test_options() {
opts.enable_debug();
assert_eq!(opts.debug, true);
}

#[test]
fn test_consistency() {
// All the RKs should be consistent.
println!("{}", RKF54::a_coeffs());
}
@@ -1 +1,3 @@
extern crate nyx;

mod propagators;
@@ -0,0 +1,68 @@
extern crate nalgebra as na;
extern crate nyx;
use self::na::DVector;

fn two_body_dynamics(_t: f64, state: DVector<f64>) -> DVector<f64> {
let mut radius = state.slice((0, 0), (3, 1)); // TODO: Change to compile time slice
let velocity = state.slice((3, 0), (3, 1));
let body_acceleration = (-398600.4 / radius.norm().powi(3)) * radius;
let test = DVector::from_iterator(6, velocity.iter().chain(body_acceleration.iter()).cloned());
println!("State:\n{}", radius, test);
DVector::from_row_slice(6, &[-2436.45, -2436.45, 6891.037, 5.088611, -5.088611, 0.0])
}

#[test]
fn geo_day_prop() {
extern crate nyx;
extern crate nalgebra as na;
use nyx::propagators::*;
use self::na::DVector;
//let init_state = Vector6::new(-2436.45, -2436.45, 6891.037, 5.088611, -5.088611, 0.0);
let init_state =
DVector::from_row_slice(6, &[-2436.45, -2436.45, 6891.037, 5.088611, -5.088611, 0.0]);
let method = RKF54::from_options(Options::with_fixed_step(1.0));
let mut prop = Propagator::new::<RKF54>(0.0, init_state, method);
prop.prop(two_body_dynamics);

/*
virtObj := CelestialObject{"virtObj", 6378.145, 149598023, 398600.4, 23.4, 0.00005, 924645.0, 0.00108248, -2.5324e-6, -1.6204e-6, 0, nil}
orbit := NewOrbitFromRV([]float64{-2436.45, -2436.45, 6891.037}, []float64{5.088611, -5.088611, 0}, virtObj)
startDT := time.Date(2017, 1, 1, 0, 0, 0, 0, time.UTC)
endDT := startDT.Add(24 * time.Hour).Add(time.Second)
NewPreciseMission(NewEmptySC("est", 0), orbit, startDT, endDT, Perturbations{}, time.Second, false, ExportConfig{}).Propagate()
expR := []float64{-5971.19544867343, 3945.58315019255, 2864.53021742433}
expV := []float64{0.049002818030, -4.185030861883, 5.848985672439}
if !floats.EqualApprox(orbit.rVec, expR, 1e-8) {
t.Fatalf("Incorrect R:\ngot: %+v\nexp: %+v", orbit.rVec, expR)
}
if !floats.EqualApprox(orbit.vVec, expV, 1e-8) {
t.Fatalf("Incorrect R:\ngot: %+v\nexp: %+v", orbit.vVec, expV)
}
R := []float64{f[0], f[1], f[2]}
V := []float64{f[3], f[4], f[5]}
tmpOrbit = NewOrbitFromRV(R, V, a.Orbit.Origin)
bodyAcc := -tmpOrbit.Origin.μ / math.Pow(Norm(R), 3)
_, _, i, Ω, _, _, _, _, u := tmpOrbit.Elements()
// Check if any impulse burn, and execute them if needed.
if maneuver, exists := a.Vehicle.Maneuvers[a.CurrentDT.Truncate(a.step)]; exists {
if !maneuver.done {
a.Vehicle.logger.Log("level", "info", "subsys", "astro", "date", a.CurrentDT, "thrust", "impulse", "v(km/s)", maneuver.Δv())
Δv[0] += maneuver.R
Δv[1] += maneuver.N
Δv[2] += maneuver.C
maneuver.done = true
a.Vehicle.Maneuvers[a.CurrentDT.Truncate(a.step)] = maneuver
}
}
Δv = Rot313Vec(-u, -i, -Ω, Δv)
// d\vec{R}/dt
fDot[0] = f[3]
fDot[1] = f[4]
fDot[2] = f[5]
// d\vec{V}/dt
fDot[3] = bodyAcc*f[0] + Δv[0]
fDot[4] = bodyAcc*f[1] + Δv[1]
fDot[5] = bodyAcc*f[2] + Δv[2]
*/
}

0 comments on commit 1dbaa88

Please sign in to comment.