blackholesfactory 0.0.4

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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

use crate::config::parameters::{C, G, SOLAR_MASS};
use std::f64::consts::PI;

#[derive(Clone, Copy)]
pub enum JetLaunchMechanism {
    BlandfordZnajek,
    BlandfordPayne,
}

pub struct RelativisticJet {
    pub mass: f64,
    pub spin: f64,
    pub magnetic_field: f64,
    pub bulk_lorentz_factor: f64,
    pub opening_angle: f64,
    pub mechanism: JetLaunchMechanism,
}

impl RelativisticJet {
    pub fn new(mass: f64, spin: f64) -> Self {
        Self {
            mass,
            spin: spin.clamp(-1.0, 1.0),
            magnetic_field: 100.0,
            bulk_lorentz_factor: 5.0,
            opening_angle: 0.1,
            mechanism: JetLaunchMechanism::BlandfordZnajek,
        }
    }

    pub fn with_magnetic_field(mut self, b: f64) -> Self {
        self.magnetic_field = b.max(1e-10);
        self
    }

    pub fn with_lorentz_factor(mut self, gamma: f64) -> Self {
        self.bulk_lorentz_factor = gamma.max(1.0);
        self
    }

    pub fn with_opening_angle(mut self, angle: f64) -> Self {
        self.opening_angle = angle.clamp(0.001, PI / 4.0);
        self
    }

    pub fn with_mechanism(mut self, mech: JetLaunchMechanism) -> Self {
        self.mechanism = mech;
        self
    }

    pub fn blandford_znajek_power(&self) -> f64 {
        let rg = G * self.mass / (C * C);
        let rh = rg + (rg * rg - (self.spin * rg).powi(2)).max(0.0).sqrt();
        let omega_h = self.spin * C / (2.0 * rh);
        let b2 = self.magnetic_field * self.magnetic_field;
        let area = 4.0 * PI * rh * rh;
        b2 * area * rh * rh * omega_h * omega_h / (6.0 * C)
    }

    pub fn kinetic_luminosity(&self) -> f64 {
        let bz = self.blandford_znajek_power();
        let eta = 1.0 - 1.0 / self.bulk_lorentz_factor;
        bz * eta
    }

    pub fn doppler_factor(&self, viewing_angle: f64) -> f64 {
        let beta = (1.0 - 1.0 / (self.bulk_lorentz_factor * self.bulk_lorentz_factor))
            .max(0.0)
            .sqrt();
        1.0 / (self.bulk_lorentz_factor * (1.0 - beta * viewing_angle.cos()))
    }

    pub fn apparent_velocity(&self, viewing_angle: f64) -> f64 {
        let beta = (1.0 - 1.0 / (self.bulk_lorentz_factor * self.bulk_lorentz_factor))
            .max(0.0)
            .sqrt();
        beta * C * viewing_angle.sin() / (1.0 - beta * viewing_angle.cos())
    }

    pub fn max_apparent_velocity(&self) -> f64 {
        let beta = (1.0 - 1.0 / (self.bulk_lorentz_factor * self.bulk_lorentz_factor))
            .max(0.0)
            .sqrt();
        let gamma = self.bulk_lorentz_factor;
        beta * gamma * C
    }

    pub fn synchrotron_cooling_time(&self, electron_gamma: f64) -> f64 {
        let sigma_t = 6.652e-29;
        let me = 9.109e-31;
        let u_b = self.magnetic_field * self.magnetic_field / (2.0 * 4.0 * PI * 1e-7);
        6.0 * PI * me * C / (sigma_t * electron_gamma * u_b)
    }

    pub fn magnetic_field_at(&self, distance: f64) -> f64 {
        let r0 = G * self.mass / (C * C);
        self.magnetic_field * (r0 / distance).max(0.0)
    }

    pub fn beam_radius_at(&self, distance: f64) -> f64 {
        distance * self.opening_angle.tan()
    }

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