suns 0.0.3

Sun celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

use sciforge::hub::domain::common::constants::{K_B, SIGMA_SB};
use sciforge::hub::domain::meteorology::radiation::stefan_boltzmann_flux;

use crate::{
    PHOTOSPHERE_THICKNESS, PROTON_MASS, SOLAR_EFFECTIVE_TEMP, SOLAR_RADIUS, SPEED_OF_LIGHT,
};

pub const PHOTOSPHERE_BASE_DENSITY: f64 = 2e-4;
pub const PHOTOSPHERE_BASE_PRESSURE: f64 = 1.25e4;
pub const PHOTOSPHERE_SCALE_HEIGHT: f64 = 1.5e5;
pub const GRANULATION_CELL_SIZE_M: f64 = 1e6;
pub const GRANULATION_VELOCITY_MS: f64 = 1000.0;
pub const SUPERGRANULATION_CELL_SIZE_M: f64 = 3e7;
pub const SUPERGRANULATION_VELOCITY_MS: f64 = 300.0;
pub const PHOTOSPHERIC_OPACITY_CM2G: f64 = 0.264;
pub const LIMB_DARKENING_COEFF_U: f64 = 0.6;

pub struct Photosphere {
    pub temperature_k: f64,
    pub density_kgm3: f64,
    pub pressure_pa: f64,
    pub thickness_m: f64,
}

impl Default for Photosphere {
    fn default() -> Self {
        Self::new()
    }
}

impl Photosphere {
    pub fn new() -> Self {
        Self {
            temperature_k: SOLAR_EFFECTIVE_TEMP,
            density_kgm3: PHOTOSPHERE_BASE_DENSITY,
            pressure_pa: PHOTOSPHERE_BASE_PRESSURE,
            thickness_m: PHOTOSPHERE_THICKNESS,
        }
    }

    pub fn surface_flux(&self) -> f64 {
        stefan_boltzmann_flux(self.temperature_k)
    }

    pub fn luminosity(&self) -> f64 {
        4.0 * std::f64::consts::PI * SOLAR_RADIUS * SOLAR_RADIUS * self.surface_flux()
    }

    pub fn opacity(&self) -> f64 {
        PHOTOSPHERIC_OPACITY_CM2G * 0.1 * self.density_kgm3
    }

    pub fn optical_depth_at_height(&self, height_m: f64) -> f64 {
        let kappa = self.opacity();
        kappa
            * self.density_kgm3
            * (-height_m / PHOTOSPHERE_SCALE_HEIGHT).exp()
            * PHOTOSPHERE_SCALE_HEIGHT
    }

    pub fn temperature_at_optical_depth(&self, tau: f64) -> f64 {
        self.temperature_k * (0.75 * (tau + 2.0 / 3.0)).powf(0.25)
    }

    pub fn limb_darkening_intensity(&self, theta: f64) -> f64 {
        1.0 - LIMB_DARKENING_COEFF_U * (1.0 - theta.cos())
    }

    pub fn density_at_height(&self, height_m: f64) -> f64 {
        self.density_kgm3 * (-height_m / PHOTOSPHERE_SCALE_HEIGHT).exp()
    }

    pub fn pressure_at_height(&self, height_m: f64) -> f64 {
        self.pressure_pa * (-height_m / PHOTOSPHERE_SCALE_HEIGHT).exp()
    }

    pub fn sound_speed(&self) -> f64 {
        let gamma = 5.0 / 3.0;
        let mu = 0.61;
        (gamma * K_B * self.temperature_k / (mu * PROTON_MASS)).sqrt()
    }

    pub fn thermal_velocity(&self) -> f64 {
        (2.0 * K_B * self.temperature_k / PROTON_MASS).sqrt()
    }

    pub fn spectral_peak_wavelength_m(&self) -> f64 {
        2.898e-3 / self.temperature_k
    }

    pub fn radiation_pressure(&self) -> f64 {
        4.0 * SIGMA_SB * self.temperature_k.powi(4) / (3.0 * SPEED_OF_LIGHT)
    }

    pub fn mean_molecular_weight(&self) -> f64 {
        let mu_ion = 1.26;
        let ionization_frac = 0.5;
        mu_ion / (1.0 + ionization_frac)
    }
}

pub fn granulation_lifetime() -> f64 {
    GRANULATION_CELL_SIZE_M / GRANULATION_VELOCITY_MS
}

pub fn supergranulation_lifetime() -> f64 {
    SUPERGRANULATION_CELL_SIZE_M / SUPERGRANULATION_VELOCITY_MS
}

pub fn effective_temperature_from_luminosity(luminosity: f64) -> f64 {
    (luminosity / (4.0 * std::f64::consts::PI * SOLAR_RADIUS * SOLAR_RADIUS * SIGMA_SB)).powf(0.25)
}

pub fn spectral_irradiance_planck(wavelength_m: f64, temperature_k: f64) -> f64 {
    let h = crate::PLANCK_H;
    let c = SPEED_OF_LIGHT;
    let k = K_B;
    let a = 2.0 * h * c * c / wavelength_m.powi(5);
    let b = (h * c / (wavelength_m * k * temperature_k)).exp() - 1.0;
    a / b
}

pub fn bolometric_correction(temperature_k: f64) -> f64 {
    let log_t = temperature_k.log10();
    if log_t < 3.7 {
        -0.4 * (log_t - 3.7)
    } else {
        -0.09 * (log_t - 3.76).powi(2)
    }
}

pub fn wilson_depression_depth(field_t: f64) -> f64 {
    let beta = crate::MU_0;
    let p_mag = field_t.powi(2) / (2.0 * beta);
    p_mag / (PHOTOSPHERE_BASE_DENSITY * crate::SOLAR_SURFACE_GRAVITY)
}