OPM — Orbit Parameter Message¶

use brahe as bh;
use bh::ccsds::OPM;
use bh::traits::DStatePropagator;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();
    bh::initialize_sw().unwrap();

    // Parse OPM — use Example1 which has spacecraft mass
    let opm = OPM::from_file("test_assets/ccsds/opm/OPMExample1.txt").unwrap();
    println!(
        "Object: {} ({})",
        opm.metadata.object_name, opm.metadata.object_id
    );
    println!("Epoch:  {}", opm.state_vector.epoch);
    println!("Frame:  {}", opm.metadata.ref_frame);

    // Extract initial conditions from OPM
    let pos = opm.state_vector.position;
    let vel = opm.state_vector.velocity;
    println!("\nInitial state:");
    println!(
        "  Position: [{:.3}, {:.3}, {:.3}] km",
        pos[0] / 1e3,
        pos[1] / 1e3,
        pos[2] / 1e3
    );
    println!(
        "  Velocity: [{:.3}, {:.3}, {:.3}] m/s",
        vel[0], vel[1], vel[2]
    );

    // Build spacecraft parameters from OPM
    let sc = opm.spacecraft_parameters.as_ref();
    let mass = sc.and_then(|s| s.mass).unwrap_or(500.0);
    let drag_area = sc.and_then(|s| s.drag_area).unwrap_or(2.0);
    let drag_coeff = sc.and_then(|s| s.drag_coeff).unwrap_or(2.2);
    let srp_area = sc.and_then(|s| s.solar_rad_area).unwrap_or(2.0);
    let srp_coeff = sc.and_then(|s| s.solar_rad_coeff).unwrap_or(1.3);
    let params = na::DVector::from_vec(vec![mass, drag_area, drag_coeff, srp_area, srp_coeff]);
    println!(
        "\nSpacecraft params: mass={}kg, Cd={}, Cr={}",
        mass, drag_coeff, srp_coeff
    );

    // The message declares its state in ITRF2000; the propagator expects GCRF,
    // and state_in_frame maps the declared frame through the reference frame router
    let state_eci = opm.state_in_frame(bh::CelestialFrame::GCRF).unwrap();
    let state_dyn = na::DVector::from_column_slice(state_eci.as_slice());

    // Initialize propagator from OPM state
    let mut prop = bh::DNumericalOrbitPropagator::builder(
        opm.state_vector.epoch,
        state_dyn,
        bh::ForceModelConfig::default(),
    )
    .params(params)
    .build()
    .unwrap();

    // Propagate for ~1 orbit period
    let r = (pos[0].powi(2) + pos[1].powi(2) + pos[2].powi(2)).sqrt();
    let period = 2.0 * std::f64::consts::PI * (r.powi(3) / bh::GM_EARTH).sqrt();
    println!("\nEstimated period: {:.0}s ({:.1} min)", period, period / 60.0);

    let target_epoch = opm.state_vector.epoch + period;
    prop.propagate_to(target_epoch).unwrap();

    // Check final state
    let final_state = prop.current_state();
    println!("\nAfter 1 orbit:");
    println!("  Epoch: {}", prop.current_epoch());
    println!(
        "  Position: [{:.3}, {:.3}, {:.3}] km",
        final_state[0] / 1e3,
        final_state[1] / 1e3,
        final_state[2] / 1e3
    );
    println!(
        "  Velocity: [{:.3}, {:.3}, {:.3}] m/s",
        final_state[3], final_state[4], final_state[5]
    );
}

use brahe as bh;
use brahe::ccsds::OPM;

