//! Orbit and attitude propagated together, with the area that follows from them.
//!
//! Run with `cargo run --example attitude`.

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

fn main() {
    let epoch = Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0);
    let frame = Frame::EARTH_J2000;
    let forces = vec![TwoBody { body: Body::EARTH }];

    // A 400 km circular orbit at the inclination of the space station. Circular, so
    // the argument of periapsis is undefined and the last angle is the argument of
    // latitude; `Keplerian::new` works that out rather than being told.
    let orbit = Keplerian::new(6.7781e6, 0.0, 51.6_f64.to_radians(), 0.4, 0.0, 0.0, frame);
    let spacecraft = Spacecraft::rectangular_box(500.0, Vector3::new(2.0, 1.0, 3.0), 2.2, 1.3);
    // A tumble about no principal axis, so all three rates change.
    let rate = Vector3::new(0.0012, 0.0028, -0.0007);

    let initial = Full
        .with_keplerian(orbit)
        .with_mass(spacecraft.mass)
        .with_angular_rate(rate);

    let period = orbit.period();
    let mut propagation = Propagation::new(epoch, initial, &forces)
        .frame(frame)
        .spacecraft(spacecraft.clone())
        .duration(Duration::from_seconds(period))
        .accuracy(Accuracy::High);
    let solution = propagation.run();

    let state = Full.state(*solution.states.last().unwrap());
    // Switches of the parameters to their shadow set, which the propagation does for
    // itself whenever the layout carries an attitude.
    let switches = jump_times(&solution.times).len();

    // The same arc without the attitude, to confirm the orbit is untouched by it.
    let mut translational = Propagation::new(epoch, Orbit.with_keplerian(orbit), &forces)
        .frame(frame)
        .duration(Duration::from_seconds(period))
        .accuracy(Accuracy::High);
    let reference = Orbit
        .state(*translational.run().states.last().unwrap())
        .cartesian(frame);

    let inertia = &spacecraft.inertia;
    let momentum = |w: Vector3<f64>| (inertia * w).norm();
    let energy = |w: Vector3<f64>| (w.transpose() * inertia * w)[0] * 0.5;

    // Drag sees the face turned into the flow, which is the velocity direction
    // taken into body components.
    let presented_area = |s: &State<f64, 14, Full>| {
        let incident = attitude::to_dcm(s.attitude()) * -s.velocity().normalize();
        spacecraft.surface.drag_coefficient_area(incident)
    };

    println!("epoch:            {epoch}");
    println!(
        "orbit:            a = {:.1} km, i = {:.1} deg ({:?})",
        orbit.semi_major_axis / 1e3,
        orbit.inclination.to_degrees(),
        orbit.degeneracy
    );
    println!("period:           {}", Duration::from_seconds(period));
    println!(
        "inertia:          diag({:.1}, {:.1}, {:.1}) kg m^2",
        inertia[(0, 0)],
        inertia[(1, 1)],
        inertia[(2, 2)]
    );
    println!(
        "steps:            {}, shadow-set switches: {switches}",
        solution.times.len() - 1 - switches
    );
    let final_rate = state.angular_rate();
    println!(
        "rate:             {:.4} deg/s about ({:.2}, {:.2}, {:.2})",
        rate.norm().to_degrees(),
        rate[0] / rate.norm(),
        rate[1] / rate.norm(),
        rate[2] / rate.norm()
    );
    println!(
        "          ->      {:.4} deg/s about ({:.2}, {:.2}, {:.2})",
        final_rate.norm().to_degrees(),
        final_rate[0] / final_rate.norm(),
        final_rate[1] / final_rate.norm(),
        final_rate[2] / final_rate.norm()
    );
    // An asymmetric body under no torque changes both its axis and its rate. What
    // it keeps is the angular momentum and the rotational energy.
    println!(
        "momentum drift:   {:.3e} relative",
        (momentum(final_rate) - momentum(rate)).abs() / momentum(rate)
    );
    println!(
        "energy drift:     {:.3e} relative",
        (energy(final_rate) - energy(rate)).abs() / energy(rate)
    );
    println!(
        "presented area:   {:.3} -> {:.3} m^2 of C_d A",
        presented_area(&initial),
        presented_area(&state)
    );
    println!(
        "orbit, against a translational propagation of the same arc: {:.3e} m",
        (state.position() - reference.position).norm()
    );
}