solarsystems 0.0.1

use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};
use solarsystems::Vec3;
use solarsystems::orbits::keplerian::synodic_period;
use solarsystems::orbits::keplerian::*;
use solarsystems::orbits::perturbations::*;
use solarsystems::orbits::transfers::*;
use std::f64::consts::PI;

const MU_SUN: f64 = G * SOLAR_MASS;
const MU_EARTH: f64 = G * 5.972e24;

fn earth_oe() -> OrbitalElements {
    OrbitalElements {
        semi_major_axis: AU,
        eccentricity: 0.0167,
        inclination: 0.0,
        longitude_ascending_node: 0.0,
        argument_periapsis: 0.0,
        true_anomaly: 0.0,
    }
}

// ── OrbitalElements ──

#[test]
fn circular_orbit_eccentricity_zero() {
    let r = AU;
    let v = (MU_SUN / r).sqrt();
    let pos = Vec3::new(r, 0.0, 0.0);
    let vel = Vec3::new(0.0, v, 0.0);
    let oe = OrbitalElements::from_state_vectors(pos, vel, MU_SUN);
    assert!(oe.eccentricity < 1e-6, "e = {}", oe.eccentricity);
    assert!((oe.semi_major_axis - AU).abs() / AU < 1e-6);
}

#[test]
fn earth_orbital_period() {
    let oe = earth_oe();
    let period_days = oe.orbital_period(MU_SUN) / 86400.0;
    assert!(
        (period_days - 365.25).abs() < 1.0,
        "T = {} days",
        period_days
    );
    let period_from_kepler = 2.0 * PI * (AU.powi(3) / MU_SUN).sqrt();
    assert!((oe.orbital_period(MU_SUN) - period_from_kepler).abs() < 1.0);
}

#[test]
fn periapsis_apoapsis_consistency() {
    let oe = earth_oe();
    let q = oe.periapsis();
    let big_q = oe.apoapsis();
    assert!((q + big_q) / 2.0 - AU < AU * 1e-6);
    assert!(big_q > q);
}

#[test]
fn state_vector_roundtrip() {
    let r = 1.5 * AU;
    let v = (MU_SUN / r).sqrt() * 0.9;
    let pos = Vec3::new(r, 0.0, 0.0);
    let vel = Vec3::new(0.0, v, 0.0);
    let oe = OrbitalElements::from_state_vectors(pos, vel, MU_SUN);
    let (pos2, vel2) = oe.to_state_vectors(MU_SUN);
    assert!((pos.x - pos2.x).abs() / pos.x < 1e-8);
    assert!((vel.y - vel2.y).abs() / vel.y < 1e-8);
}

#[test]
fn kepler_equation_convergence() {
    let ma = 1.5;
    let e = 0.5;
    let ea = solve_kepler(ma, e);
    let residual = (ea - e * ea.sin() - ma).abs();
    assert!(residual < 1e-12, "residual = {}", residual);
}

#[test]
fn kepler_equation_circular() {
    let ma = 2.0;
    let ea = solve_kepler(ma, 0.0);
    assert!((ea - ma).abs() < 1e-14);
}

#[test]
fn kepler_equation_high_eccentricity() {
    let ma = 0.5;
    let e = 0.99;
    let ea = solve_kepler(ma, e);
    let residual = (ea - e * ea.sin() - ma).abs();
    assert!(residual < 1e-10);
}

#[test]
fn synodic_period_earth_mars() {
    let t_earth = 365.25;
    let t_mars = 687.0;
    let syn = synodic_period(t_earth, t_mars);
    assert!((syn - 779.9).abs() < 5.0, "synodic = {} days", syn);
}

#[test]
fn advance_orbit_one_period() {
    let mut oe = earth_oe();
    let period = oe.orbital_period(MU_SUN);
    let nu0 = oe.true_anomaly;
    oe.advance(period, MU_SUN);
    assert!(
        (oe.true_anomaly - nu0).abs() < 0.01,
        "Should return to same anomaly"
    );
}

#[test]
fn mean_motion_consistent_with_period() {
    let oe = earth_oe();
    let n = oe.mean_motion(MU_SUN);
    let t = oe.orbital_period(MU_SUN);
    assert!((n * t - 2.0 * PI).abs() < 1e-8);
}

