sagittariusas 0.0.1

Simulation engine for Sagittarius A* — Kerr spacetime, accretion, jets, lensing, and shadow observables
Documentation
use crate::config::parameters::C;
use std::f64::consts::PI;

#[derive(Debug, Clone, Copy)]
pub enum JetLaunchMechanism {
    BlandfordZnajek,
    BlandfordPayne,
}

#[derive(Debug, Clone)]
pub struct RelativisticJet {
    pub black_hole_mass: f64,
    pub spin: f64,
    pub magnetic_field: f64,
    pub bulk_lorentz_factor: f64,
    pub half_opening_angle: f64,
    pub mechanism: JetLaunchMechanism,
}

impl RelativisticJet {
    pub fn new(bh_mass: f64, spin: f64) -> Self {
        Self {
            black_hole_mass: bh_mass,
            spin: spin.clamp(0.0, 1.0),
            magnetic_field: 100.0,
            bulk_lorentz_factor: 5.0,
            half_opening_angle: 0.1,
            mechanism: JetLaunchMechanism::BlandfordZnajek,
        }
    }

    pub fn with_magnetic_field(mut self, b: f64) -> Self {
        self.magnetic_field = b;
        self
    }

    pub fn with_lorentz_factor(mut self, gamma: f64) -> Self {
        self.bulk_lorentz_factor = gamma.max(1.0);
        self
    }

    pub fn with_opening_angle(mut self, angle: f64) -> Self {
        self.half_opening_angle = angle.clamp(0.01, PI / 4.0);
        self
    }

    pub fn jet_velocity(&self) -> f64 {
        let g = self.bulk_lorentz_factor;
        C * (1.0 - 1.0 / (g * g)).sqrt()
    }

    pub fn blandford_znajek_power(&self) -> f64 {
        let rs =
            sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(self.black_hole_mass);
        let rg = rs / 2.0;
        let r_h = rg * (1.0 + (1.0 - self.spin * self.spin).sqrt());
        let omega_h = self.spin * C / (2.0 * r_h);
        let b2 = self.magnetic_field * self.magnetic_field;
        let mu_0 = 4.0 * PI * 1e-7;
        let phi_bh = b2 * r_h * r_h / mu_0;
        phi_bh * (omega_h * r_h / C).powi(2) * C / (6.0 * PI)
    }

    pub fn kinetic_luminosity(&self) -> f64 {
        let g = self.bulk_lorentz_factor;
        self.blandford_znajek_power() * (g - 1.0) / g
    }

    pub fn doppler_factor(&self, viewing_angle: f64) -> f64 {
        let beta = self.jet_velocity() / C;
        let g = self.bulk_lorentz_factor;
        1.0 / (g * (1.0 - beta * viewing_angle.cos()))
    }

    pub fn apparent_superluminal_velocity(&self, viewing_angle: f64) -> f64 {
        let beta = self.jet_velocity() / C;
        beta * viewing_angle.sin() / (1.0 - beta * viewing_angle.cos())
    }

    pub fn max_apparent_velocity(&self) -> f64 {
        let beta = self.jet_velocity() / C;
        let g = self.bulk_lorentz_factor;
        beta * g * C
    }

    pub fn beam_radius_at(&self, distance_along_jet: f64) -> f64 {
        distance_along_jet * self.half_opening_angle.tan()
    }

    pub fn magnetic_field_at(&self, distance_along_jet: f64) -> f64 {
        let rs =
            sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(self.black_hole_mass);
        self.magnetic_field * (rs / distance_along_jet).powi(1)
    }

    pub fn synchrotron_cooling_time(&self, electron_gamma: f64) -> f64 {
        let sigma_t = 6.652e-29;
        let m_e = 9.109e-31;
        let b2 = self.magnetic_field * self.magnetic_field;
        let u_b = b2 / (2.0 * 4.0 * PI * 1e-7);
        6.0 * PI * m_e * C / (sigma_t * u_b * electron_gamma)
    }
}