use ode_solvers::dopri5::*;
use ode_solvers::*;
type State = DVector<f64>;
type Time = f64;
use std::{f64::consts::PI, fs::File, io::BufWriter, io::Write, path::Path};
fn main() {
let system = KeplerOrbit { mu: 398600.435436 };
let a: f64 = 20000.0;
let period = 2.0 * PI * (a.powi(3) / system.mu).sqrt();
let y0 = State::from(vec![
-5007.248417988539,
-1444.918140151374,
3628.534606178356,
0.717716656891,
-10.224093784269,
0.748229399696,
]);
let mut stepper = Dopri5::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10);
let res = stepper.integrate();
match res {
Ok(stats) => {
println!("{}", stats);
let path = Path::new("./outputs/kepler_orbit_dopri5_dvector.dat");
save(stepper.x_out(), stepper.y_out(), path);
println!("Results saved in: {:?}", path);
}
Err(e) => println!("An error occured: {}", e),
}
}
struct KeplerOrbit {
mu: f64,
}
impl ode_solvers::System<f64, State> for KeplerOrbit {
fn system(&self, _t: Time, y: &State, dy: &mut State) {
let r = (y[0] * y[0] + y[1] * y[1] + y[2] * y[2]).sqrt();
dy[0] = y[3];
dy[1] = y[4];
dy[2] = y[5];
dy[3] = -self.mu * y[0] / r.powi(3);
dy[4] = -self.mu * y[1] / r.powi(3);
dy[5] = -self.mu * y[2] / r.powi(3);
}
}
pub fn save(times: &Vec<Time>, states: &Vec<State>, filename: &Path) {
let file = match File::create(filename) {
Err(e) => {
println!("Could not open file. Error: {:?}", e);
return;
}
Ok(buf) => buf,
};
let mut buf = BufWriter::new(file);
for (i, state) in states.iter().enumerate() {
buf.write_fmt(format_args!("{}", times[i])).unwrap();
for val in state.iter() {
buf.write_fmt(format_args!(", {}", val)).unwrap();
}
buf.write_fmt(format_args!("\n")).unwrap();
}
if let Err(e) = buf.flush() {
println!("Could not write to file. Error: {:?}", e);
}
}