suns 0.0.1

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
159
160
161

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

use crate::{
    FLARE_CLASS_C_FLUX, FLARE_CLASS_M_FLUX,
    FLARE_CLASS_X_FLUX, MU_0,
};

pub const FLARE_DURATION_IMPULSIVE_S: f64 = 60.0;
pub const FLARE_DURATION_GRADUAL_S: f64 = 3600.0;
pub const FLARE_RISE_TIME_S: f64 = 300.0;
pub const FLARE_RECONNECTION_RATE: f64 = 0.1;
pub const NONTHERMAL_ELECTRON_FRACTION: f64 = 0.1;
pub const FLARE_LOOP_LENGTH_M: f64 = 5e7;
pub const FLARE_LOOP_FOOTPOINT_AREA_M2: f64 = 1e14;

pub struct SolarFlare {
    pub classification: FlareClass,
    pub peak_flux_wm2: f64,
    pub duration_s: f64,
    pub loop_length_m: f64,
    pub magnetic_field_t: f64,
    pub area_m2: f64,
}

#[derive(Clone, Copy, PartialEq)]
pub enum FlareClass {
    A,
    B,
    C,
    M,
    X,
}

impl FlareClass {
    pub fn base_flux(&self) -> f64 {
        match self {
            FlareClass::A => 1e-8,
            FlareClass::B => 1e-7,
            FlareClass::C => FLARE_CLASS_C_FLUX,
            FlareClass::M => FLARE_CLASS_M_FLUX,
            FlareClass::X => FLARE_CLASS_X_FLUX,
        }
    }

    pub fn label(&self) -> &'static str {
        match self {
            FlareClass::A => "A",
            FlareClass::B => "B",
            FlareClass::C => "C",
            FlareClass::M => "M",
            FlareClass::X => "X",
        }
    }
}

impl SolarFlare {
    pub fn from_class(class: FlareClass, multiplier: f64) -> Self {
        let base = class.base_flux();
        let peak = base * multiplier;
        let (duration, b, area) = match class {
            FlareClass::A | FlareClass::B => (300.0, 0.01, 1e13),
            FlareClass::C => (600.0, 0.05, 5e13),
            FlareClass::M => (1800.0, 0.1, 1e14),
            FlareClass::X => (3600.0, 0.3, 5e14),
        };
        Self {
            classification: class,
            peak_flux_wm2: peak,
            duration_s: duration,
            loop_length_m: FLARE_LOOP_LENGTH_M,
            magnetic_field_t: b,
            area_m2: area,
        }
    }

    pub fn total_energy_joules(&self) -> f64 {
        let volume = self.area_m2 * self.loop_length_m;
        self.magnetic_field_t.powi(2) * volume / (2.0 * MU_0)
    }

    pub fn peak_luminosity(&self) -> f64 {
        self.peak_flux_wm2 * 4.0 * std::f64::consts::PI * (1.496e11_f64).powi(2)
    }

    pub fn peak_temperature(&self) -> f64 {
        let e = self.total_energy_joules();
        let volume = self.area_m2 * self.loop_length_m;
        let n_e = 1e16;
        e / (3.0 * n_e * K_B * volume)
    }

    pub fn reconnection_inflow_speed(&self) -> f64 {
        let rho = 1e-11;
        let v_a = self.magnetic_field_t / (MU_0 * rho).sqrt();
        FLARE_RECONNECTION_RATE * v_a
    }

    pub fn reconnection_electric_field(&self) -> f64 {
        self.reconnection_inflow_speed() * self.magnetic_field_t
    }

    pub fn nonthermal_electron_energy(&self) -> f64 {
        NONTHERMAL_ELECTRON_FRACTION * self.total_energy_joules()
    }

    pub fn hard_xray_flux_estimate(&self) -> f64 {
        0.01 * self.peak_flux_wm2
    }

    pub fn chromospheric_evaporation_speed(&self) -> f64 {
        let energy_flux = self.total_energy_joules() / (self.area_m2 * self.duration_s);
        let rho_chromo = 1e-8;
        (2.0 * energy_flux / rho_chromo).powf(1.0 / 3.0)
    }

    pub fn confinement_time(&self) -> f64 {
        let n_e = 1e16;
        let volume = self.area_m2 * self.loop_length_m;
        3.0 * n_e * K_B * self.peak_temperature() * volume / self.peak_luminosity().max(1.0)
    }

    pub fn neupert_effect_ratio(&self) -> f64 {
        self.peak_flux_wm2 * self.duration_s
    }
}

pub fn flare_frequency_per_day(class: FlareClass) -> f64 {
    match class {
        FlareClass::A => 100.0,
        FlareClass::B => 30.0,
        FlareClass::C => 8.0,
        FlareClass::M => 1.0,
        FlareClass::X => 0.1,
    }
}

pub fn goes_class_from_flux(flux_wm2: f64) -> FlareClass {
    if flux_wm2 >= FLARE_CLASS_X_FLUX {
        FlareClass::X
    } else if flux_wm2 >= FLARE_CLASS_M_FLUX {
        FlareClass::M
    } else if flux_wm2 >= FLARE_CLASS_C_FLUX {
        FlareClass::C
    } else if flux_wm2 >= 1e-7 {
        FlareClass::B
    } else {
        FlareClass::A
    }
}

pub fn magnetic_energy_available(field_t: f64, volume_m3: f64) -> f64 {
    field_t.powi(2) * volume_m3 / (2.0 * MU_0)
}

pub fn flare_ribbon_speed(reconnection_rate: f64, alfven_speed: f64) -> f64 {
    reconnection_rate * alfven_speed
}

pub fn white_light_flare_contrast(flare_luminosity: f64, photosphere_flux: f64) -> f64 {
    flare_luminosity / (photosphere_flux * FLARE_LOOP_FOOTPOINT_AREA_M2)
}