//! Two-body propagation over one orbital period.
//!
//! Run with `cargo run --example two_body`.

use astrokit::prelude::*;
use nalgebra::Vector3;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
    let frame = Frame::EARTH_J2000;

    let initial = Cartesian::new(
        Vector3::new(
            4.263868426884883e6,
            5.146189057155391e6,
            1.1310208421331816e6,
        ),
        Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558),
        frame,
    );

    let elements = initial.keplerian()?;
    let period = Duration::from_seconds(initial.period()?);
    let forces = vec![TwoBody { body: Body::EARTH }];

    let mut propagation = Propagation::from_cartesian(epoch, initial, &forces)
        .duration(period)
        .accuracy(Accuracy::High);
    let solution = propagation.run();

    let end = Orbit
        .state(*solution.states.last().unwrap())
        .cartesian(frame);

    println!("epoch:          {epoch}");
    println!("body:           {}", initial.frame.center.name());
    println!("period:         {period}");
    println!(
        "elements:       a = {:.3} km, e = {:.6}, i = {:.3} deg ({:?})",
        elements.semi_major_axis / 1e3,
        elements.eccentricity,
        elements.inclination.to_degrees(),
        elements.degeneracy
    );
    println!(
        "altitude:       {:.3} km",
        initial.spherical_altitude() / 1e3
    );
    println!("steps:          {}", solution.states.len() - 1);
    println!(
        "position error: {:.3e} m",
        (end.position - initial.position).norm()
    );
    println!(
        "velocity error: {:.3e} m/s",
        (end.velocity - initial.velocity).norm()
    );
    println!(
        "energy drift:   {:.3e} J/kg",
        end.specific_energy() - initial.specific_energy()
    );

    // Dense output: the state at any time in the arc, not only at a step.
    let half = Orbit
        .state(solution.interpolate(period.to_seconds() / 2.0))
        .cartesian(frame);
    println!(
        "half a period:  {:.3} km altitude",
        half.spherical_altitude() / 1e3
    );

    Ok(())
}