blackholesfactory 0.0.1

use std::f64::consts::FRAC_PI_2;
use std::process;

use blackholesfactory::config::parameters::*;
use blackholesfactory::engine::accretion_disk::AccretionDisk;
use blackholesfactory::engine::evaporation::HawkingEvaporation;
use blackholesfactory::engine::jet_dynamics::RelativisticJet;
use blackholesfactory::engine::spacetime::{KerrSpacetime, SpacetimePoint};
use blackholesfactory::observables::lensing::GravitationalLens;
use blackholesfactory::observables::radiation::AccretionRadiation;
use blackholesfactory::observables::shadows::{BlackHoleShadow, eht_m87_expected_diameter_uas};
use blackholesfactory::physics::gravitation::GravitationalField;
use blackholesfactory::physics::relativity::{SchwarzschildMetric, gravitational_wave_power};
use blackholesfactory::types::intermediate_mass::IntermediateMassBlackHole;
use blackholesfactory::types::stellar_mass::StellarMassBlackHole;
use blackholesfactory::types::supermassive::SupermassiveBlackHole;

struct FactoryValidation {
    stellar_mass_solar: f64,
    stellar_strain: f64,
    intermediate_mass_solar: f64,
    imbh_qpo: f64,
    imbh_recoil: f64,
    supermassive_mass_solar: f64,
    m_sigma_dev: f64,
    jet_power_est: f64,
    can_disrupt: bool,
    shadow_deviation: f64,
    predicted_uas: f64,
    model_uas: f64,
    contour_area: f64,
    cygnus_precession: f64,
    disk_lum: f64,
    disk_eff: f64,
    peak_temp: f64,
    visc: f64,
    edd_ratio: f64,
    inner_temp: f64,
    bz_power: f64,
    kinetic_power: f64,
    max_apparent: f64,
    event_horizon: f64,
    final_photon_r: f64,
    tidal: f64,
    bondi: f64,
    einstein: f64,
    shapiro: f64,
    gw_power: f64,
    evap_time: f64,
    page_time: f64,
    peak_wavelength: f64,
    mass_remaining: f64,
    inspiral: f64,
    friction: f64,
    combined_power: f64,
    ratio_to_edd: f64,
    photon_captured: bool,
}

impl FactoryValidation {
    fn passed(&self) -> bool {
        let stellar_ok = self.stellar_mass_solar >= STELLAR_MASS_MIN_SOLAR
            && self.stellar_mass_solar <= STELLAR_MASS_MAX_SOLAR
            && self.stellar_strain > 0.0;

        let imbh_ok = self.intermediate_mass_solar >= INTERMEDIATE_MASS_MIN_SOLAR
            && self.intermediate_mass_solar <= INTERMEDIATE_MASS_MAX_SOLAR
            && self.imbh_qpo > 0.0
            && self.imbh_recoil > 0.0;

        let smbh_ok = self.supermassive_mass_solar >= SUPERMASSIVE_MASS_MIN_SOLAR
            && self.supermassive_mass_solar <= SUPERMASSIVE_MASS_MAX_SOLAR
            && self.m_sigma_dev.abs() < 1.0
            && self.jet_power_est > 0.0
            && self.can_disrupt;

        let shadow_ok = self.shadow_deviation < 0.3
            && self.predicted_uas > 0.0
            && self.model_uas > 0.0
            && self.contour_area > 0.0;

        let disk_ok = self.disk_lum > 0.0
            && self.disk_eff > 0.0
            && self.disk_eff < 0.5
            && self.peak_temp > 0.0
            && self.visc > 0.0
            && self.edd_ratio > 0.0
            && self.inner_temp > 0.0;

        let jet_ok = self.bz_power > 0.0 && self.kinetic_power > 0.0 && self.max_apparent > C;

        let spacetime_ok = self.event_horizon > 0.0 && self.final_photon_r > 0.0;

        let gravity_ok = self.cygnus_precession > 0.0 && self.tidal > 0.0 && self.bondi > 0.0;

        let lens_ok = self.einstein > 0.0 && self.shapiro > 0.0;

        let evap_ok = self.evap_time > 0.0
            && self.evap_time < f64::INFINITY
            && self.page_time > 0.0
            && self.peak_wavelength > 0.0
            && self.mass_remaining > 0.0
            && self.mass_remaining < 1.0;

        let multi_ok = self.inspiral > 0.0 && self.friction > 0.0;

        let combined_ok = self.combined_power > 0.0
            && self.ratio_to_edd > 0.0
            && self.gw_power > 0.0
            && self.photon_captured;

        stellar_ok
            && imbh_ok
            && smbh_ok
            && shadow_ok
            && disk_ok
            && jet_ok
            && spacetime_ok
            && gravity_ok
            && lens_ok
            && evap_ok
            && multi_ok
            && combined_ok
    }
}