fn main() {
    bh::initialize_eop().unwrap();

    // Parse OPM with Keplerian elements and maneuvers
    let opm = OPM::from_file("test_assets/ccsds/opm/OPMExample2.txt").unwrap();

    // Header
    println!("Format version: {}", opm.header.format_version);
    println!("Originator:     {}", opm.header.originator);
    println!("Creation date:  {}", opm.header.creation_date);

    // Metadata
    println!("\nObject name:  {}", opm.metadata.object_name);
    println!("Object ID:    {}", opm.metadata.object_id);
    println!("Center name:  {}", opm.metadata.center_name);
    println!("Ref frame:    {}", opm.metadata.ref_frame);
    println!("Time system:  {}", opm.metadata.time_system);

    // State vector (SI units: meters, m/s)
    println!("\nEpoch:    {}", opm.state_vector.epoch);
    let pos = opm.state_vector.position;
    let vel = opm.state_vector.velocity;
    println!(
        "Position: [{:.4}, {:.4}, {:.4}] km",
        pos[0] / 1e3,
        pos[1] / 1e3,
        pos[2] / 1e3
    );
    println!(
        "Velocity: [{:.8}, {:.8}, {:.8}] m/s",
        vel[0], vel[1], vel[2]
    );

    // Keplerian elements
    if let Some(ref kep) = opm.keplerian_elements {
        println!("\nKeplerian elements:");
        println!("  Semi-major axis:    {:.4} km", kep.semi_major_axis / 1e3);
        println!("  Eccentricity:       {:.9}", kep.eccentricity);
        println!("  Inclination:        {:.6} deg", kep.inclination);
        println!("  RAAN:               {:.6} deg", kep.ra_of_asc_node);
        println!(
            "  Arg of pericenter:  {:.6} deg",
            kep.arg_of_pericenter
        );
        if let Some(ta) = kep.true_anomaly {
            println!("  True anomaly:       {:.6} deg", ta);
        }
        println!("  GM:                 {:.4e} m³/s²", kep.gm.unwrap_or(0.0));
    }

    // Spacecraft parameters
    if let Some(ref sc) = opm.spacecraft_parameters {
        println!("\nMass:           {} kg", sc.mass.unwrap_or(0.0));
        println!(
            "Solar rad area: {} m²",
            sc.solar_rad_area.unwrap_or(0.0)
        );
        println!(
            "Solar rad coef: {}",
            sc.solar_rad_coeff.unwrap_or(0.0)
        );
        println!("Drag area:      {} m²", sc.drag_area.unwrap_or(0.0));
        println!("Drag coeff:     {}", sc.drag_coeff.unwrap_or(0.0));
    }

    // Maneuvers
    println!("\nManeuvers: {}", opm.maneuvers.len());
    for (i, man) in opm.maneuvers.iter().enumerate() {
        println!("\n  Maneuver {}:", i);
        println!("    Epoch ignition: {}", man.epoch_ignition);
        println!("    Duration:       {} s", man.duration);
        println!(
            "    Delta mass:     {} kg",
            man.delta_mass.map_or("None".to_string(), |m| format!("{}", m))
        );
        println!("    Ref frame:      {}", man.ref_frame);
        println!(
            "    Delta-V:        [{:.5}, {:.5}, {:.5}] m/s",
            man.dv[0], man.dv[1], man.dv[2]
        );
    }
}

