solarsystems 0.0.1

use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};
use solarsystems::gravity::interactions::*;
use solarsystems::{BodyId, CelestialBody, Vec3};

fn make_sun() -> CelestialBody {
    CelestialBody {
        id: BodyId::Sun,
        name: "Sun",
        mass: SOLAR_MASS,
        radius: 6.9634e8,
        position: Vec3::ZERO,
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 1.0e-7,
    }
}

fn make_earth() -> CelestialBody {
    CelestialBody {
        id: BodyId::Earth,
        name: "Earth",
        mass: 5.972e24,
        radius: 6.371e6,
        position: Vec3::new(AU, 0.0, 0.0),
        velocity: Vec3::new(0.0, 29784.0, 0.0),
        acceleration: Vec3::ZERO,
        j2: 1.08263e-3,
    }
}

fn make_jupiter() -> CelestialBody {
    CelestialBody {
        id: BodyId::Jupiter,
        name: "Jupiter",
        mass: 1.898e27,
        radius: 6.9911e7,
        position: Vec3::new(5.2 * AU, 0.0, 0.0),
        velocity: Vec3::new(0.0, 13070.0, 0.0),
        acceleration: Vec3::ZERO,
        j2: 0.01475,
    }
}

fn make_moon() -> CelestialBody {
    CelestialBody {
        id: BodyId::Moon,
        name: "Moon",
        mass: 7.342e22,
        radius: 1.7374e6,
        position: Vec3::new(AU + 3.844e8, 0.0, 0.0),
        velocity: Vec3::new(0.0, 30806.0, 0.0),
        acceleration: Vec3::ZERO,
        j2: 0.0,
    }
}

// ── gravitational_acceleration ──

#[test]
fn grav_accel_points_toward_source() {
    let earth = make_earth();
    let sun = make_sun();
    let acc = gravitational_acceleration(&earth, &sun);
    assert!(acc.x < 0.0);
}

#[test]
fn grav_accel_magnitude_inverse_square() {
    let sun = make_sun();
    let mut near = make_earth();
    near.position = Vec3::new(AU, 0.0, 0.0);
    let mut far = make_earth();
    far.position = Vec3::new(2.0 * AU, 0.0, 0.0);
    let acc_near = gravitational_acceleration(&near, &sun).magnitude();
    let acc_far = gravitational_acceleration(&far, &sun).magnitude();
    assert!((acc_near / acc_far - 4.0).abs() < 0.1);
}

#[test]
fn grav_accel_proportional_to_source_mass() {
    let earth = make_earth();
    let sun = make_sun();
    let mut light_sun = make_sun();
    light_sun.mass = SOLAR_MASS / 2.0;
    let acc_full = gravitational_acceleration(&earth, &sun).magnitude();
    let acc_half = gravitational_acceleration(&earth, &light_sun).magnitude();
    assert!((acc_full / acc_half - 2.0).abs() < 0.1);
}

// ── total_acceleration ──

#[test]
fn total_accel_sum_of_individuals() {
    let earth = make_earth();
    let sun = make_sun();
    let jupiter = make_jupiter();
    let bodies = vec![sun.clone(), earth.clone(), jupiter.clone()];
    let total = total_acceleration(&earth, &bodies);
    let from_sun = gravitational_acceleration(&earth, &sun);
    let from_jupiter = gravitational_acceleration(&earth, &jupiter);
    let manual_sum = from_sun + from_jupiter;
    assert!((total - manual_sum).magnitude() / total.magnitude() < 1e-6);
}

#[test]
fn total_accel_excludes_self() {
    let earth = make_earth();
    let bodies = vec![earth.clone()];
    let total = total_acceleration(&earth, &bodies);
    assert!(total.magnitude() < 1e-10);
}

// ── tidal_acceleration ──

#[test]
fn tidal_accel_nonzero() {
    let body_pos = Vec3::new(AU + 3.844e8, 0.0, 0.0);
    let perturber_pos = make_sun().position;
    let central_pos = Vec3::new(AU, 0.0, 0.0);
    let acc = tidal_acceleration(body_pos, SOLAR_MASS, perturber_pos, central_pos);
    assert!(acc.magnitude() > 0.0);
}

#[test]
fn tidal_accel_zero_at_central() {
    let central_pos = Vec3::new(AU, 0.0, 0.0);
    let perturber_pos = Vec3::ZERO;
    let acc = tidal_acceleration(central_pos, SOLAR_MASS, perturber_pos, central_pos);
    assert!(acc.magnitude() < 1e-10);
}

// ── mutual_gravitational_force ──

#[test]
fn mutual_force_positive() {
    let earth = make_earth();
    let sun = make_sun();
    let f = mutual_gravitational_force(&earth, &sun);
    assert!(f > 0.0);
}

#[test]
fn mutual_force_symmetric() {
    let earth = make_earth();
    let sun = make_sun();
    let f1 = mutual_gravitational_force(&earth, &sun);
    let f2 = mutual_gravitational_force(&sun, &earth);
    assert!((f1 - f2).abs() / f1 < 1e-10);
}

#[test]
fn mutual_force_earth_sun_known() {
    let earth = make_earth();
    let sun = make_sun();
    let f = mutual_gravitational_force(&earth, &sun);
    let expected = G * SOLAR_MASS * 5.972e24 / (AU * AU);
    assert!((f - expected).abs() / expected < 1e-4);
}

// ── tidal_force_magnitude ──

#[test]
fn tidal_force_mag_positive() {
    let f = tidal_force_magnitude(SOLAR_MASS, AU, 6.371e6);
    assert!(f > 0.0);
}

#[test]
fn tidal_force_stronger_closer() {
    let f_near = tidal_force_magnitude(SOLAR_MASS, AU, 6.371e6);
    let f_far = tidal_force_magnitude(SOLAR_MASS, 2.0 * AU, 6.371e6);
    assert!(f_near > f_far);
}

#[test]
fn tidal_force_cube_law() {
    let f1 = tidal_force_magnitude(SOLAR_MASS, AU, 6.371e6);
    let f2 = tidal_force_magnitude(SOLAR_MASS, 2.0 * AU, 6.371e6);
    assert!((f1 / f2 - 8.0).abs() < 0.5);
}

// ── j2_acceleration ──

#[test]
fn j2_accel_nonzero_for_oblate_body() {
    let earth = make_earth();
    let mut moon = make_moon();
    moon.position = Vec3::new(AU + 3.844e8, 0.0, 1e5);
    let acc = j2_acceleration(&moon, &earth);
    assert!(acc.magnitude() > 0.0);
}

#[test]
fn j2_accel_zero_j2_gives_zero() {
    let mut central = make_sun();
    central.j2 = 0.0;
    let earth = make_earth();
    let acc = j2_acceleration(&earth, &central);
    assert!(acc.magnitude() < 1e-20);
}

// ── relativistic_correction ──

#[test]
fn relativistic_correction_small() {
    let earth = make_earth();
    let sun = make_sun();
    let corr = relativistic_correction(&earth, &sun);
    let grav = gravitational_acceleration(&earth, &sun);
    assert!(corr.magnitude() / grav.magnitude() < 1e-4);
}

#[test]
fn relativistic_correction_nonzero() {
    let earth = make_earth();
    let sun = make_sun();
    let corr = relativistic_correction(&earth, &sun);
    assert!(corr.magnitude() > 0.0);
}