fn main() {
    let stellar = StellarMassBlackHole::cygnus_x1();
    let stellar_mass_solar = stellar.mass_solar();
    let stellar_strain =
        stellar.gravitational_wave_strain(10.0 * SOLAR_MASS, CYGNUS_X1_DISTANCE, 1e10);

    let imbh = IntermediateMassBlackHole::hlx1();
    let imbh_qpo = imbh.ulf_qpo_frequency();
    let imbh_recoil = imbh.merger_recoil_velocity(0.5);

    let smbh = SupermassiveBlackHole::m87();
    let m_sigma_dev = smbh.m_sigma_deviation();
    let jet_power_est = smbh.jet_power_estimate();
    let can_disrupt = smbh.can_disrupt_star(SOLAR_MASS, 6.957e8);

    let shadow = BlackHoleShadow::new(M87_MASS_SOLAR * SOLAR_MASS, M87_SPIN, 0.3);
    let predicted_uas = shadow.shadow_diameter_uas(M87_DISTANCE);
    let model_uas = eht_m87_expected_diameter_uas();
    let shadow_deviation = (predicted_uas - M87_SHADOW_UAS).abs() / M87_SHADOW_UAS;
    let contour = shadow.shadow_contour(64);
    let contour_area: f64 = contour
        .windows(2)
        .map(|w| w[0].0 * w[1].1 - w[1].0 * w[0].1)
        .sum::<f64>()
        .abs()
        * 0.5;

    let schw = SchwarzschildMetric::new(CYGNUS_X1_MASS_SOLAR * SOLAR_MASS);
    let cygnus_precession = schw.perihelion_precession(0.2 * PARSEC, 0.0);

    let mdot_m87 = 0.1 * SOLAR_MASS / (365.25 * 86400.0);
    let disk = AccretionDisk::with_spin(M87_MASS_SOLAR * SOLAR_MASS, mdot_m87, M87_SPIN);
    let disk_lum = disk.luminosity();
    let disk_eff = disk.radiative_efficiency();
    let peak_temp = disk.peak_temperature();
    let visc = disk.viscous_timescale(disk.inner_radius * 10.0);

    let accr = AccretionRadiation::new(M87_MASS_SOLAR * SOLAR_MASS, mdot_m87);
    let edd_ratio = accr.eddington_ratio();
    let inner_temp = accr.disk_inner_temperature(M87_SPIN);

    let jet = RelativisticJet::new(M87_MASS_SOLAR * SOLAR_MASS, M87_SPIN)
        .with_lorentz_factor(6.0)
        .with_magnetic_field(50.0)
        .with_opening_angle(0.05);
    let bz = jet.blandford_znajek_power();
    let kin = jet.kinetic_luminosity();
    let max_app = jet.max_apparent_velocity();

    let ks = KerrSpacetime::new(M87_MASS_SOLAR * SOLAR_MASS, M87_SPIN);
    let rh = ks.event_horizon();
    let traj = ks.trace_photon(
        SpacetimePoint::new(0.0, 50.0 * rh, FRAC_PI_2, 0.0),
        [C, 0.0, C / (50.0 * rh)],
        500,
        0.01,
    );
    let final_r = traj.last().map(|p| p.r).unwrap_or(0.0);

    let gf = GravitationalField::new(M87_MASS_SOLAR * SOLAR_MASS);
    let tidal = gf.tidal_force(rh, 1.0);
    let bondi = gf.bondi_radius(1e6);

    let lens = GravitationalLens::new(
        M87_MASS_SOLAR * SOLAR_MASS,
        M87_DISTANCE,
        2.0 * M87_DISTANCE,
    );
    let einstein = lens.einstein_radius();
    let shapiro = lens.shapiro_delay(1e15);

    let gw = gravitational_wave_power(stellar.singularity.mass, 10.0 * SOLAR_MASS, 1e10);

    let mut evap = HawkingEvaporation::new(1e12);
    let evap_time = evap.evaporation_time();
    let page = evap.page_time();
    let peak_wl = evap.peak_emission_wavelength();
    evap.evolve(1e10);
    let remaining = evap.mass_fraction_remaining();

    let inspiral = smbh.inspiral_timescale(imbh.singularity.mass, 0.1 * PARSEC);
    let friction =
        smbh.dynamical_friction_timescale(imbh.singularity.mass, PARSEC, smbh.host_sigma);

    let combined_power = disk_lum + bz + kin;

    let v = FactoryValidation {
        stellar_mass_solar,
        stellar_strain,
        intermediate_mass_solar: imbh.mass_solar(),
        imbh_qpo,
        imbh_recoil,
        supermassive_mass_solar: smbh.mass_solar(),
        m_sigma_dev,
        jet_power_est,
        can_disrupt,
        shadow_deviation,
        predicted_uas,
        model_uas,
        contour_area,
        cygnus_precession,
        disk_lum,
        disk_eff,
        peak_temp,
        visc,
        edd_ratio,
        inner_temp,
        bz_power: bz,
        kinetic_power: kin,
        max_apparent: max_app,
        event_horizon: rh,
        final_photon_r: final_r,
        tidal,
        bondi,
        einstein,
        shapiro,
        gw_power: gw,
        evap_time,
        page_time: page,
        peak_wavelength: peak_wl,
        mass_remaining: remaining,
        inspiral,
        friction,
        combined_power,
        ratio_to_edd: combined_power / accr.eddington_luminosity(),
        photon_captured: final_r < rh * 1.1,
    };

    process::exit(if v.passed() { 0 } else { 1 });
}