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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

use sciforge::hub::domain::common::constants::K_B;

use crate::{
    CHROMOSPHERE_TEMP_MAX, CHROMOSPHERE_TEMP_MIN, CHROMOSPHERE_THICKNESS, PROTON_MASS,
    SOLAR_SURFACE_GRAVITY,
};

pub const CHROMOSPHERE_BASE_DENSITY: f64 = 1e-6;
pub const CHROMOSPHERE_TOP_DENSITY: f64 = 1e-12;
pub const CHROMOSPHERE_BASE_PRESSURE: f64 = 1.0;
pub const HALPHA_WAVELENGTH_M: f64 = 6.563e-7;
pub const CA_II_K_WAVELENGTH_M: f64 = 3.934e-7;
pub const SPICULE_HEIGHT_M: f64 = 1e7;
pub const SPICULE_SPEED_MS: f64 = 2.5e4;
pub const SPICULE_LIFETIME_S: f64 = 300.0;
pub const SPICULE_WIDTH_M: f64 = 5e5;

pub struct Chromosphere {
    pub base_temperature_k: f64,
    pub top_temperature_k: f64,
    pub base_density_kgm3: f64,
    pub top_density_kgm3: f64,
    pub thickness_m: f64,
}

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

impl Chromosphere {
    pub fn new() -> Self {
        Self {
            base_temperature_k: CHROMOSPHERE_TEMP_MIN,
            top_temperature_k: CHROMOSPHERE_TEMP_MAX,
            base_density_kgm3: CHROMOSPHERE_BASE_DENSITY,
            top_density_kgm3: CHROMOSPHERE_TOP_DENSITY,
            thickness_m: CHROMOSPHERE_THICKNESS,
        }
    }

    pub fn temperature_at_height(&self, height_fraction: f64) -> f64 {
        let f = height_fraction.clamp(0.0, 1.0);
        if f < 0.8 {
            self.base_temperature_k + (6500.0 - self.base_temperature_k) * (f / 0.8)
        } else {
            6500.0 + (self.top_temperature_k - 6500.0) * ((f - 0.8) / 0.2)
        }
    }

    pub fn density_at_height(&self, height_m: f64) -> f64 {
        let scale_height = K_B * 7000.0 / (0.6 * PROTON_MASS * SOLAR_SURFACE_GRAVITY);
        self.base_density_kgm3 * (-height_m / scale_height).exp()
    }

    pub fn pressure_at_height(&self, height_m: f64) -> f64 {
        let rho = self.density_at_height(height_m);
        let t = self.temperature_at_height(height_m / self.thickness_m);
        rho * K_B * t / (0.6 * PROTON_MASS)
    }

    pub fn scale_height(&self, temperature_k: f64) -> f64 {
        K_B * temperature_k / (0.6 * PROTON_MASS * SOLAR_SURFACE_GRAVITY)
    }

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

    pub fn halpha_optical_depth(&self, column_density: f64) -> f64 {
        let sigma_ha = 1e-18;
        column_density * sigma_ha
    }

    pub fn radiative_loss_rate(&self, height_m: f64) -> f64 {
        let n_e = self.density_at_height(height_m) / PROTON_MASS;
        let t = self.temperature_at_height(height_m / self.thickness_m);
        let lambda = 1e-22 * (t / 1e4).powf(-0.5);
        n_e * n_e * lambda
    }

    pub fn acoustic_heating_flux(&self) -> f64 {
        4.5e3
    }

    pub fn mass_column_density(&self) -> f64 {
        self.base_density_kgm3 * self.scale_height(7000.0)
    }
}

pub struct Spicule {
    pub height_m: f64,
    pub velocity_ms: f64,
    pub width_m: f64,
    pub temperature_k: f64,
    pub density_kgm3: f64,
}

impl Default for Spicule {
    fn default() -> Self {
        Self::type_i()
    }
}

impl Spicule {
    pub fn type_i() -> Self {
        Self {
            height_m: SPICULE_HEIGHT_M,
            velocity_ms: SPICULE_SPEED_MS,
            width_m: SPICULE_WIDTH_M,
            temperature_k: 1.5e4,
            density_kgm3: 1e-10,
        }
    }

    pub fn type_ii() -> Self {
        Self {
            height_m: 7e6,
            velocity_ms: 1e5,
            width_m: 3e5,
            temperature_k: 8e4,
            density_kgm3: 5e-11,
        }
    }

    pub fn kinetic_energy(&self) -> f64 {
        let volume = std::f64::consts::PI * (self.width_m / 2.0).powi(2) * self.height_m;
        0.5 * self.density_kgm3 * volume * self.velocity_ms.powi(2)
    }

    pub fn lifetime(&self) -> f64 {
        self.height_m / self.velocity_ms
    }

    pub fn mass_flux(&self) -> f64 {
        self.density_kgm3 * self.velocity_ms
    }
}

pub fn number_of_spicules() -> f64 {
    3e5
}

pub fn total_spicule_mass_flux() -> f64 {
    let s = Spicule::type_i();
    let area = std::f64::consts::PI * (s.width_m / 2.0).powi(2);
    s.mass_flux() * area * number_of_spicules()
}

pub fn transition_region_temperature(height_above_chromo: f64) -> f64 {
    let t_chromo = CHROMOSPHERE_TEMP_MAX;
    let t_corona = 1e6;
    let thickness = crate::TRANSITION_REGION_THICKNESS;
    let f = (height_above_chromo / thickness).clamp(0.0, 1.0);
    t_chromo * ((t_corona / t_chromo).ln() * f).exp()
}