blackholesfactory 0.0.4

Black hole factory — create and simulate stellar, intermediate-mass, and supermassive black holes with full Kerr physics
Documentation
use crate::config::parameters::*;
use crate::physics::singularity::Singularity;

pub struct SupermassiveBlackHole {
    pub singularity: Singularity,
    pub host_sigma: f64,
}

impl SupermassiveBlackHole {
    pub fn new(mass_solar: f64, spin: f64, host_sigma: f64) -> Self {
        let clamped = mass_solar.clamp(SUPERMASSIVE_MASS_MIN_SOLAR, SUPERMASSIVE_MASS_MAX_SOLAR);
        Self {
            singularity: Singularity::new(clamped * SOLAR_MASS).with_spin(spin),
            host_sigma,
        }
    }

    pub fn from_m_sigma(sigma: f64) -> Self {
        let mass_solar = 1e8 * (sigma / 200e3).powf(4.38);
        let spin = 0.7;
        Self::new(mass_solar, spin, sigma)
    }

    pub fn mass_solar(&self) -> f64 {
        self.singularity.mass / SOLAR_MASS
    }

    pub fn m_sigma_predicted(&self) -> f64 {
        1e8 * (self.host_sigma / 200e3).powf(4.38) * SOLAR_MASS
    }

    pub fn m_sigma_deviation(&self) -> f64 {
        let predicted = self.m_sigma_predicted();
        (self.singularity.mass - predicted).abs() / predicted
    }

    pub fn eddington_luminosity(&self) -> f64 {
        eddington_luminosity(self.singularity.mass)
    }

    pub fn bondi_accretion_rate(&self, ambient_density: f64, sound_speed: f64) -> f64 {
        let r_b = 2.0 * G * self.singularity.mass / (sound_speed * sound_speed);
        4.0 * std::f64::consts::PI * ambient_density * r_b * r_b * sound_speed
    }

    pub fn jet_power_estimate(&self) -> f64 {
        let a = self.singularity.spin;
        1e38 * (a * a) * (self.singularity.mass / (1e9 * SOLAR_MASS)).powi(2)
    }

    pub fn tidal_disruption_radius(&self, star_mass: f64, star_radius: f64) -> f64 {
        star_radius * (self.singularity.mass / star_mass).cbrt()
    }

    pub fn hills_mass(&self, star_radius: f64, star_mass: f64) -> f64 {
        let rs = schwarzschild_radius(self.singularity.mass);
        star_mass * (rs / star_radius).powi(3).cbrt()
    }

    pub fn can_disrupt_star(&self, star_mass: f64, star_radius: f64) -> bool {
        self.tidal_disruption_radius(star_mass, star_radius) > self.singularity.event_horizon()
    }

    pub fn shadow_angular_diameter_uas(&self, distance: f64) -> f64 {
        let r_shadow = 5.2 * gravitational_radius(self.singularity.mass);
        let angle_rad = 2.0 * r_shadow / distance;
        angle_rad * 180.0 / std::f64::consts::PI * 3.6e9
    }

    pub fn dynamical_friction_timescale(
        &self,
        companion_mass: f64,
        separation: f64,
        sigma: f64,
    ) -> f64 {
        let ln_coulomb = (separation * sigma * sigma / (G * companion_mass))
            .ln()
            .max(1.0);
        sigma.powi(3) * separation
            / (4.0
                * std::f64::consts::PI
                * G
                * G
                * companion_mass
                * self.singularity.mass
                * ln_coulomb)
    }

    pub fn inspiral_timescale(&self, companion_mass: f64, separation: f64) -> f64 {
        let m1 = self.singularity.mass;
        5.0 / 256.0 * C.powi(5) * separation.powi(4)
            / (G.powi(3) * m1 * companion_mass * (m1 + companion_mass))
    }
}