blackholesfactory 0.0.2

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

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

pub struct GravitationalLens {
    pub lens_mass: f64,
    pub lens_distance: f64,
    pub source_distance: f64,
}

impl GravitationalLens {
    pub fn new(mass: f64, d_l: f64, d_s: f64) -> Self {
        Self {
            lens_mass: mass,
            lens_distance: d_l,
            source_distance: d_s,
        }
    }

    pub fn einstein_radius(&self) -> f64 {
        let d_ls = self.source_distance - self.lens_distance;
        if d_ls <= 0.0 {
            return 0.0;
        }
        (4.0 * G * self.lens_mass * d_ls / (C * C * self.lens_distance * self.source_distance))
            .max(0.0)
            .sqrt()
    }

    pub fn deflection_angle(&self, impact_parameter: f64) -> f64 {
        4.0 * G * self.lens_mass / (C * C * impact_parameter)
    }

    pub fn image_positions(&self, source_angle: f64) -> (f64, f64) {
        let theta_e = self.einstein_radius();
        let disc = source_angle * source_angle + 4.0 * theta_e * theta_e;
        let sqrt_disc = disc.max(0.0).sqrt();
        (
            0.5 * (source_angle + sqrt_disc),
            0.5 * (source_angle - sqrt_disc),
        )
    }

    pub fn magnification(&self, image_angle: f64) -> f64 {
        let theta_e = self.einstein_radius();
        let u = (image_angle / theta_e).powi(4);
        if (u - 1.0).abs() < 1e-30 {
            return f64::INFINITY;
        }
        u / (u - 1.0)
    }

    pub fn total_magnification(&self, source_angle: f64) -> f64 {
        let (t1, t2) = self.image_positions(source_angle);
        self.magnification(t1).abs() + self.magnification(t2).abs()
    }

    pub fn time_delay(&self, theta1: f64, theta2: f64) -> f64 {
        let theta_e = self.einstein_radius();
        let d_eff = self.lens_distance * self.source_distance
            / (self.source_distance - self.lens_distance).max(1e-30);
        d_eff / (2.0 * C)
            * ((theta1 * theta1 - theta2 * theta2)
                - 2.0 * theta_e * theta_e * (theta1.abs().ln() - theta2.abs().max(1e-30).ln()))
    }

    pub fn shapiro_delay(&self, closest_approach: f64) -> f64 {
        let rs = schwarzschild_radius(self.lens_mass);
        rs / C
            * (4.0 * self.lens_distance * self.source_distance
                / (closest_approach * closest_approach))
                .ln()
    }
}

pub fn light_bending_strong_field(mass: f64, impact_parameter: f64, n_orbits: u32) -> f64 {
    let r_ph = photon_sphere_radius(mass);
    let b_crit = r_ph / (1.0 - schwarzschild_radius(mass) / r_ph).max(1e-30).sqrt();
    let db = (impact_parameter - b_crit).abs() / b_crit;
    if db < 1e-10 {
        return f64::INFINITY;
    }
    -PI + 2.0 * n_orbits as f64 * PI + (1.0 / db).ln()
}