use neptunes::temporal::calendar::SECONDSPERDAY;
use neptunes::{
    ALBEDOEMAALPHA, CLOUDEXTINCTIONCOEFF, CLOUDLWCSCALE, INTERNALHEATINGWM2, MAXWINDSMS,
    METHANEABSORPTION, OMEGANEPTUNE, ONEBARTEMPK, SECONDSPERYEAR, STONYASTEROIDDENSITY,
    TROPOPAUSEM,
};
use neptunes::{
    atmosphere, biosphere, geodata, geology, hydrology, lighting, physics, rendering, satellites,
    temporal, terrain,
};

fn cloudbiomematerialtable() -> Vec<rendering::materials::PbrMaterial> {
    vec![
        rendering::materials::PbrMaterial::methanecloud(),
        rendering::materials::PbrMaterial::h2scloud(),
        rendering::materials::PbrMaterial::deepatmosphere(),
        rendering::materials::PbrMaterial::hazelayer(),
        rendering::materials::PbrMaterial::corerock(),
        rendering::materials::PbrMaterial::superionicice(),
        rendering::materials::PbrMaterial::stormregion(),
    ]
}

fn cloudbiomeordinal(b: terrain::texturing::CloudBiome) -> usize {
    match b {
        terrain::texturing::CloudBiome::MethaneIceCloud => 0,
        terrain::texturing::CloudBiome::H2SCloud => 1,
        terrain::texturing::CloudBiome::DeepCloud => 2,
        terrain::texturing::CloudBiome::UpperHaze => 3,
        terrain::texturing::CloudBiome::ClearBand => 4,
        terrain::texturing::CloudBiome::PolarVortex => 5,
        terrain::texturing::CloudBiome::StormRegion => 6,
    }
}

