use sciforge::hub::prelude::constants::elements::atomic_mass;
use sciforge::hub::prelude::constants::{EARTH_RADIUS, N_A};
pub struct MolecularSpecies {
pub name: &'static str,
pub symbol: &'static str,
pub molar_mass_kg_mol: f64,
pub volume_fraction: f64,
pub refractive_index_stp: f64,
pub depolarization_factor: f64,
}
pub struct AtmosphereEndpoint {
pub planet_radius_m: f64,
pub atmosphere_height_m: f64,
pub sea_level_pressure_pa: f64,
pub sea_level_temperature_k: f64,
pub sea_level_number_density_m3: f64,
pub mean_molar_mass_kg_mol: f64,
pub rayleigh_scale_height_m: f64,
pub mie_scale_height_m: f64,
pub mie_asymmetry_g: f64,
pub mie_coefficient: f64,
pub ozone_peak_altitude_m: f64,
pub ozone_column_du: f64,
pub rayleigh_coefficients_rgb: [f64; 3],
pub species: Vec<MolecularSpecies>,
pub sun_irradiance_w_m2: f64,
}
fn n2_molar() -> f64 {
2.0 * atomic_mass(7) * 1e-3
}
fn o2_molar() -> f64 {
2.0 * atomic_mass(8) * 1e-3
}
fn ar_molar() -> f64 {
atomic_mass(18) * 1e-3
}
fn co2_molar() -> f64 {
(atomic_mass(6) + 2.0 * atomic_mass(8)) * 1e-3
}
fn ne_molar() -> f64 {
atomic_mass(10) * 1e-3
}
fn he_molar() -> f64 {
atomic_mass(2) * 1e-3
}
fn ch4_molar() -> f64 {
(atomic_mass(6) + 4.0 * atomic_mass(1)) * 1e-3
}
fn h2o_molar() -> f64 {
(2.0 * atomic_mass(1) + atomic_mass(8)) * 1e-3
}
fn mean_molar_mass() -> f64 {
0.780_84 * n2_molar()
+ 0.209_46 * o2_molar()
+ 0.009_34 * ar_molar()
+ 0.000_417 * co2_molar()
+ 0.000_018_18 * ne_molar()
+ 0.000_005_24 * he_molar()
+ 0.000_001_87 * ch4_molar()
}
fn sea_level_number_density() -> f64 {
let m_air = mean_molar_mass();
let rho = 101_325.0 / (8.314_462_618 / m_air * 288.15);
N_A * rho / m_air
}
fn rayleigh_beta_rgb(n_density: f64) -> [f64; 3] {
let n_minus_1 = 2.9e-4;
let ns2 = (2.0 * n_minus_1) * (2.0 * n_minus_1);
let coeff = 8.0 * std::f64::consts::PI.powi(3) * ns2 / (3.0 * n_density);
[
coeff / (680e-9_f64).powi(4),
coeff / (550e-9_f64).powi(4),
coeff / (440e-9_f64).powi(4),
]
}
impl AtmosphereEndpoint {
pub fn earth() -> Self {
let n_density = sea_level_number_density();
let beta = rayleigh_beta_rgb(n_density);
let species = vec![
MolecularSpecies {
name: "Dinitrogen",
symbol: "N2",
molar_mass_kg_mol: n2_molar(),
volume_fraction: 0.780_84,
refractive_index_stp: 1.000_298_0,
depolarization_factor: 0.030,
},
MolecularSpecies {
name: "Dioxygen",
symbol: "O2",
molar_mass_kg_mol: o2_molar(),
volume_fraction: 0.209_46,
refractive_index_stp: 1.000_271_0,
depolarization_factor: 0.054,
},
MolecularSpecies {
name: "Argon",
symbol: "Ar",
molar_mass_kg_mol: ar_molar(),
volume_fraction: 0.009_34,
refractive_index_stp: 1.000_281_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Carbon dioxide",
symbol: "CO2",
molar_mass_kg_mol: co2_molar(),
volume_fraction: 4.17e-4,
refractive_index_stp: 1.000_450_0,
depolarization_factor: 0.075,
},
MolecularSpecies {
name: "Neon",
symbol: "Ne",
molar_mass_kg_mol: ne_molar(),
volume_fraction: 1.818e-5,
refractive_index_stp: 1.000_067_1,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Helium",
symbol: "He",
molar_mass_kg_mol: he_molar(),
volume_fraction: 5.24e-6,
refractive_index_stp: 1.000_035_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Methane",
symbol: "CH4",
molar_mass_kg_mol: ch4_molar(),
volume_fraction: 1.87e-6,
refractive_index_stp: 1.000_444_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Water vapor",
symbol: "H2O",
molar_mass_kg_mol: h2o_molar(),
volume_fraction: 0.01,
refractive_index_stp: 1.000_256_0,
depolarization_factor: 0.17,
},
];
Self {
planet_radius_m: EARTH_RADIUS,
atmosphere_height_m: 100_000.0,
sea_level_pressure_pa: 101_325.0,
sea_level_temperature_k: 288.15,
sea_level_number_density_m3: n_density,
mean_molar_mass_kg_mol: mean_molar_mass(),
rayleigh_scale_height_m: 8_500.0,
mie_scale_height_m: 1_200.0,
mie_asymmetry_g: 0.76,
mie_coefficient: 21e-6,
ozone_peak_altitude_m: 25_000.0,
ozone_column_du: 300.0,
rayleigh_coefficients_rgb: beta,
species,
sun_irradiance_w_m2: 1_361.0,
}
}
pub fn sky_luminance(&self, observer_alt: f64, cos_zenith: f64) -> f64 {
let cos_z = cos_zenith.max(0.01);
let molecular_rayleigh: f64 = self
.species
.iter()
.map(|s| {
let n_s = self.sea_level_number_density_m3 * s.volume_fraction;
let delta = s.depolarization_factor;
let king = (6.0 + 3.0 * delta) / (6.0 - 7.0 * delta);
let ns = s.refractive_index_stp - 1.0;
8.0 * std::f64::consts::PI.powi(3) * (2.0 * ns).powi(2) * king / (3.0 * n_s)
})
.sum();
let altitude_factor = (-observer_alt / self.rayleigh_scale_height_m).exp();
let mie_factor = (-observer_alt / self.mie_scale_height_m).exp() * self.mie_coefficient;
let g = self.mie_asymmetry_g;
let hg_phase = (1.0 - g * g)
/ (4.0 * std::f64::consts::PI * (1.0 + g * g - 2.0 * g * cos_z).powf(1.5));
let ozone_abs = self.ozone_column_du
* 1e-5
* (-(observer_alt - self.ozone_peak_altitude_m).powi(2) / 1e9).exp();
let scatter_550 = molecular_rayleigh / (550e-9_f64).powi(4) * altitude_factor;
let sky_r = self.rayleigh_coefficients_rgb[0] * cos_z * altitude_factor * 1e4;
let sky_g = self.rayleigh_coefficients_rgb[1] * cos_z * altitude_factor * 1e4;
let sky_b = self.rayleigh_coefficients_rgb[2] * cos_z * altitude_factor * 1e4;
sky_r * 0.3
+ sky_g * 0.59
+ sky_b * 0.11
+ scatter_550 * 1e-20
+ mie_factor * 1e-6
+ ozone_abs * 1e-8
+ hg_phase * mie_factor * 1e-4
}
pub fn barometric_density(&self) -> f64 {
self.sea_level_pressure_pa
/ (self.sea_level_temperature_k * 8.314_462_618 / self.mean_molar_mass_kg_mol)
}
pub fn shell_volume(&self) -> f64 {
4.0 / 3.0
* std::f64::consts::PI
* ((self.planet_radius_m + self.atmosphere_height_m).powi(3)
- self.planet_radius_m.powi(3))
}
}