Force Models¶

use bh::GravityModelType;
use brahe as bh;

fn main() {
    // Create a fully-configured force model
    let force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::Earth,
        // Gravity: Spherical harmonic model (EGM2008, 20x20 degree/order)
        gravity: bh::GravityConfiguration::SphericalHarmonic {
            source: bh::GravityModelSource::ModelType(GravityModelType::EGM2008_120),
            degree: 20,
            order: 20,
            parallel: bh::orbit_dynamics::ParallelMode::Auto,
        },
        // Atmospheric drag: Harris-Priester model with parameter indices
        drag: Some(bh::DragConfiguration {
            model: bh::AtmosphericModel::HarrisPriester,
            area: bh::ParameterSource::ParameterIndex(1), // Index into parameter vector
            cd: bh::ParameterSource::ParameterIndex(2),
            body: None,
        }),
        // Solar radiation pressure: Conical eclipse model
        srp: Some(bh::SolarRadiationPressureConfiguration {
            area: bh::ParameterSource::ParameterIndex(3),
            cr: bh::ParameterSource::ParameterIndex(4),
            eclipse_model: bh::EclipseModel::Conical,
            occulting_bodies: vec![bh::OccultingBody::Earth],
        }),
        // Third-body: Sun and Moon with DE440s ephemeris (the default source
        // for bare bodies converted via .into())
        third_body: Some(vec![bh::ThirdBody::Sun.into(), bh::ThirdBody::Moon.into()]),
        // General relativistic corrections
        relativity: true,
        // Spacecraft mass (can also use ParameterIndex for estimation)
        mass: Some(bh::ParameterSource::Value(1000.0)), // kg
        // ECI<->ECEF rotation precision used by gravity, drag, and density models.
        // FullEarthRotation is the default and applies the full IAU 2006/2000A
        // bias-precession-nutation + Earth rotation + polar motion. Use
        // EarthRotationOnly for a faster, slightly less accurate alternative.
        frame_transform: bh::FrameTransformationModel::FullEarthRotation,
        // Tidal corrections to the gravity field (None disables tides)
        tides: None,
    };

    println!("Gravity: {:?}", force_config.gravity);
    println!("Drag: {:?}", force_config.drag);
    println!("SRP: {:?}", force_config.srp);
    println!("Third-body: {:?}", force_config.third_body);
    println!("Relativity: {:?}", force_config.relativity);
    println!("Mass: {:?}", force_config.mass);
    println!("Frame transform: {:?}", force_config.frame_transform);
}

use brahe as bh;

fn main() {
    // Point mass gravity configuration
    // Uses only central body gravity (mu/r^2)
    // No spherical harmonics, J2, or higher-order terms
    let _force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::Earth,
        gravity: bh::GravityConfiguration::PointMass,
        drag: None,
        srp: None,
        third_body: None,
        relativity: false,
        mass: None,
        frame_transform: bh::FrameTransformationModel::default(),
        tides: None,
    };
}

use brahe as bh;
use bh::GravityModelType;

fn main() {
    // ==========================================================================
    // Packaged Gravity Models
    // ==========================================================================

    // EGM2008 - High-fidelity NGA model (120x120 max)
    let _gravity_egm2008 = bh::GravityConfiguration::SphericalHarmonic {
        source: bh::GravityModelSource::ModelType(GravityModelType::EGM2008_120),
        degree: 20,
        order: 20,
        parallel: bh::orbit_dynamics::ParallelMode::Auto,
    };

    // GGM05S - GRACE mission model (180x180 max)
    let _gravity_ggm05s = bh::GravityConfiguration::SphericalHarmonic {
        source: bh::GravityModelSource::ModelType(GravityModelType::GGM05S),
        degree: 20,
        order: 20,
        parallel: bh::orbit_dynamics::ParallelMode::Auto,
    };

    // JGM3 - Legacy model, fast computation (70x70 max)
    let _gravity_jgm3 = bh::GravityConfiguration::SphericalHarmonic {
        source: bh::GravityModelSource::ModelType(GravityModelType::JGM3),
        degree: 20,
        order: 20,
        parallel: bh::orbit_dynamics::ParallelMode::Auto,
    };

    // ==========================================================================
    // Custom Model from File
    // ==========================================================================

    // Load custom gravity model from GFC format file
    // GravityModelType::from_file validates the path exists
    let custom_model_type =
        GravityModelType::from_file("data/gravity_models/EGM2008_120.gfc").unwrap();
    let _gravity_custom = bh::GravityConfiguration::SphericalHarmonic {
        source: bh::GravityModelSource::ModelType(custom_model_type),
        degree: 20,
        order: 20,
        parallel: bh::orbit_dynamics::ParallelMode::Auto,
    };
}

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