#[test]
fn semi_latus_rectum_formula() {
    let oe = earth_oe();
    let p = oe.semi_latus_rectum();
    let expected = oe.semi_major_axis * (1.0 - oe.eccentricity * oe.eccentricity);
    assert!((p - expected).abs() < 1.0);
}

#[test]
fn radius_at_true_anomaly_periapsis_and_apoapsis() {
    let oe = earth_oe();
    let r_peri = oe.radius_at_true_anomaly(0.0);
    let r_apo = oe.radius_at_true_anomaly(PI);
    assert!((r_peri - oe.periapsis()).abs() < 1.0);
    assert!((r_apo - oe.apoapsis()).abs() < 1.0);
}

#[test]
fn velocity_at_radius_vis_viva() {
    let oe = earth_oe();
    let r = oe.periapsis();
    let v = oe.velocity_at_radius(r, MU_SUN);
    let v_expected = (MU_SUN * (2.0 / r - 1.0 / oe.semi_major_axis)).sqrt();
    assert!((v - v_expected).abs() / v_expected < 1e-8);
}

#[test]
fn eccentric_anomaly_at_periapsis() {
    let oe = earth_oe();
    let ea = oe.eccentric_anomaly();
    assert!(
        ea.abs() < 0.02,
        "E at periapsis should be near 0, got {}",
        ea
    );
}

#[test]
fn mean_anomaly_at_periapsis() {
    let oe = earth_oe();
    let ma = oe.mean_anomaly();
    assert!(
        ma.abs() < 0.02,
        "M at periapsis should be near 0, got {}",
        ma
    );
}

#[test]
fn true_anomaly_from_mean_roundtrip() {
    let ecc = 0.3;
    let original_ta = 1.5;
    let oe = OrbitalElements {
        semi_major_axis: AU,
        eccentricity: ecc,
        inclination: 0.0,
        longitude_ascending_node: 0.0,
        argument_periapsis: 0.0,
        true_anomaly: original_ta,
    };
    let ma = oe.mean_anomaly();
    let recovered_ta = OrbitalElements::true_anomaly_from_mean(ma, ecc);
    assert!(
        (recovered_ta - original_ta).abs() < 1e-8,
        "roundtrip: {} vs {}",
        recovered_ta,
        original_ta
    );
}

#[test]
fn from_body_state_uses_g_times_mass() {
    let r = AU;
    let v = (MU_SUN / r).sqrt();
    let pos = Vec3::new(r, 0.0, 0.0);
    let vel = Vec3::new(0.0, v, 0.0);
    let oe = OrbitalElements::from_body_state(pos, vel, SOLAR_MASS);
    assert!((oe.semi_major_axis - AU).abs() / AU < 1e-6);
}

#[test]
fn clone_orbital_elements() {
    let oe = earth_oe();
    let oe2 = oe.clone();
    assert!((oe.semi_major_axis - oe2.semi_major_axis).abs() < 1e-10);
    assert!((oe.eccentricity - oe2.eccentricity).abs() < 1e-10);
}

// ── Perturbations ──

#[test]
fn j2_raan_precession_sign() {
    let rate = j2_secular_precession_raan(1.08263e-3, 6.371e6, 7e6, 0.001, 0.5, MU_EARTH);
    assert!(rate < 0.0);
}

#[test]
fn j2_argp_precession_nonzero() {
    let rate = j2_secular_precession_argp(1.08263e-3, 6.371e6, 7e6, 0.001, 0.5, MU_EARTH);
    assert!(rate.abs() > 0.0);
}

#[test]
fn mercury_gr_precession_positive() {
    let a = 0.387 * AU;
    let e = 0.2056;
    let rate = mercury_gr_precession_rate(a, e, MU_SUN);
    assert!(rate > 0.0);
    let arcsec_per_century = rate * 365.25 * 86400.0 * 100.0 * 206265.0;
    assert!(
        (arcsec_per_century - 43.0).abs() < 2.0,
        "precession = {} ''/cy",
        arcsec_per_century
    );
}

#[test]
fn kozai_lidov_max_eccentricity_high_inc() {
    let e_max = kozai_lidov_max_eccentricity(80.0_f64.to_radians());
    assert!(e_max > 0.9);
}

#[test]
fn atmospheric_drag_always_negative() {
    let accel = atmospheric_drag_acceleration(1e-12, 7800.0, 2.2, 0.01);
    assert!(accel < 0.0, "Drag should decelerate");
}

