use suns::{events, interactions, layers, physics, plasma};
const DEFAULT_DT_S: f64 = 3600.0;
const STEP_COUNT: usize = 500;
struct SunDiagnostics {
elapsed_days: f64,
core_luminosity_w: f64,
surface_flux_wm2: f64,
surface_temperature_k: f64,
photon_diffusion_speed_ms: f64,
convective_velocity_ms: f64,
chromosphere_scale_height_m: f64,
radiative_zone_volume_m3: f64,
alfven_speed_corona_ms: f64,
coronal_temperature_k: f64,
coronal_plasma_beta: f64,
solar_wind_speed_ms: f64,
solar_wind_density_m3: f64,
parker_spiral_angle_deg: f64,
sunspot_number: f64,
sunspot_latitude_deg: f64,
cme_rate_per_day: f64,
flare_rate_xclass: f64,
prominence_mass_kg: f64,
tidal_force_total_n: f64,
barycenter_offset_m: f64,
magnetic_cycle_phase: f64,
tsi_wm2: f64,
checksum: f64,
}
impl Default for SunDiagnostics {
fn default() -> Self {
Self {
elapsed_days: 0.0,
core_luminosity_w: 0.0,
surface_flux_wm2: 0.0,
surface_temperature_k: 0.0,
photon_diffusion_speed_ms: 0.0,
convective_velocity_ms: 0.0,
chromosphere_scale_height_m: 0.0,
radiative_zone_volume_m3: 0.0,
alfven_speed_corona_ms: 0.0,
coronal_temperature_k: 0.0,
coronal_plasma_beta: 0.0,
solar_wind_speed_ms: 0.0,
solar_wind_density_m3: 0.0,
parker_spiral_angle_deg: 0.0,
sunspot_number: 0.0,
sunspot_latitude_deg: 0.0,
cme_rate_per_day: 0.0,
flare_rate_xclass: 0.0,
prominence_mass_kg: 0.0,
tidal_force_total_n: 0.0,
barycenter_offset_m: 0.0,
magnetic_cycle_phase: 0.0,
tsi_wm2: 0.0,
checksum: 0.0,
}
}
}
impl SunDiagnostics {
fn recompute_checksum(&mut self) {
self.checksum = self.elapsed_days
+ self.core_luminosity_w * 1e-26
+ self.surface_flux_wm2 * 1e-7
+ self.surface_temperature_k * 1e-3
+ self.photon_diffusion_speed_ms * 1e-5
+ self.convective_velocity_ms * 1e-3
+ self.chromosphere_scale_height_m * 1e-6
+ self.radiative_zone_volume_m3 * 1e-30
+ self.alfven_speed_corona_ms * 1e-6
+ self.coronal_temperature_k * 1e-6
+ self.coronal_plasma_beta
+ self.solar_wind_speed_ms * 1e-5
+ self.solar_wind_density_m3 * 1e-6
+ self.parker_spiral_angle_deg * 0.01
+ self.sunspot_number * 0.01
+ self.sunspot_latitude_deg * 0.1
+ self.cme_rate_per_day
+ self.flare_rate_xclass
+ self.prominence_mass_kg * 1e-12
+ self.tidal_force_total_n * 1e-20
+ self.barycenter_offset_m * 1e-9
+ self.magnetic_cycle_phase
+ self.tsi_wm2 * 1e-3;
}
}
struct SunSimulation {
cycle: events::sunspots::SunspotCycle,
diagnostics: SunDiagnostics,
elapsed_s: f64,
}
impl SunSimulation {
fn new() -> Self {
Self {
cycle: events::sunspots::cycle_25(),
diagnostics: SunDiagnostics::default(),
elapsed_s: 0.0,
}
}
fn tick(&mut self, dt_s: f64) {
self.elapsed_s += dt_s;
let cycle_period = suns::SOLAR_CYCLE_YEARS * suns::SECONDS_PER_YEAR;
self.cycle.phase = (self.elapsed_s % cycle_period) / cycle_period;
let photosphere = layers::photosphere::Photosphere::new();
let chromosphere = layers::chromosphere::Chromosphere::new();
let convective = layers::convective_zone::ConvectiveZone::new();
let radiative = layers::radiative_zone::RadiativeZone::new();
let corona = plasma::corona::quiet_corona();
let wind = plasma::solar_wind::SolarWindState::slow_wind_1au();
let prominence = plasma::prominences::Prominence::quiescent();
let rz_transfer = physics::energy_transport::radiative_zone_transfer();
let core_lum = physics::nuclear_fusion::core_luminosity();
let ssn = self.cycle.current_ssn();
let tsi = self.cycle.total_solar_irradiance();
let butterfly_lat = self.cycle.butterfly_latitude_deg();
let cme_rate = events::coronal_mass_ejections::cme_rate_at_cycle_phase(self.cycle.phase);
let tidal_total = interactions::tidal_forces::total_planetary_tidal_force();
let bary_offset = interactions::tidal_forces::solar_barycenter_offset();
self.diagnostics = SunDiagnostics {
elapsed_days: self.elapsed_s / 86_400.0,
core_luminosity_w: core_lum,
surface_flux_wm2: photosphere.surface_flux(),
surface_temperature_k: photosphere.temperature_k,
photon_diffusion_speed_ms: rz_transfer.photon_diffusion_speed(),
convective_velocity_ms: convective.convective_velocity(0.85),
chromosphere_scale_height_m: chromosphere
.scale_height(chromosphere.temperature_at_height(0.5)),
radiative_zone_volume_m3: radiative.volume_m3(),
alfven_speed_corona_ms: corona.alfven_speed(),
coronal_temperature_k: corona.temperature_k,
coronal_plasma_beta: corona.plasma_beta(),
solar_wind_speed_ms: wind.speed_ms,
solar_wind_density_m3: wind.density_m3,
parker_spiral_angle_deg: wind.parker_spiral_angle_deg(),
sunspot_number: ssn,
sunspot_latitude_deg: butterfly_lat,
cme_rate_per_day: cme_rate,
flare_rate_xclass: events::solar_flares::FlareCategory::XClass.frequency_per_day(),
prominence_mass_kg: prominence.mass_kg(),
tidal_force_total_n: tidal_total,
barycenter_offset_m: bary_offset,
magnetic_cycle_phase: self.cycle.phase,
tsi_wm2: tsi,
checksum: 0.0,
};
self.diagnostics.recompute_checksum();
}
}
fn main() {
let mut simulation = SunSimulation::new();
for _ in 0..STEP_COUNT {
simulation.tick(DEFAULT_DT_S);
}
}