solarsystems 0.0.1

use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};
use solarsystems::config::parameters::{
    GRAVITATIONAL_SOFTENING, J2_EARTH, J2_JUPITER, J2_MARS, J2_SATURN, J2_SUN, MU_SUN,
    default_bodies, default_bodies_with_belts, find_body, find_body_mut, system_barycenter,
    total_angular_momentum, total_energy,
};
use solarsystems::exports::state::{
    extract_3d, extract_5d, extract_6d, extract_7d, extract_8d, extract_by_id, extract_full,
    extract_full_flat, extract_pairwise_distances, extract_positions_flat,
    extract_state_vectors_flat, relative_state,
};
use solarsystems::gravity::interactions::*;
use solarsystems::gravity::nbody::*;
use solarsystems::{BodyId, CelestialBody, Vec3};

#[test]
fn gravitational_acceleration_inverse_square() {
    let body = CelestialBody {
        id: BodyId::Earth,
        name: "Earth",
        mass: 5.972e24,
        radius: 6.371e6,
        position: Vec3::new(AU, 0.0, 0.0),
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let sun = CelestialBody {
        id: BodyId::Sun,
        name: "Sun",
        mass: SOLAR_MASS,
        radius: 6.9634e8,
        position: Vec3::ZERO,
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let acc = gravitational_acceleration(&body, &sun);
    let expected = G * SOLAR_MASS / (AU * AU);
    assert!((acc.magnitude() - expected).abs() / expected < 1e-4);
    assert!(acc.x < 0.0);
}

#[test]
fn gravitational_force_earth_sun() {
    let body = CelestialBody {
        id: BodyId::Earth,
        name: "Earth",
        mass: 5.972e24,
        radius: 6.371e6,
        position: Vec3::new(AU, 0.0, 0.0),
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let sun = CelestialBody {
        id: BodyId::Sun,
        name: "Sun",
        mass: SOLAR_MASS,
        radius: 6.9634e8,
        position: Vec3::ZERO,
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let f = mutual_gravitational_force(&body, &sun);
    let expected = 3.542e22;
    assert!((f - expected).abs() / expected < 0.01);
}

#[test]
fn tidal_force_moon_on_earth() {
    let earth_radius = 6.371e6;
    let moon_mass = 7.342e22;
    let distance = 3.844e8;
    let tf = tidal_force_magnitude(moon_mass, distance, earth_radius);
    assert!(tf > 0.0);
    assert!(tf < 1e-5);
}

#[test]
fn j2_acceleration_equatorial_vs_polar() {
    let central = CelestialBody {
        id: BodyId::Earth,
        name: "Earth",
        mass: 5.972e24,
        radius: 6.371e6,
        position: Vec3::ZERO,
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 1.08263e-3,
    };
    let equatorial = CelestialBody {
        id: BodyId::Moon,
        name: "Moon",
        mass: 7.342e22,
        radius: 1.7374e6,
        position: Vec3::new(4e8, 0.0, 0.0),
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let polar = CelestialBody {
        id: BodyId::Moon,
        name: "Moon",
        mass: 7.342e22,
        radius: 1.7374e6,
        position: Vec3::new(0.0, 0.0, 4e8),
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let acc_eq = j2_acceleration(&equatorial, &central);
    let acc_po = j2_acceleration(&polar, &central);
    assert!(acc_eq.magnitude() != acc_po.magnitude());
}

#[test]
fn nbody_conserves_barycenter_momentum() {
    let mut bodies = default_bodies();
    compute_accelerations(&mut bodies);
    let (_, v0) = barycenter(&bodies);
    for _ in 0..100 {
        solarsystems::dynamics::integrator::leapfrog_kdk(&mut bodies, 3600.0);
    }
    let (_, v1) = barycenter(&bodies);
    let dv = (v1 - v0).magnitude();
    assert!(dv < 1e-4, "Barycenter velocity drift: {}", dv);
}

#[test]
fn relativistic_correction_mercury() {
    let sun = CelestialBody {
        id: BodyId::Sun,
        name: "Sun",
        mass: SOLAR_MASS,
        radius: 6.9634e8,
        position: Vec3::ZERO,
        velocity: Vec3::ZERO,
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let mercury = CelestialBody {
        id: BodyId::Mercury,
        name: "Mercury",
        mass: 3.3011e23,
        radius: 2.4397e6,
        position: Vec3::new(0.387 * AU, 0.0, 0.0),
        velocity: Vec3::new(0.0, 47_362.0, 0.0),
        acceleration: Vec3::ZERO,
        j2: 0.0,
    };
    let gr = relativistic_correction(&mercury, &sun);
    let newtonian = gravitational_acceleration(&mercury, &sun);
    assert!(gr.magnitude() / newtonian.magnitude() < 1e-6);
}

#[test]
fn total_system_energy_negative() {
    let bodies = default_bodies();
    let e = total_system_energy(&bodies);
    assert!(e < 0.0, "Bound system should have negative total energy");
}

#[test]
fn tidal_acceleration_is_differential() {
    let body_pos = Vec3::new(AU + 3.844e8, 0.0, 0.0);
    let perturber_pos = Vec3::new(5.203 * AU, 0.0, 0.0);
    let central_pos = Vec3::new(AU, 0.0, 0.0);
    let ta = tidal_acceleration(body_pos, 1.8982e27, perturber_pos, central_pos);
    assert!(ta.magnitude() > 0.0);
}

#[test]
fn total_acceleration_includes_all_bodies() {
    let bodies = default_bodies();
    let earth = find_body(&bodies, BodyId::Earth).unwrap();
    let acc = total_acceleration(earth, &bodies);
    assert!(acc.magnitude() > 0.0);
    let grav_sun_only = gravitational_acceleration(earth, find_body(&bodies, BodyId::Sun).unwrap());
    let diff = (acc.magnitude() - grav_sun_only.magnitude()).abs() / grav_sun_only.magnitude();
    assert!(
        diff < 0.01,
        "Other bodies should contribute small perturbations"
    );
}

#[test]
fn compute_accelerations_relativistic_gives_extra_correction() {
    let mut bodies1 = default_bodies();
    compute_accelerations(&mut bodies1);
    let a_newton: f64 = bodies1.iter().map(|b| b.acceleration.magnitude()).sum();

    let mut bodies2 = default_bodies();
    compute_accelerations_relativistic(&mut bodies2);
    let a_rel: f64 = bodies2.iter().map(|b| b.acceleration.magnitude()).sum();

    assert!(
        (a_newton - a_rel).abs() / a_newton < 1e-4,
        "Relativistic correction is small"
    );
    assert!(a_newton != a_rel, "Should differ slightly");
}

#[test]
fn compute_pure_gravitational_no_j2() {
    let mut bodies1 = default_bodies();
    compute_accelerations(&mut bodies1);
    let a_full: f64 = bodies1.iter().map(|b| b.acceleration.magnitude()).sum();

    let mut bodies2 = default_bodies();
    compute_pure_gravitational(&mut bodies2);
    let a_pure: f64 = bodies2.iter().map(|b| b.acceleration.magnitude()).sum();

    let diff = (a_full - a_pure).abs() / a_full;
    assert!(diff < 0.1, "Pure gravitational should be close to full");
}

#[test]
fn total_system_angular_momentum_nonzero() {
    let bodies = default_bodies();
    let am = total_system_angular_momentum(&bodies);
    assert!(am.magnitude() > 0.0);
}

#[test]
fn default_bodies_category_counts() {
    let bodies = default_bodies();
    let sun_count = bodies.iter().filter(|b| b.id == BodyId::Sun).count();
    let planets = bodies.iter().filter(|b| b.id.is_planet()).count();
    let satellites = bodies.iter().filter(|b| b.id.is_satellite()).count();
    let dwarfs = bodies.iter().filter(|b| b.id.is_dwarf_planet()).count();
    let asteroids = bodies.iter().filter(|b| b.id.is_asteroid()).count();
    let comets = bodies.iter().filter(|b| b.id.is_comet()).count();

    assert_eq!(sun_count, 1);
    assert_eq!(planets, 8);
    assert_eq!(satellites, 12);
    assert_eq!(dwarfs, 5);
    assert_eq!(asteroids, 5);
    assert_eq!(comets, 5);
    assert_eq!(
        bodies.len(),
        sun_count + planets + satellites + dwarfs + asteroids + comets
    );
}

#[test]
fn default_bodies_with_belts_has_extra() {
    let bodies = default_bodies_with_belts(50, 50);
    let belt_a = bodies.iter().filter(|b| b.id.is_belt_asteroid()).count();
    let kuiper = bodies.iter().filter(|b| b.id.is_kuiper_object()).count();
    assert_eq!(belt_a, 50);
    assert_eq!(kuiper, 50);
    assert!(bodies.len() > default_bodies().len());
}

#[test]
fn all_bodies_have_positive_mass() {
    let bodies = default_bodies();
    for body in &bodies {
        assert!(body.mass > 0.0, "{} should have positive mass", body.name);
    }
}

#[test]
fn all_heliocentric_bodies_have_nonzero_velocity() {
    let bodies = default_bodies();
    for body in &bodies {
        if body.id == BodyId::Sun {
            continue;
        }
        if body.id.is_satellite() {
            continue;
        }
        assert!(
            body.velocity.magnitude() > 0.0,
            "{} should have nonzero heliocentric velocity",
            body.name
        );
    }
}

#[test]
fn find_body_and_find_body_mut() {
    let bodies = default_bodies();
    assert!(find_body(&bodies, BodyId::Earth).is_some());
    assert!(find_body(&bodies, BodyId::Sun).is_some());
    assert!(find_body(&bodies, BodyId::BeltAsteroid(999)).is_none());

    let mut bodies_mut = default_bodies();
    if let Some(earth) = find_body_mut(&mut bodies_mut, BodyId::Earth) {
        earth.position.x += 1.0;
    }
    let earth = find_body(&bodies_mut, BodyId::Earth).unwrap();
    assert!(
        (earth.position.x - find_body(&bodies, BodyId::Earth).unwrap().position.x - 1.0).abs()
            < 1e-10
    );
}

#[test]
fn system_barycenter_near_sun() {
    let bodies = default_bodies();
    let (pos, _vel) = system_barycenter(&bodies);
    let sun = find_body(&bodies, BodyId::Sun).unwrap();
    let dist = pos.distance(&sun.position);
    assert!(dist < 2e9, "Barycenter should be near Sun, dist={}", dist);
}

#[test]
fn total_energy_and_angular_momentum_consistent() {
    let bodies = default_bodies();
    let te = total_energy(&bodies);
    let tam = total_angular_momentum(&bodies);
    let nbe = total_system_energy(&bodies);
    let nbam = total_system_angular_momentum(&bodies);
    assert!((te - nbe).abs() / te.abs() < 1e-10, "Should be consistent");
    assert!((tam.magnitude() - nbam.magnitude()).abs() / tam.magnitude() < 1e-10);
}

#[test]
fn config_constants_sane() {
    let gs = GRAVITATIONAL_SOFTENING;
    let je = J2_EARTH;
    let jj = J2_JUPITER;
    let js = J2_SATURN;
    let jm = J2_MARS;
    let jsun = J2_SUN;
    let mu = MU_SUN;
    assert!(gs > 0.0);
    assert!(je > 0.0);
    assert!(jj > 0.0);
    assert!(js > 0.0);
    assert!(jm > 0.0);
    assert!(jsun >= 0.0);
    assert!(mu > 0.0);
}

#[test]
fn body_id_all_methods() {
    let all = BodyId::all();
    let planets = BodyId::all_planets();
    let sats = BodyId::all_satellites();
    let dwarfs = BodyId::all_dwarf_planets();
    let asteroids = BodyId::all_asteroids();
    let comets = BodyId::all_comets();

    assert_eq!(planets.len(), 8);
    assert!(sats.len() >= 12);
    assert!(dwarfs.len() >= 5);
    assert!(asteroids.len() >= 5);
    assert!(comets.len() >= 5);
    assert!(
        all.len() >= planets.len() + sats.len() + dwarfs.len() + asteroids.len() + comets.len()
    );

    for p in planets {
        assert!(p.is_planet());
        assert!(!p.is_satellite());
        assert!(!p.is_dwarf_planet());
    }
    for s in sats {
        assert!(s.is_satellite());
        assert!(s.parent().is_some());
    }
    for d in dwarfs {
        assert!(d.is_dwarf_planet());
    }
    for a in asteroids {
        assert!(a.is_asteroid());
    }
    for c in comets {
        assert!(c.is_comet());
    }

    assert!(!BodyId::BeltAsteroid(0).is_planet());
    assert!(BodyId::BeltAsteroid(0).is_belt_asteroid());
    assert!(BodyId::BeltAsteroid(0).is_belt_particle());
    assert!(BodyId::KuiperObject(0).is_kuiper_object());
    assert!(BodyId::KuiperObject(0).is_belt_particle());
}

#[test]
fn celestial_body_derived_properties() {
    let bodies = default_bodies();
    let sun = find_body(&bodies, BodyId::Sun).unwrap();
    let earth = find_body(&bodies, BodyId::Earth).unwrap();
    let jupiter = find_body(&bodies, BodyId::Jupiter).unwrap();
    let earth_a = earth.position.distance(&sun.position);

    let gp = earth.gravitational_parameter();
    assert!(gp > 0.0);

    let sg = earth.surface_gravity();
    assert!((sg - 9.81).abs() < 0.5, "g = {}", sg);

    let ev = earth.escape_velocity();
    assert!(ev > 10_000.0 && ev < 12_000.0, "v_esc = {}", ev);

    let md = earth.mean_density();
    assert!(md > 4000.0 && md < 6000.0, "density = {}", md);

    let sr = sun.schwarzschild_radius();
    assert!(sr > 2000.0 && sr < 4000.0, "r_s = {}", sr);

    let hs = earth.hill_sphere(sun.mass, earth_a);
    assert!(hs > 1e9, "Hill sphere = {}", hs);

    let rl = earth.roche_limit(3344.0);
    assert!(rl > earth.radius, "Roche limit = {}", rl);

    let soi = earth.sphere_of_influence(sun.mass, earth_a);
    assert!(soi > 1e8, "SOI = {}", soi);

    let period = earth.orbital_period(sun.mass);
    let period_days = period / 86400.0;
    assert!(
        (period_days - 365.25).abs() < 2.0,
        "T = {} days",
        period_days
    );

    let ke = earth.kinetic_energy();
    assert!(ke > 0.0);

    let pe = earth.potential_energy(sun);
    assert!(pe < 0.0);

    let am = earth.angular_momentum();
    assert!(am.magnitude() > 0.0);

    let spec_am = earth.specific_angular_momentum();
    assert!(spec_am.magnitude() > 0.0);

    let spec_oe = earth.specific_orbital_energy(sun.mass);
    assert!(
        spec_oe < 0.0,
        "Bound orbit should have negative specific energy"
    );

    let vv = earth.vis_viva_speed(sun.mass, earth_a);
    let v_circ = (G * sun.mass / earth_a).sqrt();
    assert!((vv - v_circ).abs() / v_circ < 0.01, "vis-viva ~ circular");

    let jev = jupiter.escape_velocity();
    assert!(jev > 50_000.0, "Jupiter escape = {}", jev);
}

#[test]
fn vec3_extended_methods() {
    let a = Vec3::new(3.0, 4.0, 0.0);
    let b = Vec3::new(0.0, 5.0, 0.0);

    assert!((a.magnitude_squared() - 25.0).abs() < 1e-10);
    let n = a.normalize();
    assert!((n.magnitude() - 1.0).abs() < 1e-10);
    assert!((a.dot(&b) - 20.0).abs() < 1e-10);
    let cross = a.cross(&b);
    assert!((cross.z - 15.0).abs() < 1e-10);

    let angle = a.angle_between(&b);
    assert!(angle > 0.0 && angle < std::f64::consts::PI);

    let rx = a.rotate_x(std::f64::consts::FRAC_PI_2);
    assert!((rx.x - 3.0).abs() < 1e-10);
    assert!(rx.z.abs() > 3.0);

    let ry = a.rotate_y(std::f64::consts::FRAC_PI_2);
    assert!(ry.x.abs() < 1e-10);

    let rz = a.rotate_z(std::f64::consts::FRAC_PI_2);
    assert!((rz.x + 4.0).abs() < 1e-10);

    let lerp = a.lerp(&b, 0.0);
    assert!((lerp.x - a.x).abs() < 1e-10);
    let lerp1 = a.lerp(&b, 1.0);
    assert!((lerp1.x - b.x).abs() < 1e-10);
    let lerp_mid = a.lerp(&b, 0.5);
    assert!((lerp_mid.x - 1.5).abs() < 1e-10);
}

#[test]
fn extract_3d_returns_correct_count() {
    let bodies = default_bodies();
    let states = extract_3d(&bodies);
    assert_eq!(states.len(), bodies.len());
}

#[test]
fn extract_5d_has_speed_and_mass() {
    let bodies = default_bodies();
    let states = extract_5d(&bodies);
    assert_eq!(states.len(), bodies.len());
    for s in &states {
        assert!(s.mass > 0.0);
        assert!(s.speed >= 0.0);
    }
}

#[test]
fn extract_6d_has_velocities() {
    let bodies = default_bodies();
    let states = extract_6d(&bodies);
    assert_eq!(states.len(), bodies.len());
    let earth_state = &states[bodies.iter().position(|b| b.id == BodyId::Earth).unwrap()];
    assert!(earth_state.velocity.magnitude() > 0.0);
}

#[test]
fn extract_7d_has_mass() {
    let bodies = default_bodies();
    let states = extract_7d(&bodies);
    assert_eq!(states.len(), bodies.len());
    for s in &states {
        assert!(s.mass > 0.0);
    }
}

#[test]
fn extract_8d_has_radius() {
    let bodies = default_bodies();
    let states = extract_8d(&bodies);
    assert_eq!(states.len(), bodies.len());
    for s in &states {
        assert!(s.radius > 0.0);
    }
}

#[test]
fn extract_full_and_by_id() {
    let bodies = default_bodies();
    let full = extract_full(&bodies);
    assert_eq!(full.len(), bodies.len());
    let earth_full = extract_by_id(&bodies, BodyId::Earth).unwrap();
    assert!(earth_full.kinetic_energy > 0.0);
    assert!(earth_full.distance_from_origin > 0.0);
    assert!(extract_by_id(&bodies, BodyId::BeltAsteroid(999)).is_none());
}

#[test]
fn extract_flat_formats() {
    let bodies = default_bodies();
    let pos_flat = extract_positions_flat(&bodies);
    assert_eq!(pos_flat.len(), bodies.len() * 3);
    let sv_flat = extract_state_vectors_flat(&bodies);
    assert_eq!(sv_flat.len(), bodies.len() * 6);
    let full_flat = extract_full_flat(&bodies);
    assert!(full_flat.len() >= bodies.len() * 6);
}

#[test]
fn extract_pairwise_distances_symmetric() {
    let bodies = default_bodies();
    let pw = extract_pairwise_distances(&bodies);
    let n = bodies.len();
    assert_eq!(pw.len(), n * (n - 1) / 2);
    for (_, _, d) in &pw {
        assert!(*d > 0.0);
    }
}

#[test]
fn relative_state_moon_earth() {
    let bodies = default_bodies();
    let rel = relative_state(&bodies, BodyId::Moon, BodyId::Earth).unwrap();
    assert!(rel.position.magnitude() > 3e8);
    assert!(rel.position.magnitude() < 5e8);
    assert!(rel.velocity.magnitude() > 0.0);
    assert!(relative_state(&bodies, BodyId::BeltAsteroid(999), BodyId::Earth).is_none());
}