sagittariusas 0.0.3

use sagittariusas::config::parameters::*;
use sagittariusas::engine::accretion_disk::*;
use sagittariusas::engine::jet_dynamics::*;
use sagittariusas::engine::spacetime::*;
use std::f64::consts::PI;

const REL_TOL: f64 = 1e-3;

fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
    if b.abs() < 1e-30 {
        return a.abs() < tol;
    }
    ((a - b) / b).abs() < tol
}

#[test]
fn accretion_disk_inner_at_isco() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let rs = schwarzschild_radius(SGR_A_MASS);
    assert!(approx_eq(disk.inner_radius, 3.0 * rs, REL_TOL));
}

#[test]
fn temperature_zero_inside_isco() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let t = disk.temperature_profile(disk.inner_radius * 0.5);
    assert_eq!(t, 0.0);
}

#[test]
fn temperature_positive_outside_isco() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let t = disk.temperature_profile(disk.inner_radius * 2.0);
    assert!(t > 0.0);
}

#[test]
fn luminosity_positive() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    assert!(disk.luminosity() > 0.0);
}

#[test]
fn radiative_efficiency_bounded() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let eta = disk.radiative_efficiency();
    assert!(eta > 0.0 && eta < 1.0);
}

#[test]
fn radial_profile_correct_length() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let profile = disk.radial_profile(50);
    assert_eq!(profile.len(), 50);
}

#[test]
fn alpha_viscosity_clamped() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20).with_alpha(5.0);
    assert!(disk.alpha_viscosity <= 1.0);
}

#[test]
fn orbital_frequency_decreases_with_radius() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    let omega_inner = disk.orbital_frequency(disk.inner_radius * 2.0);
    let omega_outer = disk.orbital_frequency(disk.outer_radius * 0.5);
    assert!(omega_inner > omega_outer);
}

#[test]
fn total_mass_positive() {
    let disk = AccretionDisk::new(SGR_A_MASS, 1e20);
    assert!(disk.total_mass() > 0.0);
}

#[test]
fn jet_velocity_subluminal() {
    let jet = RelativisticJet::new(SGR_A_MASS, 0.9);
    assert!(jet.jet_velocity() < C);
}

#[test]
fn bz_power_positive() {
    let jet = RelativisticJet::new(SGR_A_MASS, 0.9).with_magnetic_field(100.0);
    assert!(jet.blandford_znajek_power() > 0.0);
}

#[test]
fn kinetic_luminosity_positive() {
    let jet = RelativisticJet::new(SGR_A_MASS, 0.9);
    assert!(jet.kinetic_luminosity() > 0.0);
}

#[test]
fn doppler_boosting() {
    let jet = RelativisticJet::new(SGR_A_MASS, 0.9);
    let d_toward = jet.doppler_factor(0.0);
    let d_away = jet.doppler_factor(PI);
    assert!(d_toward > 1.0);
    assert!(d_away < 1.0);
}

#[test]
fn apparent_superluminal_velocity() {
    let jet = RelativisticJet::new(SGR_A_MASS, 0.99);
    let v_app = jet.apparent_superluminal_velocity(PI / 6.0);
    assert!(v_app > 1.0);
}

#[test]
fn kerr_event_horizon_positive() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.9);
    assert!(ks.event_horizon() > 0.0);
}

#[test]
fn kerr_cauchy_less_than_event_horizon() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.5);
    assert!(ks.cauchy_horizon() < ks.event_horizon());
}

#[test]
fn kerr_ergosphere_ge_horizon_at_equator() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.9);
    let ergo = ks.ergosphere(PI / 2.0);
    assert!(ergo >= ks.event_horizon());
}

#[test]
fn kerr_metric_tensor_diagonal_dominant() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.5);
    let p = SpacetimePoint::new(0.0, 1e12, PI / 2.0, 0.0);
    let g = ks.metric_tensor(&p);
    assert!(g[0][0] < 0.0);
    assert!(g[1][1] > 0.0);
    assert!(g[2][2] > 0.0);
    assert!(g[3][3] > 0.0);
}

#[test]
fn photon_trace_terminates() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.5);
    let initial = SpacetimePoint::new(0.0, 1e12, PI / 2.0, 0.0);
    let traj = ks.trace_photon(initial, [C, 0.0, 0.01], 1000, 0.1);
    assert!(!traj.is_empty());
}

#[test]
fn photon_radial_infall_captured() {
    let ks = KerrSpacetime::new(SGR_A_MASS, SGR_A_SPIN);
    let rh = ks.event_horizon();
    let initial = SpacetimePoint::new(0.0, 3.0 * rh, PI / 2.0, 0.0);
    let traj = ks.trace_photon(initial, [-C, 0.0, 0.0], 100000, 0.001);
    let last = traj.last().unwrap();
    assert!(
        last.r < rh * 1.1,
        "photon should be captured: last r={}, rh={}, rh*1.1={}, traj.len={}",
        last.r,
        rh,
        rh * 1.1,
        traj.len()
    );
}

#[test]
fn schwarzschild_limit_horizons_degenerate() {
    let ks = KerrSpacetime::new(SGR_A_MASS, 0.0);
    let r_plus = ks.event_horizon();
    let r_minus = ks.cauchy_horizon();
    assert!(r_minus.abs() < 1e-6);
    assert!(r_plus > 0.0);
}