suns 0.0.1

Sun celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation
use sciforge::hub::domain::meteorology::radiation::stefan_boltzmann_flux;

use crate::{
    MU_0, SOLAR_EFFECTIVE_TEMP, SOLAR_LUMINOSITY,
    SOLAR_MAGNETIC_FIELD_SUNSPOT_T, SOLAR_RADIUS, SUNSPOT_PENUMBRA_TEMP, SUNSPOT_UMBRA_TEMP,
};

pub const SUNSPOT_MIN_DIAMETER_M: f64 = 2.5e6;
pub const SUNSPOT_MAX_DIAMETER_M: f64 = 5e7;
pub const SUNSPOT_LIFETIME_SCALE_DAYS: f64 = 10.0;
pub const ZURICH_NUMBER_SCALE: f64 = 10.0;
pub const WOLF_NUMBER_K: f64 = 0.6;
pub const BUTTERFLY_LATITUDE_MAX_DEG: f64 = 35.0;
pub const BUTTERFLY_LATITUDE_MIN_DEG: f64 = 5.0;
pub const MAUNDER_MINIMUM_SSN: f64 = 0.0;
pub const MODERN_MAXIMUM_SSN: f64 = 190.0;
pub const JOY_LAW_TILT_COEFF: f64 = 0.5;

pub struct Sunspot {
    pub latitude_deg: f64,
    pub longitude_deg: f64,
    pub diameter_m: f64,
    pub umbra_temperature_k: f64,
    pub penumbra_temperature_k: f64,
    pub magnetic_field_t: f64,
    pub polarity: i8,
}

impl Default for Sunspot {
    fn default() -> Self {
        Self::typical()
    }
}

impl Sunspot {
    pub fn typical() -> Self {
        Self {
            latitude_deg: 15.0,
            longitude_deg: 0.0,
            diameter_m: 2e7,
            umbra_temperature_k: SUNSPOT_UMBRA_TEMP,
            penumbra_temperature_k: SUNSPOT_PENUMBRA_TEMP,
            magnetic_field_t: SOLAR_MAGNETIC_FIELD_SUNSPOT_T,
            polarity: 1,
        }
    }

    pub fn large() -> Self {
        Self {
            latitude_deg: 12.0,
            longitude_deg: 45.0,
            diameter_m: SUNSPOT_MAX_DIAMETER_M,
            umbra_temperature_k: 3500.0,
            penumbra_temperature_k: 4300.0,
            magnetic_field_t: 0.4,
            polarity: -1,
        }
    }

    pub fn umbra_area_m2(&self) -> f64 {
        let r = self.diameter_m * 0.2;
        std::f64::consts::PI * r * r
    }

    pub fn penumbra_area_m2(&self) -> f64 {
        let r_total = self.diameter_m / 2.0;
        let r_umbra = self.diameter_m * 0.2;
        std::f64::consts::PI * (r_total * r_total - r_umbra * r_umbra)
    }

    pub fn total_area_m2(&self) -> f64 {
        let r = self.diameter_m / 2.0;
        std::f64::consts::PI * r * r
    }

    pub fn area_millionths_hemisphere(&self) -> f64 {
        let hemisphere = 2.0 * std::f64::consts::PI * SOLAR_RADIUS * SOLAR_RADIUS;
        self.total_area_m2() / hemisphere * 1e6
    }

    pub fn umbra_flux(&self) -> f64 {
        stefan_boltzmann_flux(self.umbra_temperature_k)
    }

    pub fn penumbra_flux(&self) -> f64 {
        stefan_boltzmann_flux(self.penumbra_temperature_k)
    }

    pub fn photosphere_flux(&self) -> f64 {
        stefan_boltzmann_flux(SOLAR_EFFECTIVE_TEMP)
    }