fn main() {
    // Reference an ICGEM Earth model. Use
    // `bh::datasets::icgem::list_icgem_models(ICGEMBody::Earth)` to discover
    // the full catalog. Append "-<DEGREE>" to pin a specific variant.
    let grav_type = GravityModelType::ICGEMModel {
        body: ICGEMBody::Earth,
        name: "JGM3".to_string(),
    };

    let _force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::Earth,
        gravity: bh::GravityConfiguration::SphericalHarmonic {
            source: bh::GravityModelSource::ModelType(grav_type),
            degree: 20,
            order: 20,
            parallel: bh::orbit_dynamics::ParallelMode::Auto,
        },
        drag: None,
        srp: None,
        third_body: None,
        relativity: false,
        mass: None,
        frame_transform: bh::FrameTransformationModel::default(),
        tides: None,
    };

    println!("Built ForceModelConfig with ICGEM gravity source: Earth/JGM3 (20x20)");
}

use brahe as bh;

fn main() {
    // Harris-Priester atmospheric drag configuration
    // - Valid for altitudes 100-1000 km
    // - Accounts for latitude-dependent diurnal bulge
    // - Does not require space weather data (F10.7, Ap)

    // Using parameter indices (default layout)
    let _drag_config = bh::DragConfiguration {
        model: bh::AtmosphericModel::HarrisPriester,
        area: bh::ParameterSource::ParameterIndex(1), // drag_area from params[1]
        cd: bh::ParameterSource::ParameterIndex(2),   // Cd from params[2]
        body: None, // Drag about the central body
    };
}

use brahe as bh;

fn main() {
    // Initialize space weather data provider
    bh::initialize_sw().unwrap();

    // NRLMSISE-00 atmospheric drag configuration
    // - Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar
    // - High-fidelity empirical model
    // - Valid from ground to thermospheric heights
    // - Uses space weather data (F10.7, Ap) when available
    // - More computationally expensive than Harris-Priester

    let _drag_config = bh::DragConfiguration {
        model: bh::AtmosphericModel::NRLMSISE00,
        area: bh::ParameterSource::ParameterIndex(1), // drag_area from params[1]
        cd: bh::ParameterSource::ParameterIndex(2),   // Cd from params[2]
        body: None, // Drag about the central body
    };
}

use bh::GravityModelType;
use brahe as bh;

fn main() {
    let drag_config = bh::DragConfiguration {
        model: bh::AtmosphericModel::Exponential {
            scale_height: 53000.0, // Scale height H in meters (53 km for ~300 km altitude)
            rho0: 1.225e-11,       // Reference density at h0 in kg/m^3
            h0: 300000.0,          // Reference altitude in meters (300 km)
        },
        area: bh::ParameterSource::ParameterIndex(1),
        cd: bh::ParameterSource::ParameterIndex(2),
        body: None,
    };

    // Create force model with exponential drag
    let _force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::Earth,
        gravity: bh::GravityConfiguration::SphericalHarmonic {
            source: bh::GravityModelSource::ModelType(GravityModelType::EGM2008_120),
            degree: 20,
            order: 20,
            parallel: bh::orbit_dynamics::ParallelMode::Auto,
        },
        drag: Some(drag_config),
        srp: None,
        third_body: None,
        relativity: false,
        mass: Some(bh::ParameterSource::ParameterIndex(0)),
        frame_transform: bh::FrameTransformationModel::default(),
        tides: None,
    };
}

use brahe as bh;

