use earths::temporal::calendar::SECONDSPERDAY;
use earths::{
ALBEDOEMAALPHA, CELSIUSTOKELVIN, CPSEAWATER, LVAPORIZATION, MANNINGNGRAVEL, MANNINGNROCK,
MANNINGNSAND, MANNINGNVEGETATED, OCEANSURFACEAREA, PARFRACTION, PARWTOUMOL, QUARTZDENSITY,
RVAPOR, SECONDSPERYEAR, SMITHDRAGCOEFF, STONEMERIDIONALD, SURFACEGRAVITY, TROPOPAUSEEQUATORM,
TROPOPAUSEPOLEM, VAPORPRESSURE0C,
};
use earths::{
atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
temporal, terrain,
};
fn biomematerialtable() -> Vec<rendering::materials::PbrMaterial> {
vec![
rendering::materials::PbrMaterial::deep_ocean(),
rendering::materials::PbrMaterial::desert_sand(),
rendering::materials::PbrMaterial::desert_sand(),
rendering::materials::PbrMaterial::grassland(),
rendering::materials::PbrMaterial::tropical_forest(),
rendering::materials::PbrMaterial::tropical_forest(),
rendering::materials::PbrMaterial::fresh_snow(),
rendering::materials::PbrMaterial::fresh_snow(),
rendering::materials::PbrMaterial::granite(),
rendering::materials::PbrMaterial::volcanic_lava(),
rendering::materials::PbrMaterial::glacier_ice(),
]
}
fn biomeordinal(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,
timescale: temporal::time_scale::TimeScale,
calendar: temporal::calendar::DateTime,
orbit: physics::orbit::EarthOrbit,
meananomalyrad: f64,
rotation: physics::rotation::EarthRotation,
rotationanglerad: f64,
moon: satellites::moon::MoonState,
constellation: satellites::artificial::Constellation,
climate: atmosphere::climate::ClimateState,
sstc: f64,
magma: geology::volcanism::Magma,
tectonicplate: geology::plate_tectonics::TectonicPlate,
mountain: geology::mountains::Mountain,
cumulativetectonicstrain: f64,
predatorprey: biosphere::fauna::PredatorPrey,
vegetation: biosphere::vegetation::Vegetation,
ecosystem: biosphere::ecosystems::Ecosystem,
referencenpp: f64,
basecarryingcap: f64,
glacier: hydrology::glaciers::Glacier,
river: hydrology::rivers::River,
rivernetwork: hydrology::rivers::RiverNetwork,
lake: hydrology::lakes::Lake,
impactaccumulatoryr: f64,
impactprngstate: u64,
heightmap: terrain::heightmap::Heightmap,
lodterrain: terrain::lod::LodTerrain,
erosionratems: f64,
winderosionthresh: f64,
cloudsystem: rendering::clouds::CloudSystemEndpoint,
atmorender: rendering::atmosphere_scattering::AtmosphereEndpoint,
oceanrender: rendering::ocean_rendering::OceanEndpoint,
biomeclassifier: terrain::texturing::BiomeClassifier,
materials: Vec<rendering::materials::PbrMaterial>,
terrainshader: rendering::shaders::ShaderEndpoint,
atmoshader: rendering::shaders::ShaderEndpoint,
oceanshader: rendering::shaders::ShaderEndpoint,
nightshader: rendering::shaders::ShaderEndpoint,
observer: geodata::coordinates::LatLon,
observercam: [f64; 3],
troposphere: atmosphere::layers::AtmosphericLayer,
oceanlayer: hydrology::oceans::OceanLayer,
regions: geodata::regions::RegionDatabase,
elevation: geodata::elevation::ElevationProvider,
bathymetry: geodata::bathymetry::BathymetryData,
oceanheatcap: f64,
prevclouddensity: f64,
prevsurfacealbedo: f64,
prevbiome: terrain::texturing::Biome,
}
impl Earth {
fn new() -> Self {
let dt = temporal::calendar::DateTime::new(2026, 3, 31, 12, 0, 0.0);
let initialjd = dt.tojuliandate();
let epoch = temporal::epoch::Epoch::fromjd(initialjd);
let timescale = temporal::time_scale::TimeScale::fastforward(SECONDSPERDAY);
let mut orbit = physics::orbit::EarthOrbit::new();
let daysj2000 = epoch.dayssincej2000();
let meanmotion = 2.0 * std::f64::consts::PI / orbit.orbitalperiods();
let meananomalyrad =
(daysj2000 * SECONDSPERDAY * meanmotion) % (2.0 * std::f64::consts::PI);
let e = orbit.eccentricity;
let mut eccanomaly = meananomalyrad + e * meananomalyrad.sin();
for iter in 0..15 {
let f = eccanomaly - e * eccanomaly.sin() - meananomalyrad;
let fp = 1.0 - e * eccanomaly.cos();
eccanomaly -= f / fp;
if f.abs() < 1e-12 || iter == 14 {
break;
}
}
orbit.trueanomalyrad = 2.0
* f64::atan2(
(1.0 + e).sqrt() * (eccanomaly / 2.0).sin(),
(1.0 - e).sqrt() * (eccanomaly / 2.0).cos(),
);
let rotation = physics::rotation::EarthRotation::new();
let rotationanglerad = epoch.gmstdegrees().to_radians();
let moon = satellites::moon::MoonState::new();
let mut constellation = satellites::artificial::Constellation::new("Observers");
constellation.add(satellites::artificial::ArtificialSatellite::leo(
"ISS", 420000.0, 408000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::geo(
"GEO-1", 3000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::new(
"MEO-1", 1500.0, 20200000.0, 0.01, 0.96,
));
let climate = atmosphere::climate::ClimateState::current();
let oceanlayer = hydrology::oceans::surfacemixedlayer();
let sstc = oceanlayer.meantemperaturec;
let vegetation = biosphere::vegetation::tropicalbroadleaf();
let initphoto = vegetation.photosynthesisrate(climate.co2ppm, 500.0, 25.0);
let referencenpp = vegetation.nppkgcm2yr(vegetation.canopyphotosynthesis(initphoto));
let ecosystem = biosphere::ecosystems::tropicalrainforest();
let prey = biosphere::fauna::africanelephant();
let basecarryingcap = prey.carryingcapacity;
let predatorprey = biosphere::fauna::PredatorPrey {
prey,
predator: biosphere::fauna::Population {
speciesname: "Lion",
count: 20000.0,
carryingcapacity: 50000.0,
intrinsicgrowthrate: 0.1,
bodymasskg: 190.0,
},
attackrate: 1e-7,
conversionefficiency: 0.1,
predatordeathrate: 0.15,
};
let glacier = hydrology::glaciers::antarcticicesheet();
let river = hydrology::rivers::amazon();
let rivernetwork = hydrology::rivers::amazonbasin();
let lake = hydrology::lakes::baikal();
let impactaccumulatoryr = 0.0;
let impactseed: u64 = std::env::var("EARTHIMPACTSEED")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(42);
let impactprngstate = impactseed;
let tectonicplate = geology::plate_tectonics::eurasianplate();
let magma = geology::volcanism::Magma::basaltic();
let mountain = geology::mountains::everest();
let cumulativetectonicstrain = 0.0;
let calendar = temporal::calendar::DateTime::fromjuliandate(epoch.juliandate);
let mut heightmap = terrain::heightmap::Heightmap::new(256);
heightmap.generateearthtopography();
heightmap.maxelevationm = mountain.height();
let lodterrain = terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default());
let tempc = climate.globalmeantempk - CELSIUSTOKELVIN;
let erosionmmyr = geology::erosion::chemicalweatheringrate(tempc, 1000.0)
+ geology::erosion::frostweatheringrate(100.0, 0.05)
+ geology::erosion::glacialerosionrate(
glacier.deformationvelocitymyr(),
glacier.basalshearstresspa(),
);
let erosionratems = erosionmmyr / (1000.0 * SECONDSPERYEAR);
let winderosionthresh =
geology::erosion::winderosionthresholdvelocity(0.001, *QUARTZDENSITY);
let cloudsystem = rendering::clouds::CloudSystemEndpoint::earth_default();
let atmorender = rendering::atmosphere_scattering::AtmosphereEndpoint::earth();
let oceanrender = rendering::ocean_rendering::OceanEndpoint::earth_atlantic();
let biomeclassifier = terrain::texturing::BiomeClassifier::default();
let materials = biomematerialtable();
let terrainshader = rendering::shaders::ShaderEndpoint::terrain();
let atmoshader = rendering::shaders::ShaderEndpoint::atmosphere();
let oceanshader = rendering::shaders::ShaderEndpoint::ocean();
let nightshader = rendering::shaders::ShaderEndpoint::night_lights();
let obslat: f64 = std::env::var("EARTHLAT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(48.8566);
let obslon: f64 = std::env::var("EARTHLON")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(2.3522);
let obsalt: f64 = std::env::var("EARTHALT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(35.0);
let observer = geodata::coordinates::LatLon::new(obslat, obslon, obsalt);
let observercam = observer.tocartesiansimple();
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 mixingdepth = oceanlayer.depthrangem.1 - oceanlayer.depthrangem.0;
let oceanheatcap = oceanlayer.densitykgm3() * CPSEAWATER * mixingdepth;
let inittropopause = TROPOPAUSEEQUATORM
- (TROPOPAUSEEQUATORM - TROPOPAUSEPOLEM) * observer.latdeg.to_radians().sin().abs();
let prevclouddensity = cloudsystem
.layers
.iter()
.map(|l| {
let alt_frac = (l.base_altitude_m / inittropopause).clamp(0.0, 1.0);
l.density * l.coverage * (-alt_frac).exp()
})
.sum::<f64>()
/ cloudsystem.layers.len().max(1) as f64;
let prevsurfacealbedo = climate.albedo;
let observerelev = elevation.sample(observer.latdeg, observer.londeg);
let prevbiome = biomeclassifier.classify(observerelev, observer.latdeg, 0.5);
Earth {
epoch,
timescale,
calendar,
orbit,
meananomalyrad,
rotation,
rotationanglerad,
moon,
constellation,
climate,
sstc,
magma,
tectonicplate,
mountain,
cumulativetectonicstrain,
predatorprey,
vegetation,
ecosystem,
referencenpp,
basecarryingcap,
glacier,
river,
rivernetwork,
lake,
impactaccumulatoryr,
impactprngstate,
heightmap,
lodterrain,
erosionratems,
winderosionthresh,
cloudsystem,
atmorender,
oceanrender,
biomeclassifier,
materials,
terrainshader,
atmoshader,
oceanshader,
nightshader,
observer,
observercam,
troposphere,
oceanlayer,
regions,
elevation,
bathymetry,
oceanheatcap,
prevclouddensity,
prevsurfacealbedo,
prevbiome,
}
}
fn tick(&mut self, realdts: f64) {
let pi2 = 2.0 * std::f64::consts::PI;
let simdt = self.timescale.simulationdt(realdts);
self.timescale.step(realdts);
self.epoch.advanceseconds(simdt);
let jd = self.epoch.juliandate;
let lat = self.observer.latdeg;
let lon = self.observer.londeg;
let simyears = simdt / SECONDSPERYEAR;
let meanmotion = pi2 / self.orbit.orbitalperiods();
self.meananomalyrad = (self.meananomalyrad + meanmotion * simdt) % pi2;
let m = self.meananomalyrad;
let e = self.orbit.eccentricity;
let mut eccanomaly = m + e * m.sin();
for iter in 0..15 {
let f = eccanomaly - e * eccanomaly.sin() - m;
let fp = 1.0 - e * eccanomaly.cos();
eccanomaly -= f / fp;
if f.abs() < 1e-12 || iter == 14 {
break;
}
}
self.orbit.trueanomalyrad = 2.0
* f64::atan2(
(1.0 + e).sqrt() * (eccanomaly / 2.0).sin(),
(1.0 - e).sqrt() * (eccanomaly / 2.0).cos(),
);
let rsun = self.orbit.current_radius();
self.rotationanglerad =
(self.rotationanglerad + self.rotation.angularvelocityrads * simdt) % pi2;
let jdtt = self.calendar.tojuliandatett();
let sun = lighting::solar_position::SolarPosition::compute(jdtt, lat, lon);
let coszenith = sun.elevationdeg.to_radians().sin().max(0.0);
let season = lighting::seasons::seasonat(jd, lat);
let dnc = lighting::day_night::DayNightCycle::new(jd);
let ambient = dnc.ambientlight(lat, lon);
let observeralt = self.elevation.sample(lat, lon).max(0.0);
let pressure = self.troposphere.pressureat(observeralt);
let airtempk = self.troposphere.temperatureat(observeralt);
let airdensity = self.troposphere.densityat(observeralt);
let latrad = lat.to_radians();
let deltat = 40.0 * latrad.cos().abs();
let lmeridional = 1.0e7;
let gradp = airdensity * *SURFACEGRAVITY * deltat / (airtempk * lmeridional);
let windspeed = {
let ws = atmosphere::winds::geostrophicwindspeed(gradp, lat, airdensity);
if ws.is_finite() { ws.min(100.0) } else { 0.0 }
};
let ekmandepth = {
let ed = atmosphere::winds::ekmanspiraldepth(0.1, lat);
if ed.is_finite() { ed.min(500.0) } else { 200.0 }
};
self.oceanheatcap = self.oceanlayer.densitykgm3() * CPSEAWATER * ekmandepth;
let omegaearth = self.rotation.angularvelocityrads;
let fcoriolis = 2.0 * omegaearth * latrad.sin();
let dragcoeff = SMITHDRAGCOEFF;
let windstress = airdensity * dragcoeff * windspeed * windspeed;
let surfacecurrent = if fcoriolis.abs() > 1e-8 {
(windstress / (self.oceanlayer.densitykgm3() * fcoriolis.abs() * ekmandepth.max(1.0)))
.min(2.0)
} else {
windspeed * 0.03
};
self.moon.step(simdt);
let (mx, my, mz) = self.moon.position();
let moondist = (mx * mx + my * my + mz * mz).sqrt();
let tidalheight = physics::tides::TidalForce {
bodymass: satellites::moon::LUNARMASS,
bodydistance: moondist,
}
.tidalbulgeheight()
+ physics::tides::TidalForce::fromsun().tidalbulgeheight();
self.constellation.stepall(simdt);
let tropopauseh = TROPOPAUSEEQUATORM
- (TROPOPAUSEEQUATORM - TROPOPAUSEPOLEM) * lat.to_radians().sin().abs();
let cloudsamplealt = tropopauseh * 0.25;
let clouddensity = self.cloudsystem.sample_density(cloudsamplealt);
let solarconst = self.atmorender.sun_irradiance_w_m2;
let invsq = (physics::orbit::SEMIMAJORAXIS / rsun).powi(2);
let cloud_wind_transport = self.cloudsystem.wind_transport();
let cloudthickness = self.cloudsystem.total_thickness();
let cloudopticaldepth = self.cloudsystem.optical_depth();
let cloudtrans = (-cloudopticaldepth).exp();
let cloud_advection = cloud_wind_transport * cloudthickness * 1e-10;
let direct = solarconst * invsq * coszenith * (cloudtrans - cloud_advection.min(0.01));
let scattered = solarconst * 0.05 * ambient;
let irradiance = direct + scattered;
let effectivewind = windspeed + surfacecurrent;
self.oceanrender.wind_speed_m_s = effectivewind.min(50.0);
self.oceanrender.surface_temperature_k = self.sstc + CELSIUSTOKELVIN;
let oceandepth = if self.bathymetry.isocean(lat, lon) {
self.bathymetry.sample(lat, lon).abs().max(1.0)
} else {
self.oceanrender.mean_depth_m
};
self.oceanrender.mean_depth_m = oceandepth;
let waveamplitude = self.oceanrender.wave_amplitude();
self.biomeclassifier.sealevel = tidalheight + waveamplitude * 0.01;
let sstforcing = irradiance * (1.0 - self.prevsurfacealbedo)
- hydrology::oceans::surfacelongwaveradiationwm2(self.sstc, 0.97);
let salinity_heat_factor = 1.0 - 0.0005 * (self.oceanrender.salinity_psu - 35.0);
let fresnel_r = self.oceanrender.fresnel_r0();
let barometric_scale = self.atmorender.barometric_density();
let atmo_shell_volume = self.atmorender.shell_volume();
let total_atmo_mass = barometric_scale * atmo_shell_volume * 1e-9;
let ocean_absorbed =
sstforcing * (1.0 - fresnel_r) * salinity_heat_factor * (1.0 + total_atmo_mass * 1e-25);
self.sstc += ocean_absorbed * simdt / self.oceanheatcap;
let imbalance = self.climate.energyimbalance();
let ddiff = STONEMERIDIONALD;
let tlocal = airtempk;
let tglobal = self.climate.globalmeantempk;
let declinationfactor = 1.0
+ 0.3
* (season.solardeclinationdeg / lighting::solar_position::EARTHAXIALTILTDEG).abs();
let meridionalflux = -ddiff * declinationfactor * (tlocal - tglobal);
let earthsurfacearea =
4.0 * std::f64::consts::PI * (earths::geodata::coordinates::EARTHSEMIMAJORM).powi(2);
let oceanfraction = OCEANSURFACEAREA / earthsurfacearea;
self.climate.globalmeantempk +=
(imbalance + meridionalflux * oceanfraction) * simdt / self.oceanheatcap;
let surfacetempc = airtempk - CELSIUSTOKELVIN;
let massbalance = self.glacier.massbalancemyr();
let sealevelmm = self
.glacier
.sealevelequivalentmm(self.glacier.lengthkm * 1000.0);
self.biomeclassifier.sealevel += -massbalance * simyears * 0.001;
let glaciertempc = surfacetempc.min(-5.0);
let sloperad = self.glacier.surfaceslopedeg.to_radians();
let icevelocity = hydrology::glaciers::shallowicevelocity(
self.glacier.thicknessm,
sloperad,
glaciertempc,
);
let glocal = *crate::SURFACEGRAVITY;
let overburdenpa = hydrology::glaciers::icedensity() * glocal * self.glacier.thicknessm;
let waterpressurepa = 0.8 * overburdenpa;
let glacialerosion = hydrology::glaciers::glacialbederosionrate(
icevelocity,
self.glacier.thicknessm,
waterpressurepa,
);
let networkdischarge = self.rivernetwork.routedischarge();
let outletdischarge = self.rivernetwork.outletdischarge();
let networksedimentfactor = networkdischarge.iter().fold(0.0f64, |acc, &d| acc + d * d)
/ (outletdischarge * outletdischarge).max(1.0);
let effectivedepth =
(outletdischarge / (self.river.drainageareakm2 * 1e6)).powf(0.4) * 100.0;
let hydraulicradius = effectivedepth.max(1.0) * 0.7;
let manningn = match self.prevbiome {
terrain::texturing::Biome::Desert | terrain::texturing::Biome::Beach => MANNINGNSAND,
terrain::texturing::Biome::Mountain | terrain::texturing::Biome::Volcanic => {
MANNINGNROCK
}
terrain::texturing::Biome::Forest | terrain::texturing::Biome::Taiga => {
MANNINGNVEGETATED
}
terrain::texturing::Biome::Grassland
| terrain::texturing::Biome::Tundra
| terrain::texturing::Biome::Snow
| terrain::texturing::Biome::Ocean => MANNINGNGRAVEL,
};
let rivervelocity = self.river.manningvelocity(hydraulicradius, manningn);
let froude = self
.river
.froudenumber(rivervelocity, effectivedepth.max(1.0));
let fluvialerosionfactor = if froude > 0.5 { froude } else { 0.5 };
let streampower = self.river.streampowerwperm();
let lakeesat = VAPORPRESSURE0C
* (*LVAPORIZATION / *RVAPOR * (1.0 / CELSIUSTOKELVIN - 1.0 / airtempk)).exp();
let lakerhpct = (lakeesat * 0.75 / lakeesat * 100.0).clamp(10.0, 100.0);
let lakeevap =
self.lake
.evaporationratemmday(self.sstc, surfacetempc, windspeed.max(1.0), lakerhpct);
self.sstc -= lakeevap * 0.001 * *LVAPORIZATION / self.oceanheatcap * simdt;
let impactmeanintervalyr = 1000.0;
self.impactaccumulatoryr += simyears;
let impactcratererosion = if self.impactaccumulatoryr >= 1.0 {
let yearselapsed = self.impactaccumulatoryr;
self.impactaccumulatoryr = 0.0;
self.impactprngstate ^= self.impactprngstate << 13;
self.impactprngstate ^= self.impactprngstate >> 7;
self.impactprngstate ^= self.impactprngstate << 17;
let u = (self.impactprngstate as f64) / (u64::MAX as f64);
let pimpact = 1.0 - (-yearselapsed / impactmeanintervalyr).exp();
if u < pimpact {
let impactor = physics::collisions::tunguskaequivalent();
let craterd = impactor.craterdiameterm(*QUARTZDENSITY);
let ejecta = impactor.ejectavolumem3(*QUARTZDENSITY);
let fireball = impactor.fireballradiusm();
let energymt = impactor.kineticenergymt();
let eartharea = 4.0
* std::f64::consts::PI
* (earths::geodata::coordinates::EARTHSEMIMAJORM).powi(2);
(ejecta / eartharea)
* (1.0 + craterd * 1e-6)
* (1.0 + fireball * 1e-9)
* (1.0 + energymt * 1e-12)
} else {
0.0
}
} else {
0.0
};
let par = irradiance * PARFRACTION * PARWTOUMOL;
let photo = self
.vegetation
.photosynthesisrate(self.climate.co2ppm, par, surfacetempc);
let canopy = self.vegetation.canopyphotosynthesis(photo);
let transpiration = self.vegetation.transpirationmmday(
atmosphere::heatwaves::saturationvaporpressurepa(surfacetempc) / 1000.0,
0.01,
);
let npp = self.vegetation.nppkgcm2yr(canopy);
let ecosystemnpp = self.ecosystem.totalnppgcyr();
let diversity = self.ecosystem.shannonindex();
let nppratio = (npp / self.referencenpp * diversity).clamp(0.01, 10.0);
self.predatorprey.prey.carryingcapacity = self.basecarryingcap * nppratio;
self.predatorprey.step(simyears);
self.calendar = temporal::calendar::DateTime::fromjuliandate(self.epoch.juliandate);
let platevelocity = self.tectonicplate.velocityatpoint(lat, lon);
self.cumulativetectonicstrain +=
platevelocity * simdt / geology::plate_tectonics::LITHOSPHERETHICKNESS;
let rootdepth = self.mountain.rootdepthm();
let isostaticelev = self.mountain.isostaticelevation();
let strainthreshold = 1e-4;
let earthquakeenergy = if self.cumulativetectonicstrain > strainthreshold {
let moment = self.cumulativetectonicstrain
* geology::plate_tectonics::ASTHENOSPHEREVISCOSITY
* self.tectonicplate.areakm2
* 1e6;
let quake = geology::earthquakes::Earthquake::frommoment(lat, lon, 15.0, moment);
let energy = quake.energyjoules();
let mw = quake.momentmagnitude();
let pga = quake.pgaatdistance(500.0);
self.cumulativetectonicstrain = 0.0;
let aftershockrate = geology::earthquakes::aftershockrate(mw, 1.0);
energy * pga * 1e-10 * (1.0 + aftershockrate * 0.01)
} else {
0.0
};
let geothermalinput = geology::plate_tectonics::surfaceheatflow(3.0, 0.03);
self.magma.temperaturec += (geothermalinput * 0.001 - 0.01) * simdt;
self.magma.temperaturec = self.magma.temperaturec.clamp(600.0, 1400.0);
let magmaviscosity = self.magma.viscositypas();
let mgnum = self.magma.mgnumber();
let eruptionpotential = self.magma.h2owtpercent / magmaviscosity.log10().max(1.0);
let vei = if eruptionpotential > 0.1 {
geology::volcanism::volcanicexplosivityindex(eruptionpotential * 0.01)
} else {
0
};
let streamerosion = streampower * 1e-12 * networksedimentfactor;
let effectiveerosion = if windspeed > self.winderosionthresh {
self.erosionratems * fluvialerosionfactor * 2.0
+ earthquakeenergy
+ glacialerosion
+ streamerosion
+ impactcratererosion
} else {
self.erosionratems * fluvialerosionfactor
+ earthquakeenergy
+ glacialerosion
+ streamerosion
+ impactcratererosion
};
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 -= effectiveerosion * simdt;
*cell += isostaticelev * 1e-12 * simdt;
if vei > 0 {
*cell += vei as f64 * 0.01 * simyears;
}
}
let minelev = rootdepth
* (geology::plate_tectonics::CRUSTDENSITY / geology::plate_tectonics::MANTLEDENSITY
- 1.0);
let terrainh = self.heightmap.sample(lat, lon).max(minelev * 0.001);
let volcanicfertility = mgnum * 0.01;
self.lodterrain.update([
self.observercam[0],
self.observercam[1],
self.heightmap.radiusat(lat, lon),
]);
let iscontinental = self
.regions
.pointinregion(lat, lon)
.map(|r| {
matches!(
r.regiontype,
geodata::regions::RegionType::Continent
| geodata::regions::RegionType::Country
| geodata::regions::RegionType::Island
)
})
.unwrap_or(false);
let moisturefactor = if iscontinental { 0.7 } else { 1.0 };
let esat = VAPORPRESSURE0C
* (*LVAPORIZATION / *RVAPOR * (1.0 / CELSIUSTOKELVIN - 1.0 / airtempk)).exp();
let pressureratio = pressure / atmosphere::layers::SEALEVELPRESSURE;
let eactual = esat * pressureratio * if iscontinental { 0.5 } else { 0.85 };
let relativehumidity = (eactual / esat).clamp(0.0, 1.0);
let moisture =
(relativehumidity * moisturefactor + transpiration * 0.005 + volcanicfertility)
.clamp(0.0, 1.0);
let biome = self
.biomeclassifier
.classify(terrainh, lat, moisture.clamp(0.0, 1.0));
let mat = &self.materials[biomeordinal(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 - ALBEDOEMAALPHA) + albedo * ALBEDOEMAALPHA;
let splat = self
.biomeclassifier
.splat(terrainh, lat, moisture.clamp(0.0, 1.0));
let roughness = mat.roughness as f64 * splat.weights[0].1 as f64;
for uniform in &mut self.terrainshader.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.atmoshader.uniforms {
if let rendering::shaders::UniformValue::Float(ref mut f) = uniform.value
&& uniform.name == "u_sun_intensity"
{
*f = (irradiance / solarconst) as f32;
}
}
for uniform in &mut self.oceanshader.uniforms {
if let rendering::shaders::UniformValue::Float(ref mut f) = uniform.value {
if uniform.name == "u_wind_speed" {
*f = windspeed as f32;
}
if uniform.name == "u_wave_amplitude" {
*f = waveamplitude as f32;
}
}
}
let is_night = coszenith <= 0.0;
let night_intensity = if is_night {
self.nightshader
.uniforms
.iter()
.filter_map(|u| {
if let rendering::shaders::UniformValue::Float(f) = u.value {
Some(f as f64)
} else {
None
}
})
.sum::<f64>()
} else {
0.0
};
let night_color = if is_night {
self.nightshader
.uniforms
.iter()
.find_map(|u| {
if let rendering::shaders::UniformValue::Vec3(v) = u.value {
Some(v)
} else {
None
}
})
.unwrap_or([0.0; 3])
} else {
[0.0; 3]
};
let coszen = coszenith.max(0.01);
let ocean_depth_color = self.oceanrender.depth_color(oceandepth);
let foam_active =
if self.oceanrender.wind_speed_m_s > self.oceanrender.foam_threshold_wind_m_s {
1.0
} else {
0.0
};
let skylum = self.atmorender.sky_luminance(observeralt, coszen)
+ roughness * 0.01
+ ocean_depth_color[1] * 0.001
+ foam_active * 0.001;
let mesh = terrain::mesh::TerrainMesh::fromregion(
lat - 0.5,
lat + 0.5,
lon - 0.5,
lon + 0.5,
4,
&|la, lo| self.heightmap.sample(la, lo),
);
let renderbudget = mesh.vertexcount()
+ mesh.trianglecount()
+ self.terrainshader.name.len()
+ self.atmoshader.name.len()
+ self.oceanshader.name.len()
+ self.nightshader.name.len();
if let Some(vtx) = mesh.vertices.first() {
let surfacenormal = vtx.normal;
let ndotl = surfacenormal[0] * sun.direction[0]
+ surfacenormal[1] * sun.direction[1]
+ surfacenormal[2] * sun.direction[2];
let fresnel_surface = mat.fresnel_r0();
let sss_contrib = mat.subsurface as f64 * (1.0 - ndotl.max(0.0).powi(2)) * 0.1;
let normal_detail = mat.normal_strength as f64;
let emissivestrength = mat.emissive[0] as f64 * ndotl.max(0.0)
+ skylum * 0.01
+ sss_contrib
+ fresnel_surface * 0.001
+ normal_detail * 0.0001
+ night_intensity * night_color[0] as f64 * 0.001;
self.prevsurfacealbedo = albedo + emissivestrength * 0.001;
} else {
self.prevsurfacealbedo = albedo;
}
self.prevclouddensity = clouddensity;
self.prevbiome = biome;
self.biomeclassifier.sealevel += sealevelmm * 1e-6 * simyears;
self.climate.co2ppm -= ecosystemnpp * 1e-15 * simyears;
if renderbudget < 64 {
self.lodterrain.update(self.observercam);
}
}
}
fn main() {
let mut earth = Earth::new();
let targetdt = 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 realdt = now.duration_since(last).as_secs_f64();
last = now;
earth.tick(realdt);
let elapsed = now.elapsed();
if elapsed < targetdt {
std::thread::sleep(targetdt - elapsed);
}
}
}