    pub fn luminosity_deficit(&self) -> f64 {
        let f_photo = self.photosphere_flux();
        let f_umbra = self.umbra_flux();
        let f_penumbra = self.penumbra_flux();
        (f_photo - f_umbra) * self.umbra_area_m2()
            + (f_photo - f_penumbra) * self.penumbra_area_m2()
    }

    pub fn brightness_contrast(&self) -> f64 {
        self.umbra_flux() / self.photosphere_flux()
    }

    pub fn magnetic_pressure(&self) -> f64 {
        self.magnetic_field_t.powi(2) / (2.0 * MU_0)
    }

    pub fn magnetic_energy(&self) -> f64 {
        let volume = self.total_area_m2() * 1e7;
        self.magnetic_field_t.powi(2) * volume / (2.0 * MU_0)
    }

    pub fn wilson_depression_m(&self) -> f64 {
        let p_mag = self.magnetic_pressure();
        let rho = 2e-4;
        let g = 274.0;
        p_mag / (rho * g)
    }

    pub fn evershed_flow_speed(&self) -> f64 {
        2000.0
    }

    pub fn estimated_lifetime_days(&self) -> f64 {
        let area_msh = self.area_millionths_hemisphere();
        0.1 * area_msh
    }

    pub fn joy_law_tilt_deg(&self) -> f64 {
        JOY_LAW_TILT_COEFF * self.latitude_deg.abs()
    }
}

pub struct SunspotCycle {
    pub cycle_number: u32,
    pub phase: f64,
    pub ssn_max: f64,
}

impl SunspotCycle {
    pub fn current_ssn(&self) -> f64 {
        self.ssn_max * (std::f64::consts::PI * self.phase).sin().max(0.0)
    }

    pub fn activity_level(&self) -> &'static str {
        let ssn = self.current_ssn();
        if ssn < 20.0 {
            "Minimum"
        } else if ssn < 80.0 {
            "Low"
        } else if ssn < 150.0 {
            "Moderate"
        } else {
            "Maximum"
        }
    }

    pub fn butterfly_latitude_deg(&self) -> f64 {
        let max = BUTTERFLY_LATITUDE_MAX_DEG;
        let min = BUTTERFLY_LATITUDE_MIN_DEG;
        max - (max - min) * self.phase
    }

    pub fn f107_index(&self) -> f64 {
        67.0 + 0.572 * self.current_ssn()
    }

    pub fn tsi_variation_wm2(&self) -> f64 {
        0.1 * (self.current_ssn() / 100.0)
    }

    pub fn total_solar_irradiance(&self) -> f64 {
        1361.0 + self.tsi_variation_wm2()
    }
}

pub fn wolf_number(groups: u32, individual_spots: u32) -> f64 {
    WOLF_NUMBER_K * (ZURICH_NUMBER_SCALE * groups as f64 + individual_spots as f64)
}

pub fn sunspot_number_to_area_msh(ssn: f64) -> f64 {
    16.7 * ssn
}

pub fn maunder_minimum_cycle() -> SunspotCycle {
    SunspotCycle {
        cycle_number: 0,
        phase: 0.5,
        ssn_max: MAUNDER_MINIMUM_SSN,
    }
}

pub fn cycle_25() -> SunspotCycle {
    SunspotCycle {
        cycle_number: 25,
        phase: 0.45,
        ssn_max: 150.0,
    }
}

pub fn sporer_law_latitude(phase: f64) -> f64 {
    BUTTERFLY_LATITUDE_MAX_DEG * (1.0 - phase)
}

pub fn active_region_probability(ssn: f64) -> f64 {
    1.0 - (-ssn / 100.0).exp()
}

pub fn facular_brightening(ssn: f64) -> f64 {
    0.0011 * ssn
}

pub fn total_luminosity_variation(ssn: f64) -> f64 {
    let spot_dimming = -0.0003 * ssn;
    let facular = facular_brightening(ssn);
    SOLAR_LUMINOSITY * (1.0 + spot_dimming + facular)
}