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
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

use std::f64::consts::{FRAC_PI_2, PI};
use std::process;

use sagittariusas::config::parameters::{
    C, SGR_A_DISTANCE, SGR_A_MASS, SGR_A_SPIN, SIGMA_SB, SOLAR_MASS, isco_radius,
};
use sagittariusas::engine::accretion_disk::AccretionDisk;
use sagittariusas::engine::jet_dynamics::RelativisticJet;
use sagittariusas::engine::spacetime::{KerrSpacetime, SpacetimePoint};
use sagittariusas::observables::lensing::GravitationalLens;
use sagittariusas::observables::radiation::AccretionRadiation;
use sagittariusas::observables::shadows::{BlackHoleShadow, eht_expected_diameter_uas};
use sagittariusas::physics::relativity::{SchwarzschildMetric, gravitational_wave_power};
use sagittariusas::physics::singularity::Singularity;
use sagittariusas::utils::math::log_range;

struct Validation {
    shadow_deviation: f64,
    model_deviation: f64,
    precession_deviation: f64,
    sub_eddington: bool,
    photon_captured: bool,
    jet_efficiency: f64,
    disk_luminosity: f64,
    gw_power: f64,
    einstein_radius: f64,
}

impl Validation {
    fn passed(&self) -> bool {
        self.shadow_deviation < 0.2
            && self.model_deviation < 0.3
            && self.precession_deviation < 0.2
            && self.sub_eddington
            && self.photon_captured
            && self.jet_efficiency > 0.0
            && self.jet_efficiency < 1.0
            && self.disk_luminosity > 0.0
            && self.gw_power > 0.0
            && self.einstein_radius > 0.0
    }
}

fn main() {
    let bh = Singularity::kerr(SGR_A_MASS, SGR_A_SPIN);
    if bh.is_naked() {
        process::exit(1);
    }
    let r_isco = isco_radius(SGR_A_MASS, SGR_A_SPIN);

    let shadow = BlackHoleShadow::new(SGR_A_MASS, SGR_A_SPIN, 0.785);
    let predicted_uas = shadow.shadow_diameter_uas(SGR_A_DISTANCE);
    let model_uas = eht_expected_diameter_uas();
    let shadow_deviation = (predicted_uas - 51.8).abs() / 51.8;
    let model_deviation = (predicted_uas - model_uas).abs() / model_uas;

    let schw = SchwarzschildMetric::new(SGR_A_MASS);
    let s2_a = 970.0 * 1.496e11;
    let s2_e = 0.8843;
    let s2_peri = s2_a * (1.0 - s2_e);
    let precession_arcmin = schw.perihelion_precession(s2_a, s2_e).to_degrees() * 60.0;
    let precession_deviation = (precession_arcmin - 12.1).abs() / 12.1;

    let mdot = 1e-9 * SOLAR_MASS / (365.25 * 86400.0);
    let disk = AccretionDisk::new(SGR_A_MASS, mdot);
    let accr = AccretionRadiation::new(SGR_A_MASS, mdot);
    let radii = log_range(r_isco, 500.0 * bh.gravitational_radius(), 200);
    let disk_lum: f64 = radii
        .windows(2)
        .map(|w| {
            let r = 0.5 * (w[0] + w[1]);
            let dr = w[1] - w[0];
            SIGMA_SB * disk.temperature_profile(r).powi(4) * 4.0 * PI * r * dr
        })
        .sum();

    let jet = RelativisticJet::new(SGR_A_MASS, SGR_A_SPIN);
    let bz = jet.blandford_znajek_power();
    let kin = jet.kinetic_luminosity();

    let gw = gravitational_wave_power(SGR_A_MASS, 14.0 * SOLAR_MASS, s2_peri);

    let ks = KerrSpacetime::new(SGR_A_MASS, SGR_A_SPIN);
    let rh = ks.event_horizon();
    let traj = ks.trace_photon(
        SpacetimePoint::new(0.0, 10.0 * rh, FRAC_PI_2, 0.0),
        [-C, 0.0, 0.0],
        10000,
        0.01,
    );

    let lens = GravitationalLens::new(SGR_A_MASS, SGR_A_DISTANCE, 2.0 * SGR_A_DISTANCE);

    let v = Validation {
        shadow_deviation,
        model_deviation,
        precession_deviation,
        sub_eddington: accr.eddington_ratio() < 1.0,
        photon_captured: traj.last().is_none_or(|p| p.r < rh * 1.1),
        jet_efficiency: kin / bz,
        disk_luminosity: disk_lum,
        gw_power: gw,
        einstein_radius: lens.einstein_radius(),
    };

    process::exit(if v.passed() { 0 } else { 1 });
}