use sciforge::hub::prelude::constants::N_A;
use sciforge::hub::prelude::constants::elements::atomic_mass;
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 dust_optical_depth: f64,
pub rayleigh_coefficients_rgb: [f64; 3],
pub species: Vec<MolecularSpecies>,
pub sun_irradiance_w_m2: f64,
}
fn co2_molar() -> f64 {
(atomic_mass(6) + 2.0 * atomic_mass(8)) * 1e-3
}
fn n2_molar() -> f64 {
2.0 * atomic_mass(7) * 1e-3
}
fn ar_molar() -> f64 {
atomic_mass(18) * 1e-3
}
fn o2_molar() -> f64 {
2.0 * atomic_mass(8) * 1e-3
}
fn co_molar() -> f64 {
(atomic_mass(6) + atomic_mass(8)) * 1e-3
}
fn h2o_molar() -> f64 {
(2.0 * atomic_mass(1) + atomic_mass(8)) * 1e-3
}
fn mean_molar_mass() -> f64 {
crate::CO2_FRACTION * co2_molar()
+ crate::N2_FRACTION * n2_molar()
+ crate::AR_FRACTION * ar_molar()
+ 0.001_3 * o2_molar()
+ 0.000_8 * co_molar()
+ 0.000_21 * h2o_molar()
}
fn sea_level_number_density() -> f64 {
let m_air = mean_molar_mass();
let rho = crate::SURFACE_PRESSURE_PA / (crate::R_SPECIFIC * crate::MEAN_TEMPERATURE_K);
N_A * rho / m_air
}
fn rayleigh_beta_rgb(n_density: f64) -> [f64; 3] {
let n_minus_1 = 4.49e-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 mars() -> Self {
let n_density = sea_level_number_density();
let beta = rayleigh_beta_rgb(n_density);
let species = vec![
MolecularSpecies {
name: "Carbon dioxide",
symbol: "CO2",
molar_mass_kg_mol: co2_molar(),
volume_fraction: crate::CO2_FRACTION,
refractive_index_stp: 1.000_449_0,
depolarization_factor: 0.075,
},
MolecularSpecies {
name: "Dinitrogen",
symbol: "N2",
molar_mass_kg_mol: n2_molar(),
volume_fraction: crate::N2_FRACTION,
refractive_index_stp: 1.000_298_0,
depolarization_factor: 0.030,
},
MolecularSpecies {
name: "Argon",
symbol: "Ar",
molar_mass_kg_mol: ar_molar(),
volume_fraction: crate::AR_FRACTION,
refractive_index_stp: 1.000_281_0,
depolarization_factor: 0.0,
},
MolecularSpecies {
name: "Dioxygen",
symbol: "O2",
molar_mass_kg_mol: o2_molar(),
volume_fraction: 0.001_3,
refractive_index_stp: 1.000_271_0,
depolarization_factor: 0.054,
},
MolecularSpecies {
name: "Carbon monoxide",
symbol: "CO",
molar_mass_kg_mol: co_molar(),
volume_fraction: 0.000_8,
refractive_index_stp: 1.000_338_0,
depolarization_factor: 0.030,
},
MolecularSpecies {
name: "Water vapor",
symbol: "H2O",
molar_mass_kg_mol: h2o_molar(),
volume_fraction: 0.000_21,
refractive_index_stp: 1.000_256_0,
depolarization_factor: 0.17,
},
];
Self {
planet_radius_m: crate::MARS_EQUATORIAL_RADIUS,
atmosphere_height_m: 100_000.0,
sea_level_pressure_pa: crate::SURFACE_PRESSURE_PA,
sea_level_temperature_k: crate::MEAN_TEMPERATURE_K,
sea_level_number_density_m3: n_density,
mean_molar_mass_kg_mol: mean_molar_mass(),
rayleigh_scale_height_m: crate::SCALE_HEIGHT_M,
mie_scale_height_m: 5_000.0,
mie_asymmetry_g: 0.63,
mie_coefficient: 5e-4,
dust_optical_depth: crate::MEAN_DUST_OPTICAL_DEPTH,
rayleigh_coefficients_rgb: beta,
species,
sun_irradiance_w_m2: 1_361.0 / (crate::SEMI_MAJOR_AXIS_AU * crate::SEMI_MAJOR_AXIS_AU),
}
}
pub fn rayleigh_density(&self, altitude_m: f64) -> f64 {
(-altitude_m / self.rayleigh_scale_height_m).exp()
}
pub fn dust_density(&self, altitude_m: f64) -> f64 {
self.dust_optical_depth * (-altitude_m / self.mie_scale_height_m).exp()
}
pub fn sky_color(&self, sun_elevation_deg: f64) -> [f64; 3] {
let el = sun_elevation_deg.clamp(-5.0, 90.0);
let zenith_cos = el.to_radians().sin().max(0.01);
let airmass = 1.0 / zenith_cos;
let dust_tau = self.dust_optical_depth * airmass;
let dust_extinction = (-dust_tau).exp();
let ray_r = (-self.rayleigh_coefficients_rgb[0] * airmass * 1e4).exp();
let ray_g = (-self.rayleigh_coefficients_rgb[1] * airmass * 1e4).exp();
let ray_b = (-self.rayleigh_coefficients_rgb[2] * airmass * 1e4).exp();
let ssa = 0.92;
let dust_scatter_r = ssa * 0.75 * (1.0 - dust_extinction);
let dust_scatter_g = ssa * 0.55 * (1.0 - dust_extinction);
let dust_scatter_b = ssa * 0.35 * (1.0 - dust_extinction);
let r = (dust_scatter_r + ray_r * dust_extinction * 0.01).clamp(0.0, 1.0);
let g = (dust_scatter_g + ray_g * dust_extinction * 0.01).clamp(0.0, 1.0);
let b = (dust_scatter_b + ray_b * dust_extinction * 0.01).clamp(0.0, 1.0);
if sun_elevation_deg < 0.0 {
let fade = ((sun_elevation_deg + 5.0) / 5.0).clamp(0.0, 1.0);
[r * fade, g * fade, b * fade]
} else {
[r, g, b]
}
}
pub fn direct_transmission(&self, sun_elevation_deg: f64) -> f64 {
let el = sun_elevation_deg.clamp(0.01, 90.0);
let airmass = 1.0 / el.to_radians().sin().max(0.01);
let total_tau = self.dust_optical_depth + 0.01;
(-total_tau * airmass).exp()
}
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 dust_abs = self.dust_optical_depth * (-observer_alt / self.mie_scale_height_m).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
+ dust_abs * 1e-4
+ 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))
}
}