suns 0.0.4

Sun celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation
use suns::layers::chromosphere::{
    Chromosphere, Spicule, number_of_spicules, total_spicule_mass_flux,
};
use suns::layers::convective_zone::{ConvectiveZone, differential_rotation_angular_velocity};
use suns::layers::photosphere::Photosphere;
use suns::layers::radiative_zone::{RadiativeZone, Tachocline};
use suns::physics::magnetic_fields::{
    active_region, coronal_field, differential_rotation_period_days, quiet_sun, sunspot_umbra,
};
use suns::physics::stellar_structure::{
    convective_zone, core, dynamical_timescale, gravitational_binding_energy,
    kelvin_helmholtz_timescale, luminosity_from_temperature, nuclear_timescale, radiative_zone,
    surface_gravity,
};
use suns::plasma::corona::{
    active_region_corona, coronal_hole, coronal_streamer, density_at_height,
    hydrostatic_scale_height, quiet_corona,
};
use suns::plasma::prominences::Prominence;
use suns::plasma::solar_wind::{SolarWindState, parker_critical_radius, solar_wind_mass_loss_rate};

fn main() {
    println!("=== CORONA MAPPING & SOLAR OVERVIEW ===\n");

    // --- Stellar Structure ---
    println!("--- Stellar Structure ---");
    let c = core();
    let rz = radiative_zone();
    let cz = convective_zone();
    for layer in [&c, &rz, &cz] {
        println!(
            "  {:>16}: T_mean={:.2e}K  ρ_mean={:.2e}kg/m³  cs={:.0}m/s",
            layer.name,
            layer.mean_temperature(),
            layer.mean_density(),
            layer.sound_speed()
        );
    }

    // --- Timescales ---
    println!("\n--- Fundamental Timescales ---");
    println!("  Dynamical:  {:.0} s", dynamical_timescale());
    println!(
        "  K-H:        {:.2e} yr",
        kelvin_helmholtz_timescale() / 3.15576e7
    );
    println!("  Nuclear:    {:.2e} yr", nuclear_timescale() / 3.15576e7);

    // --- Corona Regions ---
    println!("\n--- Corona Regions ---");
    let regions = [
        quiet_corona(),
        active_region_corona(),
        coronal_hole(),
        coronal_streamer(),
    ];
    for r in &regions {
        println!(
            "  {:>20}: T={:.1e}K  n_e={:.1e}m⁻³  β={:.2e}  v_A={:.0}km/s",
            r.name,
            r.temperature_k,
            r.electron_density_m3,
            r.plasma_beta(),
            r.alfven_speed() / 1e3
        );
    }

    // --- Coronal Density Profile ---
    println!("\n--- Coronal Density Profile ---");
    let base_n = 1e14;
    let temp = 1.5e6;
    let h_scale = hydrostatic_scale_height(temp);
    println!(
        "  Scale height @ {:.1e}K: {:.2e} m ({:.1} R_sun)",
        temp,
        h_scale,
        h_scale / suns::SOLAR_RADIUS
    );
    for height_rsun in [0.0, 0.1, 0.5, 1.0, 2.0, 5.0] {
        let h = height_rsun * suns::SOLAR_RADIUS;
        println!(
            "  h={:.1}R: n={:.2e} m⁻³",
            height_rsun,
            density_at_height(base_n, h, temp)
        );
    }

    // --- Magnetic Regions ---
    println!("\n--- Magnetic Regions ---");
    for m in [
        quiet_sun(),
        active_region(),
        sunspot_umbra(),
        coronal_field(),
    ] {
        println!(
            "  {:>14}: B={:.1e}T  β={:.2e}  v_A={:.0}km/s",
            m.name,
            m.field_strength_t,
            m.plasma_beta(),
            m.alfven_speed() / 1e3
        );
    }

    // --- Differential Rotation ---
    println!("\n--- Differential Rotation ---");
    for lat in [0.0, 15.0, 30.0, 45.0, 60.0, 75.0] {
        println!(
            "  lat={:>5.1}°: P={:.2}d  Ω={:.2e}rad/s",
            lat,
            differential_rotation_period_days(lat),
            differential_rotation_angular_velocity(lat)
        );
    }

    // --- Solar Layers ---
    println!("\n--- Solar Layers ---");
    let photo = Photosphere::new();
    let chromo = Chromosphere::new();
    let conv = ConvectiveZone::new();
    let rad = RadiativeZone::new();
    println!(
        "  Photosphere:  T={:.0}K  λ_peak={:.0}nm  cs={:.0}m/s",
        photo.temperature_k,
        photo.spectral_peak_wavelength_m() * 1e9,
        photo.sound_speed()
    );
    println!(
        "  Chromosphere: T_base={:.0}K  T_top={:.0}K  Δh={:.0}km",
        chromo.base_temperature_k,
        chromo.top_temperature_k,
        chromo.thickness_m / 1e3
    );
    println!(
        "  Convective:   Δr={:.0}Mm  BV(0.85)={:.2e}Hz",
        conv.thickness_m() / 1e6,
        conv.brunt_vaisala_frequency(0.85)
    );
    println!(
        "  Radiative:    Δr={:.0}Mm  t_photon={:.2e}yr",
        rad.thickness_m() / 1e6,
        rad.photon_random_walk_time() / 3.15576e7
    );

    // --- Tachocline ---
    let tach = Tachocline::new();
    println!("\n--- Tachocline ---");
    println!(
        "  Thickness: {:.0} km  Shear: {:.2e} s⁻¹/m  τ_diff: {:.2e} yr",
        tach.thickness_m / 1e3,
        tach.shear_rate(),
        tach.diffusion_time() / 3.15576e7
    );

    // --- Spicules ---
    println!("\n--- Spicules ---");
    let sp1 = Spicule::type_i();
    let sp2 = Spicule::type_ii();
    println!(
        "  Type-I: h={:.0}km  v={:.0}km/s  τ={:.0}s",
        sp1.height_m / 1e3,
        sp1.velocity_ms / 1e3,
        sp1.lifetime()
    );
    println!(
        "  Type-II: h={:.0}km  v={:.0}km/s  τ={:.0}s",
        sp2.height_m / 1e3,
        sp2.velocity_ms / 1e3,
        sp2.lifetime()
    );
    println!("  Active spicules: {:.2e}", number_of_spicules());
    println!("  Total mass flux: {:.2e} kg/s", total_spicule_mass_flux());

    // --- Prominences ---
    println!("\n--- Prominences ---");
    let pq = Prominence::quiescent();
    let pa = Prominence::active();
    println!(
        "  Quiescent: m={:.1e}kg  E_th={:.1e}J  E_mag={:.1e}J  KS={:.2e}",
        pq.mass_kg(),
        pq.thermal_energy(),
        pq.magnetic_energy(),
        pq.kippenhahn_schluter_support()
    );
    println!(
        "  Active:    m={:.1e}kg  v_erupt={:.0}km/s  τ_cool={:.0}s",
        pa.mass_kg(),
        pa.eruption_speed() / 1e3,
        pa.cooling_time()
    );

    // --- Solar Wind ---
    println!("\n--- Solar Wind ---");
    let slow = SolarWindState::slow_wind_1au();
    let fast = SolarWindState::fast_wind_1au();
    for (name, w) in [("Slow", &slow), ("Fast", &fast)] {
        println!(
            "  {} @1AU: v={:.0}km/s  n={:.1e}m⁻³  T={:.0}K  B={:.1}nT  spiral={:.1}°",
            name,
            w.speed_ms / 1e3,
            w.density_m3,
            w.temperature_k,
            w.magnetic_field_t * 1e9,
            w.parker_spiral_angle_deg()
        );
    }
    println!("  Mass loss: {:.2e} kg/s", solar_wind_mass_loss_rate());
    println!(
        "  r_critical: {:.2} R_sun",
        parker_critical_radius() / suns::SOLAR_RADIUS
    );

    // --- Global ---
    println!("\n--- Global Solar Properties ---");
    println!("  L(T_eff):  {:.4e} W", luminosity_from_temperature());
    println!("  g:         {:.1} m/s²", surface_gravity());
    println!("  E_bind:    {:.3e} J", gravitational_binding_energy());
    println!(
        "  v_esc:     {:.1} km/s",
        *suns::SOLAR_ESCAPE_VELOCITY / 1e3
    );

    println!("\n=== CORONA MAPPING COMPLETE ===");
}