#[test]
fn third_body_perturbation_nonzero() {
    let body_pos = Vec3::new(AU, 0.0, 0.0);
    let perturber_pos = Vec3::new(5.203 * AU, 0.0, 0.0);
    let acc = third_body_perturbation(body_pos, perturber_pos, 1.8982e27);
    assert!(acc.magnitude() > 0.0);
}

#[test]
fn yarkovsky_acceleration_positive() {
    let acc = yarkovsky_acceleration(0.01, 0.1, 1.0, 3000.0, AU, 3.828e26);
    assert!(acc > 0.0);
    assert!(acc < 1e-8, "Yarkovsky is very small: {}", acc);
}

#[test]
fn poynting_robertson_sign() {
    let acc = poynting_robertson_drag(3.828e26, 1e-3, 2000.0, AU, 1000.0);
    assert!(
        acc < 0.0,
        "PR drag should be negative for outward radial velocity"
    );
}

#[test]
fn lense_thirring_positive() {
    let rate = lense_thirring_precession(7e41, 7e6, MU_EARTH);
    assert!(rate > 0.0);
    assert!(rate < 1.0, "LT rate should be tiny: {}", rate);
}

#[test]
fn tidal_orbital_decay_negative() {
    let rate = tidal_orbital_decay_rate(7.342e22, 5.972e24, 6.371e6, 3.844e8, 12.0, 0.3);
    assert!(rate < 0.0, "Tidal decay should shrink orbit");
}

#[test]
fn solar_radiation_pressure_positive() {
    let accel = solar_radiation_pressure_accel(3.828e26, 0.01, 0.5, AU);
    assert!(accel > 0.0);
    assert!(accel < 1e-4);
}

// ── Transfers ──

#[test]
fn hohmann_transfer_dv_earth_to_mars() {
    let r1 = AU;
    let r2 = 1.524 * AU;
    let (dv1, dv2) = hohmann_transfer_dv(r1, r2, MU_SUN);
    assert!(dv1 > 2000.0 && dv1 < 4000.0, "dv1 = {}", dv1);
    assert!(dv2 > 1000.0 && dv2 < 3000.0, "dv2 = {}", dv2);
}

#[test]
fn hohmann_transfer_full_result() {
    let r1 = AU;
    let r2 = 1.524 * AU;
    let result = hohmann_transfer(r1, r2, MU_SUN);
    assert!(result.delta_v1 > 0.0);
    assert!(result.delta_v2 > 0.0);
    assert!((result.total_delta_v - result.delta_v1 - result.delta_v2).abs() < 1e-6);
    assert!(result.transfer_time > 0.0);
    let expected_sma = (r1 + r2) / 2.0;
    assert!((result.semi_major_axis - expected_sma).abs() / expected_sma < 1e-10);
}

#[test]
fn bielliptic_transfer_three_burns() {
    let r1 = AU;
    let r2 = 20.0 * AU;
    let r_int = 30.0 * AU;
    let result = bielliptic_transfer(r1, r2, r_int, MU_SUN);
    assert!(result.delta_v1 > 0.0);
    assert!(result.delta_v2 > 0.0);
    assert!(result.delta_v3 > 0.0);
    assert!(
        (result.total_delta_v - result.delta_v1 - result.delta_v2 - result.delta_v3).abs() < 1e-6
    );
    assert!(result.transfer_time > 0.0);
}

#[test]
fn hohmann_vs_bielliptic_threshold_positive() {
    let threshold = hohmann_vs_bielliptic_threshold(AU, MU_SUN);
    assert!(threshold > 0.0);
}

#[test]
fn launch_window_phase_angle_earth_mars() {
    let r1 = AU;
    let r2 = 1.524 * AU;
    let angle = launch_window_phase_angle(r1, r2, MU_SUN);
    assert!(
        angle.abs() < PI,
        "Phase angle should be reasonable: {}",
        angle
    );
}

#[test]
fn launch_window_next_returns_positive_time() {
    let t = launch_window_next(0.5, 0.8, 779.9 * 86400.0);
    assert!(t > 0.0);
    assert!(t < 779.9 * 86400.0);
}

#[test]
fn characteristic_energy_and_v_infinity() {
    let v_inf = 3000.0;
    let c3 = characteristic_energy(v_inf);
    assert!((c3 - v_inf * v_inf).abs() < 1e-6);
    let v_back = v_infinity_from_c3(c3);
    assert!((v_back - v_inf).abs() < 1e-6);
}