use brahe as bh;
use bh::ccsds::{CCSDSRefFrame, OPM};
use bh::events::{DTimeEvent, EventAction};
use bh::traits::DStatePropagator;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();
    bh::initialize_sw().unwrap();

    // Parse OPM with maneuvers
    let opm = OPM::from_file("test_assets/ccsds/opm/OPMExample2.txt").unwrap();
    println!("Object: {}", opm.metadata.object_name);
    println!("Epoch:  {}", opm.state_vector.epoch);
    println!("Maneuvers: {}", opm.maneuvers.len());

    // Extract initial state; the OPM declares its state in the TOD frame
    let state_eci = opm.state_in_frame(bh::CelestialFrame::GCRF).unwrap();

    // Spacecraft parameters
    let sc = opm.spacecraft_parameters.as_ref();
    let mass = sc.and_then(|s| s.mass).unwrap_or(500.0);
    let params = na::DVector::from_vec(vec![
        mass,
        sc.and_then(|s| s.drag_area).unwrap_or(10.0),
        sc.and_then(|s| s.drag_coeff).unwrap_or(2.3),
        sc.and_then(|s| s.solar_rad_area).unwrap_or(10.0),
        sc.and_then(|s| s.solar_rad_coeff).unwrap_or(1.3),
    ]);

    // Create propagator
    let mut prop = bh::DNumericalOrbitPropagator::builder(
        opm.state_vector.epoch,
        na::DVector::from_column_slice(state_eci.as_slice()),
        bh::ForceModelConfig::default(),
    )
    .params(params)
    .build()
    .unwrap();

    // The message's ignition epochs precede its state epoch, so maneuvers are
    // scheduled relative to the state epoch, preserving the spacing between them
    let first_ignition = opm.maneuvers[0].epoch_ignition;
    let scheduled_epochs: Vec<bh::Epoch> = opm
        .maneuvers
        .iter()
        .map(|man| opm.state_vector.epoch + 3600.0 + (man.epoch_ignition - first_ignition))
        .collect();

    // The frame bias between EME2000 (alias J2000) and GCRF is epoch-independent
    let r_eme2000_to_gcrf = bh::rotation_eme2000_to_gcrf();

    // Add an event detector for each maneuver
    for (i, man) in opm.maneuvers.iter().enumerate() {
        let sched_epoch = scheduled_epochs[i];
        let dv = na::Vector3::new(man.dv[0], man.dv[1], man.dv[2]);
        let dv_mag = dv.norm();
        let idx = i;

        // Inertial delta-Vs rotate into GCRF once; RTN delta-Vs depend on the
        // state at the ignition epoch and rotate inside the callback
        let is_rtn = match man.ref_frame {
            CCSDSRefFrame::J2000 | CCSDSRefFrame::EME2000 => false,
            CCSDSRefFrame::RTN => true,
            ref other => panic!("Unsupported maneuver reference frame: {}", other),
        };
        let dv_frame = if is_rtn { dv } else { r_eme2000_to_gcrf * dv };

        let callback: bh::events::DEventCallback = Box::new(
            move |_t: bh::Epoch,
                  state: &na::DVector<f64>,
                  _params: Option<&na::DVector<f64>>|
                  -> (Option<na::DVector<f64>>, Option<na::DVector<f64>>, EventAction) {
                let dv_gcrf = if is_rtn {
                    let x = na::SVector::<f64, 6>::from_column_slice(&state.as_slice()[..6]);
                    bh::rotation_rtn_to_eci(x) * dv_frame
                } else {
                    dv_frame
                };
                let mut new_state = state.clone();
                new_state[3] += dv_gcrf[0];
                new_state[4] += dv_gcrf[1];
                new_state[5] += dv_gcrf[2];
                println!("  Applied maneuver {}: |dv|={:.3} m/s", idx, dv_gcrf.norm());
                (Some(new_state), None, EventAction::Continue)
            },
        );

        let event =
            DTimeEvent::new(sched_epoch, format!("Maneuver-{}", i)).with_callback(callback);
        prop.add_event_detector(Box::new(event));
        println!(
            "  Registered maneuver {}: epoch={}, frame={}, |dv|={:.3} m/s",
            i, sched_epoch, man.ref_frame, dv_mag
        );
    }

    // Propagate past all maneuvers
    let target = *scheduled_epochs.last().unwrap() + 3600.0;
    println!("\nPropagating to {}...", target);
    prop.propagate_to(target).unwrap();

    // Report final state
    let final_state = prop.current_state();
    println!("\nFinal state at {}:", prop.current_epoch());
    println!(
        "  Position: [{:.3}, {:.3}, {:.3}] km",
        final_state[0] / 1e3,
        final_state[1] / 1e3,
        final_state[2] / 1e3
    );
    println!(
        "  Velocity: [{:.3}, {:.3}, {:.3}] m/s",
        final_state[3], final_state[4], final_state[5]
    );

    println!("\nExample completed successfully!");
}