satellitesfactory 0.0.3

use satellitesfactory::config::parameters::*;
use satellitesfactory::engine::evolution::*;
use satellitesfactory::engine::formation::*;
use satellitesfactory::engine::generator::*;
use satellitesfactory::engine::orbits::*;
use satellitesfactory::observables::astrometry::*;
use satellitesfactory::observables::photometry::*;
use satellitesfactory::observables::spectra::*;
use satellitesfactory::utils::helpers::*;

// ====================== Orbits ======================

#[test]
fn orbital_elements_from_km_deg() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 125.08, 318.15, 135.27);
    assert!((e.a - 384_400.0e3).abs() < 1.0);
    assert!((e.e - 0.054_9).abs() < 1e-8);
}

#[test]
fn orbital_elements_from_au_deg() {
    let e = OrbitalElements::from_au_deg(1.0, 0.017, 0.0, 0.0, 0.0, 0.0);
    assert!((e.a - AU).abs() < 1e6);
}

#[test]
fn orbital_elements_zero() {
    let z = OrbitalElements::zero();
    assert!(z.a == 0.0 && z.e == 0.0);
}

#[test]
fn period_moon() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    let p_days = e.period(G * (EARTH_MASS + MOON_MASS)) / 86400.0;
    assert!((p_days - 27.3).abs() < 0.5);
}

#[test]
fn mean_motion_positive() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    let n = e.mean_motion(G * EARTH_MASS);
    assert!(n > 0.0);
}

#[test]
fn periapsis_apoapsis() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    assert!(e.periapsis() < e.a);
    assert!(e.apoapsis() > e.a);
}

#[test]
fn semi_latus_rectum() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    let p = e.semi_latus_rectum();
    let expected = e.a * (1.0 - e.e * e.e);
    assert!((p - expected).abs() < 1.0);
}

#[test]
fn specific_energy_negative() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    assert!(e.specific_energy(G * EARTH_MASS) < 0.0);
}

#[test]
fn specific_angular_momentum_positive() {
    let e = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 0.0, 0.0, 0.0);
    assert!(e.specific_angular_momentum(G * EARTH_MASS) > 0.0);
}

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

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

#[test]
fn true_from_eccentric_at_zero() {
    let nu = true_from_eccentric(0.0, 0.5);
    assert!(nu.abs() < 1e-12);
}

#[test]
fn elements_to_state_radius() {
    let elem = OrbitalElements::from_km_deg(384_400.0, 0.054_9, 5.145, 125.08, 318.15, 135.27);
    let (pos, vel) = elements_to_state(&elem, EARTH_MASS + MOON_MASS);
    let r = (pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]).sqrt();
    assert!(r > 3.5e8 && r < 4.1e8);
    let v = (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]).sqrt();
    assert!(v > 900.0 && v < 1200.0);
}

#[test]
fn vis_viva_circular() {
    let mu = G * EARTH_MASS;
    let r = 384_400.0e3;
    let v = vis_viva(mu, r, r);
    let v_circ = (mu / r).sqrt();
    assert!((v - v_circ).abs() / v_circ < 1e-10);
}

#[test]
fn hohmann_positive() {
    let mu = G * EARTH_MASS;
    let (dv1, dv2) = hohmann_delta_v(mu, 6.8e6, 42.0e6);
    assert!(dv1 > 0.0 && dv2 > 0.0);
}

#[test]
fn laplace_resonance() {
    let (r12, r23) = laplace_resonance_check(1.769, 3.551, 7.155);
    assert!((r12 - 2.0).abs() < 0.02);
    assert!((r23 - 2.0).abs() < 0.02);
}

// ====================== Formation ======================

#[test]
fn formation_disk_mass() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    assert!(fm.disk_mass() > 0.0);
    assert!(fm.disk_mass() < JUPITER_MASS);
}

#[test]
fn formation_smaller_disk() {
    let jov = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    let sat = FormationModel::new(
        SATURN_MASS,
        SATURN_RADIUS,
        0.01,
        25.0 * SATURN_RADIUS,
        250.0,
    );
    assert!(sat.disk_mass() < jov.disk_mass());
}

