sagittariusas 0.0.1

Simulation engine for Sagittarius A* — Kerr spacetime, accretion, jets, lensing, and shadow observables
Documentation
use crate::config::parameters::{C, G};
use sciforge::hub::prelude::astronomy::stellar;
use sciforge::hub::prelude::physics::relativity;
use std::f64::consts::PI;

#[derive(Debug, Clone, Copy)]
pub struct SchwarzschildMetric {
    pub mass: f64,
    pub rs: f64,
}

impl SchwarzschildMetric {
    pub fn new(mass: f64) -> Self {
        let rs = stellar::schwarzschild_radius(mass);
        Self { mass, rs }
    }

    pub fn g_tt(&self, r: f64) -> f64 {
        -(1.0 - self.rs / r)
    }

    pub fn g_rr(&self, r: f64) -> f64 {
        1.0 / (1.0 - self.rs / r)
    }

    pub fn proper_distance(&self, r1: f64, r2: f64, steps: usize) -> f64 {
        let f = |r: f64| -> f64 { self.g_rr(r).sqrt() };
        sciforge::hub::prelude::maths::integration::composite_simpson(f, r1, r2, steps)
    }

    pub fn proper_time_ratio(&self, r: f64) -> f64 {
        (1.0 - self.rs / r).sqrt()
    }

    pub fn radial_freefall_velocity(&self, r: f64) -> f64 {
        C * (self.rs / r).sqrt()
    }

    pub fn perihelion_precession(&self, semi_major: f64, eccentricity: f64) -> f64 {
        let l = semi_major * (1.0 - eccentricity * eccentricity);
        6.0 * PI * G * self.mass / (C * C * l)
    }
}

#[derive(Debug, Clone, Copy)]
pub struct KerrMetric {
    pub mass: f64,
    pub a: f64,
    pub rs: f64,
}

impl KerrMetric {
    pub fn new(mass: f64, spin: f64) -> Self {
        let rs = stellar::schwarzschild_radius(mass);
        let rg = rs / 2.0;
        let a = spin * rg;
        Self { mass, a, rs }
    }

    pub fn sigma(&self, r: f64, theta: f64) -> f64 {
        r * r + self.a * self.a * theta.cos().powi(2)
    }

    pub fn delta(&self, r: f64) -> f64 {
        r * r - self.rs * r + self.a * self.a
    }

    pub fn event_horizon_outer(&self) -> f64 {
        let rg = self.rs / 2.0;
        rg + (rg * rg - self.a * self.a).sqrt()
    }

    pub fn event_horizon_inner(&self) -> f64 {
        let rg = self.rs / 2.0;
        rg - (rg * rg - self.a * self.a).sqrt()
    }

    pub fn ergosphere_radius(&self, theta: f64) -> f64 {
        let rg = self.rs / 2.0;
        rg + (rg * rg - self.a * self.a * theta.cos().powi(2)).sqrt()
    }

    pub fn frame_dragging_rate(&self, r: f64, theta: f64) -> f64 {
        let sig = self.sigma(r, theta);
        self.rs * r * self.a * C / (sig * sig + self.rs * r * self.a * self.a * theta.sin().powi(2))
    }
}

pub fn lorentz_factor(v: f64) -> f64 {
    relativity::gamma_factor(v)
}

pub fn relativistic_energy(rest_mass: f64, v: f64) -> f64 {
    relativity::relativistic_energy(rest_mass, v)
}

pub fn relativistic_momentum(rest_mass: f64, v: f64) -> f64 {
    relativity::relativistic_momentum(rest_mass, v)
}

pub fn relativistic_doppler(v: f64, angle: f64) -> f64 {
    relativity::relativistic_doppler(1.0, v, angle)
}

pub fn time_dilation(v: f64, proper_time: f64) -> f64 {
    relativity::time_dilation(proper_time, v)
}

pub fn length_contraction(v: f64, proper_length: f64) -> f64 {
    relativity::length_contraction(proper_length, v)
}

pub fn gravitational_wave_power(m1: f64, m2: f64, r: f64) -> f64 {
    let mu = m1 * m2 / (m1 + m2);
    let m_total = m1 + m2;
    let prefactor = 32.0 / 5.0 * G.powi(4) / C.powi(5);
    prefactor * mu * mu * m_total * m_total * m_total / r.powi(5)
}