use sciforge::hub::domain::common::constants::{G, K_B, SIGMA_SB, SOLAR_MASS};
use crate::{
PROTON_MASS, RADIATIVE_ZONE_OUTER_FRAC, SOLAR_LUMINOSITY, SOLAR_RADIUS, SPEED_OF_LIGHT,
};
pub const RADIATIVE_ZONE_BASE_TEMP: f64 = 7.0e6;
pub const RADIATIVE_ZONE_TOP_TEMP: f64 = 2.0e6;
pub const RADIATIVE_ZONE_BASE_DENSITY: f64 = 2.0e4;
pub const RADIATIVE_ZONE_TOP_DENSITY: f64 = 200.0;
pub const RADIATIVE_ZONE_INNER_FRAC: f64 = 0.25;
pub const PHOTON_MEAN_FREE_PATH_RZ: f64 = 2e-2;
pub const PHOTON_ESCAPE_TIME_YEARS: f64 = 1.7e5;
pub const TACHOCLINE_INNER_FRAC: f64 = 0.69;
pub const TACHOCLINE_OUTER_FRAC: f64 = 0.71;
pub struct RadiativeZone {
pub inner_radius_frac: f64,
pub outer_radius_frac: f64,
pub base_temperature_k: f64,
pub top_temperature_k: f64,
pub base_density: f64,
pub top_density: f64,
}
impl Default for RadiativeZone {
fn default() -> Self {
Self::new()
}
}
impl RadiativeZone {
pub fn new() -> Self {
Self {
inner_radius_frac: RADIATIVE_ZONE_INNER_FRAC,
outer_radius_frac: RADIATIVE_ZONE_OUTER_FRAC,
base_temperature_k: RADIATIVE_ZONE_BASE_TEMP,
top_temperature_k: RADIATIVE_ZONE_TOP_TEMP,
base_density: RADIATIVE_ZONE_BASE_DENSITY,
top_density: RADIATIVE_ZONE_TOP_DENSITY,
}
}
pub fn thickness_m(&self) -> f64 {
(self.outer_radius_frac - self.inner_radius_frac) * SOLAR_RADIUS
}
pub fn volume_m3(&self) -> f64 {
let r_out = self.outer_radius_frac * SOLAR_RADIUS;
let r_in = self.inner_radius_frac * SOLAR_RADIUS;
(4.0 / 3.0) * std::f64::consts::PI * (r_out.powi(3) - r_in.powi(3))
}
pub fn temperature_at_frac(&self, radius_frac: f64) -> f64 {
let f = ((radius_frac - self.inner_radius_frac)
/ (self.outer_radius_frac - self.inner_radius_frac))
.clamp(0.0, 1.0);
self.base_temperature_k * (1.0 - f) + self.top_temperature_k * f
}
pub fn density_at_frac(&self, radius_frac: f64) -> f64 {
let f = ((radius_frac - self.inner_radius_frac)
/ (self.outer_radius_frac - self.inner_radius_frac))
.clamp(0.0, 1.0);
(self.base_density.ln() * (1.0 - f) + self.top_density.ln() * f).exp()
}
pub fn pressure_at_frac(&self, radius_frac: f64) -> f64 {
let rho = self.density_at_frac(radius_frac);
let t = self.temperature_at_frac(radius_frac);
rho * K_B * t / (0.61 * PROTON_MASS)
}
pub fn opacity_at_frac(&self, radius_frac: f64) -> f64 {
let rho = self.density_at_frac(radius_frac);
let t = self.temperature_at_frac(radius_frac);
4.34e25 * rho * t.powf(-3.5)
}
pub fn mean_free_path_at_frac(&self, radius_frac: f64) -> f64 {
let kappa = self.opacity_at_frac(radius_frac);
let rho = self.density_at_frac(radius_frac);
1.0 / (kappa * 0.1 * rho).max(1e-30)
}
pub fn radiative_temperature_gradient(&self, radius_frac: f64) -> f64 {
let kappa = self.opacity_at_frac(radius_frac);
let rho = self.density_at_frac(radius_frac);
let t = self.temperature_at_frac(radius_frac);
let r = radius_frac * SOLAR_RADIUS;
3.0 * kappa * 0.1 * rho * SOLAR_LUMINOSITY
/ (16.0
* std::f64::consts::PI
* SPEED_OF_LIGHT
* 4.0
* SIGMA_SB
* t.powi(3)
* r.powi(2))
}
pub fn photon_diffusion_speed(&self, radius_frac: f64) -> f64 {
let l = self.mean_free_path_at_frac(radius_frac);
SPEED_OF_LIGHT * l / (3.0 * self.thickness_m())
}
pub fn photon_random_walk_time(&self) -> f64 {
let l = PHOTON_MEAN_FREE_PATH_RZ;
let d = self.thickness_m();
let n_steps = (d / l).powi(2);
n_steps * l / SPEED_OF_LIGHT
}
pub fn luminosity_at_frac(&self, radius_frac: f64) -> f64 {
let frac = radius_frac.clamp(self.inner_radius_frac, self.outer_radius_frac);
let t = (frac - self.inner_radius_frac) / (self.outer_radius_frac - self.inner_radius_frac);
SOLAR_LUMINOSITY * (1.0 + 0.0 * t)
}
pub fn radiation_pressure_at_frac(&self, radius_frac: f64) -> f64 {
let t = self.temperature_at_frac(radius_frac);
4.0 * SIGMA_SB * t.powi(4) / (3.0 * SPEED_OF_LIGHT)
}
pub fn gas_to_radiation_pressure_ratio(&self, radius_frac: f64) -> f64 {
let p_gas = self.pressure_at_frac(radius_frac);
let p_rad = self.radiation_pressure_at_frac(radius_frac);
p_gas / p_rad.max(1e-30)
}
pub fn brunt_vaisala_frequency(&self, radius_frac: f64) -> f64 {
let g = G * SOLAR_MASS / ((radius_frac * SOLAR_RADIUS).powi(2));
let hp = K_B * self.temperature_at_frac(radius_frac) / (0.61 * PROTON_MASS * g);
let grad_ad = 0.4;
let grad_rad = self.radiative_temperature_gradient(radius_frac);
let n2 = g / hp * (grad_ad - grad_rad);
if n2 > 0.0 { n2.sqrt() } else { 0.0 }
}
}
pub struct Tachocline {
pub inner_radius_frac: f64,
pub outer_radius_frac: f64,
pub thickness_m: f64,
}
impl Default for Tachocline {
fn default() -> Self {
Self::new()
}
}
impl Tachocline {
pub fn new() -> Self {
let inner = TACHOCLINE_INNER_FRAC;
let outer = TACHOCLINE_OUTER_FRAC;
Self {
inner_radius_frac: inner,
outer_radius_frac: outer,
thickness_m: (outer - inner) * SOLAR_RADIUS,
}
}
pub fn shear_rate(&self) -> f64 {
let omega_inner = 2.0 * std::f64::consts::PI / (25.0 * crate::SECONDS_PER_DAY);
let omega_outer = 2.0 * std::f64::consts::PI / (27.0 * crate::SECONDS_PER_DAY);
(omega_inner - omega_outer).abs() / self.thickness_m
}
pub fn magnetic_field_amplification_rate(&self) -> f64 {
self.shear_rate()
}
pub fn diffusion_time(&self) -> f64 {
let eta = 1e3;
self.thickness_m.powi(2) / eta
}
pub fn mean_radius_m(&self) -> f64 {
(self.inner_radius_frac + self.outer_radius_frac) / 2.0 * SOLAR_RADIUS
}
}
pub fn gravity_wave_period(radius_frac: f64, n_bv: f64) -> f64 {
if n_bv <= 0.0 {
return f64::INFINITY;
}
let local_r = radius_frac * SOLAR_RADIUS;
2.0 * std::f64::consts::PI * local_r / (local_r * n_bv)
}
pub fn photon_escape_time_estimate() -> f64 {
PHOTON_ESCAPE_TIME_YEARS * crate::SECONDS_PER_YEAR
}