fn main() {
    // Solar Radiation Pressure configuration
    // Parameters:
    // - area: Cross-sectional area facing the Sun (m^2)
    // - cr: Coefficient of reflectivity (1.0=absorbing to 2.0=perfectly reflecting)
    // - eclipse_model: How to handle Earth's shadow

    // Option 1: No eclipse model (always illuminated)
    // Fast but inaccurate during eclipse periods
    let _srp_no_eclipse = bh::SolarRadiationPressureConfiguration {
        area: bh::ParameterSource::ParameterIndex(3), // srp_area from params[3]
        cr: bh::ParameterSource::ParameterIndex(4),   // Cr from params[4]
        eclipse_model: bh::EclipseModel::None,
        occulting_bodies: vec![bh::OccultingBody::Earth],
    };

    // Option 2: Cylindrical shadow model
    // Simple and fast, sharp shadow boundary (no penumbra)
    let _srp_cylindrical = bh::SolarRadiationPressureConfiguration {
        area: bh::ParameterSource::ParameterIndex(3),
        cr: bh::ParameterSource::ParameterIndex(4),
        eclipse_model: bh::EclipseModel::Cylindrical,
        occulting_bodies: vec![bh::OccultingBody::Earth],
    };

    // Option 3: Conical shadow model (recommended)
    // Accounts for penumbra and umbra regions
    let _srp_conical = bh::SolarRadiationPressureConfiguration {
        area: bh::ParameterSource::ParameterIndex(3),
        cr: bh::ParameterSource::ParameterIndex(4),
        eclipse_model: bh::EclipseModel::Conical,
        occulting_bodies: vec![bh::OccultingBody::Earth],
    };
}

use bh::GravityModelType;
use brahe as bh;

fn main() {
    // Third-body perturbations configuration: one entry per perturbing body,
    // each carrying its own ephemeris source and gravity model (point-mass
    // by default).

    // Option 1: Low-precision analytical ephemerides
    // Fast but less accurate (~km level errors for Sun/Moon)
    // Only Sun and Moon are available
    let _third_bodies_low = vec![
        bh::ThirdBodyConfiguration {
            ephemeris_source: bh::EphemerisSource::LowPrecision,
            ..bh::ThirdBody::Sun.into()
        },
        bh::ThirdBodyConfiguration {
            ephemeris_source: bh::EphemerisSource::LowPrecision,
            ..bh::ThirdBody::Moon.into()
        },
    ];

    // Option 2: DE440s high-precision ephemerides (recommended)
    // Uses JPL Development Ephemeris 440 (small bodies version)
    // ~m level accuracy, valid 1550-2650 CE
    // All planets available, ~17 MB file. DE440s is the default source, so
    // bare bodies convert directly into point-mass entries.
    let third_bodies_de440s: Vec<bh::ThirdBodyConfiguration> =
        vec![bh::ThirdBody::Sun.into(), bh::ThirdBody::Moon.into()];

    // Option 3: DE440 full-precision ephemerides
    // Highest accuracy (~mm level), valid 13200 BCE-17191 CE
    // All planets available, ~114 MB file
    let _third_bodies_de440 = vec![
        bh::ThirdBodyConfiguration {
            ephemeris_source: bh::EphemerisSource::DE440,
            ..bh::ThirdBody::Sun.into()
        },
        bh::ThirdBodyConfiguration {
            ephemeris_source: bh::EphemerisSource::DE440,
            ..bh::ThirdBody::Moon.into()
        },
    ];

    // Option 4: Include all major planets (high-fidelity). The *Barycenter
    // variants use the planetary-system barycenters with system GMs — the
    // classical third-body formulation, resolvable from the DE kernel alone.
    let _third_bodies_all_planets: Vec<bh::ThirdBodyConfiguration> = vec![
        bh::ThirdBody::Sun.into(),
        bh::ThirdBody::Moon.into(),
        bh::ThirdBody::Mercury.into(),
        bh::ThirdBody::Venus.into(),
        bh::ThirdBody::MarsBarycenter.into(),
        bh::ThirdBody::JupiterBarycenter.into(),
        bh::ThirdBody::SaturnBarycenter.into(),
        bh::ThirdBody::UranusBarycenter.into(),
        bh::ThirdBody::NeptuneBarycenter.into(),
    ];

    // Create force model with Sun/Moon perturbations (common case)
    let _force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::Earth,
        gravity: bh::GravityConfiguration::SphericalHarmonic {
            source: bh::GravityModelSource::ModelType(GravityModelType::EGM2008_120),
            degree: 20,
            order: 20,
            parallel: bh::orbit_dynamics::ParallelMode::Auto,
        },
        drag: None,
        srp: None,
        third_body: Some(third_bodies_de440s),
        relativity: false,
        mass: None,
        frame_transform: bh::FrameTransformationModel::default(),
        tides: None,
    };
}

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

