blackholesfactory 0.0.1

Black hole factory — create and simulate stellar, intermediate-mass, and supermassive black holes with full Kerr physics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

use crate::config::parameters::*;
use std::f64::consts::PI;

pub struct BlackHoleShadow {
    pub mass: f64,
    pub spin: f64,
    pub inclination: f64,
}

impl BlackHoleShadow {
    pub fn new(mass: f64, spin: f64, inclination: f64) -> Self {
        Self {
            mass,
            spin: spin.clamp(-1.0, 1.0),
            inclination,
        }
    }

    pub fn photon_ring_radius(&self) -> f64 {
        photon_sphere_radius(self.mass)
    }

    pub fn critical_impact_parameter(&self) -> f64 {
        let rph = self.photon_ring_radius();
        let rs = schwarzschild_radius(self.mass);
        rph / (1.0 - rs / rph).max(1e-30).sqrt()
    }

    pub fn shadow_radius(&self) -> f64 {
        let rg = gravitational_radius(self.mass);
        let a = self.spin;
        let base = 3.0 * 3.0_f64.sqrt() * rg;
        let spin_correction = 1.0 - 0.07 * a * a - 0.04 * a * a * self.inclination.cos().powi(2);
        base * spin_correction
    }

    pub fn shadow_diameter_uas(&self, distance: f64) -> f64 {
        let r_shadow = self.shadow_radius();
        let angle = 2.0 * r_shadow / distance;
        angle * 180.0 / PI * 3.6e9
    }

    pub fn asymmetry_parameter(&self) -> f64 {
        self.spin.abs() * self.inclination.sin()
    }

    pub fn shadow_area(&self) -> f64 {
        let r = self.shadow_radius();
        PI * r * r * (1.0 + 0.5 * self.asymmetry_parameter().powi(2))
    }

    pub fn shadow_contour(&self, n_points: usize) -> Vec<(f64, f64)> {
        let r = self.shadow_radius();
        let asym = self.asymmetry_parameter();
        (0..n_points)
            .map(|i| {
                let phi = 2.0 * PI * i as f64 / n_points as f64;
                let r_local = r * (1.0 + 0.1 * asym * phi.cos());
                (r_local * phi.cos(), r_local * phi.sin())
            })
            .collect()
    }
}

pub fn eht_m87_expected_diameter_uas() -> f64 {
    let mass = M87_MASS_SOLAR * SOLAR_MASS;
    let shadow = BlackHoleShadow::new(mass, M87_SPIN, 0.3);
    shadow.shadow_diameter_uas(M87_DISTANCE)
}