use sciforge::hub::domain::meteorology::radiation::{solar_constant, stefan_boltzmann_flux};
use sciforge::hub::prelude::constants::{ATM_TO_PASCAL, EARTH_MASS, EARTH_RADIUS, G, R_GAS};
pub const SEA_LEVEL_PRESSURE: f64 = ATM_TO_PASCAL;
pub const SEA_LEVEL_TEMPERATURE: f64 = 288.15;
pub fn dry_air_molar_mass() -> f64 {
crate::m_dry_air()
}
pub const LAPSE_RATE_TROPOSPHERE: f64 = 0.0065;
pub struct AtmosphericLayer {
pub name: &'static str,
pub base_altitude_m: f64,
pub top_altitude_m: f64,
pub base_temperature_k: f64,
pub lapse_rate_k_per_m: f64,
}
impl AtmosphericLayer {
pub fn temperature_at(&self, altitude_m: f64) -> f64 {
let h = (altitude_m - self.base_altitude_m)
.max(0.0)
.min(self.top_altitude_m - self.base_altitude_m);
self.base_temperature_k - self.lapse_rate_k_per_m * h
}
pub fn pressure_at(&self, altitude_m: f64) -> f64 {
let t = self.temperature_at(altitude_m);
let h = altitude_m - self.base_altitude_m;
let g0 = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
if self.lapse_rate_k_per_m.abs() < 1e-10 {
let p_base = barometric_pressure(self.base_altitude_m);
p_base * (-g0 * dry_air_molar_mass() * h / (R_GAS * self.base_temperature_k)).exp()
} else {
let p_base = barometric_pressure(self.base_altitude_m);
p_base
* (t / self.base_temperature_k)
.powf(g0 * dry_air_molar_mass() / (R_GAS * self.lapse_rate_k_per_m))
}
}
pub fn density_at(&self, altitude_m: f64) -> f64 {
let p = self.pressure_at(altitude_m);
let t = self.temperature_at(altitude_m);
p * dry_air_molar_mass() / (R_GAS * t)
}
}
pub fn troposphere() -> AtmosphericLayer {
AtmosphericLayer {
name: "Troposphere",
base_altitude_m: 0.0,
top_altitude_m: 12_000.0,
base_temperature_k: 288.15,
lapse_rate_k_per_m: 0.0065,
}
}
pub fn stratosphere() -> AtmosphericLayer {
AtmosphericLayer {
name: "Stratosphere",
base_altitude_m: 12_000.0,
top_altitude_m: 50_000.0,
base_temperature_k: 216.65,
lapse_rate_k_per_m: -0.001,
}
}
pub fn mesosphere() -> AtmosphericLayer {
AtmosphericLayer {
name: "Mesosphere",
base_altitude_m: 50_000.0,
top_altitude_m: 85_000.0,
base_temperature_k: 270.65,
lapse_rate_k_per_m: 0.0028,
}
}
pub fn thermosphere() -> AtmosphericLayer {
AtmosphericLayer {
name: "Thermosphere",
base_altitude_m: 85_000.0,
top_altitude_m: 600_000.0,
base_temperature_k: 186.87,
lapse_rate_k_per_m: -0.0017,
}
}
pub fn barometric_pressure(altitude_m: f64) -> f64 {
let g0 = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
if altitude_m <= 12_000.0 {
let t = SEA_LEVEL_TEMPERATURE - LAPSE_RATE_TROPOSPHERE * altitude_m;
SEA_LEVEL_PRESSURE
* (t / SEA_LEVEL_TEMPERATURE)
.powf(g0 * dry_air_molar_mass() / (R_GAS * LAPSE_RATE_TROPOSPHERE))
} else {
let p_tropo = barometric_pressure(12_000.0);
let t_tropo = 216.65;
p_tropo * (-g0 * dry_air_molar_mass() * (altitude_m - 12_000.0) / (R_GAS * t_tropo)).exp()
}
}
pub fn scale_height() -> f64 {
let g0 = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
R_GAS * SEA_LEVEL_TEMPERATURE / (dry_air_molar_mass() * g0)
}
pub fn mean_solar_irradiance() -> f64 {
solar_constant()
}
pub fn effective_temperature_k() -> f64 {
let s = solar_constant();
let albedo = 0.306;
stefan_boltzmann_flux(((s * (1.0 - albedo)) / 4.0).powf(0.25))
}