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

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

pub struct HawkingEvaporation {
    pub initial_mass: f64,
    pub current_mass: f64,
}

impl HawkingEvaporation {
    pub fn new(mass: f64) -> Self {
        Self {
            initial_mass: mass,
            current_mass: mass,
        }
    }

    pub fn temperature(&self) -> f64 {
        hawking_temperature(self.current_mass)
    }

    pub fn luminosity(&self) -> f64 {
        let t = self.temperature();
        let rs = schwarzschild_radius(self.current_mass);
        4.0 * PI * rs * rs * SIGMA_SB * t.powi(4)
    }

    pub fn evaporation_time(&self) -> f64 {
        evaporation_time(self.current_mass)
    }

    pub fn mass_loss_rate(&self) -> f64 {
        self.luminosity() / (C * C)
    }

    pub fn evolve(&mut self, dt: f64) {
        let dm = self.mass_loss_rate() * dt;
        self.current_mass = (self.current_mass - dm).max(0.0);
    }

    pub fn evolve_to_fraction(&mut self, fraction: f64) -> Vec<(f64, f64)> {
        let target = self.initial_mass * fraction.clamp(0.0, 1.0);
        let mut history = Vec::new();
        let mut t = 0.0;
        let dt_base = self.evaporation_time() * 1e-4;
        history.push((t, self.current_mass));

        while self.current_mass > target {
            let dt = dt_base.min(self.current_mass / self.mass_loss_rate().max(1e-100) * 0.01);
            self.evolve(dt);
            t += dt;
            history.push((t, self.current_mass));
            if self.current_mass <= 0.0 {
                break;
            }
        }

        history
    }

    pub fn entropy(&self) -> f64 {
        let rh = schwarzschild_radius(self.current_mass) / 2.0 + {
            let rg = gravitational_radius(self.current_mass);
            (rg * rg).sqrt()
        };
        let area = 4.0 * PI * rh * rh;
        K_B * C * C * C * area / (4.0 * HBAR * G)
    }

    pub fn page_time(&self) -> f64 {
        self.evaporation_time() / 3.0
    }

    pub fn scrambling_time(&self) -> f64 {
        let rh = schwarzschild_radius(self.current_mass);
        rh / C * (self.entropy() / K_B).ln().max(1.0)
    }

    pub fn peak_emission_wavelength(&self) -> f64 {
        2.898e-3 / self.temperature()
    }

    pub fn mass_fraction_remaining(&self) -> f64 {
        self.current_mass / self.initial_mass
    }
}