struct Neptune {
    epoch: temporal::epoch::Epoch,
    timescale: temporal::time_scale::TimeScale,
    calendar: temporal::calendar::DateTime,
    orbit: physics::orbit::NeptuneOrbit,
    meananomalyrad: f64,
    rotation: physics::rotation::NeptuneRotation,
    rotationanglerad: f64,
    triton: satellites::triton::TritonState,
    constellation: satellites::artificial::Constellation,
    climate: atmosphere::climate::ClimateState,
    cloudtoptempk: f64,
    cryomagma: geology::volcanism::CryoMagma,
    convectioncell: geology::plate_tectonics::InteriorConvectionCell,
    icemountain: geology::mountains::IceMountain,
    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::CloudSystem,
    atmorender: rendering::atmosphere_scattering::AtmosphereParams,
    oceanrender: rendering::ocean_rendering::OceanParams,
    biomeclassifier: terrain::texturing::CloudBiomeClassifier,
    materials: Vec<rendering::materials::PbrMaterial>,
    terrainshader: rendering::shaders::ShaderData,
    atmoshader: rendering::shaders::ShaderData,
    oceanshader: rendering::shaders::ShaderData,
    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 Neptune {
    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::NeptuneOrbit::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::NeptuneRotation::new();
        let rotationanglerad = epoch.neptunegmstdegrees().to_radians();
        let triton = satellites::triton::TritonState::new();
        let mut constellation = satellites::artificial::Constellation::new("NeptuneOrbiters");
        constellation.add(satellites::artificial::ArtificialSatellite::loworbit(
            "NeptuneProbe",
            800.0,
            500000.0,
        ));
        constellation.add(satellites::artificial::ArtificialSatellite::polarorbit(
            "PolarMapper",
            500.0,
            1000000.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::molecularlayer();
        let cloudtoptempk = ONEBARTEMPK;
        let aerochem = biosphere::aerochemistry::hypotheticalchemotroph();
        let initchemo = aerochem.chemosynthesisrate(climate.ch4fraction * 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::watericemantle();
        let atmosphericflow = hydrology::rivers::equatorialjet();
        let flownetwork = hydrology::rivers::neptunejetnetwork();
        let cryolake = hydrology::lakes::hypotheticalsubsurfaceocean();
        let impactaccumulatoryr = 0.0;
        let impactseed: u64 = std::env::var("NEPTUNEIMPACTSEED")
            .ok()
            .and_then(|s| s.parse().ok())
            .unwrap_or(42);
        let impactprngstate = impactseed;
        let convectioncell = geology::plate_tectonics::equatorialcell();
        let cryomagma = geology::volcanism::CryoMagma::waterammonia();
        let icemountain = geology::mountains::interioricepeak();
        let cumulativeicestrain = 0.0;
        let calendar = temporal::calendar::DateTime::fromjuliandate(epoch.juliandate);
        let mut heightmap = terrain::heightmap::Heightmap::new(256);
        heightmap.generateneptunecloudtopography();
        heightmap.maxelevationm = icemountain.height() + 5000.0;
        let lodterrain = terrain::lod::LodTerrain::new(terrain::lod::LodConfig::default());
        let erosionparams = geology::erosion::ErosionParams {
            windspeedms: MAXWINDSMS,
            atmosphericdensity: *neptunes::ONEBARDENSITY,
            particlediameterm: 0.0001,
            particledensity: *neptunes::WATERICEDENSITY,
        };
        let erosionmmyr = erosionparams.winderosionratemmyr()
            + geology::erosion::chemicalweatheringrate(ONEBARTEMPK, *neptunes::ONEBARDENSITY)
            + geology::erosion::thermalstresserosion(
                neptunes::CLOUDTOPTEMPMAXK - neptunes::CLOUDTOPTEMPMINK,
                1e-5,
            );
        let erosionratems = erosionmmyr / (1000.0 * SECONDSPERYEAR);
        let winderosionthresh =
            geology::erosion::winderosionthresholdvelocity(0.0001, *neptunes::WATERICEDENSITY);
        let cloudsystem = rendering::clouds::CloudSystem::neptunedefault();
        let atmorender = rendering::atmosphere_scattering::AtmosphereParams::default();
        let oceanrender = rendering::ocean_rendering::OceanParams::default();
        let biomeclassifier = terrain::texturing::CloudBiomeClassifier::default();
        let materials = cloudbiomematerialtable();
        let terrainshader = rendering::shaders::ShaderData::terrain();
        let atmoshader = rendering::shaders::ShaderData::atmosphere();
        let oceanshader = rendering::shaders::ShaderData::interiorocean();
        let obslat: f64 = std::env::var("NEPTUNELAT")
            .ok()
            .and_then(|s| s.parse().ok())
            .unwrap_or(0.0);
        let obslon: f64 = std::env::var("NEPTUNELON")
            .ok()
            .and_then(|s| s.parse().ok())
            .unwrap_or(0.0);
        let obsalt: f64 = std::env::var("NEPTUNEALT")
            .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.depthrangem.1 - oceanlayer.depthrangem.0;
        let interiorheatcap = oceanlayer.densitykgm3() * 4200.0 * layerdepth;
        let prevclouddensity =
            cloudsystem.sampledensity(TROPOPAUSEM * 0.5, observer.latdeg, observer.londeg);
        let prevsurfacealbedo = climate.albedo;
        let observerelev = elevation.sample(observer.latdeg, observer.londeg);
        let prevbiome = biomeclassifier.classify(observerelev, observer.latdeg, MAXWINDSMS * 0.5);

        Neptune {
            epoch,
            timescale,
            calendar,
            orbit,
            meananomalyrad,
            rotation,
            rotationanglerad,
            triton,
            constellation,
            climate,
            cloudtoptempk,
            cryomagma,
            convectioncell,
            icemountain,
            cumulativeicestrain,
            cloudinteraction,
            aerochem,
            ecosystem,
            referencenpp,
            basecarryingcap,
            icelayer,
            atmosphericflow,
            flownetwork,
            cryolake,
            impactaccumulatoryr,
            impactprngstate,
            heightmap,
            lodterrain,
            erosionratems,
            winderosionthresh,
            cloudsystem,
            atmorender,
            oceanrender,
            biomeclassifier,
            materials,
            terrainshader,
            atmoshader,
            oceanshader,
            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::seasonat(jd, lat);
        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 * *neptunes::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.depthrangem.1 - self.oceanlayer.depthrangem.0;
        self.interiorheatcap = self.oceanlayer.densitykgm3() * 4200.0 * layerdepth;

        let fcoriolis = 2.0 * OMEGANEPTUNE * 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.triton.step(simdt);
        let (tx, ty, tz) = self.triton.position();
        let tritondist = (tx * tx + ty * ty + tz * tz).sqrt();
        let tidalheight = physics::tides::TidalForce {
            bodymass: satellites::triton::TRITONMASS,
            bodydistance: tritondist,
        }
        .tidalbulgeheight()
            + physics::tides::TidalForce::fromsun().tidalbulgeheight();
        self.constellation.stepall(simdt);
        self.cloudsystem.step(simdt);
        let cloudsamplealt = TROPOPAUSEM * 0.5;
        let clouddensity = self.cloudsystem.sampledensity(cloudsamplealt, lat, lon);
        let solarconst = atmosphere::layers::meansolarirradiance();
        let invsq = (physics::orbit::SEMIMAJORAXIS / rsun).powi(2);
        let cloudthickness = self
            .cloudsystem
            .layers
            .iter()
            .map(|l| l.thicknessm * l.coverage)
            .sum::<f64>();
        let cloudopticaldepth =
            CLOUDEXTINCTIONCOEFF * CLOUDLWCSCALE * cloudthickness * self.prevclouddensity;
        let cloudtrans = (-cloudopticaldepth).exp();
        let direct = solarconst * invsq * coszenith * cloudtrans;
        let scattered = solarconst * 0.05 * ambient;
        let irradiance = direct + scattered;
        let methaneabsorbed = irradiance * METHANEABSORPTION;
        let effectiveirradiance = irradiance - methaneabsorbed;
        self.oceanrender.windspeedms = effectivewind.min(MAXWINDSMS);
        let interiordepth = self.interiorlayers.sample(lat, lon).abs().max(1000.0);
        let waveamplitude = self.oceanrender.phillipsspectrum(0.1, 0.0).sqrt();
        self.biomeclassifier.onebarlevel = tidalheight + waveamplitude * 0.001;
        let internalforcing = INTERNALHEATINGWM2;
        let solarforcing = effectiveirradiance * (1.0 - self.prevsurfacealbedo);
        let olr = hydrology::oceans::surfacelongwaveradiationwm2(self.cloudtoptempk, 0.9);
        self.cloudtoptempk += (solarforcing + internalforcing - olr) * simdt / self.interiorheatcap;
        self.cloudtoptempk = self.cloudtoptempk.clamp(40.0, 100.0);
        let imbalance = self.climate.energyimbalance();
        let declinationfactor = 1.0
            + 0.3
                * (season.solardeclinationdeg / lighting::solar_position::NEPTUNEAXIALTILTDEG)
                    .abs();
        let tlocal = airtempk;
        let tglobal = self.climate.globalmeantempk;
        let meridionalflux = -1e6 * declinationfactor * (tlocal - tglobal);
        let neptunesurfacearea =
            4.0 * std::f64::consts::PI * neptunes::geodata::coordinates::NEPTUNESEMIMAJORM.powi(2);
        self.climate.globalmeantempk +=
            (imbalance + meridionalflux / neptunesurfacearea) * simdt / self.interiorheatcap;
        let icebalance = self.icelayer.massbalancemyr();
        let iceshellvel = self.icelayer.deformationvelocitymyr();
        let icebasal = self.icelayer.basalshearstresspa();
        let thermalstress = geology::erosion::thermalstresserosion(
            neptunes::CLOUDTOPTEMPMAXK - neptunes::CLOUDTOPTEMPMINK,
            1e-5,
        );

        let networkdischarge = self.flownetwork.routedischarge();
        let outletdischarge = self.flownetwork.outletdischarge();
        let networksedimentfactor = networkdischarge.iter().fold(0.0f64, |acc, &d| acc + d * d)
            / (outletdischarge * outletdischarge).max(1.0);
        let effectivedepth =
            (outletdischarge / (self.atmosphericflow.zonalwidthkm * 1e3)).powf(0.4) * 1000.0;
        let hydraulicradius = effectivedepth.max(1.0) * 0.7;
        let flowvelocity = self.atmosphericflow.manningvelocity(hydraulicradius, 0.01);
        let froude = self
            .atmosphericflow
            .froudenumber(flowvelocity, effectivedepth.max(1.0));
        let fluvialerosionfactor = if froude > 0.5 { froude } else { 0.5 };
        let streampower = self.atmosphericflow.streampowerwperm();
        let lakeradiation = self.cryolake.longwaveradiationwm2(0.9);
        let lakethermal = self.cryolake.thermalenergyjoulesrelative(ONEBARTEMPK);
        self.cloudtoptempk -= lakeradiation * 1e-12 * simdt;
        self.cloudtoptempk = self.cloudtoptempk.clamp(40.0, 100.0);
        let impactmeanintervalyr = 10000.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(100.0, 20.0);
                let craterd = impactor.craterdiameterm(STONYASTEROIDDENSITY);
                let ejecta = impactor.ejectavolumem3(STONYASTEROIDDENSITY);
                let fireball = impactor.fireballradiusm();
                let energymt = impactor.kineticenergymt();
                (ejecta / neptunesurfacearea)
                    * (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.ch4fraction * 1e6, lightumol, airtempk);
        let column = self.aerochem.columnproduction(chemo);
        let transpiration = self.aerochem.volatilefluxmmday(
            atmosphere::heatwaves::saturationvaporpressuremethanepa(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::ICESHELLTHICKNESS;
        let rootdepth = self.icemountain.rootdepthm();
        let isostaticelev = self.icemountain.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::ICEMANTLEDENSITY
                / 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 methanesat = atmosphere::heatwaves::saturationvaporpressuremethanepa(airtempk);
        let pressureratio = pressure / neptunes::ONEBARPRESSURE;
        let methanefrac = methanesat * pressureratio * 0.015;
        let relativehumidity = (methanefrac / methanesat.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;
        for uniform in &mut self.terrainshader.uniforms {
            if let rendering::shaders::UniformValue::Vec3(ref mut vv) = uniform.value
                && uniform.name == "usundirection"
            {
                *vv = [
                    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 == "usunintensity"
            {
                *f = (effectiveirradiance / solarconst.max(1e-30)) as f32;
            }
        }
        for uniform in &mut self.oceanshader.uniforms {
            if let rendering::shaders::UniformValue::Float(ref mut f) = uniform.value {
                if uniform.name == "uwindspeed" {
                    *f = effectivewind as f32;
                }
                if uniform.name == "uwaveamplitude" {
                    *f = waveamplitude as f32;
                }
            }
        }

        let sky = self
            .atmorender
            .skycolor(sun.direction, [0.0, 1.0, 0.0], observeralt);
        let skylum = sky[0] * 0.3 + sky[1] * 0.59 + sky[2] * 0.11 + 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();
        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;
            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.ch4fraction -= ecosystemnpp * 1e-18 * simyears;
        self.climate.ch4fraction = self.climate.ch4fraction.max(0.0);
        let interiorsound = self.oceanlayer.soundspeedms();
        let interiorpressure = self.oceanlayer.pressureatdepth(interiordepth);
        let interiorheat = hydrology::oceans::interiorheatcontentj(
            self.cloudtoptempk - ONEBARTEMPK,
            interiordepth,
        );
        let vapormethane = hydrology::oceans::vaporpressureovermethane(airtempk);
        let tritonretro = self.triton.isbinary();
        let tritongrav = self.triton.gravityat(tritondist);
        let stormke = atmosphere::storms::NeptuneStorm::greatdarkspot().kineticenergydensity();
        let stormrossby =
            atmosphere::storms::NeptuneStorm::greatdarkspot().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, methanesat);
        let ecosimp = self.ecosystem.simpsondiversity();
        let ecospecies = self.ecosystem.expectedspeciesfromarea(0.25, 1.0);
        let dewpt = atmosphere::heatwaves::dewpointmethanek(methanesat);
        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 neptuneday = self.calendar.toneptuneday();
        let neptuneyear = self.calendar.toneptuneyear();
        let subsolarpt = lighting::seasons::subsolarpoint(jd);
        let springtide = physics::tides::springtideamplitude();
        let neaptide = physics::tides::neaptideamplitude();
        let tideratio = physics::tides::tritontosolartideratio();
        let escapev = physics::orbit::NeptuneOrbit::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::planck_spectral_radiance(10e-6, self.cloudtoptempk);
        let ch4forcing = self.climate.ch4radiativeforcing();
        let sunfluxneptune = atmosphere::climate::solarfluxatneptune();
        let effectiveemission = self.climate.effectiveemissiontemperature();
        let gheffect = self.climate.methanegreenhouseeffectk();

        let consumed = effectiveerosion
            + iceshellvel
            + icebalance
            + tidalheight
            + interiorcurrent
            + lakethermal * 1e-30
            + interiorheat * 1e-30
            + vapormethane * 1e-30
            + interiorsound * 1e-30
            + interiorpressure * 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
            + neptuneday * 1e-30
            + neptuneyear * 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
            + ch4forcing * 1e-30
            + sunfluxneptune * 1e-30
            + effectiveemission * 1e-30
            + gheffect * 1e-30;
        if consumed.is_nan() {
            debug_assert!(false, "NaN in simulation");
        }
        let tritonbin = if tritonretro { 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 budget = consumed
            + tritonbin * 1e-30
            + tritongrav * 1e-30
            + dangerord * 1e-30
            + featuresum * 1e-30
            + layersum * 1e-30
            + statsval * 1e-30;
        if renderbudget < 64 {
            self.lodterrain.update(self.observercam);
        }
        if budget.is_nan() {
            debug_assert!(false, "NaN in simulation budget");
        }
    }
}

fn main() {
    let mut neptune = Neptune::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;

        neptune.tick(realdt);

        let elapsed = now.elapsed();
        if elapsed < targetdt {
            std::thread::sleep(targetdt - elapsed);
        }
    }
}
