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
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

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

#[derive(Clone, Copy)]
pub enum DiskModel {
    ShakuraSunyaev,
    Adaf,
    Slim,
}

pub struct AccretionDisk {
    pub mass: f64,
    pub mdot: f64,
    pub spin: f64,
    pub inner_radius: f64,
    pub outer_radius: f64,
    pub alpha: f64,
    pub model: DiskModel,
}

impl AccretionDisk {
    pub fn new(mass: f64, mdot: f64) -> Self {
        let r_isco = isco_radius(mass, 0.0);
        Self {
            mass,
            mdot,
            spin: 0.0,
            inner_radius: r_isco,
            outer_radius: r_isco * 1e4,
            alpha: 0.1,
            model: DiskModel::ShakuraSunyaev,
        }
    }

    pub fn with_spin(mass: f64, mdot: f64, spin: f64) -> Self {
        let r_isco = isco_radius(mass, spin);
        Self {
            mass,
            mdot,
            spin,
            inner_radius: r_isco,
            outer_radius: r_isco * 1e4,
            alpha: 0.1,
            model: DiskModel::ShakuraSunyaev,
        }
    }

    pub fn with_model(mut self, model: DiskModel) -> Self {
        self.model = model;
        self
    }

    pub fn with_alpha(mut self, alpha: f64) -> Self {
        self.alpha = alpha.clamp(0.001, 1.0);
        self
    }

    pub fn with_outer_radius(mut self, r_out: f64) -> Self {
        self.outer_radius = r_out.max(self.inner_radius * 2.0);
        self
    }

    pub fn temperature_profile(&self, r: f64) -> f64 {
        if r < self.inner_radius {
            return 0.0;
        }
        let factor = 3.0 * G * self.mass * self.mdot / (8.0 * PI * SIGMA_SB * r.powi(3))
            * (1.0 - (self.inner_radius / r).sqrt());
        let base = factor.max(0.0).powf(0.25);
        match self.model {
            DiskModel::ShakuraSunyaev => base,
            DiskModel::Adaf => base * 0.3 * (r / self.inner_radius).powf(-0.5),
            DiskModel::Slim => {
                let edd_ratio = self.mdot * C * C * 0.1 / eddington_luminosity(self.mass);
                base * (1.0 + edd_ratio).ln().max(1.0).powf(0.25)
            }
        }
    }

    pub fn radial_profile(&self, n_points: usize) -> Vec<(f64, f64)> {
        let log_min = self.inner_radius.ln();
        let log_max = self.outer_radius.ln();
        (0..n_points)
            .map(|i| {
                let r =
                    (log_min + (log_max - log_min) * i as f64 / (n_points - 1).max(1) as f64).exp();
                (r, self.temperature_profile(r))
            })
            .collect()
    }

    pub fn luminosity(&self) -> f64 {
        let n = 500;
        let log_min = self.inner_radius.ln();
        let log_max = self.outer_radius.ln();
        let d_log = (log_max - log_min) / n as f64;
        (0..n)
            .map(|i| {
                let r = (log_min + (i as f64 + 0.5) * d_log).exp();
                let t = self.temperature_profile(r);
                2.0 * PI * r * SIGMA_SB * t.powi(4) * 2.0 * r * d_log
            })
            .sum()
    }

    pub fn radiative_efficiency(&self) -> f64 {
        let total = self.mdot * C * C;
        if total < 1e-30 {
            return 0.0;
        }
        (self.luminosity() / total).clamp(0.0, 0.42)
    }

    pub fn peak_temperature(&self) -> f64 {
        let r_peak = self.inner_radius * (49.0 / 36.0);
        self.temperature_profile(r_peak)
    }

    pub fn viscous_timescale(&self, r: f64) -> f64 {
        let h_over_r = 0.01;
        let omega = (G * self.mass / (r * r * r)).sqrt();
        1.0 / (self.alpha * h_over_r * h_over_r * omega)
    }

    pub fn total_mass(&self) -> f64 {
        let sigma_0 = self.mdot / (3.0 * PI * self.alpha * 0.01 * 0.01);
        let n = 200;
        let log_min = self.inner_radius.ln();
        let log_max = self.outer_radius.ln();
        let d_log = (log_max - log_min) / n as f64;
        (0..n)
            .map(|i| {
                let r = (log_min + (i as f64 + 0.5) * d_log).exp();
                let omega = (G * self.mass / (r * r * r)).sqrt();
                let sigma = sigma_0 / omega;
                2.0 * PI * r * sigma * r * d_log
            })
            .sum()
    }
}