Gravity Models¶

use brahe as bh;
use bh::GravityModelType;
use bh::datasets::icgem::ICGEMBody;

fn main() {
    // Earth — JGM3 (~70x70), small and stable
    let earth_type = GravityModelType::ICGEMModel {
        body: ICGEMBody::Earth,
        name: "JGM3".to_string(),
    };
    let earth_model = bh::GravityModel::from_model_type(&earth_type).unwrap();
    println!(
        "Loaded {} ({}x{}, GM={:.6e} m^3/s^2)",
        earth_model.model_name, earth_model.n_max, earth_model.m_max, earth_model.gm
    );

    // Inspect a coefficient — normalized C_{2,0} (J2 term) from the Earth model
    let (c20, _s20) = earth_model.get(2, 0).unwrap();
    println!("\nEarth normalized C(2,0) = {c20:.6e}");
}

//! nalgebra = "0.33"
//! ```

//! Compute point-mass gravitational acceleration for an Earth satellite

#[allow(unused_imports)]
use brahe as bh;
use nalgebra as na;

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

    // Define satellite position in ECI frame (LEO satellite at 500 km altitude)
    // Using Keplerian elements and converting to Cartesian
    let a = bh::constants::R_EARTH + 500e3; // Semi-major axis (m)
    let e = 0.001;                          // Eccentricity
    let i = 97.8;                           // Inclination (deg)
    let raan = 0.0;                         // RAAN (deg)
    let argp = 0.0;                         // Argument of perigee (deg)
    let nu = 0.0;                           // True anomaly (deg)

    // Convert to Cartesian state
    let oe = na::SVector::<f64, 6>::new(a, e, i, raan, argp, nu);
    let state = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let r_sat = na::Vector3::new(state[0], state[1], state[2]); // Position vector (m)

    println!("Satellite position (ECI, m):");
    println!("  x = {:.3}", r_sat[0]);
    println!("  y = {:.3}", r_sat[1]);
    println!("  z = {:.3}", r_sat[2]);

    // Compute point-mass gravitational acceleration
    // For Earth-centered case, central body is at origin
    let r_earth = na::Vector3::<f64>::zeros();
    let accel = bh::orbit_dynamics::accel_point_mass_gravity(r_sat, r_earth, bh::constants::GM_EARTH);

    println!("\nPoint-mass gravity acceleration (m/s²):");
    println!("  ax = {:.6}", accel[0]);
    println!("  ay = {:.6}", accel[1]);
    println!("  az = {:.6}", accel[2]);

    // Compute magnitude
    let accel_mag = accel.norm();
    println!("\nAcceleration magnitude: {:.6} m/s²", accel_mag);

    // Compare to theoretical value: GM/r²
    let r_mag = r_sat.norm();
    let accel_theoretical = bh::constants::GM_EARTH / (r_mag * r_mag);
    println!("Theoretical magnitude: {:.6} m/s²", accel_theoretical);

}

//! nalgebra = "0.33"
//! ```

//! Compute spherical harmonic gravitational acceleration for an Earth satellite

#[allow(unused_imports)]
use brahe as bh;
use nalgebra as na;

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

    // Create an epoch for frame transformations
    let epoch = bh::Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh::TimeSystem::UTC);

    // Define satellite position in ECI frame (LEO satellite at 500 km altitude)
    let a = bh::constants::R_EARTH + 500e3; // Semi-major axis (m)
    let e = 0.001;                          // Eccentricity
    let i = 97.8;                           // Inclination (deg)
    let raan = 45.0;                        // RAAN (deg)
    let argp = 30.0;                        // Argument of perigee (deg)
    let nu = 60.0;                          // True anomaly (deg)

    // Convert to Cartesian state
    let oe = na::SVector::<f64, 6>::new(a, e, i, raan, argp, nu);
    let state_eci = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let r_eci = na::Vector3::new(state_eci[0], state_eci[1], state_eci[2]); // Position vector (m)

    println!("Satellite position (ECI, m):");
    println!("  x = {:.3}", r_eci[0]);
    println!("  y = {:.3}", r_eci[1]);
    println!("  z = {:.3}", r_eci[2]);

    // Load gravity model (GGM05S - degree/order 180)
    let gravity_model = bh::orbit_dynamics::GravityModel::from_model_type(
        &bh::orbit_dynamics::GravityModelType::GGM05S
    ).unwrap();
    println!("\nGravity model: GGM05S (max degree {}, max order {})",
             gravity_model.n_max, gravity_model.m_max);

    // For spherical harmonics, we need the ECI to body-fixed rotation matrix
    // This rotates from ECI (inertial) to ECEF (Earth-fixed) frame
    let r_eci_ecef = bh::rotation_eci_to_ecef(epoch);

    // Compute spherical harmonic acceleration (degree 10, order 10)
    let n_max = 10;
    let m_max = 10;
    let accel_sh = bh::orbit_dynamics::accel_gravity_spherical_harmonics(
        r_eci, r_eci_ecef, &gravity_model, n_max, m_max
    );

    println!("\nSpherical harmonic acceleration (degree {}, order {}):", n_max, m_max);
    println!("  ax = {:.9} m/s²", accel_sh[0]);
    println!("  ay = {:.9} m/s²", accel_sh[1]);
    println!("  az = {:.9} m/s²", accel_sh[2]);

    // Compute point-mass for comparison
    let accel_pm = bh::orbit_dynamics::accel_point_mass_gravity(
        r_eci, na::Vector3::<f64>::zeros(), bh::constants::GM_EARTH
    );

    println!("\nPoint-mass acceleration:");
    println!("  ax = {:.9} m/s²", accel_pm[0]);
    println!("  ay = {:.9} m/s²", accel_pm[1]);
    println!("  az = {:.9} m/s²", accel_pm[2]);

    // Compute difference (perturbation due to non-spherical Earth)
    let accel_pert = accel_sh - accel_pm;

    println!("\nPerturbation (spherical harmonics - point mass):");
    println!("  Δax = {:.9} m/s²", accel_pert[0]);
    println!("  Δay = {:.9} m/s²", accel_pert[1]);
    println!("  Δaz = {:.9} m/s²", accel_pert[2]);
    println!("  Magnitude: {:.9} m/s²", accel_pert.norm());

}