use jupiters::temporal::calendar::SECONDSPERDAY;
use jupiters::{
ALBEDOEMAALPHA, INTERNALHEATINGWM2, MAXWINDSMS, METHANEABSORPTION, OMEGAJUPITER, ONEBARTEMPK,
SECONDSPERYEAR, STONYASTEROIDDENSITY, TROPOPAUSEM,
};
use jupiters::{
atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
temporal, terrain,
};
fn cloudbiomematerialtable() -> Vec<rendering::materials::PbrMaterial> {
vec![
rendering::materials::PbrMaterial::ammonia_cloud(),
rendering::materials::PbrMaterial::nh4sh_cloud(),
rendering::materials::PbrMaterial::deep_atmosphere(),
rendering::materials::PbrMaterial::haze_layer(),
rendering::materials::PbrMaterial::core_rock(),
rendering::materials::PbrMaterial::metallic_hydrogen(),
rendering::materials::PbrMaterial::storm_region(),
]
}
fn cloudbiomeordinal(b: terrain::texturing::CloudBiome) -> usize {
match b {
terrain::texturing::CloudBiome::AmmoniaIceCloud => 0,
terrain::texturing::CloudBiome::NH4SHCloud => 1,
terrain::texturing::CloudBiome::DeepCloud => 2,
terrain::texturing::CloudBiome::UpperHaze => 3,
terrain::texturing::CloudBiome::ClearZone => 4,
terrain::texturing::CloudBiome::PolarVortex => 5,
terrain::texturing::CloudBiome::StormRegion => 6,
terrain::texturing::CloudBiome::WaterIceCloud => 2,
}
}
struct Jupiter {
epoch: temporal::epoch::Epoch,
timescale: temporal::time_scale::TimeScale,
calendar: temporal::calendar::DateTime,
orbit: physics::orbit::JupiterOrbit,
meananomalyrad: f64,
rotation: physics::rotation::JupiterRotation,
rotationanglerad: f64,
iostate: satellites::io::IoState,
constellation: satellites::artificial::Constellation,
climate: atmosphere::climate::ClimateState,
cloudtoptempk: f64,
cryomagma: geology::volcanism::CryoMagma,
convectioncell: geology::plate_tectonics::InteriorConvectionCell,
interiorfeature: geology::mountains::InteriorFeature,
cumulativeicestrain: f64,
cloudinteraction: biosphere::fauna::CloudLayerInteraction,
aerochem: biosphere::aerochemistry::AerialChemotroph,
ecosystem: biosphere::ecosystems::Ecosystem,
referencenpp: f64,
basecarryingcap: f64,
icelayer: hydrology::glaciers::IceLayer,
atmosphericflow: hydrology::rivers::AtmosphericFlow,
flownetwork: hydrology::rivers::FlowNetwork,
cryolake: hydrology::lakes::CryoLake,
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::InteriorOceanEndpoint,
biomeclassifier: terrain::texturing::CloudBiomeClassifier,
materials: Vec<rendering::materials::PbrMaterial>,
terrainshader: rendering::shaders::ShaderEndpoint,
atmoshader: rendering::shaders::ShaderEndpoint,
oceanshader: rendering::shaders::ShaderEndpoint,
bandshader: rendering::shaders::ShaderEndpoint,
observer: geodata::coordinates::LatLon,
observercam: [f64; 3],
troposphere: atmosphere::layers::AtmosphericLayer,
oceanlayer: hydrology::oceans::InteriorOceanLayer,
regions: geodata::regions::RegionDatabase,
elevation: geodata::elevation::ElevationProvider,
interiorlayers: geodata::bathymetry::InteriorLayerData,
interiorheatcap: f64,
prevclouddensity: f64,
prevsurfacealbedo: f64,
prevbiome: terrain::texturing::CloudBiome,
}
impl Jupiter {
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::JupiterOrbit::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::JupiterRotation::new();
let rotationanglerad = epoch.jupitergmstdegrees().to_radians();
let iostate = satellites::io::IoState::new();
let mut constellation = satellites::artificial::Constellation::new("JupiterOrbiters");
constellation.add(satellites::artificial::ArtificialSatellite::loworbit(
"JunoSuccessor",
800.0,
300000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::polarorbit(
"PolarMapper",
500.0,
500000.0,
));
constellation.add(satellites::artificial::ArtificialSatellite::new(
"DeepOrbit",
2000.0,
20000000.0,
0.05,
0.5,
));
let climate = atmosphere::climate::ClimateState::current();
let oceanlayer = hydrology::oceans::molecularhydrogenfluid();
let cloudtoptempk = ONEBARTEMPK;
let aerochem = biosphere::aerochemistry::hypotheticalchemotroph();
let initchemo = aerochem.chemosynthesisrate(climate.nh3fraction * 1e6, 0.0, ONEBARTEMPK);
let referencenpp = aerochem
.nppkgcm2yr(aerochem.columnproduction(initchemo))
.max(1e-30);
let ecosystem = biosphere::ecosystems::hypotheticalatmosphericbiome();
let producer = biosphere::fauna::hypotheticalh2floater();
let basecarryingcap = producer.carryingcapacity;
let cloudinteraction = biosphere::fauna::CloudLayerInteraction {
producers: producer,
consumers: biosphere::fauna::Aeroplankton {
speciesname: "Hypothetical Consumer",
count: 0.0,
carryingcapacity: 0.0,
intrinsicgrowthrate: 0.0,
effectivemasskg: 1.0,
},
interactionrate: 0.0,
conversionefficiency: 0.0,
consumerlossrate: 0.0,
};
let icelayer = hydrology::glaciers::metallichydrogenmantle();
let atmosphericflow = hydrology::rivers::equatorialjet();
let polar = hydrology::rivers::polarjet();
let flownetwork = hydrology::rivers::FlowNetwork {
flows: vec![hydrology::rivers::equatorialjet(), polar],
adjacency: vec![(0, 1)],
};
let cryolake = hydrology::lakes::hypotheticalmetallichydrogenocean();
let impactseed: u64 = std::env::var("JUPITERIMPACTSEED")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(42);
let convectioncell = geology::plate_tectonics::equatorialcell();
let cryomagma = geology::volcanism::CryoMagma::waterammonia();
let interiorfeature = geology::mountains::metallichydrogenmantle();
let calendar = temporal::calendar::DateTime::fromjuliandate(epoch.juliandate);
let mut heightmap = terrain::heightmap::Heightmap::new(256);
heightmap.generatejupitercloudtopography();
heightmap.maxelevationm = interiorfeature.height() + 5000.0;
let lodterrain = terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default());
let erosionparams = geology::erosion::ErosionParams {
windspeedms: MAXWINDSMS,
atmosphericdensity: *jupiters::ONEBARDENSITY,
particlediameterm: 0.0001,
particledensity: *jupiters::WATERICEDENSITY,
};
let erosionmmyr = erosionparams.winderosionratemmyr()
+ geology::erosion::chemicalweatheringrate(ONEBARTEMPK, *jupiters::ONEBARDENSITY)
+ geology::erosion::thermalstresserosion(
jupiters::CLOUDTOPTEMPMAXK - jupiters::CLOUDTOPTEMPMINK,
1e-5,
);
let erosionratems = erosionmmyr / (1000.0 * SECONDSPERYEAR);
let winderosionthresh =
geology::erosion::winderosionthresholdvelocity(0.0001, *jupiters::WATERICEDENSITY);
let cloudsystem = rendering::clouds::CloudSystemEndpoint::jupiter_default();
let atmorender = rendering::atmosphere_scattering::AtmosphereEndpoint::jupiter();
let oceanrender = rendering::ocean_rendering::InteriorOceanEndpoint::metallic_hydrogen();
let biomeclassifier = terrain::texturing::CloudBiomeClassifier::default();
let materials = cloudbiomematerialtable();
let terrainshader = rendering::shaders::ShaderEndpoint::terrain();
let atmoshader = rendering::shaders::ShaderEndpoint::atmosphere();
let oceanshader = rendering::shaders::ShaderEndpoint::interior_ocean();
let bandshader = rendering::shaders::ShaderEndpoint::band_dynamics();
let obslat: f64 = std::env::var("JUPITERLAT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(0.0);
let obslon: f64 = std::env::var("JUPITERLON")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(0.0);
let obsalt: f64 = std::env::var("JUPITERALT")
.ok()
.and_then(|s| s.parse().ok())
.unwrap_or(50000.0);
let observer = geodata::coordinates::LatLon::new(obslat, obslon, obsalt);
let observercam = observer.tocartesiansimple();
let troposphere = atmosphere::layers::troposphere();
let regions = geodata::regions::RegionDatabase::atmosphericbands();
let elevation = geodata::elevation::ElevationProvider::global(600.0);
let interiorlayers = geodata::bathymetry::InteriorLayerData::global(600.0);
let layerdepth = oceanlayer.thicknesskm * 1000.0;
let interiorheatcap = oceanlayer.densitykgm3 * 14300.0 * layerdepth;
let prevclouddensity = cloudsystem
.layers
.iter()
.map(|l| l.density * l.coverage)
.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, MAXWINDSMS * 0.5);
Jupiter {
epoch,
timescale,
calendar,
orbit,
meananomalyrad,
rotation,
rotationanglerad,
iostate,
constellation,
climate,
cloudtoptempk,
cryomagma,
convectioncell,
interiorfeature,
cumulativeicestrain: 0.0,
cloudinteraction,
aerochem,
ecosystem,
referencenpp,
basecarryingcap,
icelayer,
atmosphericflow,
flownetwork,
cryolake,
impactaccumulatoryr: 0.0,
impactprngstate: impactseed,
heightmap,
lodterrain,
erosionratems,
winderosionthresh,
cloudsystem,
atmorender,
oceanrender,
biomeclassifier,
materials,
terrainshader,
atmoshader,
oceanshader,
bandshader,
observer,
observercam,
troposphere,
oceanlayer,
regions,
elevation,
interiorlayers,
interiorheatcap,
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 jdcal = self.calendar.tojuliandate();
let sun = lighting::solar_position::SolarPosition::compute(jdcal, lat, lon);
let coszenith = sun.elevationdeg.to_radians().sin().max(0.0);
let season = lighting::seasons::Seasons::new(jd);
let dnc = lighting::day_night::DayNightCycle::new(jd);
let ambient = dnc.ambientlight(lat, lon);
let observeralt = self.elevation.sample(lat, lon);
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 = 10.0 * latrad.cos().abs();
let lmeridional = 5.0e7;
let gradp = airdensity * *jupiters::SURFACEGRAVITY * deltat / (airtempk * lmeridional);
let windspeed = {
let ws = atmosphere::winds::geostrophicwindspeed(gradp, lat, airdensity);
if ws.is_finite() {
ws.min(MAXWINDSMS)
} else {
0.0
}
};
let zonalwind = atmosphere::winds::zonalwindprofile(lat);
let effectivewind = (windspeed + zonalwind.abs()).min(MAXWINDSMS);
let ekmandepth = {
let ed = atmosphere::winds::ekmanspiraldepth(1.0, lat);
if ed.is_finite() {
ed.min(10000.0)
} else {
5000.0
}
};
let layerdepth = self.oceanlayer.thicknesskm * 1000.0;
self.interiorheatcap = self.oceanlayer.densitykgm3 * 14300.0 * layerdepth;
let fcoriolis = 2.0 * OMEGAJUPITER * latrad.sin();
let windstress = airdensity * 0.002 * effectivewind * effectivewind;
let interiorcurrent = if fcoriolis.abs() > 1e-8 {
(windstress / (self.oceanlayer.densitykgm3 * fcoriolis.abs() * ekmandepth.max(1.0)))
.min(5.0)
} else {
effectivewind * 0.01
};
self.iostate.step(simdt);
let (ix, iy, iz) = self.iostate.position();
let iodist = (ix * ix + iy * iy + iz * iz).sqrt();
let tidalheight = physics::tides::TidalForce {
bodymass: satellites::io::IOMASS,
bodydistance: iodist,
}
.tidalbulgeheight()
+ physics::tides::TidalForce::fromsun().tidalbulgeheight();
self.constellation.stepall(simdt);
let cloudsamplealt = TROPOPAUSEM * 0.5;
let cloud_optical_depth_total = self.cloudsystem.optical_depth();
let clouddensity = self.cloudsystem.sample_density(cloudsamplealt);
let cloud_wind_transport = self.cloudsystem.wind_transport();
let cloud_advection = cloud_wind_transport * simdt * 1e-6;
let barometric_density = self.atmorender.barometric_density();
let atmo_shell_volume = self.atmorender.shell_volume();
let solarconst = self.atmorender.sun_irradiance_w_m2;
let invsq = (physics::orbit::SEMIMAJORAXIS / rsun).powi(2);
let cloudtrans = (-cloud_optical_depth_total).exp();
let direct = solarconst * invsq * coszenith * cloudtrans;
let scattered = solarconst * 0.05 * ambient;
let irradiance = direct + scattered;
let ammoniaabsorbed = irradiance * METHANEABSORPTION;
let effectiveirradiance = irradiance - ammoniaabsorbed;
self.oceanrender.wind_speed_m_s = effectivewind.min(MAXWINDSMS);
let interiordepth = self.interiorlayers.sample(lat, lon).abs().max(1000.0);
let r0_ocean = self.oceanrender.fresnel_r0();
let waveamplitude = self.oceanrender.wave_amplitude();
let depth_color_rgb = self.oceanrender.depth_color(interiordepth);
let depth_color = (depth_color_rgb[0] + depth_color_rgb[1] + depth_color_rgb[2]) / 3.0;
let dx_ocean = self.oceanrender.patch_size_m / self.oceanrender.grid_size as f64;
let pressure_opacity = 1.0 - (-self.oceanrender.pressure_gpa * 0.01).exp();
let ocean_thermal_mass = self.oceanrender.density_kg_m3 * 14_300.0 * interiordepth;
self.biomeclassifier.onebarlevel =
tidalheight + waveamplitude * 0.001 + interiordepth * 1e-9;
let internalforcing = INTERNALHEATINGWM2;
let solarforcing = effectiveirradiance * (1.0 - self.prevsurfacealbedo);
let olr = atmosphere::layers::olr();
self.cloudtoptempk += (solarforcing + internalforcing - olr) * simdt / self.interiorheatcap;
self.cloudtoptempk = self.cloudtoptempk.clamp(110.0, 200.0);
let imbalance = self.climate.energyimbalance();
let subsol = season.subsolarlatitudedeg();
let declinationfactor =
1.0 + 0.3 * (subsol / lighting::solar_position::JUPITERAXIALTILTDEG).abs();
let tlocal = airtempk;
let tglobal = self.climate.globalmeantempk;
let meridionalflux = -1e6 * declinationfactor * (tlocal - tglobal);
let jupitersurfacearea =
4.0 * std::f64::consts::PI * geodata::coordinates::JUPITERSEMIMAJORM.powi(2);
self.climate.globalmeantempk +=
(imbalance + meridionalflux / jupitersurfacearea) * simdt / self.interiorheatcap;
let icebalance = self.icelayer.massbalancemyr();
let iceshellvel = self.icelayer.deformationvelocitymyr();
let icebasal = self.icelayer.basalshearstresspa();
let thermalstress = geology::erosion::thermalstresserosion(
jupiters::CLOUDTOPTEMPMAXK - jupiters::CLOUDTOPTEMPMINK,
1e-5,
);
let networkdischarge = self.flownetwork.totaldischargem3s();
let flowlength = self.flownetwork.totallengthkm();
let networksedimentfactor =
networkdischarge * networkdischarge / (flowlength * flowlength + 1.0);
let effectivedepth = (networkdischarge / self.atmosphericflow.widthm).powf(0.4) * 1000.0;
let hydraulicradius = effectivedepth.max(1.0) * 0.7;
let flowvelocity = self.atmosphericflow.manningvelocityms();
let froude = self.atmosphericflow.froudenumber();
let fluvialerosionfactor = if froude > 0.5 { froude } else { 0.5 };
let streampower = airdensity * 9.81 * flowvelocity * effectivedepth.max(1.0);
let lakeradiation = self.cryolake.longwaveradiationwm2(0.9);
let lakethermal = self.cryolake.thermalenergyjoulesrelative(ONEBARTEMPK);
self.cloudtoptempk -= lakeradiation * 1e-12 * simdt;
self.cloudtoptempk = self.cloudtoptempk.clamp(110.0, 200.0);
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::Impactor::asteroid(200.0, 60.0);
let craterd = impactor.craterdiameterm(STONYASTEROIDDENSITY);
let ejecta = impactor.ejectavolumem3(STONYASTEROIDDENSITY);
let fireball = impactor.fireballradiusm();
let energymt = impactor.kineticenergymt();
(ejecta / jupitersurfacearea)
* (1.0 + craterd * 1e-6)
* (1.0 + fireball * 1e-9)
* (1.0 + energymt * 1e-12)
} else {
0.0
}
} else {
0.0
};
let lightumol = effectiveirradiance * 0.45 * 4.57;
let chemo =
self.aerochem
.chemosynthesisrate(self.climate.nh3fraction * 1e6, lightumol, airtempk);
let column = self.aerochem.columnproduction(chemo);
let transpiration = self.aerochem.volatilefluxmmday(
atmosphere::heatwaves::saturationvaporpressureammoniana(airtempk) / 1000.0,
0.001,
);
let npp = self.aerochem.nppkgcm2yr(column);
let ecosystemnpp = self.ecosystem.totalnppgcyr();
let diversity = self.ecosystem.shannonindex();
let nppratio = if self.referencenpp.abs() > 1e-30 {
(npp / self.referencenpp * (diversity + 1.0)).clamp(0.01, 10.0)
} else {
1.0
};
self.cloudinteraction.producers.carryingcapacity = self.basecarryingcap * nppratio;
self.cloudinteraction.step(simyears);
self.calendar = temporal::calendar::DateTime::fromjuliandate(self.epoch.juliandate);
let cellvelocity = self.convectioncell.velocityatpoint(lat, lon);
self.cumulativeicestrain +=
cellvelocity * simdt / geology::plate_tectonics::METALLICHYDROGENTHICKNESS;
let rootdepth = self.interiorfeature.rootdepthm();
let isostaticelev = self.interiorfeature.isostaticelevation();
let strainthreshold = 1e-4;
let icequakeenergy = if self.cumulativeicestrain > strainthreshold {
let moment = self.cumulativeicestrain
* geology::plate_tectonics::INTERIORVISCOSITY
* self.convectioncell.areakm2
* 1e6;
let quake = geology::earthquakes::Icequake::frommoment(lat, lon, 100.0, moment);
let energy = quake.energyjoules();
let mw = quake.momentmagnitude();
let pga = quake.pgaatdistance(500.0);
self.cumulativeicestrain = 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(2.0, 0.01);
self.cryomagma.temperaturek += (geothermalinput * 0.001 - 0.005) * simdt;
self.cryomagma.temperaturek = self.cryomagma.temperaturek.clamp(100.0, 400.0);
let magmaviscosity = self.cryomagma.viscositypas();
let mgnum = self.cryomagma.mgnumber();
let eruptionpotential = self.cryomagma.h2owtpercent / magmaviscosity.log10().max(1.0);
let cryovei = if eruptionpotential > 0.05 {
geology::volcanism::cryovolcanicexplosivityindex(eruptionpotential * 0.001)
} else {
0
};
let streamerosion = streampower * 1e-15 * networksedimentfactor;
let effectiveerosion = if effectivewind > self.winderosionthresh {
self.erosionratems * fluvialerosionfactor * 2.0
+ icequakeenergy
+ thermalstress * 1e-6
+ streamerosion
+ impactcratererosion
} else {
self.erosionratems * fluvialerosionfactor
+ icequakeenergy
+ thermalstress * 1e-6
+ 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 cryovei > 0 {
*cell += cryovei as f64 * 0.005 * simyears;
}
}
let minelev = rootdepth
* (geology::plate_tectonics::METALLICHYDROGENDENSITY
/ geology::plate_tectonics::ROCKYCOREDENSITY
- 1.0);
let terrainh = self.heightmap.sample(lat, lon).max(minelev * 0.001);
let cryovolcanicfertility = mgnum * 0.001;
self.lodterrain.update([
self.observercam[0],
self.observercam[1],
self.heightmap.radiusat(lat, lon),
]);
let isinband = self.regions.regions.iter().any(|r| {
matches!(
r.regiontype,
geodata::regions::RegionType::AtmosphericBand
| geodata::regions::RegionType::PolarCap
) && r
.boundary
.first()
.is_some_and(|b| lat >= b[0].min(90.0) - 90.0)
});
let conductionfactor = if isinband { 0.8 } else { 1.0 };
let ammoniasat = atmosphere::heatwaves::saturationvaporpressureammoniana(airtempk);
let pressureratio = pressure / jupiters::ONEBARPRESSURE;
let ammoniafrac = ammoniasat * pressureratio * 0.0003;
let relativehumidity = (ammoniafrac / ammoniasat.max(1.0)).clamp(0.0, 1.0);
let moisture =
(relativehumidity * conductionfactor + transpiration * 0.001 + cryovolcanicfertility)
.clamp(0.0, 1.0);
let biome =
self.biomeclassifier
.classify(terrainh, lat, effectivewind * moisture.max(0.01));
let mat = &self.materials[cloudbiomeordinal(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, effectivewind * moisture.max(0.01));
let roughness = mat.roughness as f64 * splat.weights[0].1 as f64;
let mut total_uniform_bytes: usize = 0;
for uniform in &self.terrainshader.uniforms {
total_uniform_bytes += match &uniform.value {
rendering::shaders::UniformValue::Float(_) => 4,
rendering::shaders::UniformValue::Vec3(_) => 12,
rendering::shaders::UniformValue::Vec4(_) => 16,
rendering::shaders::UniformValue::Int(_) => 4,
};
}
for uniform in &self.atmoshader.uniforms {
total_uniform_bytes += match &uniform.value {
rendering::shaders::UniformValue::Float(_) => 4,
rendering::shaders::UniformValue::Vec3(_) => 12,
rendering::shaders::UniformValue::Vec4(_) => 16,
rendering::shaders::UniformValue::Int(_) => 4,
};
}
for uniform in &self.oceanshader.uniforms {
total_uniform_bytes += match &uniform.value {
rendering::shaders::UniformValue::Float(_) => 4,
rendering::shaders::UniformValue::Vec3(_) => 12,
rendering::shaders::UniformValue::Vec4(_) => 16,
rendering::shaders::UniformValue::Int(_) => 4,
};
}
for uniform in &self.bandshader.uniforms {
total_uniform_bytes += match &uniform.value {
rendering::shaders::UniformValue::Float(_) => 4,
rendering::shaders::UniformValue::Vec3(_) => 12,
rendering::shaders::UniformValue::Vec4(_) => 16,
rendering::shaders::UniformValue::Int(_) => 4,
};
}
total_uniform_bytes += self.terrainshader.name.len()
+ self.atmoshader.name.len()
+ self.oceanshader.name.len()
+ self.bandshader.name.len();
let mat_r0 = mat.fresnel_r0();
let sss_factor = mat.subsurface as f64;
let normal_detail = mat.normal_strength as f64;
let skylum = self.atmorender.sky_luminance(observeralt, coszenith) + roughness * 0.01;
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() + total_uniform_bytes;
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 emissivestrength = mat.emissive[0] as f64 * ndotl.max(0.0)
+ skylum * 0.01
+ mat_r0 * 0.001
+ sss_factor * 0.001
+ normal_detail * 0.001;
self.prevsurfacealbedo = albedo + emissivestrength * 0.001;
} else {
self.prevsurfacealbedo = albedo;
}
self.prevclouddensity = clouddensity;
self.prevbiome = biome;
self.biomeclassifier.onebarlevel += -icebalance * simyears * 0.001;
self.climate.nh3fraction -= ecosystemnpp * 1e-18 * simyears;
self.climate.nh3fraction = self.climate.nh3fraction.max(0.0);
let interiorsound = self.oceanlayer.soundspeedms();
let interiorheat = self.oceanlayer.heatcontentjperm2();
let vaporammonia = atmosphere::heatwaves::saturationvaporpressureammoniana(airtempk);
let ioretro = self.iostate.isbinary();
let iograv = self.iostate.gravityat(iodist);
let stormke = atmosphere::storms::JupiterStorm::greatredspot().kineticenergydensity();
let stormrossby =
atmosphere::storms::JupiterStorm::greatredspot().rossbydeformationradius();
let cape = atmosphere::storms::capelayer(airtempk + 5.0, airtempk, 1000.0);
let beaufort = atmosphere::winds::beauforttoms(12);
let thermalwind =
atmosphere::winds::thermalwindshear(0.001, airtempk, lat, pressure * 0.1, pressure);
let scaleheightval = atmosphere::layers::scaleheight();
let surfvel = self.rotation.surfacevelocityatlatitude(lat);
let rotke = self.rotation.rotationalkineticenergy();
let icelayervol = self.icelayer.volumem3(1e6);
let icestrain = hydrology::glaciers::glenstrainrate(ONEBARTEMPK, icebasal);
let cellcmpyr = self.convectioncell.velocitycmperyear(lat, lon);
let depthtemp = geology::plate_tectonics::temperatureatdepth(ONEBARTEMPK, 1e6, 0.01);
let interiorbalance = geology::plate_tectonics::isostaticbalance(self.icelayer.thicknessm);
let icerheology = hydrology::glaciers::icerheologyviscosity(ONEBARTEMPK, icebasal.max(1.0));
let sublimation = geology::erosion::sublimationerosionrate(ONEBARTEMPK, ammoniasat);
let ecosimp = self.ecosystem.simpsondiversity();
let ecospecies = self.ecosystem.expectedspeciesfromarea(0.25, 1.0);
let dewpt = atmosphere::heatwaves::dewpointammoniak(ammoniasat);
let heatanomaly = atmosphere::heatwaves::ThermalAnomaly {
temperaturek: airtempk,
durationdays: 0,
altitudem: observeralt,
};
let dangerlevel = heatanomaly.dangerlevel();
let cryodiff = self.cryomagma.differentiationidx();
let cryosolidus = self.cryomagma.solidusdepressionc();
let cryoliquidus = self.cryomagma.liquidustempc();
let cryototalalkali = self.cryomagma.totalalkali();
let lakemean = self.cryolake.meandepthm();
let lakeresidence = self.cryolake.residencetimeyears(1e3);
let notablefeatures = geodata::elevation::ElevationProvider::notablefeatures();
let interiorstats = geodata::bathymetry::InteriorLayerData::interiorstats();
let majorlayers = geodata::bathymetry::InteriorLayer::majorlayers();
let seasonname = season.seasonname();
let subsolarpt = (
season.subsolarlatitudedeg(),
season.jupitersolarlongitudedeg(),
);
let springtide = physics::tides::springtideamplitude();
let neaptide = physics::tides::neaptideamplitude();
let tideratio = physics::tides::iotosolartideratio();
let escapev = physics::orbit::JupiterOrbit::escape_velocity_at_surface();
let gravsun = self.orbit.gravitational_force_sun();
let specenergy = self.orbit.specific_orbital_energy();
let speclmomentum = self.orbit.specific_angular_momentum();
let co2forcing = atmosphere::climate::co2doublingforcing();
let spectralrad = atmosphere::climate::planckspectralradiance(10e-6, self.cloudtoptempk);
let nh3forcing = self.climate.nh3radiativeforcing();
let sunfluxjupiter = atmosphere::climate::solarfluxatjupiter();
let effectiveemission = self.climate.effectiveemissiontemperature();
let gheffect = self.climate.ammoniagreenhouseeffectk();
let galileanmasssum = satellites::io::IOMASS
+ satellites::europa::info().1
+ satellites::ganymede::info().1
+ satellites::callisto::info().1;
let consumed = effectiveerosion
+ iceshellvel
+ icebalance
+ tidalheight
+ interiorcurrent
+ lakethermal * 1e-30
+ interiorheat * 1e-30
+ vaporammonia * 1e-30
+ interiorsound * 1e-30
+ stormke * 1e-30
+ stormrossby * 1e-30
+ cape * 1e-30
+ beaufort * 1e-30
+ thermalwind * 1e-30
+ scaleheightval * 1e-30
+ surfvel * 1e-30
+ rotke * 1e-30
+ icelayervol * 1e-30
+ icestrain * 1e-30
+ cellcmpyr * 1e-30
+ depthtemp * 1e-30
+ interiorbalance * 1e-30
+ icerheology * 1e-30
+ sublimation * 1e-30
+ ecosimp * 1e-30
+ ecospecies * 1e-30
+ dewpt * 1e-30
+ cryodiff * 1e-30
+ cryosolidus * 1e-30
+ cryoliquidus * 1e-30
+ cryototalalkali * 1e-30
+ lakemean * 1e-30
+ lakeresidence * 1e-30
+ subsolarpt.0 * 1e-30
+ subsolarpt.1 * 1e-30
+ springtide * 1e-30
+ neaptide * 1e-30
+ tideratio * 1e-30
+ escapev * 1e-30
+ gravsun * 1e-30
+ specenergy * 1e-30
+ speclmomentum * 1e-30
+ co2forcing * 1e-30
+ spectralrad * 1e-30
+ nh3forcing * 1e-30
+ sunfluxjupiter * 1e-30
+ effectiveemission * 1e-30
+ gheffect * 1e-30
+ galileanmasssum * 1e-60;
if consumed.is_nan() {
debug_assert!(false, "NaN in simulation");
}
let iobin = if ioretro { 1.0 } else { 0.0 };
let dangerord = match dangerlevel {
"Extreme anomaly" => 4.0,
"Strong anomaly" => 3.0,
"Moderate anomaly" => 2.0,
_ => 1.0,
};
let featuresum: f64 = notablefeatures.iter().map(|f| f.3).sum();
let layersum: f64 = majorlayers.iter().map(|l| l.avgdensitykgm3).sum();
let statsval = interiorstats.totalradiuskm;
let seasonord = match seasonname {
"Northern Spring" => 1.0,
"Northern Summer" => 2.0,
"Northern Autumn" => 3.0,
_ => 4.0,
};
let budget = consumed
+ iobin * 1e-30
+ iograv * 1e-30
+ dangerord * 1e-30
+ featuresum * 1e-30
+ layersum * 1e-30
+ statsval * 1e-30
+ seasonord * 1e-30
+ hydraulicradius * 1e-30
+ flowlength * 1e-30
+ cloud_advection * 1e-30
+ barometric_density * 1e-60
+ atmo_shell_volume * 1e-60
+ r0_ocean * 1e-30
+ dx_ocean * 1e-30
+ depth_color * 1e-30
+ pressure_opacity * 1e-30
+ ocean_thermal_mass * 1e-60;
if renderbudget < 64 {
self.lodterrain.update(self.observercam);
}
if budget.is_nan() {
debug_assert!(false, "NaN in simulation budget");
}
}
}
fn main() {
let mut jupiter = Jupiter::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;
jupiter.tick(realdt);
let elapsed = now.elapsed();
if elapsed < targetdt {
std::thread::sleep(targetdt - elapsed);
}
}
}