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

use crate::config::parameters::*;
use crate::physics::singularity::Singularity;

pub struct BinaryBlackHole {
    pub primary: Singularity,
    pub secondary: Singularity,
    pub separation: f64,
    pub eccentricity: f64,
}

impl BinaryBlackHole {
    pub fn new(mass_1: f64, mass_2: f64, separation: f64) -> Self {
        Self {
            primary: Singularity::new(mass_1),
            secondary: Singularity::new(mass_2),
            separation,
            eccentricity: 0.0,
        }
    }

    pub fn with_spins(mut self, spin_1: f64, spin_2: f64) -> Self {
        self.primary = self.primary.with_spin(spin_1);
        self.secondary = self.secondary.with_spin(spin_2);
        self
    }

    pub fn with_eccentricity(mut self, e: f64) -> Self {
        self.eccentricity = e.clamp(0.0, 0.999);
        self
    }

    pub fn total_mass(&self) -> f64 {
        self.primary.mass + self.secondary.mass
    }

    pub fn mass_ratio(&self) -> f64 {
        self.secondary.mass / self.primary.mass
    }

    pub fn symmetric_mass_ratio(&self) -> f64 {
        let q = self.mass_ratio();
        q / (1.0 + q).powi(2)
    }

    pub fn chirp_mass(&self) -> f64 {
        let m1 = self.primary.mass;
        let m2 = self.secondary.mass;
        (m1 * m2).powf(3.0 / 5.0) / (m1 + m2).powf(1.0 / 5.0)
    }

    pub fn orbital_period(&self) -> f64 {
        2.0 * std::f64::consts::PI * (self.separation.powi(3) / (G * self.total_mass())).sqrt()
    }

    pub fn gravitational_wave_frequency(&self) -> f64 {
        2.0 / self.orbital_period()
    }

    pub fn merger_timescale(&self) -> f64 {
        let m1 = self.primary.mass;
        let m2 = self.secondary.mass;
        let a = self.separation;
        let e = self.eccentricity;
        let circ_factor = (1.0 - e * e).powf(3.5);
        5.0 / 256.0 * C.powi(5) * a.powi(4) * circ_factor / (G.powi(3) * m1 * m2 * (m1 + m2))
    }

    pub fn gravitational_wave_strain(&self, distance: f64) -> f64 {
        let mc = self.chirp_mass();
        let f = self.gravitational_wave_frequency();
        4.0 * G.powf(5.0 / 3.0) * mc * (std::f64::consts::PI * f).powf(2.0 / 3.0)
            / (C.powi(4) * distance)
    }

    pub fn gravitational_wave_luminosity(&self) -> f64 {
        let m1 = self.primary.mass;
        let m2 = self.secondary.mass;
        let a = self.separation;
        32.0 / 5.0 * G.powi(4) * (m1 * m2).powi(2) * (m1 + m2) / (C.powi(5) * a.powi(5))
    }

    pub fn final_spin_estimate(&self) -> f64 {
        let eta = self.symmetric_mass_ratio();
        (2.0 * 3.0_f64.sqrt() * eta + 0.5 * eta * eta).min(0.998)
    }

    pub fn final_mass_estimate(&self) -> f64 {
        let eta = self.symmetric_mass_ratio();
        let radiated_fraction = 0.2 * eta * (1.0 - 0.63 * eta);
        self.total_mass() * (1.0 - radiated_fraction)
    }

    pub fn energy_radiated(&self) -> f64 {
        (self.total_mass() - self.final_mass_estimate()) * C * C
    }

    pub fn orbital_velocity(&self) -> f64 {
        (G * self.total_mass() / self.separation).sqrt()
    }

    pub fn innermost_stable_separation(&self) -> f64 {
        isco_radius(self.total_mass(), self.final_spin_estimate())
    }
}