//! EMB-centered cislunar propagation with Earth-attributed force models.
//! Earth contributes an 8x8 spherical-harmonic field as a third body and
//! NRLMSISE-00 drag evaluated at the object's Earth-relative state, so a
//! trajectory passing through LEO altitudes keeps Earth-fidelity forces
//! while the integration state stays Earth-Moon-barycenter-centered.

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

fn main() {
    // Initialize EOP and space weather data (required for NRLMSISE-00)
    bh::initialize_eop().unwrap();
    bh::initialize_sw().unwrap();

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

    // EMB-centered force model: no central gravity term (the barycenter has
    // no mass of its own); Earth carries a spherical-harmonic field and the
    // atmosphere; the Moon and Sun are point-mass perturbers.
    let force_config = bh::ForceModelConfig {
        central_body: bh::CentralBody::EMB,
        gravity: bh::GravityConfiguration::Zero,
        drag: Some(bh::DragConfiguration {
            model: bh::AtmosphericModel::NRLMSISE00,
            area: bh::ParameterSource::Value(10.0),
            cd: bh::ParameterSource::Value(2.2),
            // Attribute the drag to Earth: density and relative wind are
            // evaluated at the object's state relative to Earth.
            body: Some(bh::CentralBody::Earth),
        }),
        srp: None,
        third_body: Some(vec![
            bh::ThirdBodyConfiguration {
                body: bh::ThirdBody::Earth,
                ephemeris_source: bh::EphemerisSource::DE440s,
                gravity: bh::GravityConfiguration::SphericalHarmonic {
                    source: bh::GravityModelSource::default(),
                    degree: 8,
                    order: 8,
                    parallel: bh::orbit_dynamics::ParallelMode::Auto,
                },
            },
            bh::ThirdBody::Moon.into(),
            bh::ThirdBody::Sun.into(),
        ]),
        relativity: false,
        mass: Some(bh::ParameterSource::Value(1000.0)),
        frame_transform: bh::FrameTransformationModel::default(),
        tides: None,
    };
    force_config.validate().unwrap();

    // Start from a 500 km Earth orbit, re-expressed about the EMB via the
    // ECI->EMBI frame translation.
    let oe = na::SVector::<f64, 6>::new(bh::R_EARTH + 500e3, 0.001, 51.6, 15.0, 30.0, 45.0);
    let x_earth = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);
    let x_emb = bh::frames::state_eci_to_emb(epoch, x_earth);

    let mut prop = bh::DNumericalOrbitPropagator::new(
        epoch,
        na::DVector::from_column_slice(x_emb.as_slice()),
        bh::NumericalPropagationConfig::default(),
        force_config,
        None,
        None,
        None,
        None,
    )
    .unwrap();

    // Propagate for one day
    let epoch_end = epoch + 86400.0;
    prop.propagate_to(epoch_end);

    let x_final = prop.current_state();
    println!("Initial EMB-centered state: {:?}", x_emb.as_slice());
    println!(
        "Final EMB-centered state:   [{:.3}, {:.3}, {:.3}, {:.6}, {:.6}, {:.6}]",
        x_final[0], x_final[1], x_final[2], x_final[3], x_final[4], x_final[5]
    );

    // Re-express the final state about Earth for reference
    let x_final_eci = bh::frames::state_emb_to_eci(
        epoch_end,
        na::SVector::<f64, 6>::from_iterator(x_final.iter().cloned()),
    );
    let altitude =
        (x_final_eci[0].powi(2) + x_final_eci[1].powi(2) + x_final_eci[2].powi(2)).sqrt()
            - bh::R_EARTH;
    println!("Final altitude above Earth: {:.1} km", altitude / 1e3);
}