sagittariusas 0.0.1

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

#[derive(Debug, Clone, Copy)]
pub struct ThermalRadiation {
    pub temperature: f64,
}

impl ThermalRadiation {
    pub fn new(temperature: f64) -> Self {
        Self { temperature }
    }

    pub fn planck_spectral_radiance(&self, frequency: f64) -> f64 {
        let a = 2.0 * H * frequency.powi(3) / (C * C);
        let exponent = H * frequency / (K_B * self.temperature);
        if exponent > 500.0 {
            return 0.0;
        }
        a / (exponent.exp() - 1.0)
    }

    pub fn wien_peak_frequency(&self) -> f64 {
        2.821 * K_B * self.temperature / H
    }

    pub fn luminosity(&self, area: f64) -> f64 {
        SIGMA_SB * self.temperature.powi(4) * area
    }
}

#[derive(Debug, Clone, Copy)]
pub struct AccretionRadiation {
    pub mass: f64,
    pub accretion_rate: f64,
    pub efficiency: f64,
}

impl AccretionRadiation {
    pub fn new(mass: f64, accretion_rate: f64) -> Self {
        Self {
            mass,
            accretion_rate,
            efficiency: 0.1,
        }
    }

    pub fn with_efficiency(mut self, eta: f64) -> Self {
        self.efficiency = eta.clamp(0.0, 0.42);
        self
    }

    pub fn luminosity(&self) -> f64 {
        self.efficiency * self.accretion_rate * C * C
    }

    pub fn eddington_luminosity(&self) -> f64 {
        crate::config::parameters::eddington_luminosity(self.mass)
    }

    pub fn eddington_ratio(&self) -> f64 {
        self.luminosity() / self.eddington_luminosity()
    }

    pub fn disk_temperature_at(&self, r: f64) -> f64 {
        let rs = sciforge::hub::prelude::astronomy::stellar::schwarzschild_radius(self.mass);
        let r_isco = 3.0 * rs;
        if r < r_isco {
            return 0.0;
        }
        let factor = 3.0 * G * self.mass * self.accretion_rate / (8.0 * PI * SIGMA_SB * r.powi(3));
        let inner_correction = 1.0 - (r_isco / r).sqrt();
        (factor * inner_correction).powf(0.25)
    }

    pub fn peak_temperature(&self) -> f64 {
        let rs = 2.0 * G * self.mass / (C * C);
        let r_peak = 4.08 * rs;
        self.disk_temperature_at(r_peak)
    }

    pub fn characteristic_frequency(&self) -> f64 {
        2.821 * K_B * self.peak_temperature() / H
    }
}

pub fn synchrotron_power(charge: f64, velocity: f64, b_field: f64, mass: f64) -> f64 {
    let beta = velocity / C;
    let gamma = 1.0 / (1.0 - beta * beta).sqrt();
    let r = gamma * mass * velocity / (charge * b_field);
    sciforge::hub::prelude::physics::relativity::synchrotron_power(charge, mass, gamma, r)
}

pub fn bremsstrahlung_emissivity(n_e: f64, temperature: f64, z: f64) -> f64 {
    1.4e-27 * z * z * n_e * n_e * temperature.sqrt()
}