Propagation Around Other Central Bodies¶

use brahe as bh;
use bh::traits::{DOrbitStateProvider, DStatePropagator};
use nalgebra as na;

fn main() {
    // Initialize EOP data and the DE440s planetary ephemeris used for
    // third-body perturbations and ECI <-> MCI frame conversions.
    bh::initialize_eop().unwrap();
    bh::load_common_spice_kernels().unwrap();

    // Initial epoch
    let epoch = bh::Epoch::from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh::TimeSystem::UTC);

    // Circular low Mars orbit at 400 km altitude, expressed directly in the
    // Mars-Centered Inertial (MCI) frame the propagator integrates in.
    let a = bh::R_MARS + 400e3;
    let v = bh::periapsis_velocity(a, 0.0, bh::GM_MARS);
    let state = na::DVector::from_vec(vec![a, 0.0, 0.0, 0.0, v, 0.0]);

    // Spacecraft parameters, indexed per mars_default()'s ParameterSource
    // assignments: [mass, drag_area, Cd, srp_area, Cr].
    let params = na::DVector::from_vec(vec![500.0, 2.0, 2.2, 2.0, 1.3]);

    let force_config = bh::ForceModelConfig::mars_default();
    force_config.validate().unwrap();

    let mut prop = bh::DNumericalOrbitPropagator::new(
        epoch,
        state,
        bh::NumericalPropagationConfig::default(),
        force_config,
        Some(params),
        None, // No additional dynamics
        None, // No control input
        None, // No initial covariance
    )
    .unwrap();

    // Propagate for 6 hours
    let final_epoch = epoch + 6.0 * 3600.0;
    prop.propagate_to(final_epoch);

    // state_in_frame routes the propagator's native MCI state through the
    // reference frame router into any other supported frame. MCMF is the
    // Mars-fixed, IAU/WGCCRE body-fixed frame.
    let x0_mcmf = prop.state_in_frame(bh::ReferenceFrame::MCMF, epoch).unwrap();
    let xf_mcmf = prop
        .state_in_frame(bh::ReferenceFrame::MCMF, final_epoch)
        .unwrap();

    println!("Initial epoch: {}", epoch);
    println!("Final epoch:   {}", final_epoch);
    println!("\nInitial state (MCMF, Mars-fixed):");
    println!(
        "  Position (km): [{:.3}, {:.3}, {:.3}]",
        x0_mcmf[0] / 1e3,
        x0_mcmf[1] / 1e3,
        x0_mcmf[2] / 1e3
    );
    println!(
        "  Velocity (m/s): [{:.3}, {:.3}, {:.3}]",
        x0_mcmf[3], x0_mcmf[4], x0_mcmf[5]
    );
    println!("\nFinal state (MCMF, Mars-fixed):");
    println!(
        "  Position (km): [{:.3}, {:.3}, {:.3}]",
        xf_mcmf[0] / 1e3,
        xf_mcmf[1] / 1e3,
        xf_mcmf[2] / 1e3
    );
    println!(
        "  Velocity (m/s): [{:.3}, {:.3}, {:.3}]",
        xf_mcmf[3], xf_mcmf[4], xf_mcmf[5]
    );

    // Validate propagation completed and the orbit remains bound to Mars
    assert_eq!(prop.current_epoch(), final_epoch);
    let r_final = xf_mcmf.fixed_rows::<3>(0).norm();
    assert!(r_final > bh::R_MARS && r_final < bh::R_MARS + 1000e3);
    println!("\nExample validated successfully!");
}

use brahe as bh;
use nalgebra as na;

fn main() {
    // An uncatalogued body (e.g. a newly observed asteroid): self-assign a
    // unique negative NAIF ID for its center, mirroring NAIF's convention.
    const CENTER: i32 = -20001;

    // The frame key is an arbitrary integer handle. It is unrelated to NAIF
    // IDs; it only names the callback registered below so BodyFixedCustom can
    // look it up at evaluation time.
    const KEY: u32 = 1042;

    let t0 = bh::Epoch::from_date(2024, 3, 1, bh::TimeSystem::TDB);
    let rate = 2.0e-4; // spin rate (rad/s), ~8.7 h rotation period

    // The rotation callback is a proper function: it receives an Epoch and
    // returns the 3x3 ICRF -> body-fixed rotation matrix at that instant.
    // Here the body spins uniformly about its z-axis.
    let rotation = move |epc: bh::Epoch| {
        let theta = rate * (epc - t0);
        let (s, c) = theta.sin_cos();
        na::Matrix3::new(c, s, 0.0, -s, c, 0.0, 0.0, 0.0, 1.0)
    };

    // The optional omega callback is also a function of the epoch, returning
    // the frame's angular velocity vector (rad/s) for the exact velocity
    // transport term. Passing None falls back to numeric differentiation of
    // the rotation callback.
    let omega = move |_epc: bh::Epoch| na::Vector3::new(0.0, 0.0, rate);

    bh::register_custom_frame(KEY, rotation, Some(Box::new(omega)));

    let inertial = bh::ReferenceFrame::BodyCenteredICRF(CENTER);
    let fixed = bh::ReferenceFrame::BodyFixedCustom {
        center: CENTER,
        key: KEY,
    };

    // Convert an inertial state about the body into its body-fixed frame.
    // Both frames share the same center, so no ephemeris kernel is needed.
    let epc = t0 + 600.0;
    let x_inertial = na::SVector::<f64, 6>::new(7.0e5, -2.0e5, 3.0e5, 10.0, 25.0, -5.0);
    let x_fixed = bh::state_frame_to_frame(inertial, fixed, epc, x_inertial).unwrap();

    println!("Epoch: {}", epc);
    println!("Inertial state:   {:.3?}", x_inertial.as_slice());
    println!("Body-fixed state: {:.3?}", x_fixed.as_slice());

    // Round trip back to the inertial frame recovers the input.
    let x_back = bh::state_frame_to_frame(fixed, inertial, epc, x_fixed).unwrap();
    assert!((x_back - x_inertial).norm() < 1e-6);

    // A point co-rotating with the body is stationary in the body-fixed frame.
    let r_surf = rotation(epc).transpose() * na::Vector3::new(1.0e3, 0.0, 0.0);
    let v_surf = omega(epc).cross(&r_surf);
    let x_surf = na::SVector::<f64, 6>::new(
        r_surf[0], r_surf[1], r_surf[2], v_surf[0], v_surf[1], v_surf[2],
    );
    let x_surf_fixed = bh::state_frame_to_frame(inertial, fixed, epc, x_surf).unwrap();
    assert!(x_surf_fixed.fixed_rows::<3>(3).norm() < 1e-9);

    bh::unregister_custom_frame(KEY);
    println!("\nExample validated successfully!");
}