#[test]
fn hyperbolic_excess_velocity_departure() {
    let v_dep = 32000.0;
    let v_circ = 29780.0;
    let v_excess = hyperbolic_excess_velocity(v_dep, v_circ);
    let expected = (v_dep * v_dep - v_circ * v_circ).sqrt();
    assert!((v_excess - expected).abs() < 1e-6);
}

#[test]
fn gravity_assist_delta_v_positive() {
    let v_inf = 5000.0;
    let dv = gravity_assist_delta_v(v_inf, 5.972e24, 6.571e6);
    assert!(dv > 0.0);
    assert!(dv < 2.0 * v_inf);
}

#[test]
fn gravity_assist_max_deflection_bounded() {
    let deflection = gravity_assist_max_deflection(5000.0, 5.972e24, 6.571e6);
    assert!(deflection > 0.0);
    assert!(deflection < PI);
}

#[test]
fn oberth_effect_gains_energy() {
    let v_peri = 10000.0;
    let dv_applied = 1000.0;
    let effective_dv = oberth_effect_dv(v_peri, dv_applied);
    assert!(effective_dv > dv_applied, "Oberth effect should amplify dv");
}

#[test]
fn plane_change_dv_90_degrees() {
    let v = 7800.0;
    let dv = plane_change_dv(v, PI / 2.0);
    let expected = 2.0 * v * (PI / 4.0).sin();
    assert!((dv - expected).abs() < 1e-6);
}

#[test]
fn combined_plane_change_dv_zero_angle() {
    let v1 = 7000.0;
    let v2 = 8000.0;
    let dv = combined_plane_change_dv(v1, v2, 0.0);
    assert!(
        (dv - (v2 - v1).abs()).abs() < 1e-6,
        "0 angle = simple difference"
    );
}

#[test]
fn lambert_tof_estimate_positive() {
    let r1 = Vec3::new(AU, 0.0, 0.0);
    let r2 = Vec3::new(0.0, 1.524 * AU, 0.0);
    let tof = lambert_tof_estimate(r1, r2, MU_SUN, true);
    assert!(tof > 0.0);
}

#[test]
fn vis_viva_circular_orbit() {
    let v = vis_viva(MU_SUN, AU, AU);
    let v_circ = (MU_SUN / AU).sqrt();
    assert!((v - v_circ).abs() / v_circ < 1e-8);
}

#[test]
fn circular_velocity_earth() {
    let v = circular_velocity(MU_SUN, AU);
    assert!(v > 29000.0 && v < 30500.0, "v_circ = {}", v);
}

#[test]
fn escape_velocity_from_earth_orbit() {
    let v_esc = escape_velocity(MU_SUN, AU);
    let v_circ = circular_velocity(MU_SUN, AU);
    assert!((v_esc / v_circ - 2.0_f64.sqrt()).abs() < 1e-8);
}

#[test]
fn flight_path_angle_circular() {
    let r = AU;
    let v = (MU_SUN / r).sqrt();
    let h = r * v;
    let fpa = flight_path_angle(r, v, h);
    assert!(fpa.abs() < 1e-8, "FPA for circular = 0, got {}", fpa);
}

#[test]
fn delta_v_budget_sums_absolute() {
    let maneuvers = [2900.0, -1500.0, 500.0, -200.0];
    let budget = delta_v_budget(&maneuvers);
    let expected: f64 = maneuvers.iter().map(|dv| dv.abs()).sum();
    assert!((budget - expected).abs() < 1e-10);
}

#[test]
fn tsiolkovsky_mass_ratio_positive() {
    let mr = tsiolkovsky_mass_ratio(3000.0, 300.0);
    assert!(mr > 1.0, "Mass ratio must be > 1, got {}", mr);
    let mr2 = tsiolkovsky_mass_ratio(6000.0, 300.0);
    assert!(mr2 > mr, "Doubling dv should increase mass ratio");
}

#[test]
fn burn_time_proportional_to_mass() {
    let t1 = burn_time(3000.0, 100_000.0, 10_000.0);
    let t2 = burn_time(3000.0, 100_000.0, 20_000.0);
    assert!(
        (t2 / t1 - 2.0).abs() < 1e-10,
        "Should scale linearly with mass"
    );
    assert!(t1 > 0.0);
}