use earth::temporal::calendar::SECONDS_PER_DAY;
use earth::{
ALBEDO_EMA_ALPHA, CELSIUS_TO_KELVIN, CLOUD_EXTINCTION_COEFF, CLOUD_LWC_SCALE, CP_SEAWATER,
L_VAPORIZATION, MANNING_N_GRAVEL, MANNING_N_ROCK, MANNING_N_SAND, MANNING_N_VEGETATED,
OCEAN_SURFACE_AREA, PAR_FRACTION, PAR_W_TO_UMOL, QUARTZ_DENSITY, R_VAPOR, SECONDS_PER_YEAR,
SMITH_DRAG_COEFF, STONE_MERIDIONAL_D, SURFACE_GRAVITY, TROPOPAUSE_EQUATOR_M, TROPOPAUSE_POLE_M,
VAPOR_PRESSURE_0C,
};
use earth::{
atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
temporal, terrain,
};
fn biome_material_table() -> Vec<rendering::materials::PbrMaterial> {
vec![
rendering::materials::PbrMaterial::ocean(),
rendering::materials::PbrMaterial::desert(),
rendering::materials::PbrMaterial::desert(),
rendering::materials::PbrMaterial::grassland(),
rendering::materials::PbrMaterial::forest(),
rendering::materials::PbrMaterial::forest(),
rendering::materials::PbrMaterial::snow(),
rendering::materials::PbrMaterial::snow(),
rendering::materials::PbrMaterial::rock(),
rendering::materials::PbrMaterial::volcanic(),
rendering::materials::PbrMaterial::ice(),
]
}
fn biome_ordinal(b: terrain::texturing::Biome) -> usize {
match b {
terrain::texturing::Biome::Ocean => 0,
terrain::texturing::Biome::Beach => 1,
terrain::texturing::Biome::Desert => 2,
terrain::texturing::Biome::Grassland => 3,
terrain::texturing::Biome::Forest => 4,
terrain::texturing::Biome::Taiga => 5,
terrain::texturing::Biome::Tundra => 6,
terrain::texturing::Biome::Snow => 7,
terrain::texturing::Biome::Mountain => 8,
terrain::texturing::Biome::Volcanic => 9,
}
}
struct Earth {
epoch: temporal::epoch::Epoch,
time_scale: temporal::time_scale::TimeScale,
calendar: temporal::calendar::DateTime,
orbit: physics::orbit::EarthOrbit,
mean_anomaly_rad: f64,
rotation: physics::rotation::EarthRotation,
rotation_angle_rad: f64,
moon: satellites::moon::MoonState,
constellation: satellites::artificial::Constellation,
climate: atmosphere::climate::ClimateState,
sst_c: f64,
magma: geology::volcanism::Magma,
tectonic_plate: geology::plate_tectonics::TectonicPlate,
mountain: geology::mountains::Mountain,
cumulative_tectonic_strain: f64,
predator_prey: biosphere::fauna::PredatorPrey,
vegetation: biosphere::vegetation::Vegetation,
ecosystem: biosphere::ecosystems::Ecosystem,
reference_npp: f64,
base_carrying_cap: f64,
glacier: hydrology::glaciers::Glacier,
river: hydrology::rivers::River,
river_network: hydrology::rivers::RiverNetwork,
lake: hydrology::lakes::Lake,
impact_accumulator_yr: f64,
impact_prng_state: u64,
heightmap: terrain::heightmap::Heightmap,
lod_terrain: terrain::lod::LodTerrain,
erosion_rate_m_s: f64,
wind_erosion_thresh: f64,
cloud_system: rendering::clouds::CloudSystem,
atmo_render: rendering::atmosphere_scattering::AtmosphereParams,
ocean_render: rendering::ocean_rendering::OceanParams,
biome_classifier: terrain::texturing::BiomeClassifier,
materials: Vec<rendering::materials::PbrMaterial>,
terrain_shader: rendering::shaders::ShaderData,
atmo_shader: rendering::shaders::ShaderData,
ocean_shader: rendering::shaders::ShaderData,
observer: geodata::coordinates::LatLon,
observer_cam: [f64; 3],
troposphere: atmosphere::layers::AtmosphericLayer,
ocean_layer: hydrology::oceans::OceanLayer,
regions: geodata::regions::RegionDatabase,
elevation: geodata::elevation::ElevationProvider,
bathymetry: geodata::bathymetry::BathymetryData,
ocean_heat_cap: f64,
prev_cloud_density: f64,
prev_surface_albedo: f64,
prev_biome: terrain::texturing::Biome,
}
impl Earth {
fn new() -> Self {
let dt = temporal::calendar::DateTime::new(2026, 3, 31, 12, 0, 0.0);
let initial_jd = dt.to_julian_date();
let epoch = temporal::epoch::Epoch::from_jd(initial_jd);
let time_scale = temporal::time_scale::TimeScale::fast_forward(SECONDS_PER_DAY);
let mut orbit = physics::orbit::EarthOrbit::new();
let days_j2000 = epoch.days_since_j2000();
let mean_motion = 2.0 * std::f64::consts::PI / orbit.orbital_period_s();
let mean_anomaly_rad =
(days_j2000 * SECONDS_PER_DAY * mean_motion) % (2.0 * std::f64::consts::PI);
let e = orbit.eccentricity;
let mut ecc_anomaly = mean_anomaly_rad + e * mean_anomaly_rad.sin();
for _ in 0..15 {
let f = ecc_anomaly - e * ecc_anomaly.sin() - mean_anomaly_rad;
let fp = 1.0 - e * ecc_anomaly.cos();
ecc_anomaly -= f / fp;
}
orbit.true_anomaly_rad = 2.0
* f64::atan2(
(1.0 + e).sqrt() * (ecc_anomaly / 2.0).sin(),
(1.0 - e).sqrt() * (ecc_anomaly / 2.0).cos(),
);
let rotation = physics::rotation::EarthRotation::new();
let rotation_angle_rad = epoch.gmst_degrees().to_radians();
let moon = satellites::moon::MoonState::new();
let mut constellation = satellites::artificial::Constellation::new("Observers");
constellation.add(satellites::artificial::ArtificialSatellite::leo(
"ISS", 420_000.0, 408_000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::geo(
"GEO-1", 3000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::new(
"MEO-1",
1500.0,
20_200_000.0,
0.01,
0.96,
));
let climate = atmosphere::climate::ClimateState::current();
let ocean_layer = hydrology::oceans::surface_mixed_layer();
let sst_c = ocean_layer.mean_temperature_c;
let vegetation = biosphere::vegetation::tropical_broadleaf();
let init_photo = vegetation.photosynthesis_rate(climate.co2_ppm, 500.0, 25.0);
let reference_npp = vegetation.npp_kgc_m2_yr(vegetation.canopy_photosynthesis(init_photo));
let ecosystem = biosphere::ecosystems::tropical_rainforest();
let prey = biosphere::fauna::african_elephant();
let base_carrying_cap = prey.carrying_capacity;
let predator_prey = biosphere::fauna::PredatorPrey {
prey,
predator: biosphere::fauna::Population {
species_name: "Lion",
count: 20_000.0,
carrying_capacity: 50_000.0,
intrinsic_growth_rate: 0.1,
body_mass_kg: 190.0,
},
attack_rate: 1e-7,
conversion_efficiency: 0.1,
predator_death_rate: 0.15,
};
let glacier = hydrology::glaciers::antarctic_ice_sheet();
let river = hydrology::rivers::amazon();
let river_network = hydrology::rivers::amazon_basin();
let lake = hydrology::lakes::baikal();
let impact_accumulator_yr = 0.0;
let impact_seed: u64 = std::env::var("EARTH_IMPACT_SEED")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(42);
let impact_prng_state = impact_seed;
let tectonic_plate = geology::plate_tectonics::eurasian_plate();
let magma = geology::volcanism::Magma::basaltic();
let mountain = geology::mountains::everest();
let cumulative_tectonic_strain = 0.0;
let calendar = temporal::calendar::DateTime::from_julian_date(epoch.julian_date);
let mut heightmap = terrain::heightmap::Heightmap::new(256);
heightmap.generate_earth_topography();
heightmap.max_elevation_m = mountain.height();
let lod_terrain = terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default());
let temp_c = climate.global_mean_temp_k - CELSIUS_TO_KELVIN;
let erosion_mm_yr = geology::erosion::chemical_weathering_rate(temp_c, 1000.0)
+ geology::erosion::frost_weathering_rate(100.0, 0.05)
+ geology::erosion::glacial_erosion_rate(
glacier.deformation_velocity_m_yr(),
glacier.basal_shear_stress_pa(),
);
let erosion_rate_m_s = erosion_mm_yr / (1000.0 * SECONDS_PER_YEAR);
let wind_erosion_thresh =
geology::erosion::wind_erosion_threshold_velocity(0.001, *QUARTZ_DENSITY);
let cloud_system = rendering::clouds::CloudSystem::earth_default();
let atmo_render = rendering::atmosphere_scattering::AtmosphereParams::default();
let ocean_render = rendering::ocean_rendering::OceanParams::default();
let biome_classifier = terrain::texturing::BiomeClassifier::default();
let materials = biome_material_table();
let terrain_shader = rendering::shaders::ShaderData::terrain();
let atmo_shader = rendering::shaders::ShaderData::atmosphere();
let ocean_shader = rendering::shaders::ShaderData::ocean();
let obs_lat: f64 = std::env::var("EARTH_LAT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(48.8566);
let obs_lon: f64 = std::env::var("EARTH_LON")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(2.3522);
let obs_alt: f64 = std::env::var("EARTH_ALT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(35.0);
let observer = geodata::coordinates::LatLon::new(obs_lat, obs_lon, obs_alt);
let observer_cam = observer.to_cartesian_simple();
let troposphere = atmosphere::layers::troposphere();
let regions = geodata::regions::RegionDatabase::continents();
let elevation = geodata::elevation::ElevationProvider::global(30.0);
let bathymetry = geodata::bathymetry::BathymetryData::global(60.0);
let mixing_depth = ocean_layer.depth_range_m.1 - ocean_layer.depth_range_m.0;
let ocean_heat_cap = ocean_layer.density_kg_m3() * CP_SEAWATER * mixing_depth;
let init_tropopause = TROPOPAUSE_EQUATOR_M
- (TROPOPAUSE_EQUATOR_M - TROPOPAUSE_POLE_M)
* observer.lat_deg.to_radians().sin().abs();
let prev_cloud_density =
cloud_system.sample_density(init_tropopause * 0.25, observer.lat_deg, observer.lon_deg);
let prev_surface_albedo = climate.albedo;
let observer_elev = elevation.sample(observer.lat_deg, observer.lon_deg);
let prev_biome = biome_classifier.classify(observer_elev, observer.lat_deg, 0.5);
Earth {
epoch,
time_scale,
calendar,
orbit,
mean_anomaly_rad,
rotation,
rotation_angle_rad,
moon,
constellation,
climate,
sst_c,
magma,
tectonic_plate,
mountain,
cumulative_tectonic_strain,
predator_prey,
vegetation,
ecosystem,
reference_npp,
base_carrying_cap,
glacier,
river,
river_network,
lake,
impact_accumulator_yr,
impact_prng_state,
heightmap,
lod_terrain,
erosion_rate_m_s,
wind_erosion_thresh,
cloud_system,
atmo_render,
ocean_render,
biome_classifier,
materials,
terrain_shader,
atmo_shader,
ocean_shader,
observer,
observer_cam,
troposphere,
ocean_layer,
regions,
elevation,
bathymetry,
ocean_heat_cap,
prev_cloud_density,
prev_surface_albedo,
prev_biome,
}
}
fn tick(&mut self, real_dt_s: f64) {
let pi2 = 2.0 * std::f64::consts::PI;
let sim_dt = self.time_scale.simulation_dt(real_dt_s);
self.time_scale.step(real_dt_s);
self.epoch.advance_seconds(sim_dt);
let jd = self.epoch.julian_date;
let lat = self.observer.lat_deg;
let lon = self.observer.lon_deg;
let sim_years = sim_dt / SECONDS_PER_YEAR;
let mean_motion = pi2 / self.orbit.orbital_period_s();
self.mean_anomaly_rad = (self.mean_anomaly_rad + mean_motion * sim_dt) % pi2;
let m = self.mean_anomaly_rad;
let e = self.orbit.eccentricity;
let mut ecc_anomaly = m + e * m.sin();
for _ in 0..15 {
let f = ecc_anomaly - e * ecc_anomaly.sin() - m;
let fp = 1.0 - e * ecc_anomaly.cos();
ecc_anomaly -= f / fp;
}
self.orbit.true_anomaly_rad = 2.0
* f64::atan2(
(1.0 + e).sqrt() * (ecc_anomaly / 2.0).sin(),
(1.0 - e).sqrt() * (ecc_anomaly / 2.0).cos(),
);
let r_sun = self.orbit.current_radius();
self.rotation_angle_rad =
(self.rotation_angle_rad + self.rotation.angular_velocity_rad_s * sim_dt) % pi2;
let jd_tt = self.calendar.to_julian_date_tt();
let sun = lighting::solar_position::SolarPosition::compute(jd_tt, lat, lon);
let cos_zenith = sun.elevation_deg.to_radians().sin().max(0.0);
let season = lighting::seasons::season_at(jd, lat);
let dnc = lighting::day_night::DayNightCycle::new(jd);
let ambient = dnc.ambient_light(lat, lon);
let observer_alt = self.elevation.sample(lat, lon).max(0.0);
let pressure = self.troposphere.pressure_at(observer_alt);
let air_temp_k = self.troposphere.temperature_at(observer_alt);
let air_density = self.troposphere.density_at(observer_alt);
let lat_rad = lat.to_radians();
let delta_t = 40.0 * lat_rad.cos().abs();
let l_meridional = 1.0e7;
let grad_p = air_density * *SURFACE_GRAVITY * delta_t / (air_temp_k * l_meridional);
let wind_speed = {
let ws = atmosphere::winds::geostrophic_wind_speed(grad_p, lat, air_density);
if ws.is_finite() { ws.min(100.0) } else { 0.0 }
};
let ekman_depth = {
let ed = atmosphere::winds::ekman_spiral_depth(0.1, lat);
if ed.is_finite() { ed.min(500.0) } else { 200.0 }
};
self.ocean_heat_cap = self.ocean_layer.density_kg_m3() * CP_SEAWATER * ekman_depth;
let omega_earth = self.rotation.angular_velocity_rad_s;
let f_coriolis = 2.0 * omega_earth * lat_rad.sin();
let drag_coeff = SMITH_DRAG_COEFF;
let wind_stress = air_density * drag_coeff * wind_speed * wind_speed;
let surface_current = if f_coriolis.abs() > 1e-8 {
(wind_stress
/ (self.ocean_layer.density_kg_m3() * f_coriolis.abs() * ekman_depth.max(1.0)))
.min(2.0)
} else {
wind_speed * 0.03
};
self.moon.step(sim_dt);
let (mx, my, mz) = self.moon.position();
let moon_dist = (mx * mx + my * my + mz * mz).sqrt();
let tidal_height = physics::tides::TidalForce {
body_mass: satellites::moon::LUNAR_MASS,
body_distance: moon_dist,
}
.tidal_bulge_height()
+ physics::tides::TidalForce::from_sun().tidal_bulge_height();
self.constellation.step_all(sim_dt);
self.cloud_system.step(sim_dt);
let tropopause_h = TROPOPAUSE_EQUATOR_M
- (TROPOPAUSE_EQUATOR_M - TROPOPAUSE_POLE_M) * lat.to_radians().sin().abs();
let cloud_sample_alt = tropopause_h * 0.25;
let cloud_density = self.cloud_system.sample_density(cloud_sample_alt, lat, lon);
let solar_const = atmosphere::layers::mean_solar_irradiance();
let inv_sq = (physics::orbit::SEMI_MAJOR_AXIS / r_sun).powi(2);
let cloud_thickness = self
.cloud_system
.layers
.iter()
.map(|l| l.thickness_m * l.coverage)
.sum::<f64>();
let cloud_optical_depth =
CLOUD_EXTINCTION_COEFF * CLOUD_LWC_SCALE * cloud_thickness * self.prev_cloud_density;
let cloud_trans = (-cloud_optical_depth).exp();
let direct = solar_const * inv_sq * cos_zenith * cloud_trans;
let scattered = solar_const * 0.05 * ambient;
let irradiance = direct + scattered;
let effective_wind = wind_speed + surface_current;
self.ocean_render.wind_speed_ms = effective_wind.min(50.0);
let ocean_depth = if self.bathymetry.is_ocean(lat, lon) {
self.bathymetry.sample(lat, lon).abs().max(1.0)
} else {
self.ocean_render.depth_m
};
self.ocean_render.depth_m = ocean_depth;
let wave_amplitude = self.ocean_render.phillips_spectrum(0.1, 0.0).sqrt();
self.biome_classifier.sea_level = tidal_height + wave_amplitude * 0.01;
let sst_forcing = irradiance * (1.0 - self.prev_surface_albedo)
- hydrology::oceans::surface_longwave_radiation_w_m2(self.sst_c, 0.97);
self.sst_c += sst_forcing * sim_dt / (self.ocean_heat_cap);
let imbalance = self.climate.energy_imbalance();
let d_diff = STONE_MERIDIONAL_D;
let t_local = air_temp_k;
let t_global = self.climate.global_mean_temp_k;
let declination_factor = 1.0
+ 0.3
* (season.solar_declination_deg / lighting::solar_position::EARTH_AXIAL_TILT_DEG)
.abs();
let meridional_flux = -d_diff * declination_factor * (t_local - t_global);
let earth_surface_area =
4.0 * std::f64::consts::PI * (earth::geodata::coordinates::EARTH_SEMI_MAJOR_M).powi(2);
let ocean_fraction = OCEAN_SURFACE_AREA / earth_surface_area;
self.climate.global_mean_temp_k +=
(imbalance + meridional_flux * ocean_fraction) * sim_dt / self.ocean_heat_cap;
let surface_temp_c = air_temp_k - CELSIUS_TO_KELVIN;
let mass_balance = self.glacier.mass_balance_m_yr();
let sea_level_mm = self
.glacier
.sea_level_equivalent_mm(self.glacier.length_km * 1000.0);
self.biome_classifier.sea_level += -mass_balance * sim_years * 0.001;
let glacier_temp_c = surface_temp_c.min(-5.0);
let slope_rad = self.glacier.surface_slope_deg.to_radians();
let ice_velocity = hydrology::glaciers::shallow_ice_velocity(
self.glacier.thickness_m,
slope_rad,
glacier_temp_c,
);
let g_local = *crate::SURFACE_GRAVITY;
let overburden_pa = hydrology::glaciers::ice_density() * g_local * self.glacier.thickness_m;
let water_pressure_pa = 0.8 * overburden_pa;
let glacial_erosion = hydrology::glaciers::glacial_bed_erosion_rate(
ice_velocity,
self.glacier.thickness_m,
water_pressure_pa,
);
let network_discharge = self.river_network.route_discharge();
let outlet_discharge = self.river_network.outlet_discharge();
let network_sediment_factor = network_discharge
.iter()
.fold(0.0_f64, |acc, &d| acc + d * d)
/ (outlet_discharge * outlet_discharge).max(1.0);
let effective_depth =
(outlet_discharge / (self.river.drainage_area_km2 * 1e6)).powf(0.4) * 100.0;
let hydraulic_radius = effective_depth.max(1.0) * 0.7;
let manning_n = match self.prev_biome {
terrain::texturing::Biome::Desert | terrain::texturing::Biome::Beach => MANNING_N_SAND,
terrain::texturing::Biome::Mountain | terrain::texturing::Biome::Volcanic => {
MANNING_N_ROCK
}
terrain::texturing::Biome::Forest | terrain::texturing::Biome::Taiga => {
MANNING_N_VEGETATED
}
_ => MANNING_N_GRAVEL,
};
let river_velocity = self.river.manning_velocity(hydraulic_radius, manning_n);
let froude = self
.river
.froude_number(river_velocity, effective_depth.max(1.0));
let fluvial_erosion_factor = if froude > 0.5 { froude } else { 0.5 };
let stream_power = self.river.stream_power_w_per_m();
let lake_e_sat = VAPOR_PRESSURE_0C
* (*L_VAPORIZATION / *R_VAPOR * (1.0 / CELSIUS_TO_KELVIN - 1.0 / air_temp_k)).exp();
let lake_rh_pct = (lake_e_sat * 0.75 / lake_e_sat * 100.0).clamp(10.0, 100.0);
let lake_evap = self.lake.evaporation_rate_mm_day(
self.sst_c,
surface_temp_c,
wind_speed.max(1.0),
lake_rh_pct,
);
self.sst_c -= lake_evap * 0.001 * *L_VAPORIZATION / self.ocean_heat_cap * sim_dt;
let impact_mean_interval_yr = 1000.0;
self.impact_accumulator_yr += sim_years;
let impact_crater_erosion = if self.impact_accumulator_yr >= 1.0 {
let years_elapsed = self.impact_accumulator_yr;
self.impact_accumulator_yr = 0.0;
self.impact_prng_state ^= self.impact_prng_state << 13;
self.impact_prng_state ^= self.impact_prng_state >> 7;
self.impact_prng_state ^= self.impact_prng_state << 17;
let u = (self.impact_prng_state as f64) / (u64::MAX as f64);
let p_impact = 1.0 - (-years_elapsed / impact_mean_interval_yr).exp();
if u < p_impact {
let impactor = physics::collisions::tunguska_equivalent();
let crater_d = impactor.crater_diameter_m(*QUARTZ_DENSITY);
let ejecta = impactor.ejecta_volume_m3(*QUARTZ_DENSITY);
let fireball = impactor.fireball_radius_m();
let energy_mt = impactor.kinetic_energy_mt();
let earth_area = 4.0
* std::f64::consts::PI
* (earth::geodata::coordinates::EARTH_SEMI_MAJOR_M).powi(2);
(ejecta / earth_area)
* (1.0 + crater_d * 1e-6)
* (1.0 + fireball * 1e-9)
* (1.0 + energy_mt * 1e-12)
} else {
0.0
}
} else {
0.0
};
let par = irradiance * PAR_FRACTION * PAR_W_TO_UMOL;
let photo = self
.vegetation
.photosynthesis_rate(self.climate.co2_ppm, par, surface_temp_c);
let canopy = self.vegetation.canopy_photosynthesis(photo);
let transpiration = self.vegetation.transpiration_mm_day(
atmosphere::heatwaves::saturation_vapor_pressure_pa(surface_temp_c) / 1000.0,
0.01,
);
let npp = self.vegetation.npp_kgc_m2_yr(canopy);
let ecosystem_npp = self.ecosystem.total_npp_gc_yr();
let diversity = self.ecosystem.shannon_index();
let npp_ratio = (npp / self.reference_npp * diversity).clamp(0.01, 10.0);
self.predator_prey.prey.carrying_capacity = self.base_carrying_cap * npp_ratio;
self.predator_prey.step(sim_years);
self.calendar = temporal::calendar::DateTime::from_julian_date(self.epoch.julian_date);
let plate_velocity = self.tectonic_plate.velocity_at_point(lat, lon);
self.cumulative_tectonic_strain +=
plate_velocity * sim_dt / geology::plate_tectonics::LITHOSPHERE_THICKNESS;
let root_depth = self.mountain.root_depth_m();
let isostatic_elev = self.mountain.isostatic_elevation();
let strain_threshold = 1e-4;
let earthquake_energy = if self.cumulative_tectonic_strain > strain_threshold {
let moment = self.cumulative_tectonic_strain
* geology::plate_tectonics::ASTHENOSPHERE_VISCOSITY
* self.tectonic_plate.area_km2
* 1e6;
let quake = geology::earthquakes::Earthquake::from_moment(lat, lon, 15.0, moment);
let energy = quake.energy_joules();
let mw = quake.moment_magnitude();
let pga = quake.pga_at_distance(500.0);
self.cumulative_tectonic_strain = 0.0;
let aftershock_rate = geology::earthquakes::aftershock_rate(mw, 1.0);
energy * pga * 1e-10 * (1.0 + aftershock_rate * 0.01)
} else {
0.0
};
let geothermal_input = geology::plate_tectonics::surface_heat_flow(3.0, 0.03);
self.magma.temperature_c += (geothermal_input * 0.001 - 0.01) * sim_dt;
self.magma.temperature_c = self.magma.temperature_c.clamp(600.0, 1400.0);
let magma_viscosity = self.magma.viscosity_pa_s();
let mg_num = self.magma.mg_number();
let eruption_potential = self.magma.h2o_wt_percent / magma_viscosity.log10().max(1.0);
let vei = if eruption_potential > 0.1 {
geology::volcanism::volcanic_explosivity_index(eruption_potential * 0.01)
} else {
0
};
let stream_erosion = stream_power * 1e-12 * network_sediment_factor;
let effective_erosion = if wind_speed > self.wind_erosion_thresh {
self.erosion_rate_m_s * fluvial_erosion_factor * 2.0
+ earthquake_energy
+ glacial_erosion
+ stream_erosion
+ impact_crater_erosion
} else {
self.erosion_rate_m_s * fluvial_erosion_factor
+ earthquake_energy
+ glacial_erosion
+ stream_erosion
+ impact_crater_erosion
};
let u = ((lon + 180.0) / 360.0).clamp(0.0, 1.0);
let v = ((lat + 90.0) / 180.0).clamp(0.0, 1.0);
let hx = (u * (self.heightmap.resolution - 1) as f64) as usize;
let hy = (v * (self.heightmap.resolution - 1) as f64) as usize;
let hidx = hy * self.heightmap.resolution as usize + hx;
if let Some(cell) = self.heightmap.data.get_mut(hidx) {
*cell -= effective_erosion * sim_dt;
*cell += isostatic_elev * 1e-12 * sim_dt;
if vei > 0 {
*cell += vei as f64 * 0.01 * sim_years;
}
}
let min_elev = root_depth
* (geology::plate_tectonics::CRUST_DENSITY / geology::plate_tectonics::MANTLE_DENSITY
- 1.0);
let terrain_h = self.heightmap.sample(lat, lon).max(min_elev * 0.001);
let volcanic_fertility = mg_num * 0.01;
self.lod_terrain.update([
self.observer_cam[0],
self.observer_cam[1],
self.heightmap.radius_at(lat, lon),
]);
let is_continental = self
.regions
.point_in_region(lat, lon)
.map(|r| {
matches!(
r.region_type,
geodata::regions::RegionType::Continent
| geodata::regions::RegionType::Country
| geodata::regions::RegionType::Island
)
})
.unwrap_or(false);
let moisture_factor = if is_continental { 0.7 } else { 1.0 };
let e_sat = VAPOR_PRESSURE_0C
* (*L_VAPORIZATION / *R_VAPOR * (1.0 / CELSIUS_TO_KELVIN - 1.0 / air_temp_k)).exp();
let pressure_ratio = pressure / atmosphere::layers::SEA_LEVEL_PRESSURE;
let e_actual = e_sat * pressure_ratio * if is_continental { 0.5 } else { 0.85 };
let relative_humidity = (e_actual / e_sat).clamp(0.0, 1.0);
let moisture =
(relative_humidity * moisture_factor + transpiration * 0.005 + volcanic_fertility)
.clamp(0.0, 1.0);
let biome = self
.biome_classifier
.classify(terrain_h, lat, moisture.clamp(0.0, 1.0));
let mat = &self.materials[biome_ordinal(biome)];
let albedo =
mat.albedo[0] as f64 * 0.3 + mat.albedo[1] as f64 * 0.59 + mat.albedo[2] as f64 * 0.11;
self.climate.albedo =
self.climate.albedo * (1.0 - ALBEDO_EMA_ALPHA) + albedo * ALBEDO_EMA_ALPHA;
let splat = self
.biome_classifier
.splat(terrain_h, lat, moisture.clamp(0.0, 1.0));
let roughness = mat.roughness as f64 * splat.weights[0].1 as f64;
for uniform in &mut self.terrain_shader.uniforms {
if let rendering::shaders::UniformValue::Vec3(ref mut v) = uniform.value
&& uniform.name == "u_sun_direction"
{
*v = [
sun.direction[0] as f32,
sun.direction[1] as f32,
sun.direction[2] as f32,
];
}
}
for uniform in &mut self.atmo_shader.uniforms {
if let rendering::shaders::UniformValue::Float(ref mut f) = uniform.value
&& uniform.name == "u_sun_intensity"
{
*f = (irradiance / solar_const) as f32;
}
}
for uniform in &mut self.ocean_shader.uniforms {
if let rendering::shaders::UniformValue::Float(ref mut f) = uniform.value {
if uniform.name == "u_wind_speed" {
*f = wind_speed as f32;
}
if uniform.name == "u_wave_amplitude" {
*f = wave_amplitude as f32;
}
}
}
let sky = self
.atmo_render
.sky_color(sun.direction, [0.0, 1.0, 0.0], observer_alt);
let sky_lum = sky[0] * 0.3 + sky[1] * 0.59 + sky[2] * 0.11 + roughness * 0.01;
let mesh = terrain::mesh::TerrainMesh::from_region(
lat - 0.5,
lat + 0.5,
lon - 0.5,
lon + 0.5,
4,
&|la, lo| self.heightmap.sample(la, lo),
);
let render_budget = mesh.vertex_count() + mesh.triangle_count();
if let Some(vtx) = mesh.vertices.first() {
let surface_normal = vtx.normal;
let n_dot_l = surface_normal[0] * sun.direction[0]
+ surface_normal[1] * sun.direction[1]
+ surface_normal[2] * sun.direction[2];
let emissive_strength = mat.emissive[0] as f64 * n_dot_l.max(0.0) + sky_lum * 0.01;
self.prev_surface_albedo = albedo + emissive_strength * 0.001;
} else {
self.prev_surface_albedo = albedo;
}
self.prev_cloud_density = cloud_density;
self.prev_biome = biome;
self.biome_classifier.sea_level += sea_level_mm * 1e-6 * sim_years;
self.climate.co2_ppm -= ecosystem_npp * 1e-15 * sim_years;
if render_budget < 64 {
self.lod_terrain.update(self.observer_cam);
}
}
}
fn main() {
let mut earth = Earth::new();
let target_dt = std::time::Duration::from_secs_f64(1.0 / 60.0);
let mut last = std::time::Instant::now();
loop {
let now = std::time::Instant::now();
let real_dt = now.duration_since(last).as_secs_f64();
last = now;
earth.tick(real_dt);
let elapsed = now.elapsed();
if elapsed < target_dt {
std::thread::sleep(target_dt - elapsed);
}
}
}