use sciforge::hub::domain::common::constants::{K_B, SOLAR_MASS};
use crate::{
CNO_CYCLE_ENERGY_MEV, MEV_TO_JOULE, PP_CHAIN_ENERGY_MEV, PROTON_MASS, SOLAR_CORE_DENSITY,
SOLAR_CORE_TEMP, SOLAR_LUMINOSITY, SOLAR_RADIUS, SPEED_OF_LIGHT,
};
pub const PP1_BRANCH_FRACTION: f64 = 0.83;
pub const PP2_BRANCH_FRACTION: f64 = 0.16;
pub const PP3_BRANCH_FRACTION: f64 = 0.01;
pub const CNO_FRACTION_SOLAR: f64 = 0.017;
pub const CORE_HYDROGEN_MASS_FRAC: f64 = 0.34;
pub const CORE_VOLUME_FRAC: f64 = 0.015625;
pub const DEUTERIUM_BINDING_ENERGY_MEV: f64 = 2.224;
pub const HE3_BINDING_ENERGY_MEV: f64 = 7.718;
pub const HE4_BINDING_ENERGY_MEV: f64 = 28.296;
pub const NEUTRINO_LOSS_PP1_MEV: f64 = 0.267;
pub const NEUTRINO_LOSS_PP2_MEV: f64 = 0.860;
pub const NEUTRINO_LOSS_PP3_MEV: f64 = 7.2;
pub const NEUTRINO_FLUX_EARTH: f64 = 6.5e14;
pub struct FusionReaction {
pub name: &'static str,
pub energy_mev: f64,
pub neutrino_loss_mev: f64,
pub temperature_sensitivity: f64,
pub branch_fraction: f64,
}
impl FusionReaction {
pub fn net_energy_mev(&self) -> f64 {
self.energy_mev - self.neutrino_loss_mev
}
pub fn net_energy_joules(&self) -> f64 {
self.net_energy_mev() * MEV_TO_JOULE
}
pub fn energy_per_kg(&self) -> f64 {
let reactions_per_kg = 1.0 / (4.0 * PROTON_MASS);
reactions_per_kg * self.net_energy_joules()
}
pub fn reaction_rate_factor(&self, temperature_k: f64) -> f64 {
let t9 = temperature_k / 1e9;
t9.powf(self.temperature_sensitivity)
}
pub fn luminosity_contribution(&self) -> f64 {
SOLAR_LUMINOSITY * self.branch_fraction
}
}
pub fn pp1_chain() -> FusionReaction {
FusionReaction {
name: "pp-I chain",
energy_mev: PP_CHAIN_ENERGY_MEV,
neutrino_loss_mev: 2.0 * NEUTRINO_LOSS_PP1_MEV,
temperature_sensitivity: 4.0,
branch_fraction: PP1_BRANCH_FRACTION,
}
}
pub fn pp2_chain() -> FusionReaction {
FusionReaction {
name: "pp-II chain",
energy_mev: PP_CHAIN_ENERGY_MEV,
neutrino_loss_mev: NEUTRINO_LOSS_PP2_MEV,
temperature_sensitivity: 4.0,
branch_fraction: PP2_BRANCH_FRACTION,
}
}
pub fn pp3_chain() -> FusionReaction {
FusionReaction {
name: "pp-III chain",
energy_mev: PP_CHAIN_ENERGY_MEV,
neutrino_loss_mev: NEUTRINO_LOSS_PP3_MEV,
temperature_sensitivity: 4.0,
branch_fraction: PP3_BRANCH_FRACTION,
}
}
pub fn cno_cycle() -> FusionReaction {
FusionReaction {
name: "CNO cycle",
energy_mev: CNO_CYCLE_ENERGY_MEV,
neutrino_loss_mev: 1.71,
temperature_sensitivity: 16.0,
branch_fraction: CNO_FRACTION_SOLAR,
}
}
pub fn pp_chain_rate(temperature_k: f64, density: f64, hydrogen_frac: f64) -> f64 {
let t9 = temperature_k / 1e9;
let t9_13 = t9.powf(1.0 / 3.0);
let screening = 1.0 + 0.0123 * (density / 1e5).powf(1.0 / 3.0) / t9;
let s11 = 4.01e-22;
let gamow = (-3.380 / t9_13).exp();
let rate = s11 * gamow * screening / (t9_13 * t9_13);
density * hydrogen_frac * hydrogen_frac * rate
}
pub fn cno_rate(temperature_k: f64, density: f64, hydrogen_frac: f64, z_cno: f64) -> f64 {
let t9 = temperature_k / 1e9;
let t9_13 = t9.powf(1.0 / 3.0);
let gamow = (-15.228 / t9_13).exp();
let rate = 8.67e-17 * gamow / (t9_13 * t9_13);
density * hydrogen_frac * z_cno * rate
}
pub fn total_energy_generation_rate(temperature_k: f64, density: f64) -> f64 {
let x_h = CORE_HYDROGEN_MASS_FRAC;
let z_cno = 0.0169;
let pp = pp_chain_rate(temperature_k, density, x_h);
let cno = cno_rate(temperature_k, density, x_h, z_cno);
let e_pp = PP_CHAIN_ENERGY_MEV * MEV_TO_JOULE;
let e_cno = CNO_CYCLE_ENERGY_MEV * MEV_TO_JOULE;
pp * e_pp / (4.0 * PROTON_MASS) + cno * e_cno / (4.0 * PROTON_MASS)
}
pub fn core_luminosity() -> f64 {
let core_radius = 0.25 * SOLAR_RADIUS;
let volume = (4.0 / 3.0) * std::f64::consts::PI * core_radius.powi(3);
let epsilon = total_energy_generation_rate(SOLAR_CORE_TEMP, SOLAR_CORE_DENSITY);
epsilon * volume
}
pub fn mass_to_energy_rate() -> f64 {
SOLAR_LUMINOSITY / (SPEED_OF_LIGHT * SPEED_OF_LIGHT)
}
pub fn hydrogen_consumption_rate() -> f64 {
mass_to_energy_rate() / 0.007
}
pub fn main_sequence_lifetime_years() -> f64 {
let fuel_mass = SOLAR_MASS * 0.10 * CORE_HYDROGEN_MASS_FRAC;
let rate = hydrogen_consumption_rate();
fuel_mass / (rate * crate::SECONDS_PER_YEAR)
}
pub fn gamow_peak_energy(temperature_k: f64) -> f64 {
let e_g = 493.0 * 1e3 * crate::ELECTRON_CHARGE;
let kt = K_B * temperature_k;
(e_g * kt * kt / 4.0).powf(1.0 / 3.0)
}
pub fn gamow_window_width(temperature_k: f64) -> f64 {
let e0 = gamow_peak_energy(temperature_k);
let kt = K_B * temperature_k;
4.0 * (e0 * kt / 3.0).sqrt()
}
pub fn coulomb_barrier_mev(z1: f64, z2: f64, r_fm: f64) -> f64 {
let e2 = 1.44;
z1 * z2 * e2 / r_fm
}
pub fn tunneling_probability(energy_mev: f64, barrier_mev: f64) -> f64 {
if energy_mev >= barrier_mev {
return 1.0;
}
let eta = std::f64::consts::PI * (barrier_mev / energy_mev - 1.0).sqrt();
(-2.0 * eta).exp()
}
pub fn neutrino_luminosity() -> f64 {
let pp1 = pp1_chain();
let pp2 = pp2_chain();
let pp3 = pp3_chain();
let cno = cno_cycle();
let reactions_per_second = SOLAR_LUMINOSITY / (PP_CHAIN_ENERGY_MEV * MEV_TO_JOULE);
pp1.branch_fraction * 2.0 * NEUTRINO_LOSS_PP1_MEV * MEV_TO_JOULE * reactions_per_second
+ pp2.branch_fraction * NEUTRINO_LOSS_PP2_MEV * MEV_TO_JOULE * reactions_per_second
+ pp3.branch_fraction * NEUTRINO_LOSS_PP3_MEV * MEV_TO_JOULE * reactions_per_second
+ cno.branch_fraction * cno.neutrino_loss_mev * MEV_TO_JOULE * reactions_per_second
}