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

use sciforge::hub::domain::common::constants::{G, K_B, SOLAR_MASS};

use crate::{
    CONVECTIVE_ZONE_BASE_FRAC, PROTON_MASS, SECONDS_PER_DAY, SOLAR_EFFECTIVE_TEMP,
    SOLAR_OMEGA_EQUATORIAL, SOLAR_RADIUS,
};

pub const CONVECTIVE_ZONE_BASE_TEMP: f64 = 2.0e6;
pub const CONVECTIVE_ZONE_TOP_TEMP: f64 = 5778.0;
pub const CONVECTIVE_ZONE_BASE_DENSITY: f64 = 200.0;
pub const CONVECTIVE_ZONE_TOP_DENSITY: f64 = 2e-4;
pub const MIXING_LENGTH_ALPHA: f64 = 1.7;
pub const CONVECTIVE_VELOCITY_DEEP: f64 = 100.0;
pub const CONVECTIVE_VELOCITY_SURFACE: f64 = 2000.0;
pub const CONVECTIVE_TURNOVER_TIME_S: f64 = 2.59e6;
pub const ROSSBY_NUMBER_SOLAR: f64 = 1.96;

pub struct ConvectiveZone {
    pub inner_radius_frac: f64,
    pub outer_radius_frac: f64,
    pub base_temperature_k: f64,
    pub top_temperature_k: f64,
    pub base_density: f64,
    pub top_density: f64,
}

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

impl ConvectiveZone {
    pub fn new() -> Self {
        Self {
            inner_radius_frac: CONVECTIVE_ZONE_BASE_FRAC,
            outer_radius_frac: 1.0,
            base_temperature_k: CONVECTIVE_ZONE_BASE_TEMP,
            top_temperature_k: CONVECTIVE_ZONE_TOP_TEMP,
            base_density: CONVECTIVE_ZONE_BASE_DENSITY,
            top_density: CONVECTIVE_ZONE_TOP_DENSITY,
        }
    }

    pub fn thickness_m(&self) -> f64 {
        (self.outer_radius_frac - self.inner_radius_frac) * SOLAR_RADIUS
    }

    pub fn volume_m3(&self) -> f64 {
        let r_out = self.outer_radius_frac * SOLAR_RADIUS;
        let r_in = self.inner_radius_frac * SOLAR_RADIUS;
        (4.0 / 3.0) * std::f64::consts::PI * (r_out.powi(3) - r_in.powi(3))
    }

    pub fn temperature_at_frac(&self, radius_frac: f64) -> f64 {
        let f = ((radius_frac - self.inner_radius_frac)
            / (self.outer_radius_frac - self.inner_radius_frac))
            .clamp(0.0, 1.0);
        self.base_temperature_k * (1.0 - f) + self.top_temperature_k * f
    }

    pub fn density_at_frac(&self, radius_frac: f64) -> f64 {
        let f = ((radius_frac - self.inner_radius_frac)
            / (self.outer_radius_frac - self.inner_radius_frac))
            .clamp(0.0, 1.0);
        (self.base_density.ln() * (1.0 - f) + self.top_density.ln() * f).exp()
    }

    pub fn pressure_scale_height(&self, radius_frac: f64) -> f64 {
        let t = self.temperature_at_frac(radius_frac);
        let r = radius_frac * SOLAR_RADIUS;
        let g = G * SOLAR_MASS / (r * r);
        K_B * t / (0.61 * PROTON_MASS * g)
    }

    pub fn mixing_length(&self, radius_frac: f64) -> f64 {
        MIXING_LENGTH_ALPHA * self.pressure_scale_height(radius_frac)
    }

    pub fn convective_velocity(&self, radius_frac: f64) -> f64 {
        let f = ((radius_frac - self.inner_radius_frac)
            / (self.outer_radius_frac - self.inner_radius_frac))
            .clamp(0.0, 1.0);
        CONVECTIVE_VELOCITY_DEEP * (1.0 - f) + CONVECTIVE_VELOCITY_SURFACE * f
    }

    pub fn convective_flux(&self, radius_frac: f64) -> f64 {
        let rho = self.density_at_frac(radius_frac);
        let v = self.convective_velocity(radius_frac);
        let cp = 5.0 * K_B / (2.0 * 0.61 * PROTON_MASS);
        let dt = 100.0;
        rho * cp * v * dt
    }

    pub fn turnover_time(&self, radius_frac: f64) -> f64 {
        let l = self.mixing_length(radius_frac);
        let v = self.convective_velocity(radius_frac);
        l / v.max(1e-10)
    }

    pub fn rossby_number(&self, radius_frac: f64) -> f64 {
        let tau = self.turnover_time(radius_frac);
        let omega = *SOLAR_OMEGA_EQUATORIAL;
        2.0 * std::f64::consts::PI / (omega * tau)
    }

    pub fn superadiabaticity(&self) -> f64 {
        1e-7
    }

    pub fn brunt_vaisala_frequency(&self, radius_frac: f64) -> f64 {
        let g = G * SOLAR_MASS / ((radius_frac * SOLAR_RADIUS).powi(2));
        let hp = self.pressure_scale_height(radius_frac);
        let delta = self.superadiabaticity();
        (g * delta / hp).sqrt()
    }

    pub fn total_kinetic_energy(&self) -> f64 {
        let mean_v = (CONVECTIVE_VELOCITY_DEEP + CONVECTIVE_VELOCITY_SURFACE) / 2.0;
        let mean_rho = (self.base_density + self.top_density) / 2.0;
        0.5 * mean_rho * mean_v * mean_v * self.volume_m3()
    }
}

pub fn differential_rotation_angular_velocity(latitude_deg: f64) -> f64 {
    let a = 14.713;
    let b = -2.396;
    let c = -1.787;
    let sin_lat = latitude_deg.to_radians().sin();
    let sin2 = sin_lat * sin_lat;
    (a + b * sin2 + c * sin2 * sin2) * std::f64::consts::PI / 180.0 / SECONDS_PER_DAY
}

pub fn meridional_flow_speed(latitude_deg: f64) -> f64 {
    let max_speed = 15.0;
    max_speed * (2.0 * latitude_deg.to_radians()).sin()
}

pub fn convective_overshoot_depth() -> f64 {
    0.05 * SOLAR_RADIUS
}

pub fn entropy_rain_speed(temperature_deficit: f64) -> f64 {
    let g = G * SOLAR_MASS / (SOLAR_RADIUS * SOLAR_RADIUS);
    let hp = 1.5e5;
    (g * temperature_deficit * hp / SOLAR_EFFECTIVE_TEMP).sqrt()
}