Search code examples
rustintegration

Rust ode_solvers get all the states throughout an integration (for plotting)


I am trying to learn to use the ode_solvers crate in rust, the documentation and (more informative) github docs can be found at:

https://docs.rs/ode_solvers/latest/ode_solvers/

https://srenevey.github.io/ode-solvers/index.html

I have been able to run simulations. However, I cannot work out how to get the values given by the steps in order to plot the path. E.g. one of the examples given in their documentation is a kepler orbit https://srenevey.github.io/ode-solvers/examples/kepler_orbit.html where they show a plot of the orbit but as far as I can see they do not show how it was plotted and I cannot work out how extract the values that the integration gives.

Please let me know how to extract the final and intermediate values so I can plot the path. Thanks,

The code for the kepler orbit example is:

use ode_solvers::{self, dop853::*};

type State = ode_solvers::Vector6<f64>;
type Time = f64;

struct KeplerOrbit {
    mu: f64,
}

impl ode_solvers::System<Time, State> for KeplerOrbit {
    // Equations of motion of the system
    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);
    }
}

fn main() {
    let system = KeplerOrbit {mu: 398600.435436};

    let a: f64 = 20000.0;
    let period = 2.0 * std::f64::consts::PI * (a.powi(3) / system.mu).sqrt();
    let y0 = State::new(-5007.248417988539, -1444.918140151374, 3628.534606178356, 0.717716656891, -10.224093784269, 0.748229399696);

    let mut stepper = Dop853::new(system, 0.0, 5.0 * period, 60.0, y0, 1.0e-10, 1.0e-10);
    let res = stepper.integrate();

    // Handle result
    match res {
        Ok(stats) => {
            println!("stats: {:?}", stats)
            
        },
        Err(_) => println!("An error occured."),
    }
}

Solution

  • Your system variable has x_out and y_out methods which you can use after integration to get the values for each integration step.