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

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

#[derive(Debug, Clone, Copy)]
pub enum DiskModel {
    ShakuraSunyaev,
    Adaf,
    Slim,
}

#[derive(Debug, Clone)]
pub struct AccretionDisk {
    pub black_hole_mass: f64,
    pub accretion_rate: f64,
    pub alpha_viscosity: f64,
    pub inner_radius: f64,
    pub outer_radius: f64,
    pub model: DiskModel,
}

impl AccretionDisk {
    pub fn new(bh_mass: f64, accretion_rate: f64) -> Self {
        let rs = sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(bh_mass);
        Self {
            black_hole_mass: bh_mass,
            accretion_rate,
            alpha_viscosity: 0.1,
            inner_radius: 3.0 * rs,
            outer_radius: 1000.0 * rs,
            model: DiskModel::ShakuraSunyaev,
        }
    }

    pub fn with_alpha(mut self, alpha: f64) -> Self {
        self.alpha_viscosity = alpha.clamp(0.001, 1.0);
        self
    }

    pub fn with_model(mut self, model: DiskModel) -> Self {
        self.model = model;
        self
    }

    pub fn with_outer_radius(mut self, r: f64) -> Self {
        self.outer_radius = r;
        self
    }

    pub fn temperature_profile(&self, r: f64) -> f64 {
        if r < self.inner_radius {
            return 0.0;
        }
        let factor = 3.0 * G * self.black_hole_mass * self.accretion_rate
            / (8.0 * PI * SIGMA_SB * r.powi(3));
        let correction = 1.0 - (self.inner_radius / r).sqrt();
        let f_rad = match self.model {
            DiskModel::ShakuraSunyaev => 1.0,
            DiskModel::Adaf => self.alpha_viscosity * 0.1,
            DiskModel::Slim => 0.5,
        };
        (f_rad * factor * correction).powf(0.25)
    }

    pub fn surface_density(&self, r: f64) -> f64 {
        if r < self.inner_radius {
            return 0.0;
        }
        let correction = 1.0 - (self.inner_radius / r).sqrt();
        self.accretion_rate * correction / (3.0 * PI * self.viscosity_at(r))
    }

    pub fn viscosity_at(&self, r: f64) -> f64 {
        let cs = self.sound_speed(r);
        let h = self.scale_height(r);
        self.alpha_viscosity * cs * h
    }

    pub fn sound_speed(&self, r: f64) -> f64 {
        let t = self.temperature_profile(r);
        let mu = 0.6;
        let m_p = 1.672_621e-27;
        (sciforge::hub::prelude::constants::K_B * t / (mu * m_p)).sqrt()
    }

    pub fn scale_height(&self, r: f64) -> f64 {
        let cs = self.sound_speed(r);
        let omega = self.orbital_frequency(r);
        if omega < 1e-30 {
            return 0.0;
        }
        cs / omega
    }

    pub fn orbital_frequency(&self, r: f64) -> f64 {
        (G * self.black_hole_mass / (r * r * r)).sqrt()
    }

    pub fn luminosity(&self) -> f64 {
        G * self.black_hole_mass * self.accretion_rate / (2.0 * self.inner_radius)
    }

    pub fn radiative_efficiency(&self) -> f64 {
        let rs =
            sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(self.black_hole_mass);
        rs / (4.0 * self.inner_radius)
    }

    pub fn radial_profile(&self, n_points: usize) -> Vec<(f64, f64)> {
        let log_rin = self.inner_radius.ln();
        let log_rout = self.outer_radius.ln();
        let dlog = (log_rout - log_rin) / n_points as f64;
        (0..n_points)
            .map(|i| {
                let r = (log_rin + (i as f64 + 0.5) * dlog).exp();
                (r, self.temperature_profile(r))
            })
            .collect()
    }

    pub fn total_mass(&self) -> f64 {
        let n = 500;
        let dr = (self.outer_radius - self.inner_radius) / n as f64;
        (0..n)
            .map(|i| {
                let r = self.inner_radius + (i as f64 + 0.5) * dr;
                2.0 * PI * r * self.surface_density(r) * dr
            })
            .sum()
    }
}