#[test]
fn ice_line_outside_planet() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    assert!(fm.ice_line_radius() > fm.parent_radius);
}

#[test]
fn isolation_mass_positive() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    let iso = fm.isolation_mass(10.0 * JUPITER_RADIUS);
    assert!(iso > 0.0);
}

#[test]
fn maximum_satellite_mass_positive() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    assert!(fm.maximum_satellite_mass() > 0.0);
}

#[test]
fn disk_temperature_decreasing() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    let t1 = fm.disk_temperature(5.0 * JUPITER_RADIUS);
    let t2 = fm.disk_temperature(20.0 * JUPITER_RADIUS);
    assert!(t1 > t2);
}

#[test]
fn disk_surface_density_decreasing() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    let s1 = fm.disk_surface_density(5.0 * JUPITER_RADIUS);
    let s2 = fm.disk_surface_density(20.0 * JUPITER_RADIUS);
    assert!(s1 > s2);
}

#[test]
fn formation_timescale_positive() {
    let fm = FormationModel::new(
        JUPITER_MASS,
        JUPITER_RADIUS,
        0.02,
        30.0 * JUPITER_RADIUS,
        400.0,
    );
    let t = fm.formation_timescale(10.0 * JUPITER_RADIUS);
    assert!(t > 0.0 && t.is_finite());
}

#[test]
fn capture_probability_zero_fast() {
    let p = FormationModel::capture_probability(2000.0, 1000.0, 0.5);
    assert!(p == 0.0);
}

#[test]
fn capture_probability_bounded() {
    let p = FormationModel::capture_probability(200.0, 1000.0, 0.5);
    assert!(p > 0.0 && p <= 1.0);
}

#[test]
fn impact_ejecta_zero_slow() {
    let f = FormationModel::impact_ejecta_fraction(1000.0, 2400.0);
    assert!(f == 0.0);
}

#[test]
fn ring_spreading_timescale_positive() {
    let t = FormationModel::ring_spreading_timescale(100_000.0e3, 1e6);
    assert!(t > 0.0);
}

// ====================== Evolution ======================

#[test]
fn moon_tidal_migration() {
    let evo = SatelliteEvolution::new(
        MOON_MASS,
        MOON_RADIUS,
        384_400.0e3,
        0.054_9,
        EARTH_MASS,
        EARTH_RADIUS,
    );
    let rate = evo.tidal_migration_rate(27.0);
    assert!(rate > 0.0);
}

#[test]
fn io_orbital_decay() {
    let evo = SatelliteEvolution::new(
        IO_MASS,
        IO_RADIUS,
        421_700.0e3,
        0.004_1,
        JUPITER_MASS,
        JUPITER_RADIUS,
    );
    let t = evo.orbital_decay_timescale(100.0);
    assert!(t > 0.0);
}

#[test]
fn triton_evolution() {
    let evo = SatelliteEvolution::new(
        TRITON_MASS,
        TRITON_RADIUS,
        354_759.0e3,
        0.000_016,
        NEPTUNE_MASS,
        NEPTUNE_RADIUS,
    );
    let t = evo.orbital_decay_timescale(12_000.0);
    assert!(t > 0.0);
}

#[test]
fn circularization_timescale() {
    let evo = SatelliteEvolution::new(
        IO_MASS,
        IO_RADIUS,
        421_700.0e3,
        0.004_1,
        JUPITER_MASS,
        JUPITER_RADIUS,
    );
    let t = evo.circularization_timescale(100.0, 0.015);
    assert!(t > 0.0 && t.is_finite());
}

#[test]
fn mean_motion_resonance_ratio() {
    let evo = SatelliteEvolution::new(
        IO_MASS,
        IO_RADIUS,
        421_700.0e3,
        0.004_1,
        JUPITER_MASS,
        JUPITER_RADIUS,
    );
    let ratio = evo.mean_motion_resonance_period_ratio(671_100.0e3);
    assert!(ratio > 0.0 && ratio.is_finite());
}

