use marss::{
CO2_FROST_POINT_K, GEOTHERMAL_GRADIENT, GLOBAL_STORM_OPTICAL_DEPTH, MARS_YEAR_S,
MEAN_DUST_OPTICAL_DEPTH, MEAN_TEMPERATURE_K, POLAR_MIN_TEMP_K, SUBSOLAR_MAX_TEMP_K,
};
use marss::{
atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
temporal, terrain,
};
const TARGET_FPS: f64 = 60.0;
const REAL_DT_S: f64 = 1.0 / TARGET_FPS;
const FRAME_COUNT: usize = 100_000;
fn biome_material(biome: terrain::texturing::MarsBiome) -> rendering::materials::PbrMaterial {
match biome {
terrain::texturing::MarsBiome::DarkDune => rendering::materials::dark_dune(),
terrain::texturing::MarsBiome::BrightDust => rendering::materials::bright_dust(),
terrain::texturing::MarsBiome::Regolith => rendering::materials::regolith(),
terrain::texturing::MarsBiome::Basalt => rendering::materials::basalt(),
terrain::texturing::MarsBiome::SulfateSalt => rendering::materials::sulfate_salt(),
terrain::texturing::MarsBiome::PolarIce => rendering::materials::water_ice(),
terrain::texturing::MarsBiome::Co2Frost => rendering::materials::co2_ice(),
terrain::texturing::MarsBiome::VolcanicSummit => rendering::materials::volcanic_summit(),
terrain::texturing::MarsBiome::ImpactEjecta => rendering::materials::basalt(),
}
}
struct MarsDiagnostics {
ls_deg: f64,
local_time_h: f64,
solar_irradiance_w_m2: f64,
surface_pressure_pa: f64,
cloudless_transmission: f64,
dust_optical_depth: f64,
ambient_light: f64,
daylight_state: lighting::day_night::DaylightState,
season: lighting::seasons::SeasonalState,
sky_zenith: [f64; 3],
sky_sunset: [f64; 3],
scatter_sky: [f64; 3],
current_biome: terrain::texturing::MarsBiome,
material_roughness: f32,
phobos_distance_m: f64,
deimos_distance_m: f64,
solar_tide_bulge_m: f64,
phobos_tide_bulge_m: f64,
quake_probability: f64,
marsquake_triggered: bool,
erosion_flux: f64,
thermal_stress_mpa: f64,
chemical_weathering: f64,
north_polar_gel_m: f64,
south_polar_gel_m: f64,
paleolake_volume_km3: f64,
subsurface_water_km3: f64,
best_habitability_score: f64,
best_habitability_zone: &'static str,
volcanic_coverage: f64,
volcanic_age_myr: f64,
rayleigh_number: f64,
heat_flow_w_m2: f64,
surface_temp_k: f64,
co2_frost_deposited: bool,
thermal_inertia: f64,
leaf_count: usize,
orbiters_mean_period_h: f64,
checksum: f64,
}
impl Default for MarsDiagnostics {
fn default() -> Self {
Self {
ls_deg: 0.0,
local_time_h: 12.0,
solar_irradiance_w_m2: 0.0,
surface_pressure_pa: 0.0,
cloudless_transmission: 0.0,
dust_optical_depth: MEAN_DUST_OPTICAL_DEPTH,
ambient_light: 0.0,
daylight_state: lighting::day_night::DaylightState::Night,
season: lighting::seasons::season_at(0.0),
sky_zenith: [0.0; 3],
sky_sunset: [0.0; 3],
scatter_sky: [0.0; 3],
current_biome: terrain::texturing::MarsBiome::Regolith,
material_roughness: rendering::materials::regolith().roughness,
phobos_distance_m: 0.0,
deimos_distance_m: 0.0,
solar_tide_bulge_m: 0.0,
phobos_tide_bulge_m: 0.0,
quake_probability: 0.0,
marsquake_triggered: false,
erosion_flux: 0.0,
thermal_stress_mpa: 0.0,
chemical_weathering: 0.0,
north_polar_gel_m: 0.0,
south_polar_gel_m: 0.0,
paleolake_volume_km3: 0.0,
subsurface_water_km3: 0.0,
best_habitability_score: 0.0,
best_habitability_zone: "",
volcanic_coverage: 0.0,
volcanic_age_myr: 0.0,
rayleigh_number: 0.0,
heat_flow_w_m2: 0.0,
surface_temp_k: MEAN_TEMPERATURE_K,
co2_frost_deposited: false,
thermal_inertia: atmosphere::heatwaves::global_mean_thermal_inertia(),
leaf_count: 0,
orbiters_mean_period_h: 0.0,
checksum: 0.0,
}
}
}
impl MarsDiagnostics {
fn state_checksum(&self) -> f64 {
let daylight = match self.daylight_state {
lighting::day_night::DaylightState::Day => 1.0,
lighting::day_night::DaylightState::CivilTwilight => 0.75,
lighting::day_night::DaylightState::NauticalTwilight => 0.5,
lighting::day_night::DaylightState::AstronomicalTwilight => 0.25,
lighting::day_night::DaylightState::Night => 0.0,
};
let north_season = match self.season.season_north {
lighting::seasons::Season::NorthernSpring => 0.1,
lighting::seasons::Season::NorthernSummer => 0.2,
lighting::seasons::Season::NorthernAutumn => 0.3,
lighting::seasons::Season::NorthernWinter => 0.4,
};
let south_season = match self.season.season_south {
lighting::seasons::Season::NorthernSpring => 0.1,
lighting::seasons::Season::NorthernSummer => 0.2,
lighting::seasons::Season::NorthernAutumn => 0.3,
lighting::seasons::Season::NorthernWinter => 0.4,
};
let biome = match self.current_biome {
terrain::texturing::MarsBiome::DarkDune => 1.0,
terrain::texturing::MarsBiome::BrightDust => 2.0,
terrain::texturing::MarsBiome::Regolith => 3.0,
terrain::texturing::MarsBiome::Basalt => 4.0,
terrain::texturing::MarsBiome::SulfateSalt => 5.0,
terrain::texturing::MarsBiome::PolarIce => 6.0,
terrain::texturing::MarsBiome::Co2Frost => 7.0,
terrain::texturing::MarsBiome::VolcanicSummit => 8.0,
terrain::texturing::MarsBiome::ImpactEjecta => 9.0,
};
self.ls_deg
+ self.local_time_h
+ self.solar_irradiance_w_m2 * 1e-3
+ self.surface_pressure_pa * 1e-3
+ self.cloudless_transmission
+ self.dust_optical_depth
+ self.ambient_light
+ daylight
+ north_season
+ south_season
+ self.season.ls_deg * 1e-2
+ self.season.solar_declination_deg.abs() * 1e-2
+ self.season.sub_solar_latitude_deg.abs() * 1e-2
+ self.sky_zenith.iter().sum::<f64>()
+ self.sky_sunset.iter().sum::<f64>()
+ self.scatter_sky.iter().sum::<f64>()
+ biome
+ self.material_roughness as f64
+ self.phobos_distance_m * 1e-6
+ self.deimos_distance_m * 1e-6
+ self.solar_tide_bulge_m * 1e6
+ self.phobos_tide_bulge_m * 1e6
+ self.quake_probability * 1e6
+ if self.marsquake_triggered { 1.0 } else { 0.0 }
+ self.erosion_flux
+ self.thermal_stress_mpa
+ self.chemical_weathering
+ self.north_polar_gel_m
+ self.south_polar_gel_m
+ self.paleolake_volume_km3 * 1e-3
+ self.subsurface_water_km3 * 1e-6
+ self.best_habitability_score
+ self.best_habitability_zone.len() as f64 * 1e-2
+ self.volcanic_coverage
+ self.volcanic_age_myr * 1e-2
+ self.rayleigh_number * 1e-24
+ self.heat_flow_w_m2
+ self.surface_temp_k * 1e-2
+ if self.co2_frost_deposited { 1.0 } else { 0.0 }
+ self.thermal_inertia * 1e-2
+ self.leaf_count as f64 * 1e-3
+ self.orbiters_mean_period_h * 1e-2
}
}
struct MarsSimulation {
epoch: temporal::epoch::MarsEpoch,
time_scale: temporal::time_scale::TimeScale,
orbit: physics::orbit::MarsOrbit,
rotation: physics::rotation::MarsRotation,
rotation_angle_rad: f64,
phobos: satellites::phobos::PhobosState,
deimos: satellites::deimos::DeimosState,
orbiters: Vec<satellites::artificial::MarsSatellite>,
active_lander: satellites::artificial::Lander,
heightmap: terrain::heightmap::Heightmap,
lod: terrain::lod::LodTerrain,
biome_classifier: terrain::texturing::BiomeClassifier,
climate: atmosphere::climate::MarsClimateState,
atmosphere_params: rendering::atmosphere_scattering::AtmosphereEndpoint,
dust_layer: rendering::dust_rendering::DustSystemEndpoint,
current_material: rendering::materials::PbrMaterial,
terrain_shader: rendering::shaders::ShaderEndpoint,
atmosphere_shader: rendering::shaders::ShaderEndpoint,
dust_shader: rendering::shaders::ShaderEndpoint,
observer: geodata::coordinates::LatLon,
observer_cartesian: geodata::coordinates::Cartesian,
diagnostics: MarsDiagnostics,
}
impl MarsSimulation {
fn new() -> Self {
let heightmap = terrain::heightmap::Heightmap::generate(360, 180);
let active_lander = satellites::artificial::perseverance();
let observer_alt =
heightmap.sample(active_lander.latitude_deg, active_lander.longitude_deg);
let observer = geodata::coordinates::LatLon::new(
active_lander.latitude_deg,
active_lander.longitude_deg,
observer_alt,
);
let observer_cartesian = observer.to_cartesian();
let orbiters = vec![
satellites::artificial::mro(),
satellites::artificial::maven(),
satellites::artificial::odyssey(),
satellites::artificial::tgo(),
satellites::artificial::mars_express(),
];
Self {
epoch: temporal::epoch::MarsEpoch::j2000(),
time_scale: temporal::time_scale::TimeScale::new(),
orbit: physics::orbit::MarsOrbit::new(),
rotation: physics::rotation::MarsRotation::new(),
rotation_angle_rad: 0.0,
phobos: satellites::phobos::PhobosState::new(),
deimos: satellites::deimos::DeimosState::new(),
orbiters,
active_lander,
heightmap,
lod: terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default()),
biome_classifier: terrain::texturing::BiomeClassifier::default(),
climate: atmosphere::climate::MarsClimateState::current(),
atmosphere_params: rendering::atmosphere_scattering::AtmosphereEndpoint::mars(),
dust_layer: rendering::dust_rendering::DustSystemEndpoint::mars_background(),
current_material: rendering::materials::regolith(),
terrain_shader: rendering::shaders::ShaderEndpoint::terrain(),
atmosphere_shader: rendering::shaders::ShaderEndpoint::atmosphere(),
dust_shader: rendering::shaders::ShaderEndpoint::dust_storm(),
observer,
observer_cartesian,
diagnostics: MarsDiagnostics::default(),
}
}
fn tick(&mut self, frame: u64, real_dt_s: f64) {
self.time_scale.step(real_dt_s);
let sim_dt = self.time_scale.simulation_dt(real_dt_s);
if sim_dt == 0.0 {
return;
}
self.epoch.advance_seconds(sim_dt);
self.orbit.step(sim_dt);
self.rotation_angle_rad = (self.rotation_angle_rad
+ self.rotation.angular_velocity_rad_s * sim_dt)
% (2.0 * std::f64::consts::PI);
self.phobos.step(sim_dt);
self.deimos.step(sim_dt);
let observer_alt = self.heightmap.sample(
self.active_lander.latitude_deg,
self.active_lander.longitude_deg,
);
self.observer = geodata::coordinates::LatLon::new(
self.active_lander.latitude_deg,
self.active_lander.longitude_deg,
observer_alt,
);
self.observer_cartesian = self.observer.to_cartesian();
let sol_date = temporal::calendar::MarsSolDate::from_julian_date(self.epoch.julian_date);
let ls_deg = sol_date.ls_deg();
let season = lighting::seasons::season_at(ls_deg);
let local_time_h = ((self.rotation_angle_rad / (2.0 * std::f64::consts::PI)) * 24.0
+ self.observer.lon_deg / 360.0 * 24.0
+ 12.0)
.rem_euclid(24.0);
let day_night = lighting::day_night::DayNightCycle::new(ls_deg);
let daylight_state = day_night.state_at(self.observer.lat_deg, local_time_h);
let ambient_light = day_night.ambient_light(self.observer.lat_deg, local_time_h);
let sun = lighting::solar_position::SolarPosition::compute(
ls_deg,
local_time_h,
self.observer.lat_deg,
self.observer.lon_deg,
);
let is_dust_season = lighting::seasons::is_dust_storm_season(ls_deg);
let storm_factor = if is_dust_season {
((ls_deg - 180.0) / 180.0 * std::f64::consts::PI)
.sin()
.max(0.0)
} else {
0.0
};
self.dust_layer = if storm_factor > 0.55 {
rendering::dust_rendering::DustSystemEndpoint::mars_global_storm()
} else {
rendering::dust_rendering::DustSystemEndpoint::mars_background()
};
self.dust_layer.optical_depth = MEAN_DUST_OPTICAL_DEPTH
+ storm_factor * (GLOBAL_STORM_OPTICAL_DEPTH - MEAN_DUST_OPTICAL_DEPTH);
self.atmosphere_params.dust_optical_depth = self.dust_layer.optical_depth;
let thermal_cycle = if self.observer.lat_deg.abs() >= 70.0 {
match season.season_north {
lighting::seasons::Season::NorthernSummer => {
atmosphere::heatwaves::ThermalCycle::polar_summer()
}
lighting::seasons::Season::NorthernWinter => {
atmosphere::heatwaves::ThermalCycle::polar_winter()
}
_ => atmosphere::heatwaves::ThermalCycle::polar_summer(),
}
} else {
match season.season_north {
lighting::seasons::Season::NorthernWinter => {
atmosphere::heatwaves::ThermalCycle::equatorial_winter()
}
_ => atmosphere::heatwaves::ThermalCycle::equatorial_summer(),
}
};
let solar_irradiance = self.orbit.solar_irradiance();
let direct_transmission = self
.atmosphere_params
.direct_transmission(sun.elevation_deg);
let surface_solar_fraction = self.dust_layer.surface_solar_fraction() * direct_transmission;
let solar_heating = if sun.is_above_horizon() {
solar_irradiance
* surface_solar_fraction
* sun.elevation_deg.to_radians().sin().max(0.0)
* (1.0 - self.climate.albedo)
} else {
0.0
};
let radiative_loss = 0.95
* sciforge::hub::domain::common::constants::SIGMA_SB
* self.diagnostics.surface_temp_k.powi(4);
let heat_capacity = thermal_cycle.thermal_inertia_units().powi(2) / 1000.0;
let diurnal_target = thermal_cycle.temperature_at_phase(local_time_h / 24.0);
let relaxation = (diurnal_target - self.diagnostics.surface_temp_k) * 0.02;
let dt_temp =
(solar_heating - radiative_loss) / heat_capacity.max(1.0) * sim_dt + relaxation;
let mut surface_temp_k = (self.diagnostics.surface_temp_k + dt_temp)
.clamp(POLAR_MIN_TEMP_K, SUBSOLAR_MAX_TEMP_K);
let co2_frost_deposited = surface_temp_k < CO2_FROST_POINT_K;
if co2_frost_deposited {
surface_temp_k = CO2_FROST_POINT_K;
}
let pressure_var = self.climate.seasonal_pressure_variation_pa();
let current_pressure = atmosphere::layers::barometric_pressure(observer_alt.max(0.0))
+ pressure_var * (ls_deg.to_radians().cos() * 0.5 + 0.5);
let mean_wind = atmosphere::winds::mean_surface_wind_speed();
let thermal_tide =
atmosphere::winds::thermal_tide_amplitude(thermal_cycle.temperature_swing_k());
let wind_speed = mean_wind + thermal_tide.min(25.0);
let aeolian = geology::erosion::AeolianErosion::typical();
let erosion_flux = aeolian.transport_rate(wind_speed);
let chemical_weathering = geology::erosion::chemical_weathering_rate(
surface_temp_k,
hydrology::oceans::RecurringSlopeLineae::count_observed() as f64 / 10_000.0,
);
let solar_tide = physics::tides::SolarTide::at_distance(self.orbit.current_radius());
let phobos_tide = physics::tides::PhobosTide::at_mean_distance();
let interior = geology::plate_tectonics::MarsInterior::new();
let volcanism = geology::volcanism::MarsVolcanism::new();
let quake_rate = geology::earthquakes::annual_frequency(3.0);
let quake_probability = quake_rate * sim_dt / MARS_YEAR_S;
let seed = frame
.wrapping_mul(6364136223846793005)
.wrapping_add(1442695040888963407);
let rng = (seed >> 33) as f64 / (1u64 << 31) as f64;
let marsquake_triggered = rng < quake_probability;
let north_polar = hydrology::glaciers::north_polar_cap();
let south_polar = hydrology::glaciers::south_polar_cap();
let paleolake = hydrology::lakes::jezero_crater_lake();
let hab_zones = biosphere::ecosystems::habitability_zones();
let best_hab = hab_zones
.iter()
.max_by(|a, b| {
a.habitability_score()
.partial_cmp(&b.habitability_score())
.unwrap()
})
.unwrap();
let thermal_inertia = thermal_cycle.thermal_inertia_units();
let biome =
self.biome_classifier
.classify(self.observer.lat_deg, observer_alt, thermal_inertia);
self.current_material = biome_material(biome);
let camera = (
self.observer_cartesian.x,
self.observer_cartesian.y,
self.observer_cartesian.z,
);
self.lod.update(camera);
self.terrain_shader.uniforms[0].value = rendering::shaders::UniformValue::Vec3([
sun.direction[0] as f32,
sun.direction[1] as f32,
sun.direction[2] as f32,
]);
self.terrain_shader.uniforms[1].value = rendering::shaders::UniformValue::Vec3([
camera.0 as f32,
camera.1 as f32,
camera.2 as f32,
]);
self.terrain_shader.uniforms[4].value =
rendering::shaders::UniformValue::Float(self.dust_layer.optical_depth as f32);
self.atmosphere_shader.uniforms[3].value =
rendering::shaders::UniformValue::Float(solar_irradiance as f32);
self.dust_shader.uniforms[0].value =
rendering::shaders::UniformValue::Float(self.dust_layer.optical_depth as f32);
self.dust_shader.uniforms[3].value =
rendering::shaders::UniformValue::Float(storm_factor as f32);
let orbiters_mean_period_h = self
.orbiters
.iter()
.map(|orbiter| orbiter.orbital_period_s())
.sum::<f64>()
/ self.orbiters.len() as f64
/ 3600.0;
let sky = rendering::sky::SkyState::new(sun.elevation_deg, self.dust_layer.optical_depth);
self.diagnostics = MarsDiagnostics {
ls_deg,
local_time_h,
solar_irradiance_w_m2: solar_irradiance,
surface_pressure_pa: current_pressure,
cloudless_transmission: direct_transmission,
dust_optical_depth: self.dust_layer.optical_depth,
ambient_light,
daylight_state,
season,
sky_zenith: sky.zenith_color(),
sky_sunset: sky.sunset_color(),
scatter_sky: self.atmosphere_params.sky_color(sun.elevation_deg),
current_biome: biome,
material_roughness: self.current_material.roughness,
phobos_distance_m: self.phobos.distance(),
deimos_distance_m: self.deimos.distance(),
solar_tide_bulge_m: solar_tide.tidal_bulge_height(),
phobos_tide_bulge_m: phobos_tide.tidal_bulge_height(),
quake_probability,
marsquake_triggered,
erosion_flux,
thermal_stress_mpa: thermal_cycle.thermal_stress_mpa(),
chemical_weathering,
north_polar_gel_m: north_polar.global_equivalent_layer_m(),
south_polar_gel_m: south_polar.global_equivalent_layer_m(),
paleolake_volume_km3: paleolake.volume_km3(),
subsurface_water_km3: hydrology::oceans::subsurface_water_estimate_km3(),
best_habitability_score: best_hab.habitability_score(),
best_habitability_zone: best_hab.name,
volcanic_coverage: volcanism.total_volcanic_coverage(),
volcanic_age_myr: volcanism.last_eruption_estimate_myr(),
rayleigh_number: interior.mantle_rayleigh_number(),
heat_flow_w_m2: geology::plate_tectonics::surface_heat_flow_w_m2(
2.0,
GEOTHERMAL_GRADIENT,
),
surface_temp_k,
co2_frost_deposited,
thermal_inertia,
leaf_count: self.lod.leaf_count(),
orbiters_mean_period_h,
checksum: 0.0,
};
self.diagnostics.checksum = self.diagnostics.state_checksum();
}
}
fn main() {
let mut simulation = MarsSimulation::new();
for frame in 0..FRAME_COUNT as u64 {
simulation.tick(frame, REAL_DT_S);
}
}