sciforge_lib/physics/nuclear/
stellar.rs

1use super::nuclide::binding_energy_per_nucleon_semf;
2use super::processes::ProcessType;
3use crate::constants::{
4    C, ELECTRON_MASS_KG, G, H, K_B, PROTON_MASS_KG, SIGMA_SB, SOLAR_LUMINOSITY, SOLAR_MASS,
5    SOLAR_RADIUS, TROPICAL_YEAR,
6};
7
8pub struct StellarCore {
9    pub temperature_k: f64,
10    pub density_kg_m3: f64,
11    pub hydrogen_fraction: f64,
12    pub helium_fraction: f64,
13    pub metal_fraction: f64,
14    pub mass_solar: f64,
15}
16
17impl StellarCore {
18    pub fn new(mass_solar: f64, temperature_k: f64, density_kg_m3: f64) -> Self {
19        Self {
20            temperature_k,
21            density_kg_m3,
22            hydrogen_fraction: 0.70,
23            helium_fraction: 0.28,
24            metal_fraction: 0.02,
25            mass_solar,
26        }
27    }
28
29    pub fn luminosity_solar(&self) -> f64 {
30        self.mass_solar.powf(3.5)
31    }
32
33    pub fn main_sequence_lifetime_years(&self) -> f64 {
34        1e10 / self.mass_solar.powf(2.5)
35    }
36
37    pub fn active_processes(&self) -> Vec<ProcessType> {
38        super::processes::active_processes(self.temperature_k)
39    }
40
41    pub fn dominant_process(&self) -> Option<ProcessType> {
42        if self.temperature_k >= 2.7e9 {
43            Some(ProcessType::SiliconBurning)
44        } else if self.temperature_k >= 1.5e9 {
45            Some(ProcessType::OxygenBurning)
46        } else if self.temperature_k >= 1.2e9 {
47            Some(ProcessType::NeonBurning)
48        } else if self.temperature_k >= 5e8 {
49            Some(ProcessType::CarbonBurning)
50        } else if self.temperature_k >= 1e8 {
51            Some(ProcessType::TripleAlpha)
52        } else if self.temperature_k >= 15e6 {
53            Some(ProcessType::CNOCycle)
54        } else if self.temperature_k >= 4e6 {
55            Some(ProcessType::PPChain)
56        } else {
57            None
58        }
59    }
60
61    pub fn evolve_step(&mut self, dt_years: f64) {
62        let rate = match self.dominant_process() {
63            Some(ProcessType::PPChain) => 1e-11,
64            Some(ProcessType::CNOCycle) => 5e-11,
65            Some(ProcessType::TripleAlpha) => 1e-9,
66            Some(ProcessType::CarbonBurning) => 1e-7,
67            _ => 0.0,
68        };
69        let dh = rate * dt_years;
70        if self.hydrogen_fraction > dh {
71            self.hydrogen_fraction -= dh;
72            self.helium_fraction += dh * 0.99;
73            self.metal_fraction += dh * 0.01;
74        } else if self.helium_fraction > dh {
75            self.helium_fraction -= dh;
76            self.metal_fraction += dh;
77            self.temperature_k *= 1.0 + 1e-12 * dt_years;
78        }
79    }
80}
81
82pub fn chandrasekhar_limit() -> f64 {
83    1.4
84}
85pub fn tolman_oppenheimer_volkoff_limit() -> f64 {
86    2.17
87}
88pub fn neutron_drip_density() -> f64 {
89    4e14
90}
91
92pub fn iron_peak_binding_energy() -> f64 {
93    binding_energy_per_nucleon_semf(26, 56)
94}
95
96pub fn eddington_luminosity_solar(mass_solar: f64) -> f64 {
97    3.2e4 * mass_solar
98}
99
100pub fn kelvin_helmholtz_timescale_years(
101    mass_solar: f64,
102    radius_solar: f64,
103    luminosity_solar: f64,
104) -> f64 {
105    let m = mass_solar * SOLAR_MASS;
106    let r = radius_solar * SOLAR_RADIUS;
107    let l = luminosity_solar * SOLAR_LUMINOSITY;
108    G * m * m / (2.0 * r * l) / (365.25 * 86400.0)
109}
110
111pub fn jeans_mass_solar(temperature_k: f64, density_kg_m3: f64) -> f64 {
112    let cs = (K_B * temperature_k / PROTON_MASS_KG).sqrt();
113    let mj = (std::f64::consts::PI / 6.0) * cs.powi(3) / (G.powf(1.5) * density_kg_m3.sqrt());
114    mj / SOLAR_MASS
115}
116
117pub fn schwarzschild_radius_km(mass_solar: f64) -> f64 {
118    2.953 * mass_solar
119}
120
121pub fn nuclear_timescale_years(mass_solar: f64, luminosity_solar: f64, efficiency: f64) -> f64 {
122    let e = efficiency * mass_solar * SOLAR_MASS * C * C;
123    let l = luminosity_solar * SOLAR_LUMINOSITY;
124    e / l / TROPICAL_YEAR
125}
126
127pub fn core_collapse_min_mass_solar() -> f64 {
128    8.0
129}
130
131pub fn white_dwarf_radius_km(mass_solar: f64) -> f64 {
132    let r_earth_km = 6371.0;
133    r_earth_km * (chandrasekhar_limit() / mass_solar).powf(1.0 / 3.0)
134}
135
136pub fn electron_degeneracy_pressure(density_kg_m3: f64) -> f64 {
137    let n_e = density_kg_m3 / PROTON_MASS_KG;
138    (H * H / (5.0 * ELECTRON_MASS_KG))
139        * (3.0 / (8.0 * std::f64::consts::PI)).powf(2.0 / 3.0)
140        * n_e.powf(5.0 / 3.0)
141}
142
143pub fn neutron_star_radius_km(mass_solar: f64) -> f64 {
144    10.0 * (1.4 / mass_solar).powf(1.0 / 3.0)
145}
146
147pub fn luminosity_radius_temperature(radius_solar: f64, temperature_k: f64) -> f64 {
148    let r = radius_solar * SOLAR_RADIUS;
149    4.0 * std::f64::consts::PI * r * r * SIGMA_SB * temperature_k.powi(4) / SOLAR_LUMINOSITY
150}
151
152pub fn stellar_wind_mass_loss(luminosity_solar: f64, escape_velocity_km_s: f64) -> f64 {
153    let l = luminosity_solar * SOLAR_LUMINOSITY;
154    let v = escape_velocity_km_s * 1e3;
155    let m_sun_per_yr = SOLAR_MASS / TROPICAL_YEAR;
156    l / (v * v) / m_sun_per_yr
157}
sciforge_lib/physics/nuclear/stellar.rs

sciforge_lib/physics/nuclear/
stellar.rs