#[test]
fn sputtering_rate_positive() {
    let evo = SatelliteEvolution::new(
        IO_MASS,
        IO_RADIUS,
        421_700.0e3,
        0.004_1,
        JUPITER_MASS,
        JUPITER_RADIUS,
    );
    let rate = evo.sputtering_rate(1e12);
    assert!(rate > 0.0);
}

#[test]
fn secular_acceleration_finite() {
    let a = secular_acceleration(1e-10, 384_400.0e3);
    assert!(a.is_finite());
}

#[test]
fn kozai_lidov_bounded() {
    let emax = kozai_lidov_max_eccentricity(0.7);
    assert!((0.0..1.0).contains(&emax));
}

#[test]
fn yarkovsky_drift_positive() {
    let drift = yarkovsky_drift(5_000.0, 1500.0, 3.828e26, AU);
    assert!(drift > 0.0);
}

// ====================== Generator ======================

#[test]
fn generate_system_basic() {
    let config = SystemGeneratorConfig::new(
        "Jupiter",
        JUPITER_MASS,
        JUPITER_RADIUS,
        2.0 * JUPITER_RADIUS,
        30.0 * JUPITER_RADIUS,
    )
    .with_counts(4, 2, 3, 0, 0);
    let system = generate_system(&config);
    assert_eq!(system.len(), 9);
}

#[test]
fn generate_system_all_types() {
    let config = SystemGeneratorConfig::new(
        "Saturn",
        SATURN_MASS,
        SATURN_RADIUS,
        1.5 * SATURN_RADIUS,
        25.0 * SATURN_RADIUS,
    )
    .with_counts(2, 1, 1, 1, 1)
    .with_exotic(1, 1, 1);
    let system = generate_system(&config);
    assert_eq!(system.len(), 9);
}

#[test]
fn generated_types_correct() {
    let config = SystemGeneratorConfig::new(
        "Jupiter",
        JUPITER_MASS,
        JUPITER_RADIUS,
        2.0 * JUPITER_RADIUS,
        30.0 * JUPITER_RADIUS,
    )
    .with_counts(2, 1, 1, 1, 1)
    .with_exotic(1, 1, 1);
    let system = generate_system(&config);
    let mut reg = 0;
    let mut cap = 0;
    let mut ring = 0;
    let mut volc = 0;
    let mut icy = 0;
    let mut ocean = 0;
    let mut exo = 0;
    let mut co = 0;
    for s in &system {
        match s {
            GeneratedSatellite::Regular(_) => reg += 1,
            GeneratedSatellite::Captured(_) => cap += 1,
            GeneratedSatellite::Ring(_) => ring += 1,
            GeneratedSatellite::Volcanic(_) => volc += 1,
            GeneratedSatellite::Icy(_) => icy += 1,
            GeneratedSatellite::SubsurfaceOcean(_) => ocean += 1,
            GeneratedSatellite::Exo(_) => exo += 1,
            GeneratedSatellite::CoOrbital(_) => co += 1,
        }
    }
    assert_eq!(reg, 2);
    assert_eq!(cap, 1);
    assert_eq!(ring, 1);
    assert_eq!(volc, 1);
    assert_eq!(icy, 1);
    assert_eq!(ocean, 1);
    assert_eq!(exo, 1);
    assert_eq!(co, 1);
}

#[test]
fn generate_empty_system() {
    let config = SystemGeneratorConfig::new("Test", 1e27, 7e7, 1e8, 1e10);
    let system = generate_system(&config);
    assert!(system.is_empty());
}

#[test]
fn generate_with_seed() {
    let c1 = SystemGeneratorConfig::new("A", JUPITER_MASS, JUPITER_RADIUS, 1e8, 1e10)
        .with_counts(3, 0, 0, 0, 0)
        .with_seed(99);
    let c2 = SystemGeneratorConfig::new("A", JUPITER_MASS, JUPITER_RADIUS, 1e8, 1e10)
        .with_counts(3, 0, 0, 0, 0)
        .with_seed(99);
    let s1 = generate_system(&c1);
    let s2 = generate_system(&c2);
    assert_eq!(s1.len(), s2.len());
}

// ====================== Photometry ======================

