sagittariusas 0.0.3

Simulation engine for Sagittarius A* — Kerr spacetime, accretion, jets, lensing, and shadow observables
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

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

#[derive(Debug, Clone, Copy)]
pub struct GravitationalLens {
    pub mass: f64,
    pub distance_lens: f64,
    pub distance_source: f64,
}

impl GravitationalLens {
    pub fn new(mass: f64, distance_lens: f64, distance_source: f64) -> Self {
        Self {
            mass,
            distance_lens,
            distance_source,
        }
    }

    pub fn einstein_radius(&self) -> f64 {
        let d_ls = self.distance_source - self.distance_lens;
        let theta_e_sq =
            4.0 * G * self.mass * d_ls / (C * C * self.distance_lens * self.distance_source);
        theta_e_sq.sqrt()
    }

    pub fn einstein_radius_physical(&self) -> f64 {
        self.einstein_radius() * self.distance_lens
    }

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

    pub fn magnification(&self, beta: f64) -> (f64, f64) {
        let theta_e = self.einstein_radius();
        if theta_e < 1e-30 {
            return (1.0, 0.0);
        }
        let u = beta / theta_e;
        let u2 = u * u;
        let sqrt_term = (u2 + 4.0).sqrt();
        let mu_plus = 0.5 * (u2 + 2.0) / (u * sqrt_term) + 0.5;
        let mu_minus = 0.5 * (u2 + 2.0) / (u * sqrt_term) - 0.5;
        (mu_plus, mu_minus)
    }

    pub fn total_magnification(&self, beta: f64) -> f64 {
        let (mu_p, mu_m) = self.magnification(beta);
        mu_p.abs() + mu_m.abs()
    }

    pub fn image_positions(&self, beta: f64) -> (f64, f64) {
        let theta_e = self.einstein_radius();
        let u = beta / theta_e;
        let sqrt_term = (u * u + 4.0).sqrt();
        let theta_plus = 0.5 * (u + sqrt_term) * theta_e;
        let theta_minus = 0.5 * (u - sqrt_term) * theta_e;
        (theta_plus, theta_minus)
    }

    pub fn time_delay(&self, theta_1: f64, theta_2: f64) -> f64 {
        let theta_e = self.einstein_radius();
        let factor = 4.0 * G * self.mass / (C * C * C)
            * (1.0 + self.distance_lens / (self.distance_source - self.distance_lens));
        let psi_1 = 0.5 * (theta_1 / theta_e).powi(2) - (theta_1 / theta_e).abs().ln();
        let psi_2 = 0.5 * (theta_2 / theta_e).powi(2) - (theta_2 / theta_e).abs().ln();
        factor * (psi_1 - psi_2)
    }
}

pub fn shapiro_delay(mass: f64, closest_approach: f64, distance: f64) -> f64 {
    4.0 * G * mass / (C * C * C) * (2.0 * distance / closest_approach).ln()
}

pub fn light_bending_strong_field(mass: f64, impact_parameter: f64, n_orbits: u32) -> f64 {
    let rs = sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(mass);
    let b_crit = rs * 3.0_f64.sqrt() * 1.5;
    let base_angle = 4.0 * G * mass / (C * C * impact_parameter);
    base_angle + 2.0 * PI * n_orbits as f64 * (1.0 - impact_parameter / b_crit).abs()
}