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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

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

pub struct SchwarzschildMetric {
    pub mass: f64,
    rs: f64,
}

impl SchwarzschildMetric {
    pub fn new(mass: f64) -> Self {
        Self {
            mass,
            rs: 2.0 * G * mass / (C * C),
        }
    }

    pub fn g_tt(&self, r: f64) -> f64 {
        -(1.0 - self.rs / r)
    }

    pub fn g_rr(&self, r: f64) -> f64 {
        1.0 / (1.0 - self.rs / r)
    }

    pub fn proper_time_ratio(&self, r: f64) -> f64 {
        (1.0 - self.rs / r).max(0.0).sqrt()
    }

    pub fn proper_distance(&self, r1: f64, r2: f64) -> f64 {
        let n = 1000;
        let dr = (r2 - r1) / n as f64;
        (0..n)
            .map(|i| {
                let r = r1 + (i as f64 + 0.5) * dr;
                (1.0 / (1.0 - self.rs / r)).max(0.0).sqrt() * dr
            })
            .sum()
    }

    pub fn radial_freefall_velocity(&self, r: f64) -> f64 {
        (self.rs / r).min(1.0).sqrt() * C
    }

    pub fn perihelion_precession(&self, semi_major: f64, eccentricity: f64) -> f64 {
        let p = semi_major * (1.0 - eccentricity * eccentricity);
        6.0 * PI * G * self.mass / (C * C * p)
    }

    pub fn gravitational_redshift(&self, r_emit: f64, r_obs: f64) -> f64 {
        let z_emit = (1.0 - self.rs / r_emit).max(1e-30).sqrt();
        let z_obs = (1.0 - self.rs / r_obs).max(1e-30).sqrt();
        z_obs / z_emit - 1.0
    }
}

pub struct KerrMetric {
    pub mass: f64,
    pub spin: f64,
    rg: f64,
    a: f64,
}

impl KerrMetric {
    pub fn new(mass: f64, spin: f64) -> Self {
        let rg = G * mass / (C * C);
        Self {
            mass,
            spin,
            rg,
            a: spin * rg,
        }
    }

    fn sigma(&self, r: f64, theta: f64) -> f64 {
        r * r + self.a * self.a * theta.cos().powi(2)
    }

    pub fn delta(&self, r: f64) -> f64 {
        r * r - 2.0 * self.rg * r + self.a * self.a
    }

    pub fn g_tt(&self, r: f64, theta: f64) -> f64 {
        let sigma = self.sigma(r, theta);
        -(1.0 - 2.0 * self.rg * r / sigma)
    }

    pub fn g_t_phi(&self, r: f64, theta: f64) -> f64 {
        let sigma = self.sigma(r, theta);
        -2.0 * self.rg * r * self.a * theta.sin().powi(2) / sigma
    }

    pub fn frame_dragging_rate(&self, r: f64, theta: f64) -> f64 {
        let sigma = self.sigma(r, theta);
        let a2 = self.a * self.a;
        2.0 * self.rg * r * self.a * C
            / (sigma * (r * r + a2) + 2.0 * self.rg * r * a2 * theta.sin().powi(2))
    }
}

pub fn lorentz_factor(v: f64) -> f64 {
    1.0 / (1.0 - (v / C).powi(2)).max(1e-30).sqrt()
}

pub fn relativistic_energy(mass: f64, v: f64) -> f64 {
    lorentz_factor(v) * mass * C * C
}

pub fn relativistic_momentum(mass: f64, v: f64) -> f64 {
    lorentz_factor(v) * mass * v
}

pub fn relativistic_doppler(frequency: f64, v: f64, angle: f64) -> f64 {
    let beta = v / C;
    let gamma = lorentz_factor(v);
    frequency / (gamma * (1.0 - beta * angle.cos()))
}

pub fn time_dilation(proper_time: f64, v: f64) -> f64 {
    proper_time * lorentz_factor(v)
}

pub fn length_contraction(proper_length: f64, v: f64) -> f64 {
    proper_length / lorentz_factor(v)
}

pub fn gravitational_wave_power(m1: f64, m2: f64, separation: f64) -> f64 {
    32.0 / 5.0 * G.powi(4) * (m1 * m2).powi(2) * (m1 + m2) / (C.powi(5) * separation.powi(5))
}