#[test]
fn equilibrium_temperature_jupiter() {
    let t = equilibrium_temperature(3.828e26, 5.2 * AU, 0.34);
    assert!(t > 100.0 && t < 200.0);
}

#[test]
fn transit_depth_small() {
    let td = transit_depth_satellite(EUROPA_RADIUS, 6.957e8);
    assert!(td > 0.0 && td < 1e-4);
}

#[test]
fn geometric_albedo_identity() {
    let ga = geometric_albedo(0.5, 1.0);
    assert!((ga - 0.5).abs() < 1e-12);
}

#[test]
fn opposition_surge_at_zero() {
    let s = opposition_surge(0.0, 0.5);
    assert!((s - 1.5).abs() < 1e-10);
}

#[test]
fn thermal_emission_positive() {
    let f = thermal_emission_flux(100.0, EUROPA_RADIUS, AU);
    assert!(f > 0.0);
}

// ====================== Spectra ======================

#[test]
fn band_depth_thirty_percent() {
    let bd = band_depth(1.0, 0.7);
    assert!((bd - 0.3).abs() < 0.01);
}

#[test]
fn water_ice_index_positive() {
    let wii = water_ice_index(0.8, 0.4);
    assert!(wii > 1.5);
}

#[test]
fn composition_water_ice() {
    let c = composition_from_albedo_and_density(0.99, 1200.0);
    assert_eq!(c, SurfaceComposition::WaterIce);
}

#[test]
fn composition_carbonaceous() {
    let c = composition_from_albedo_and_density(0.05, 2000.0);
    assert_eq!(c, SurfaceComposition::Carbonaceous);
}

// ====================== Astrometry ======================

#[test]
fn angular_separation_positive() {
    let ang = angular_separation_arcsec(10.0 * 3.086e16, 5.2 * AU);
    assert!(ang > 0.0);
}

#[test]
fn occultation_duration_positive() {
    let occ = occultation_duration(EUROPA_RADIUS, 13_740.0, JUPITER_RADIUS, 0.0);
    assert!(occ > 0.0);
}

#[test]
fn light_travel_time_au() {
    let ltt = light_travel_time(AU);
    assert!((ltt - 499.0).abs() < 1.0);
}

#[test]
fn astrometric_wobble_positive() {
    let w = astrometric_wobble(GANYMEDE_MASS, JUPITER_MASS, 1_070_400.0e3, 10.0 * 3.086e16);
    assert!(w > 0.0);
}

#[test]
fn mutual_event_close_orbit() {
    let prob = mutual_event_probability(IO_RADIUS, EUROPA_RADIUS, 250_000.0e3, 0.001);
    assert!(prob == 1.0);
}

// ====================== Utils ======================

#[test]
fn lerp_midpoint() {
    assert!((lerp(0.0, 10.0, 0.5) - 5.0).abs() < 1e-12);
}

#[test]
fn lerp_endpoints() {
    assert!(lerp(0.0, 10.0, 0.0).abs() < 1e-12);
    assert!((lerp(0.0, 10.0, 1.0) - 10.0).abs() < 1e-12);
}

#[test]
fn simpson_integrate_x_squared() {
    let result = simpson_integrate(|x| x * x, 0.0, 1.0, 100);
    assert!((result - 1.0 / 3.0).abs() < 1e-6);
}

#[test]
fn log_range_endpoints() {
    let r = log_range(1.0, 1000.0, 4);
    assert_eq!(r.len(), 4);
    assert!((r[0] - 1.0).abs() < 0.01);
    assert!((r[3] - 1000.0).abs() < 1.0);
}

#[test]
fn format_si_giga() {
    let s = format_si(1.5e9);
    assert!(s.contains("G"));
}

#[test]
fn clamp_within() {
    assert!((clamp(0.5, 0.0, 1.0) - 0.5).abs() < 1e-12);
}

#[test]
fn clamp_below() {
    assert!(clamp(-1.0, 0.0, 1.0).abs() < 1e-12);
}

#[test]
fn clamp_above() {
    assert!((clamp(2.0, 0.0, 1.0) - 1.0).abs() < 1e-12);
}