use sciforge::hub::domain::common::constants::{G, K_B, SOLAR_MASS};
use crate::{
CORE_RADIUS_FRAC, MEV_TO_JOULE, PP_CHAIN_ENERGY_MEV, PROTON_MASS,
SOLAR_LUMINOSITY, SOLAR_RADIUS,
};
pub const CORE_CENTER_TEMP_K: f64 = 1.57e7;
pub const CORE_EDGE_TEMP_K: f64 = 7.0e6;
pub const CORE_CENTER_DENSITY: f64 = 1.622e5;
pub const CORE_EDGE_DENSITY: f64 = 2.0e4;
pub const CORE_CENTER_PRESSURE_PA: f64 = 2.477e16;
pub const CORE_EDGE_PRESSURE_PA: f64 = 3.0e15;
pub const CORE_HYDROGEN_FRACTION: f64 = 0.34;
pub const CORE_HELIUM_FRACTION: f64 = 0.64;
pub const CORE_METAL_FRACTION: f64 = 0.02;
pub const CORE_MEAN_MOLECULAR_WEIGHT: f64 = 0.83;
pub const CORE_ROTATION_PERIOD_DAYS: f64 = 29.0;
pub const CORE_NEUTRINO_LUMINOSITY_FRACTION: f64 = 0.02;
pub const PP_CHAIN_FRACTION: f64 = 0.991;
pub const CNO_FRACTION: f64 = 0.009;
pub struct SolarCore {
pub outer_radius_frac: f64,
pub center_temperature_k: f64,
pub edge_temperature_k: f64,
pub center_density: f64,
pub edge_density: f64,
pub center_pressure_pa: f64,
pub hydrogen_fraction: f64,
pub helium_fraction: f64,
}
impl SolarCore {
pub fn new() -> Self {
Self {
outer_radius_frac: CORE_RADIUS_FRAC,
center_temperature_k: CORE_CENTER_TEMP_K,
edge_temperature_k: CORE_EDGE_TEMP_K,
center_density: CORE_CENTER_DENSITY,
edge_density: CORE_EDGE_DENSITY,
center_pressure_pa: CORE_CENTER_PRESSURE_PA,
hydrogen_fraction: CORE_HYDROGEN_FRACTION,
helium_fraction: CORE_HELIUM_FRACTION,
}
}
pub fn radius_m(&self) -> f64 {
self.outer_radius_frac * SOLAR_RADIUS
}
pub fn volume_m3(&self) -> f64 {
(4.0 / 3.0) * std::f64::consts::PI * self.radius_m().powi(3)
}
pub fn mass_kg(&self) -> f64 {
self.mean_density() * self.volume_m3()
}
pub fn mass_fraction(&self) -> f64 {
self.mass_kg() / SOLAR_MASS
}
pub fn mean_temperature(&self) -> f64 {
(self.center_temperature_k + self.edge_temperature_k) / 2.0
}
pub fn mean_density(&self) -> f64 {
(self.center_density + self.edge_density) / 2.0
}
pub fn temperature_at_frac(&self, frac: f64) -> f64 {
let f = frac.clamp(0.0, 1.0);
self.center_temperature_k - (self.center_temperature_k - self.edge_temperature_k) * f * f
}
pub fn density_at_frac(&self, frac: f64) -> f64 {
let f = frac.clamp(0.0, 1.0);
self.center_density - (self.center_density - self.edge_density) * f * f
}
pub fn pressure_at_frac(&self, frac: f64) -> f64 {
let f = frac.clamp(0.0, 1.0);
self.center_pressure_pa * (1.0 - f * f)
}
pub fn sound_speed_at_frac(&self, frac: f64) -> f64 {
let gamma = 5.0 / 3.0;
let t = self.temperature_at_frac(frac);
let mu = CORE_MEAN_MOLECULAR_WEIGHT;
(gamma * K_B * t / (mu * PROTON_MASS)).sqrt()
}
pub fn thermal_energy(&self) -> f64 {
1.5 * K_B * self.mean_temperature() * self.volume_m3() * self.mean_density()
/ (CORE_MEAN_MOLECULAR_WEIGHT * PROTON_MASS)
}
pub fn gravitational_energy(&self) -> f64 {
let m = self.mass_kg();
3.0 * G * m * m / (5.0 * self.radius_m())
}
pub fn energy_generation_rate_center(&self) -> f64 {
let rho = self.center_density;
let t = self.center_temperature_k;
let x = self.hydrogen_fraction;
let t6 = t / 1e6;
1.08e-12 * rho * x * x * t6.powf(4.0)
}
pub fn luminosity_generated(&self) -> f64 {
let epsilon_mean = self.energy_generation_rate_center() * 0.5;
epsilon_mean * self.volume_m3()
}
pub fn luminosity_fraction(&self) -> f64 {
self.luminosity_generated() / SOLAR_LUMINOSITY
}
pub fn neutrino_luminosity(&self) -> f64 {
CORE_NEUTRINO_LUMINOSITY_FRACTION * SOLAR_LUMINOSITY
}
pub fn nuclear_timescale_s(&self) -> f64 {
let m_hydrogen = self.mass_kg() * self.hydrogen_fraction;
let e_per_kg = PP_CHAIN_ENERGY_MEV * MEV_TO_JOULE / (4.0 * PROTON_MASS);
let fuel_energy = m_hydrogen * e_per_kg * 0.007;
fuel_energy / SOLAR_LUMINOSITY
}
pub fn plasma_beta(&self) -> f64 {
let n = self.center_density / (CORE_MEAN_MOLECULAR_WEIGHT * PROTON_MASS);
let p_thermal = n * K_B * self.center_temperature_k;
let b_core = 1.0;
let p_mag = b_core * b_core / (2.0 * crate::MU_0);
p_thermal / p_mag
}
pub fn debye_length_center(&self) -> f64 {
let n_e = self.center_density / (2.0 * PROTON_MASS);
(crate::EPSILON_0 * K_B * self.center_temperature_k
/ (n_e * crate::ELECTRON_CHARGE * crate::ELECTRON_CHARGE))
.sqrt()
}
pub fn rotation_angular_velocity(&self) -> f64 {
2.0 * std::f64::consts::PI / (CORE_ROTATION_PERIOD_DAYS * 86400.0)
}
pub fn mean_free_path(&self) -> f64 {
let kappa = 1.0;
1.0 / (kappa * 0.1 * self.center_density)
}
pub fn photon_diffusion_time_to_edge(&self) -> f64 {
let mfp = self.mean_free_path();
let r = self.radius_m();
let n_steps = (r / mfp).powi(2);
n_steps * mfp / crate::SPEED_OF_LIGHT
}
pub fn pp_to_cno_ratio(&self) -> f64 {
PP_CHAIN_FRACTION / CNO_FRACTION
}
}
impl Default for SolarCore {
fn default() -> Self {
Self::new()
}
}
pub fn core_temperature_profile(n_points: usize) -> Vec<(f64, f64)> {
let core = SolarCore::new();
(0..n_points)
.map(|i| {
let frac = i as f64 / (n_points - 1) as f64;
(frac * core.radius_m(), core.temperature_at_frac(frac))
})
.collect()
}
pub fn core_density_profile(n_points: usize) -> Vec<(f64, f64)> {
let core = SolarCore::new();
(0..n_points)
.map(|i| {
let frac = i as f64 / (n_points - 1) as f64;
(frac * core.radius_m(), core.density_at_frac(frac))
})
.collect()
}
pub fn hydrogen_depletion_rate() -> f64 {
let e_per_reaction = PP_CHAIN_ENERGY_MEV * MEV_TO_JOULE;
let reactions_per_second = SOLAR_LUMINOSITY / e_per_reaction;
reactions_per_second * 4.0 * PROTON_MASS
}
pub fn helium_production_rate() -> f64 {
let alpha_mass = 4.0 * PROTON_MASS * 0.993;
let e_per_reaction = PP_CHAIN_ENERGY_MEV * MEV_TO_JOULE;
let reactions_per_second = SOLAR_LUMINOSITY / e_per_reaction;
reactions_per_second * alpha_mass
}
pub fn core_convective_stability(temperature_gradient: f64) -> bool {
let adiabatic_gradient = 0.4;
temperature_gradient < adiabatic_gradient
}