use sciforge::hub::domain::common::constants::K_B;
use crate::{
CORONA_TEMP, ELECTRON_CHARGE, ELECTRON_MASS, EPSILON_0, MU_0,
PROTON_MASS, SOLAR_LUMINOSITY, SOLAR_RADIUS,
};
pub const CORONA_BASE_DENSITY: f64 = 1e-12;
pub const CORONA_ELECTRON_DENSITY: f64 = 1e15;
pub const CORONA_INNER_RADIUS_RSUN: f64 = 1.01;
pub const CORONA_OUTER_RADIUS_RSUN: f64 = 3.0;
pub const CORONA_HEATING_RATE: f64 = 3e-4;
pub const CORONAL_LOOP_LENGTH: f64 = 1e8;
pub const CORONAL_HOLE_DENSITY_FRACTION: f64 = 0.3;
pub const STREAMER_DENSITY_FACTOR: f64 = 3.0;
pub struct CoronalRegion {
pub name: &'static str,
pub temperature_k: f64,
pub electron_density_m3: f64,
pub magnetic_field_t: f64,
pub scale_height_m: f64,
}
impl CoronalRegion {
pub fn thermal_pressure(&self) -> f64 {
2.0 * self.electron_density_m3 * K_B * self.temperature_k
}
pub fn magnetic_pressure(&self) -> f64 {
self.magnetic_field_t.powi(2) / (2.0 * MU_0)
}
pub fn plasma_beta(&self) -> f64 {
self.thermal_pressure() / self.magnetic_pressure().max(1e-30)
}
pub fn alfven_speed(&self) -> f64 {
let rho = self.electron_density_m3 * PROTON_MASS;
self.magnetic_field_t / (MU_0 * rho).sqrt()
}
pub fn sound_speed(&self) -> f64 {
let gamma = 5.0 / 3.0;
let mu = 0.6;
(gamma * K_B * self.temperature_k / (mu * PROTON_MASS)).sqrt()
}
pub fn debye_length(&self) -> f64 {
(EPSILON_0 * K_B * self.temperature_k
/ (self.electron_density_m3 * ELECTRON_CHARGE * ELECTRON_CHARGE))
.sqrt()
}
pub fn plasma_frequency(&self) -> f64 {
(self.electron_density_m3 * ELECTRON_CHARGE * ELECTRON_CHARGE
/ (EPSILON_0 * ELECTRON_MASS))
.sqrt()
}
pub fn collision_frequency(&self) -> f64 {
let lambda_d = self.debye_length();
let ln_lambda = (12.0 * std::f64::consts::PI * self.electron_density_m3 * lambda_d.powi(3))
.ln()
.max(1.0);
let v_th = (K_B * self.temperature_k / ELECTRON_MASS).sqrt();
self.electron_density_m3 * ELECTRON_CHARGE.powi(4) * ln_lambda
/ (4.0
* std::f64::consts::PI
* EPSILON_0
* EPSILON_0
* ELECTRON_MASS
* ELECTRON_MASS
* v_th.powi(3))
}
pub fn radiative_loss_rate(&self) -> f64 {
let n_e = self.electron_density_m3;
let t = self.temperature_k;
let lambda = if t < 1e5 {
1e-30
} else if t < 1e7 {
1e-22 * (t / 1e6).powf(-0.7)
} else {
1e-23 * (t / 1e7).sqrt()
};
n_e * n_e * lambda
}
pub fn thermal_conduction_flux(&self, length: f64) -> f64 {
let kappa_0 = 9.2e-12;
kappa_0 * self.temperature_k.powf(3.5) / length
}
pub fn emission_measure(&self, volume_m3: f64) -> f64 {
self.electron_density_m3.powi(2) * volume_m3
}
}
pub fn quiet_corona() -> CoronalRegion {
CoronalRegion {
name: "Quiet Corona",
temperature_k: CORONA_TEMP,
electron_density_m3: CORONA_ELECTRON_DENSITY,
magnetic_field_t: 1e-4,
scale_height_m: 5e7,
}
}
pub fn active_region_corona() -> CoronalRegion {
CoronalRegion {
name: "Active Region Corona",
temperature_k: 3.0e6,
electron_density_m3: 5e15,
magnetic_field_t: 1e-2,
scale_height_m: 1e8,
}
}
pub fn coronal_hole() -> CoronalRegion {
CoronalRegion {
name: "Coronal Hole",
temperature_k: 8e5,
electron_density_m3: CORONA_ELECTRON_DENSITY * CORONAL_HOLE_DENSITY_FRACTION,
magnetic_field_t: 5e-4,
scale_height_m: 7e7,
}
}
pub fn coronal_streamer() -> CoronalRegion {
CoronalRegion {
name: "Coronal Streamer",
temperature_k: 1.5e6,
electron_density_m3: CORONA_ELECTRON_DENSITY * STREAMER_DENSITY_FACTOR,
magnetic_field_t: 2e-4,
scale_height_m: 1e8,
}
}
pub fn density_at_height(base_density: f64, height_m: f64, temperature_k: f64) -> f64 {
let mu = 0.6;
let h = K_B * temperature_k / (mu * PROTON_MASS * 274.0);
base_density * (-height_m / h).exp()
}
pub fn coronal_x_ray_luminosity() -> f64 {
SOLAR_LUMINOSITY * 1e-6
}
pub fn coronal_euv_flux(distance_au: f64) -> f64 {
let flux_1au = 4.5e-3;
flux_1au / (distance_au * distance_au)
}
pub fn hydrostatic_scale_height(temperature_k: f64) -> f64 {
let mu = 0.6;
K_B * temperature_k / (mu * PROTON_MASS * 274.0)
}
pub fn coronal_loop_temperature(length: f64, heating_rate: f64) -> f64 {
let n_e = 1e15;
let kappa_0 = 9.2e-12;
let t_rtv = (7.0 * heating_rate * length * length / (8.0 * kappa_0)).powf(2.0 / 7.0);
let p = 2.0 * n_e * K_B * t_rtv;
p / (2.0 * n_e * K_B)
}
pub fn coronal_loop_density(temperature: f64, length: f64) -> f64 {
let p = 2.0 * K_B * temperature / length;
p / (2.0 * K_B * temperature)
}
pub fn thompson_scattering_brightness(electron_density: f64, distance_rsun: f64) -> f64 {
let sigma_t = 6.6524e-29;
let r = distance_rsun * SOLAR_RADIUS;
sigma_t * electron_density * SOLAR_LUMINOSITY / (4.0 * std::f64::consts::PI * r * r)
}
pub fn type_iii_burst_frequency(electron_density: f64) -> f64 {
let omega_pe = (electron_density * ELECTRON_CHARGE * ELECTRON_CHARGE
/ (EPSILON_0 * ELECTRON_MASS))
.sqrt();
omega_pe / (2.0 * std::f64::consts::PI)
}