commit acc0e5eec9347330aa96f14f0f00fe3c7d9773a9 Author: Jonas Reith Date: Tue Jun 23 15:08:25 2026 +0200 Initial commit: fanworgen planet sim C++/raylib semi-realistic fantasy/sci-fi planet generator on a fixed icosphere grid (Eulerian: properties flow over fixed cells). World-creation stages: tectonics, continental drift & erosion, hydrology (rivers/lakes), climate (temperature + orographic precipitation), and biome classification. Engine in src/sim (raylib-free, headless-testable), viewer in src/render. See CLAUDE.md and docs/ (design-notes.md, fauna-flora-plan.md = next step). Co-Authored-By: Claude Opus 4.8 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6a03212 --- /dev/null +++ b/.gitignore @@ -0,0 +1,21 @@ +# Build output (CMake; also holds the FetchContent'd raylib source/objects) +/build/ + +# Runtime-generated, machine-specific +planet.cfg +planet.save +screenshot*.png + +# Claude Code local (machine-specific) settings +.claude/settings.local.json + +# Object/artifact cruft +*.o +*.obj +*.a + +# Editor / OS +.vscode/ +.idea/ +*.swp +.DS_Store diff --git a/BUILD.md b/BUILD.md new file mode 100644 index 0000000..d9e4976 --- /dev/null +++ b/BUILD.md @@ -0,0 +1,168 @@ +# Planet Sim - World Creation + +Fixed icosphere geometry; properties (elevation, plate, age, climate, biome) flow +over it. World creation runs as continuous, overlapping stages on a geological clock: +tectonics (boundary stress forms mountains/trenches) → continental drift & erosion → +hydrology (rivers/lakes) → climate (temperature/precipitation) → biomes. Initial +terrain forming is a generator (a couple hundred paced ticks, ~3 s, auto-pauses at +equilibrium), then drift/erosion/etc. continue. Press R to reseed, SPACE to pause. +(The eventual goal is a separate slow real-time "Live World" weather/life mode.) + +## Build (Nobara / Linux) + +Dependencies (raylib build needs these dev headers): + + sudo dnf install cmake gcc-c++ mesa-libGL-devel libX11-devel \ + libXrandr-devel libXinerama-devel libXcursor-devel libXi-devel \ + wayland-devel libxkbcommon-devel + +Then: + + cmake -B build -DCMAKE_BUILD_TYPE=Release + cmake --build build -j + ./build/planetsim + +raylib 5.5 is fetched automatically via CMake FetchContent. + +OpenMP parallelizes the tectonic step (auto-detected by CMake; ships with +gcc, no extra package needed). It's optional -- without it the sim still runs, +serial. The step is memory-bandwidth-bound, so `OMP_NUM_THREADS=4` to `8` gives +the full ~2.8x speedup; the default uses all cores for no extra gain: + + OMP_NUM_THREADS=6 ./build/planetsim + +## Controls + + LMB drag orbit camera (in the 3D view) + wheel zoom + hover show cell info (works in the 3D globe and the 2D map) + click open tile detail panel (subtiles grid, hoverable) + overlay + C close the detail panel + 1 .. 7 color by elevation / plate / age / crust type / biome / temperature / precipitation + B toggle plate borders (on by default) + D toggle per-plate drift arrows + P labels (on by default) + G toggle lat/lon graticule (+ degree numbers on the 2D map edges) + J toggle rivers (Phase 2.5 hydrology) + H start/stop Phase 2.5 (hydrology: rivers, lakes, fluvial erosion) + SPACE pause while forming / re-evolve once settled (or the on-screen button) + [ / ] drift speed (My per real second, Phase 2/3) + S single tectonic tick + F fast-forward Phase-1 forming to settled (instant) + R reseed planet (restart forming) + + / - subdivision level (detail), 1..7 + F5 / F9 save / load full state (planet.save) + F12 screenshot to screenshot.png + F2 reload planet.cfg (validated) and regenerate + +CLI flags (applied before the first load/generate): + + --seed N override the config seed + --config PATH use an alternate config file instead of planet.cfg + +## Files (written in the working directory) + + planet.cfg human-editable key=value config of every PlanetConfig parameter; + auto-created on first run, reload live with F2. Range-checked on + load; an invalid file reverts to safe defaults (not overwritten). + planet.save binary snapshot (versioned, currently v6): seed + config + full planet + state; F5 writes it, F9 reloads and resumes deterministically. As of v6 + the config is stored as a self-describing key=value block (like + planet.cfg), so adding/removing config fields no longer breaks saves + (unknown keys ignored, missing keys default). v6 cannot load pre-v6 + saves (one-time break) -- regenerate them. + +## Config (planet.cfg, or PlanetConfig defaults in code) + + radius 6.371e6 m (Earth) -- free to change + subdivisions 5 (~10k cells, ~223 km/cell) + plateCount 12 + seaLevel 0 m + axialTilt 23.44 deg obliquity: leans the 3D globe + spin axis (visual; seasons later) + +Tectonic tuning (also PlanetConfig): relief builds gradually toward an +isostatic equilibrium instead of saturating to the clamp. + + continentBase 300 m resting elevation of continental crust + oceanBase -6000 m deep abyssal floor (oldest oceanic crust) + upliftGain 1.3e5 m/tick per unit convergence stress + beltWidth 3 cell-rings a mountain belt spreads inland + relax 0.02 isostatic relaxation toward base, per tick + +Orogeny + seafloor (PlanetConfig): tall, persistent mountains and age-based +ocean depth. The orogeny boosts (collision/arc/persistence) are drift-only -- +Phase 1 forms as before; tall mountains grow during Phase-2 drift. + + collisionFactor 1.8 continent-continent uplift (Himalaya); raise = taller + arcFactor 1.4 continental subduction-arc uplift (Andes) + isostaticPersist 0.85 how strongly high crust resists relax (0..<1) + rootScale 2500 m above continentBase where persistence saturates + ridgeDepth -2500 m shallow elevation of brand-new crust at a ridge + seafloorSubsidence 280 m per sqrt(My): seafloor deepens with crustal age + seafloorSeedAge 80 My initial oceanic age spread at generation + +Hydrology (PlanetConfig, Phase 2.5): macro drainage network + fluvial erosion. +Drift keeps running at a finer timestep during hydrology. (Config keys keep the +`phase3*` names for save/cfg compatibility; the UI labels this "Phase 2.5".) + + phase3AfterMy 300 My drift before the "Start Phase 2.5?" prompt + phase3DtScale 0.2 hydrology timestep = cflDtMy()*this (finer = slower drift) + rainfall 1.0 uniform precip per cell (drainage-area unit) + riverThreshold 50 discharge above which a cell is a river + riverIncision 0.02 K in stream-power incision K*Q^m*S^n*dt (raise = carve faster) + riverDischargeExp 0.5 m (discharge exponent) + riverSlopeExp 1.0 n (slope exponent) + riverTransport 0.1 transport capacity coefficient (cap = this*Q*S) + depFrac 0.25 fraction of excess sediment deposited per cell + +Biomes (PlanetConfig, Phase 3): per-cell biome classification thresholds. Temperature += biomeEquatorTemp - biomePoleDrop*(|lat|/90)^biomeLatExp - biomeElevLapse*elevAbove. + + biomeEquatorTemp 30 C temperature at the equator, sea level + biomePoleDrop 58 C equator->pole temperature drop + biomeLatExp 1.3 >1 keeps mid-latitudes temperate (cold concentrates at poles) + biomeElevLapse 0.006 C lost per metre above sea level + biomeIceTemp -9.5 C below -> Ice (polar caps + snowcaps); raise = bigger caps + biomeTundraTemp 2 C below (and above ice) -> Tundra/Taiga + biomeTaigaTemp 10 C cool + wet -> boreal forest + biomeSavannaTemp 22 C warm + moderate moisture -> savanna + biomeMountainElev 3000 m above sea level -> Mountains + biomeHillsElev 1200 m above sea level -> Hills + biomeBeachBand 60 m above sea level + adjacent ocean -> Beach + biomeLowlandElev 500 m wetlands only below this elevation + biomeWetlandMoist 0.72 moisture above this (low lowland) -> Wetland + biomeDesertMoist 0.28 moisture below this -> Desert + biomeGrassMoist 0.50 moisture below this -> Grassland/Savanna, else Forest + biomeTaigaMoist 0.40 cool + above this -> Taiga (else Tundra) + biomeLakeMinDepth 20 m filled-basin depth above sea level counting as a Lake + +Climate (PlanetConfig, Phase 3): temperature uses the biome* temp params above; +precipitation advects ocean moisture along zonal winds (windward rain, leeward rain +shadow, dry interiors) then diffuses it. Color modes 6 (temperature) / 7 (precipitation). + + climateOceanMoisture 1.0 moisture air carries leaving the ocean (source) + climateRainEfficiency 0.5 fraction of available moisture*belt that rains per cell + climateOrographic 3.0 extra rain per unit normalized upslope (windward) + climateOroRefHeight 500 m upslope that counts as one orographic unit + climateContinentality 0.05 fractional moisture lost per inland cell (dries interiors) + climateWindPasses 50 moisture-advection iterations (steady state) + climateMoistureSmooth 12 precipitation diffusion passes (raise = smoother, more grass/forest) + +## Headless logic test (no display) + + g++ -std=c++17 -O2 -Isrc/sim test_logic.cpp src/sim/IcoSphere.cpp \ + src/sim/Planet.cpp src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp \ + src/sim/PlanetErosion.cpp src/sim/PlanetHydrology.cpp \ + src/sim/PlanetBiomes.cpp src/sim/PlanetIO.cpp \ + -o /tmp/t && /tmp/t + +Verifies geometry, plate assignment, gradual non-saturating relief and +determinism. Run after changing Planet::step(). + +## Notes for later phases + +- Cell has a `subgrid` shared_ptr hook (null in phase 1) for the future + fine-resolution per-cell mesh (phases 4/5: civilization/culture). +- `neighbors` adjacency already built -> reuse for diffusion, climate, + cross-cell-boundary interaction. +- elevation is continuous meters (double); climate in phase 3 can read it + as a smooth function, not coarse bands. diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..56c64d7 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,520 @@ +# CLAUDE.md + +Guidance for Claude Code when working in this repository. + +## Project + +A C++/raylib simulation/game for creating semi-realistic fantasy & sci-fi +worlds/planets. Built in phases. Between phases the user can edit the world +or trigger events (meteor impact, magical events, sci-fi terraforming). + +**Core principle:** the planet geometry is FIXED. Cells (vertices/faces) never +move. Only their *properties* flow over the fixed grid (Eulerian, not +Lagrangian). This keeps the data structure stable across all phases and makes +erosion, climate and drift far simpler to implement later. + +## Roadmap + +The project has two big arcs. **World Creation** builds a plausible planet through a +set of geological/environmental stages that **overlap and run together on a geological +clock** (My) — they are not strict sequential "phases" (the code keeps `phase*` names +internally for save/config compatibility, but think of them as continuous stages). Once +a world is "done", the long-term goal is a separate **Live World** mode that runs the +finished planet at a *much* slower, real-time-ish clock (hours/days/weeks/months) with +dynamic weather and life. + +**World Creation (geological clock, mostly done):** +- **Tectonics & landmass** *(done)* — icosphere geometry, plate assignment, boundary + stress forming mountains (convergent) and trenches/rifts (divergent). The initial + "forming" pass settles to isostatic equilibrium, then plate motion continues. +- **Continental drift & erosion** *(done)* — real plate motion (cm/yr, My), the plate + lifecycle (fission, stalemate kick, spreading-born plates, merging), diffusive + erosion + a sea-level controller (~30% land), taller persistent mountains + (collision + arc + isostatic persistence), seafloor aging→depth. `planet.cfg` + + `planet.save`. See [[phase2-design-direction]]. +- **Hydrology** *(done)* — rivers, lakes and fluvial erosion as a macro drainage network + (depression-fill→lakes, steepest-descent→rivers, mass-conserving stream-power + incision). Runs at a finer timestep alongside continuing drift. See [[phase3-hydrology]]. +- **Climate** *(done)* — continuous per-cell temperature + orographic precipitation + (`Planet::computeClimate`); runs from forming onward (a live "base climate" that + updates as terrain changes). Color modes `6`/`7`. +- **Biomes** *(done)* — per-cell `Cell.biome` (13 biomes incl. polar Ice) from elevation + + the climate fields (`Planet::classifyBiomes`), color mode `5`, saved per cell. +- **Fauna & flora** *(next, planned — see `docs/fauna-flora-plan.md`)* — derived + carrying-capacity densities: vegetation (flora) + a herbivore/carnivore food chain (fauna), + computed from the climate fields each tick (the living/evolving ecosystem is reserved for + Live World). Other follow-ups: feed precipitation into hydrology rainfall; seasons (obliquity). + +> Durable design context (module layout, save format, climate/biome model, conventions) +> lives in **`docs/design-notes.md`** — important because Claude's auto-memory does not +> travel with the repo. + +**Live World (future):** with planet creation finished, run the world at a slow real-time +scale with dynamic **weather** (clouds, rain, storms, fronts), day/night, and living +**ecosystems / civilization** evolving in real time. This is a separate large effort; +the fixed-grid Eulerian model + the climate fields are the groundwork for it. + +## Current state + +Working and verified (logic tested headless): +- Icosphere level 5 → 10242 cells, ~223 km/cell, topologically correct + (exactly 12 degree-5 vertices, rest degree-6). +- Plate flood-fill, drift as rotation on the sphere, boundary-stress uplift. +- Relief builds *gradually* toward an isostatic equilibrium (no clamp-rail + saturation): after ~40 ticks graded mountain belts (>2000 m) and deep + subduction trenches (<-6000 m), 0% of cells pinned to the clamp. +- `test_logic.cpp` asserts geometry, plate assignment, non-saturation and + graded relief; run it after any `Planet::step()` change (see below). +- raylib render: 3D globe (top) + 2D Equal Earth map (bottom strip), orbit + camera, 4 color modes. Phase 1 is a **generator**: it paces tectonic ticks + toward isostatic equilibrium (watchable, ~3 s), renders live, and + **auto-pauses** once the terrain settles (max per-tick change < 2 m). No + background thread — that's Phase-2 (continuous drift/erosion) groundwork. +- **Phase 2 increment 1 (drift):** after forming settles, the app switches to + continuous **drift mode** on a My clock. Each plate has a real speed (1..20 + cm/yr); `cflDtMy()` sets the timestep so the fastest plate advances ~half a + cell/step. `advect(dt)` moves plate membership + carried crust (plateId, + elevation, oceanic, geoAge) by an accumulation scheme: a boundary cell builds + `drift` (signed convergence distance) with its dominant other-plate neighbor; + +1 cell -> overrun (take that crust; but oceanic can't overrun a buoyant + continent -> it subducts under), -1 cell -> new young ridge crust (spreading). + Crust type now lives on the **cell** (`Cell.oceanic`), not the plate. The HUD + shows elapsed My; `[`/`]` set My/sec; auto-settle still gates Phase 1 -> 2. +- **Phase 2 increment 1.5 (plate lifecycle):** without this the count collapsed + to 1-3 and the world froze (a giant continental plate that can't shed cells). + Every `splitCheckEvery`(10) advect iterations (gated by `driftIter`) `advect()` + runs a Wilson-cycle lifecycle: **fission** (a plate over `splitFraction`(20%) + of cells splits along a random great circle through its centroid; prob ramps + `0.05+0.05*(pct-20)`, the new half gets a random drift — self-regulating, so + the count equilibrates ~8-12), **stalemate kick** (a plate whose size barely + changed over a window gets a new direction + speed boost to break deadlocks), + and **spreading plates** (rift cells become a `Plate.baby` young-ridge strip; + `coalesceBabyPlates()` merges connected blobs and dissolves tiny noise ones; + a strip past `babyPromoteFrac`(0.7%) is promoted to a real plate with random + drift + `volcanicLandFrac` of its interior turned to volcanic-island land). + The hard land clamp became a **soft band** (`landBand`) so volcanic land + persists; `acquirePlate()` reuses dead plate slots to keep `plates` bounded. + Stats panel separates real plates from "Young ridges: N strips"; plate-color + mode tints baby cells a uniform ridge grey. Verified by a headless soak. +- **Phase 2 increment 1.6 (calmer + cleaner plate map):** kicks were too twitchy + ("strange border movement") so a plate must now stall for `stalemateWindows`(4) + consecutive check windows (`sStaleStreak`) before a kick. `deleteEnclosedPlates()` + absorbs any plate ringed by a single other plate into it (on a mutual pair only + the smaller). `fuseMiniPlates()` lets >= `fuseMinPlates`(3) clustered mini plates + (non-baby, < `miniPlateCells`) fuse into the largest member and steal one ring of + cells from their largest big neighbour (terrane amalgamation). Borders now draw + in two colors: real plate borders yellow, young spreading-ridge borders red + (3D + 2D, both via `B`). Soak: mini plates ~0, no enclosed slivers, ~8-12 plates. +- **Phase 2 increment 2 (erosion + sea level):** `Planet::erode(dtMy)` runs each + drift step after advect+step — a mass-conserving, slope-weighted downhill + sediment transport (one double-buffered gather pass): the higher cell of each + edge gives material to the lower, faster above sea level (`erosionLandRate`) + than below (`erosionSeaRate`). Highs wear down toward an uplift<->erosion + equilibrium; sediment fills basins and builds coastal shelves/deltas. No river + carving (sub-cell at 223 km -> subgrid later). `adjustSeaLevel()` (every + `seaLevelEvery` erode calls) eases `cfg.seaLevel` toward the percentile + elevation leaving `landFractionTarget` (30%) of cells above water — percentile + targeting (`nth_element`) because a proportional nudge oscillates across the + flat continental-base elevation. Two "land" notions coexist: crust type (plate + buoyancy) vs geographic (elev > seaLevel). Stats headline + water:land are now + geographic, with separate "Sea level" and "Crust %" lines. Headless: land -> 30% + all seeds, erosion conserves sum(elevation), deterministic. +- **Phase 2 increment 3 (gradual sea level + config/save files):** the sea-level + controller is now gradual — checked every `seaLevelEvery` (100) erode calls, it + nudges sea level by a fixed `seaLevelStep` (100 m) when outside a `seaLevelTol` + (2%) deadband, and only if the nudge reduces the error (so it rests near a flat + "cliff" instead of oscillating). A human-editable **`planet.cfg`** (key=value) + holds all PlanetConfig params (auto-created on first run, `F2` reloads + + regenerates); `loadConfig`/`saveConfig` share one `CONFIG_FIELDS` X-macro. A + **`planet.save`** binary holds seed + config + full planet state + (`Planet::writeState`/`readState`; geometry rebuilt from subdivisions via + `buildGeometry()`); `F5` saves, `F9` loads and resumes (deterministic + continuation, verified headless). File I/O lives in Planet (raylib-free). +- **Phase 2 polish (config validation, QoL, sim cleanups):** `validateConfig()` + range-checks every field (+ the cross-rule `oceanBase < continentBase`) on load + and `F2`; an invalid `planet.cfg` reverts to safe defaults without overwriting + it. Save bumped to **version 2** (now also persists the `[`/`]` drift rate; + version-gated reads accept older saves). Viewer adds a **crust-type color mode** + (`4`: continental warm brown / oceanic deep blue), `F` **fast-forward** (runs + `step()` to settled instantly), `F12` screenshot, `--seed`/`--config` CLI flags, + a `P` label per plate on the drift arrows, and a min subdivision level of 1 + (level 0 disallowed). Two sim refinements: `step()` re-anchors boundary (source) + cells to their original stress each dilation ring so adjacent belts don't + cross-inflate (sharper peaks; land flank cells 509->499); and the soft land-band + rift/accrete nudge now reads a frozen snapshot of pre-nudge `oceanic` so flipping + one cell can't cascade along vertex-index chains into linear "snakes". The + earlier drift-direction "inward" fix was reverted (it worsened snaking). +- **Phase 2 increment 4 (taller mountains + seafloor aging):** mountains used to + cap ~5000 m because collisions were under-weighted and `relax` snapped uplifted + crust back to `continentBase` once a migrating front passed. Now `step()` adds a + real **continent-continent collision** factor (`colliding[]` = continental cell + facing continental, `cfg.collisionFactor`) and an **Andes-class continental arc** + factor (`cfg.arcFactor`), and high continental crust gets **isostatic + persistence**: `relaxEff = relax*(1 - isostaticPersist*clamp((elev-continentBase) + /rootScale,0,1))`, so thick ranges stand and become *erosion-limited* (by the + drift-loop `erode()`) instead of relaxing away. These three boosts are **gated on + the `Planet::drifting` flag — active only in Phase-2 drift, OFF during Phase-1 + forming** (which keeps the original mild factors + full relax). This is + deliberate: Phase-1 forming runs `step()` with no erosion, so if the strong + uplift + weak relax were active there it would never settle (uplift never + balanced) and would rail the clamp — gating to drift, where `erode()` runs every + tick, avoids both. main.cpp sets `planet.drifting=true` when forming settles (and + in `loadGame` from the saved phase), `false` on reseed/regen. **Seafloor aging->depth:** + oceanic crust subsides with `geoAge` via `oceanicBase(age) = + max(oceanBase, ridgeDepth - seafloorSubsidence*sqrt(age))` (half-space cooling); + `oceanBase` is now the **deep abyssal floor** (-6000 m), `ridgeDepth` the shallow + young value (-2500 m), and `seedInitialRelief()` seeds an oceanic age spread + (`seafloorSeedAge`) so the starting seafloor already has ridge->abyss variety. + Headless: ~half of seeds produce >7000 m ranges that persist (the rest are + legitimately low-relief ocean worlds), <2% pinned to the clamp on all seeds, + older seafloor markedly deeper, deterministic. Tune the new knobs in `planet.cfg`. +- **Phase 3 increment 1 (hydrology: rivers, lakes, fluvial erosion):** + `Planet::hydrology(dtMy)` = `routeFlow()` then mass-conserving fluvial erosion, on + the fixed grid (Eulerian, raylib-free). `routeFlow()` does priority-flood + depression-filling (epsilon tilt so flats drain; ocean cells are outlets) → + `sFill`/`sLakeDepth` (a cell with `lakeDepth>0` above sea level is a **lake**), + steepest-descent over the filled surface → `sFlowTo`, and flow accumulation in + descending-fill order → `sDischarge` (rivers = `discharge>riverThreshold`). The + erosion pass walks the network upstream→downstream carrying a sediment load: + stream-power incision `K*Q^m*S^n*dt` where under capacity, deposition where over + (`cap=riverTransport*Q*S`) — filling lakes, building deltas at mouths, depositing + the remainder at ocean sinks so **sum(elevation) is conserved**. Lakes/rivers are + **derived from elevation each tick** (no new saved per-cell field; only erosion + writes back to `elevation`). Orchestration (main.cpp): after `phase3AfterMy` + drift-My the sim **pauses and prompts** ("Continue Phase 2" / "Start Phase 3"); + `H` toggles Phase 3 manually. In Phase 3 the drift loop keeps running + (advect/step/erode) but at a **finer dt** (`cflDtMy()*phase3DtScale`) plus + `hydrology(dt)` — drift never stops, just resolves finer. Lakes shade inland-water + blue (`recolor`); rivers draw as a `centroid→downstream` line network (3D + 2D, + two widths, `J` toggles). Save bumped to **v3** (+ a `phase3` header flag). + Headless: discharge grows downstream and all land rainfall reaches the sinks, + mass conserved to ~1e-15, ~10 lake systems + rivers persist, deterministic. + (Note: the hydrology phase is now framed in the UI as **Phase 2.5**; the + internal `phase3*` names are unchanged.) +- **Phase 3 increment 1 (biomes — classify + color):** `Planet::classifyBiomes()` + (src/sim/PlanetBiomes.cpp, raylib-free) writes a per-cell `Cell.biome` (enum + `Biome`, 13 entries: Ocean, **Ice**, Lake, Beach, Wetland, Grassland, Savanna, + Desert, Forest, Taiga, Tundra, Hills, Mountains). A first **rule-based** pass with + no real climate yet: temperature = warm-equator curve (super-linear in latitude so + cold concentrates at the poles) minus an elevation lapse; moisture = latitudinal + rainfall belts (wet equator/mid-lat, dry subtropics→deserts) + river discharge + + coastal proximity; classified first-match (ice→ocean→lake→beach→mountains→hills→ + lowland-by-temp/moisture). **Polar ice caps** fall out of the temperature test (it + also snow-caps high peaks). All thresholds are **tunable in `planet.cfg`** (the + `biome*` PlanetConfig fields — `biomeIceTemp`, `biomeMountainElev`, the moisture + cutoffs, etc.; only the latitudinal rainfall-belt curve shape stays a fixed helper). + The biome is **saved per cell** — save bumped to **v4** (per-cell biome byte appended + after `invader`; `readState(is, hasBiome)` reads it for v4, reclassifies for v3, so v3 + saves remain loadable). Rendered as color mode `5` (`biomeColor`, src/render/Colors.cpp); + re-run each `refreshView`. `lakeColor` changed to bright turquoise so lakes read clearly + vs ocean (the "blue speckle" near a clicked cell is just the subgrid detail overlay, + where ±250 m value-noise dips below sea level near coasts — not lakes). Headless: + every cell valid, both poles Ice, deep equatorial water Ocean, ≥4 land biomes present, + deterministic, v4 round-trips biome, v3 loads + reclassifies. +- **Phase 3 polish (smaller caps + axial tilt + grid labels):** ice caps trimmed a few + points (~17%→~14% of cells) by lowering `ICE_TEMP` in PlanetBiomes.cpp. Added a + planetary **`axialTilt`** (obliquity, default 23.44°, in PlanetConfig/`planet.cfg`): + the 3D globe + a drawn **spin-axis rod** (through the poles, red/blue pole caps) lean + by it via an `rlRotatef` about world Z wrapping all 3D content in `renderGlobe3D`; + picking un-rotates the world hit dir by −tilt (`rotateZ`, src/render/Picking.cpp) and + 3D plate labels rotate by +tilt so everything stays consistent (the picking sphere is + rotation-invariant). Biomes/2D map are unchanged (tilt is visual + groundwork for + seasons). The graticule (`G`) now shows **lat/lon degree numbers on the 2D map** edges + (`drawGraticuleLabels2D`, plain "60N"/"120W" — the default font has no `°`). `axialTilt` + was added to PlanetConfig (`planet.cfg`). Headless: ice ~14%, biome + `axialTilt` + round-trip, deterministic. +- **Phase 3 polish (biome knobs in config + future-proof save):** the 17 biome + classification thresholds moved from constants into PlanetConfig `biome*` fields + (tunable in `planet.cfg`, `F2`). To stop config additions from breaking saves each + time, the save now stores config as a **self-describing key=value text block** (save + **v6**) parsed like `planet.cfg` (`writeConfigFields`/`parseConfigStream` shared); + doubles written at `precision(17)` round-trip exactly. Adding/removing config fields no + longer breaks saves; v6 just can't load pre-v6 saves (one-time break). Headless: cfg + (incl. non-default `biome*`) round-trips exactly, unknown/missing keys handled. +- **Phase 3 increment 2 (climate model):** `Planet::computeClimate()` + (src/sim/PlanetClimate.cpp, raylib-free, derived/not saved) builds two continuous + per-cell fields. **Temperature** `sTemp` (°C) = the latitude curve (`biome*` temp + params) − elevation lapse. **Precipitation** `sPrecip`: prevailing winds are zonal by + band (tropics/polar easterly, mid-lat westerly); ocean cells are a moisture source and + each land cell takes its **upwind** neighbour's moisture, **rains out** more on windward + upslopes (orographic) and loses a multiplicative fraction per cell (continentality), so + leeward + deep-interior cells dry out. The raw field is near-binary (saturated where the + wind hits the sea, ~0 elsewhere), so it's **diffused** `climateMoistureSmooth` passes to + create wet→dry transition zones, then normalized to `sMoist` (0..1, **median land → + 0.5**, robust to orographic spikes). `classifyBiomes()` now reads `sTemp`/`sMoist` + (dropping the old latitude+discharge+coast hack) → rain-shadow/interior **deserts** + a + varied, per-world biome spread; wetlands now require adjacency to water (ocean/lake). + Color modes `6` (temperature, blue→red) / `7` (precipitation, dry→wet). `computeClimate()` + runs before `classifyBiomes()` in `generate()` and `refreshView()`. New `climate*` config + knobs (planet.cfg). Headless: equator warm/poles cold, lapse, coastal wetter than interior, + deserts present, deterministic. +- **UI polish (full cell info + view label + framing):** the cell-info panel + (`cellInfo`, src/render/Panels.cpp) now shows everything per cell — crust type, **biome** + (`biomeName`), **temperature** + **precipitation %**, and **river/lake** when hydrology is + on — in addition to the existing cell#, lat/lon, elevation, plate, geoAge. The active + color mode is shown **top-center** of the globe ("Biome view", etc., via `colorModeName`), + updating with `1`–`7`. The HUD title/status and the hydrology prompt were reworded to drop + the rigid "Phase N" labels (now "World Creation: forming / drift & erosion / hydrology"); + internal `phase*` names are unchanged. Render/text only — no sim/save/config change. +- Mouse hover (in either view) shows per-cell info. Clicking a tile opens a + right-side detail panel: tile info header + the tile's subgrid drawn as a + flat hoverable grid of subtiles (neighbor-owned subtiles dimmed). A high-res + subgrid patch is also overlaid on the globe for context; the hovered subtile + is marked on the globe. `C` closes the panel. + +## Architecture + +The code is split into a **raylib-free engine** (`src/sim`, testable headless) +and a **raylib viewer** (`src/render`); `src/main.cpp` is a ~10-line entry point. +`Planet` is one class implemented across several `.cpp` files (one per phase/ +concern, all sharing `Planet.hpp`); the viewer is one `Viewer` struct whose state ++ methods are likewise spread across a few render files. CMake adds both folders +to the include path, so includes stay flat (`#include "Planet.hpp"`, `"Viewer.hpp"`). + +``` +src/ + main.cpp entry point: build Viewer, init(argc,argv), run() + sim/ (raylib-free engine -- testable headless) + Vec3.hpp double-precision 3D vector math + IcoSphere.* geodesic icosphere: fixed vertices + neighbor adjacency + Projection.hpp Equal Earth equal-area projection (forward + Newton inverse), + dir<->lon/lat helpers. Header-only, raylib-free, testable. + PlanetTypes.hpp Cell / Plate / SubGrid / PlanetConfig data structures + Planet.hpp the Planet class declaration + config-file func decls + Planet.cpp generation, geometry, plate seeding, shared helpers, subgrid + PlanetTectonics.cpp step() (Phase 1/2 stress->uplift->relax) + PlanetDrift.cpp cflDtMy/advect + plate lifecycle (fission/kick/baby/fuse) + PlanetErosion.cpp erode() + adjustSeaLevel() + PlanetHydrology.cpp routeFlow/computeHydrology/hydrology (Phase 3) + PlanetIO.cpp config file (text) + binary save/load + render/ (raylib viewer) + Colors.* cell color modes (elevation/plate/age/crust/lake) + Map2D.* Equal Earth 2D map: positions + projection/draw helpers + Overlays.* borders, drift arrows, rivers, graticule, segments, subgrids + Picking.* mouse ray / sphere hit / nearest-cell / angle helpers + Panels.* right-column UI: detail panel, hover info, world stats + Viewer.{hpp,cpp} Viewer struct: all state + setup + sim orchestration + ViewerInput.cpp handleInput(): camera, hover picking, click, keys + ViewerRender.cpp renderGlobe3D / renderMap2D / renderPanels / renderHUD / renderPrompt +CMakeLists.txt fetches raylib 5.5 via FetchContent; lists src/sim + src/render +BUILD.md dependencies + build/run + controls +``` + +### Key data structures (src/sim/PlanetTypes.hpp) + +- `Cell` — `unit` (fixed sphere direction), `elevation` (continuous meters, + double, NOT quantized), `plateId`, `geoAge`, `neighbors`, and a + `std::shared_ptr subgrid` HOOK (still null on the cell; the viewer + builds subgrids on demand via `Planet::makeSubGrid(cell,res)`). +- `SubGrid`/`SubCell` — a res*res patch around one macro cell; elevation is an + inverse-distance blend of that cell + neighbors plus value-noise detail, and + each subcell records `nearestMacro`. Phase-4 preview, generated on click. +- `Plate` — `type` (Oceanic/Continental), `driftAxis`, `driftSpeed`. +- `PlanetConfig` — `radius` (default 6.371e6 m, Earth), `subdivisions`, + `seaLevel`, `plateCount`, `seed`. + +### Resolution notes (important, was discussed in design) + +- **Lateral** resolution = cell spacing = `sqrt(4*pi*R^2 / N)`. Coarse on + purpose (~223 km at level 5). Only mountain-range scale, not single hills. +- **Vertical** (elevation) resolution is effectively unlimited: it is a + `double` in meters. 100 m steps or finer are free. Mount-Everest / + Mariana-Trench range is exactly representable. +- Climate (phase 3) reads elevation as a smooth continuous function, so no + resolution is lost by coarse height bands. +- For dense detail later (population/culture), generate a per-cell **subgrid** + on demand and mark which edge sub-cells border which neighbor macro-cell, + so cross-boundary interaction (rain shadow, transition zones) works. + +## Build & run + +```bash +cmake -B build -DCMAKE_BUILD_TYPE=Release +cmake --build build -j +./build/planetsim +``` + +Target OS is Nobara Linux (KDE/Wayland, Intel Arc A770). Dependency install +line is in BUILD.md (`dnf install cmake gcc-c++ mesa-libGL-devel ...`). +raylib 5.5 is fetched automatically — do not vendor it. + +### Quick headless logic test (no display needed) + +```bash +g++ -std=c++17 -O2 -Isrc/sim test_logic.cpp src/sim/IcoSphere.cpp \ + src/sim/Planet.cpp src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp \ + src/sim/PlanetErosion.cpp src/sim/PlanetHydrology.cpp \ + src/sim/PlanetBiomes.cpp src/sim/PlanetClimate.cpp src/sim/PlanetIO.cpp \ + -o /tmp/t && /tmp/t +``` + +Use this to verify tectonics after changing `Planet::step()` without launching +the window (the engine lives in `src/sim` and is raylib-free, so it links without +any render code). The `#pragma omp` lines in `step()` are ignored without +`-fopenmp`, so this serial build is correct; add `-fopenmp` to benchmark the +threaded path. + +### Performance / parallelism + +`Planet::step()` is data-parallel: each pass writes only its own cell index +(double-buffered where it reads a field it writes), so the OpenMP `parallel for` +loops are **bit-identical for any thread count** (determinism intact). It is +memory-bandwidth-bound, so the speedup tops out ~3x (sub-7: ~5.1 -> ~1.5 ms) +around 4-8 threads regardless of core count. `step()` also reuses persistent +scratch buffers (the `s*` members) so it allocates nothing per tick. Small grids +stay serial via `if(n > 20000)`. Control threads with `OMP_NUM_THREADS` (4-8 is the sweet spot; +the default uses all cores for no extra gain). OpenMP is optional and +auto-detected by CMake. For a bigger leap (or 1M+ cells) the next step is a GPU +compute-shader port -- a Phase-2 effort. + +## Controls + +LMB drag orbit · wheel zoom · hover for cell info (3D or map) · +click a tile to open its detail panel (subtiles) · `C` close panel · +drag the 2D map to pan it east/west · `1`..`7` color by +elevation/plate/age/crust-type/biome/temperature/precipitation (active mode shown +top-center of the globe) · +`B` plate borders · `D` drift vectors · `G` lat/lon grid · `J` rivers (Phase 3, +all in 3D + 2D) · `SPACE` or on-screen button pause · `[`/`]` drift speed (My/sec) · +`S` single tick · `F` fast-forward Phase-1 forming to settled · +`H` toggle Phase 3 (hydrology) · `R` reseed · +`+`/`-` subdivision level (1..7) · `F5` save (`planet.save`) · `F9` load · +`F12` screenshot (`screenshot.png`) · `F2` reload `planet.cfg` + regenerate. + +Phase 3: after `phase3AfterMy` simulated years a modal prompt asks **Continue +Phase 2** / **Start Phase 3**; `H` starts/stops it manually. In Phase 3 drift +keeps running at a finer timestep (`cflDtMy()*phase3DtScale`) while rivers, lakes +and fluvial erosion evolve. + +CLI: `--seed N` overrides `cfg.seed`; `--config PATH` uses an alternate config +file (both applied before the initial load/generate). + +Two files in the working dir: `planet.cfg` (human-editable key=value of every +PlanetConfig param, auto-created on first run, reload with `F2`) and +`planet.save` (binary: seed + config + full planet state, written/read by +`Planet::writeState`/`readState`, resumes deterministically). Config is +range-checked by `validateConfig()` on load/`F2`; an invalid file reverts to safe +defaults (without overwriting your `planet.cfg`) and shows a status message. The +save header is versioned (currently **6**; v2 adds the `[`/`]` drift rate, v3 a +`phase3` flag, v4 a per-cell biome byte, v6 stores config as a **self-describing +key=value text block** instead of a raw POD dump); newer-than-supported is rejected. +**As of v6, adding/removing PlanetConfig fields no longer breaks saves** — the saved +config is parsed like `planet.cfg` (unknown keys ignored, missing keys keep defaults), +written at `precision(17)` so doubles round-trip exactly. (v6 cannot load pre-v6 saves — +a one-time break; a length guard makes that fail gracefully.) `writeConfigFields`/ +`parseConfigStream` in PlanetIO.cpp are shared by `planet.cfg` and the save. + +Layout (1920x1080): left column 70% wide = 3D globe (top, 60% h, RenderTexture +1344x648) + 2D Equal Earth map (bottom, 40% h); right column 30% wide = cell info +(top 50% h) + subareas (bottom 50% h). 3D hover uses a custom camera ray with the +3D viewport (1344x648 at the origin -- GetScreenToWorldRay assumes the full +screen, wrong here); 2D hover uses EqualEarth::inverse (minus the `mapLon` pan). +Borders (`B`), drift vectors (`D`) and a lat/lon graticule (`G`) draw in BOTH the +3D globe and the 2D map; drag the 2D map left/right to pan longitude (`mapLon`). +The graticule (`G`) also draws **lat/lon degree numbers along the 2D map edges**. The +3D globe (and its spin-axis rod) lean by `cfg.axialTilt` (an `rlRotatef` about world Z +around all 3D content in `renderGlobe3D`); 3D picking un-rotates the world hit dir by +−tilt and plate labels rotate by +tilt (`rotateZ`, src/render/Picking.cpp) to stay in sync. +Borders draw in two colors: real plate boundaries yellow, young spreading-ridge +(baby-plate) boundaries red — both via `buildBorders` filling two segment lists. + +Drift arrows: one cyan arrow per plate at its centroid along the local drift +velocity `omega x r`; length scales with speed. Rebuilt per world-gen. Each plate +is also tagged with a `P` label drawn just above its centroid — projected by +hand in 3D (the project's own `Vec3` camera-basis math, matching `BeginMode3D`'s +fovy/aspect, behind-camera points skipped) and via Equal Earth in the 2D map. + +**Sim model (`Viewer::stepSim`, src/render/Viewer.cpp):** single-threaded. Each +frame (while not paused/settled) +it advances ~`formRate` (55) ticks/second of `planet.step()` (Phase 1 = pure +tectonic forming, no erosion — `planet.drifting` is false) and rebuilds the view +live, so you watch the terrain rise. `step()` returns the max per-tick +elevation change; after `settleNeed` (3) consecutive ticks below `settleThresh` +(2 m) it sets `settled` and stops stepping — no idle CPU. Re-evolve / reseed / +subdivision-change clear `settled` and restart the forming pass. Rendering reads +the live `cells` directly (safe: nothing else mutates them). Phase 2's +*continuous* sim will want a worker thread + a render snapshot (so the 60 fps +render never races the sim) — that's the natural place to reintroduce threading. +Plate borders are traced once per world-gen as a dual contour through boundary +triangles (plates are fixed in phase 1). + +## Tuning knobs (expect to adjust these first) + +- `elevExagg` (src/render/Viewer.hpp) — visual elevation exaggeration; without it the + sphere looks smooth. +- `axialTilt` (PlanetConfig, `planet.cfg`) — obliquity in degrees (default 23.44). Leans + the 3D globe + spin-axis rod; groundwork for seasons. Editable + `F2` to apply. +- Biome thresholds (`biome*` in PlanetConfig / `planet.cfg`) — `biomeIceTemp` (lower = + smaller polar/snow caps), `biomeMountainElev`/`biomeHillsElev`, the moisture cutoffs + (`biomeDesertMoist`/`biomeGrassMoist`/`biomeWetlandMoist`), and the temperature model + (`biomeEquatorTemp`/`biomePoleDrop`/`biomeLatExp`/`biomeElevLapse`, shared with climate). + Editable + `F2`. +- Climate (`climate*` in PlanetConfig / `planet.cfg`) — precipitation model: + `climateOrographic` (windward-rain strength), `climateRainEfficiency`, + `climateContinentality` (inland drying), `climateMoistureSmooth` (diffusion passes → + wet/dry transition zones; raise for smoother, more grassland/forest), `climateOceanMoisture`, + `climateOroRefHeight`, `climateWindPasses`. Temperature uses the `biome*` temp params. +- `upliftGain` (PlanetConfig) — m/tick per unit convergence stress; main + knob for how fast/high relief builds. +- `relax` (PlanetConfig) — isostatic relaxation toward base elevation. Peaks + asymptote at `base + perTickUplift/relax`, so raising it lowers/flattens + the equilibrium; lowering it makes relief taller and slower to settle. +- `beltWidth` (PlanetConfig) — how many cell-rings a mountain belt spreads + inland (belt width / flank extent). +- `trenchFactor` (Planet::step) — depth multiplier for subduction trenches. +- `driftSpeed` (Planet::assignPlates) — how fast boundary stress builds. +- **Plate lifecycle (PlanetConfig, Phase 2 inc. 1.5):** `splitFraction` (0.20, + plate share that may rift), `splitProbBase`/`splitProbSlope` (fission prob + ramp), `splitCheckEvery` (10, iters between periodic checks), `stalemateEps` / + `stalemateBoost` (deadlock detection + kick), `stalemateWindows` (4, consecutive + stuck windows before a kick — raise to calm border jitter), `babyMinCells` (4, + baby blobs smaller than this dissolve as noise), `babyPromoteFrac` (0.7%, + ridge-strip size to become a real plate), `volcanicLandFrac` / `volcanicElev` + (land grown on promotion), `landBand` (0.10, soft land-conservation band), + `miniPlateCells` (50, below = "mini") + `fuseMinPlates` (3, cluster size to fuse + + steal). Tune these to control how lively / fragmented the plate map stays. +- **Erosion + sea level (PlanetConfig, Phase 2 inc. 2):** `erosionLandRate` (0.08) + / `erosionSeaRate` (0.02) — fraction/My worn off above / below sea level (raise + for faster, flatter terrain); `landFractionTarget` (0.30) — geographic land + goal; `seaLevelStep` (100 m) — fixed nudge per adjustment + `seaLevelTol` (0.02) + — deadband where sea level rests (raise step or shrink tol for a tighter 30/70, + but a big step can overshoot a flat "cliff"); `seaLevelEvery` (100) — iterations + between sea-level updates (higher = more gradual). All editable in `planet.cfg`. +- **Orogeny — taller mountains (PlanetConfig, Phase 2 inc. 4):** `collisionFactor` + (1.8, continent-continent uplift, Himalaya) and `arcFactor` (1.4, continental + subduction-arc uplift, Andes) — raise either to make ranges taller/more reliable + across seeds; `isostaticPersist` (0.85, how strongly high crust resists relax — + toward 1 = ranges barely erode and continents stay high; lower = more dynamic + rise/erode) saturating at `rootScale` (2500 m above `continentBase`). These three + boosts are **drift-only** (gated on `Planet::drifting`); in Phase-2 drift the + per-tick `erode()` is what limits their height, so they don't rail the clamp. The + hard clamp `[-11000, 9000]` in `step()` is the Everest-class cap; a few peak + cells may sit there during drift (<2% — fine). +- **Seafloor aging->depth (PlanetConfig, Phase 2 inc. 4):** `oceanBase` is now the + **deep abyssal floor** (-6000 m, not a flat ocean base) and `ridgeDepth` the + shallow young value (-2500 m); `seafloorSubsidence` (280 m per sqrt(My)) sets how + fast oceanic crust deepens with `geoAge` (half-space cooling), and + `seafloorSeedAge` (80 My) is the initial oceanic age spread at generation so the + starting seafloor already has ridge->abyss variety. +- **Hydrology (PlanetConfig, Phase 3 inc. 1):** `phase3AfterMy` (300) — drift-My + before the Phase-3 prompt; `phase3DtScale` (0.2) — Phase-3 timestep = + `cflDtMy()*this` (smaller = finer carving, slower drift per step); `rainfall` + (1.0) — uniform precip per cell (drainage-area unit; orographic precip is a later + climate add-on); `riverThreshold` (50) — discharge above which a cell is a river + (also the river-render threshold); `riverIncision` K (0.02) + `riverDischargeExp` + m (0.5) + `riverSlopeExp` n (1.0) — stream-power incision `K*Q^m*S^n*dt` (raise K + for faster valley carving); `riverTransport` (0.1) — transport capacity + `cap=this*Q*S` and `depFrac` (0.25) — deposition rate of excess load (raise both + for more deltas / faster lake infill). All editable in `planet.cfg`. + +## Conventions + +- All code, identifiers, comments and filenames in **English**. +- Before generating code, summarize the approach and ask whether to proceed. +- Prefer **inline code in chat** over attachments; give exact file contents + rather than long explanations. +- Direct, concrete answers. Metric system throughout. +- C++17, multi-file CMake. Keep simulation logic free of raylib so it stays + testable headless: engine in `src/sim` (no raylib), all rendering in `src/render`. +- Geometry stays fixed — never make cells move; add new per-cell properties + and flow them over the existing grid + neighbor adjacency. diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f5dc2e3 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,54 @@ +cmake_minimum_required(VERSION 3.16) +project(PlanetSim CXX) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +include(FetchContent) +FetchContent_Declare( + raylib + GIT_REPOSITORY https://github.com/raysan5/raylib.git + GIT_TAG 5.5 +) +# Build raylib as a static lib, no examples. +set(BUILD_EXAMPLES OFF CACHE BOOL "" FORCE) +set(BUILD_GAMES OFF CACHE BOOL "" FORCE) +FetchContent_MakeAvailable(raylib) + +add_executable(planetsim + src/main.cpp + # Engine (raylib-free, testable headless) -- src/sim + src/sim/IcoSphere.cpp + src/sim/Planet.cpp + src/sim/PlanetTectonics.cpp + src/sim/PlanetDrift.cpp + src/sim/PlanetErosion.cpp + src/sim/PlanetHydrology.cpp + src/sim/PlanetBiomes.cpp + src/sim/PlanetClimate.cpp + src/sim/PlanetIO.cpp + # Viewer (raylib) -- src/render + src/render/Colors.cpp + src/render/Map2D.cpp + src/render/Overlays.cpp + src/render/Picking.cpp + src/render/Panels.cpp + src/render/Viewer.cpp + src/render/ViewerInput.cpp + src/render/ViewerRender.cpp +) +# Flat includes ("Planet.hpp", "Viewer.hpp", ...) resolve across both folders. +target_include_directories(planetsim PRIVATE src/sim src/render) +target_link_libraries(planetsim PRIVATE raylib) + +# OpenMP parallelizes the per-cell passes in Planet::step(). Optional: without +# it the #pragma omp lines are ignored and the sim runs (correctly) serial. +find_package(OpenMP) +if(OpenMP_CXX_FOUND) + target_link_libraries(planetsim PRIVATE OpenMP::OpenMP_CXX) +endif() + +# Linux system libs raylib needs at link time. +if(UNIX AND NOT APPLE) + target_link_libraries(planetsim PRIVATE m pthread dl) +endif() diff --git a/docs/design-notes.md b/docs/design-notes.md new file mode 100644 index 0000000..f39eeea --- /dev/null +++ b/docs/design-notes.md @@ -0,0 +1,94 @@ +# Design notes (durable context) + +These are the non-obvious decisions/conventions that were previously only in Claude's +auto-memory (which lives under `~/.claude/` and does **not** travel with the repo). Captured +here so the context survives a move to another machine/server. `CLAUDE.md` has the +authoritative current-state changelog; this is the "why / where things live" summary. + +## Framing: World Creation → Live World + +The roadmap is no longer rigid numbered "phases". **World Creation** is a set of +continuous, overlapping stages on a geological clock (My): tectonics → continental drift & +erosion → hydrology → climate → biomes → (fauna & flora, next). The long-term goal is a +separate **Live World** mode that runs the *finished* planet at a much slower real-time +clock (hours/days/weeks/months) with dynamic weather (clouds, rain, storms) and living +ecosystems/civilization. **Internal code still uses `phase*` names** (`Planet::drifting`, +the `phase3` flag, `phase3AfterMy`/`phase3DtScale` config keys) for save/config +compatibility — only display strings and docs use the new framing. + +## Code module layout + +Split into a raylib-free **engine** (`src/sim/`, headless-testable) and a raylib **viewer** +(`src/render/`); `src/main.cpp` is a ~10-line entry point. CMake adds both dirs to the +include path, so includes stay flat (`#include "Planet.hpp"`, `"Viewer.hpp"`). + +`Planet` is **one class implemented across several .cpp files** (all share `Planet.hpp`): +- `PlanetTypes.hpp` — `Cell` / `Plate` / `SubGrid` / `Biome` enum / `PlanetConfig`. +- `Planet.cpp` — generation, geometry, plate seeding, RNG + shared helpers, subgrid, min/max. +- `PlanetTectonics.cpp` — `step()` (stress→uplift→relax; orogeny boosts gated on `drifting`). +- `PlanetDrift.cpp` — `cflDtMy`/`advect` + plate lifecycle (fission/kick/baby/fuse/enclosed). +- `PlanetErosion.cpp` — `erode` + `adjustSeaLevel`. +- `PlanetHydrology.cpp` — `routeFlow`/`computeHydrology`/`hydrology` (depression-fill→lakes, + steepest-descent→rivers, mass-conserving stream-power incision). +- `PlanetClimate.cpp` — `computeClimate()` (temperature + orographic precipitation). +- `PlanetBiomes.cpp` — `classifyBiomes()` (per-cell `Cell.biome` from elevation + climate). +- `PlanetIO.cpp` — text config + binary save/load. +- (`PlanetBiosphere.cpp` — fauna/flora, planned; see `fauna-flora-plan.md`.) + +The viewer is one `Viewer` struct: `Viewer.{hpp,cpp}` (state + setup + sim orchestration), +`ViewerInput.cpp` (camera/picking/keys), `ViewerRender.cpp` (globe/map/panels/HUD/prompt), +plus topical helpers `Colors` / `Map2D` / `Overlays` / `Picking` / `Panels`. + +Per-tick order in `Viewer::refreshView()`: `computeHydrology()` (if hydrology on) → +`computeClimate()` → `classifyBiomes()` → (`computeBiosphere()` when added) → `recolor()`. + +## Core principle (do not violate) + +Geometry is **fixed** — cells (icosphere vertices) never move. Only per-cell *properties* +flow over the fixed grid + neighbor adjacency (Eulerian). New phenomena = new per-cell +fields flowed over the grid, never moving cells. + +## Axial tilt render convention (non-obvious) + +The 3D globe is rendered leaned by `cfg.axialTilt` via `rlRotatef(tilt,0,0,1)` wrapping all +3D content in `renderGlobe3D`. Because that rotation isn't in the data, anything mapping +between world and model space must compensate with `rotateZ(v, ±tilt)` (src/render/ +Picking.cpp): 3D picking un-rotates the ray hit by `−tilt` before `nearestCell`; 3D plate +labels rotate by `+tilt` before projecting. The 2D map + biome/climate are tilt-independent. + +## Save format (v6) — self-describing config + +`planet.save` stores `PlanetConfig` as a **self-describing key=value text block** (not a raw +POD dump), parsed like `planet.cfg` (`writeConfigFields`/`parseConfigStream` shared in +PlanetIO.cpp), written at `precision(17)` so doubles round-trip exactly. Consequence: +**adding/removing PlanetConfig fields no longer breaks saves** (unknown keys ignored, missing +keys keep defaults). v6 cannot load pre-v6 saves (one-time break; a length guard fails it +gracefully). Per-cell `Cell.biome` is saved (a byte appended after `invader`). + +## Climate + biome model (derived, not saved) + +`computeClimate()` builds two derived per-cell fields: +- **Temperature** (°C) = latitude curve (`biomeEquatorTemp/PoleDrop/LatExp`, super-linear so + cold concentrates at poles) − `biomeElevLapse` × elevation. +- **Precipitation**: zonal prevailing winds (easterly tropics/poles, westerly mid-lat); ocean + cells are a moisture source; each land cell takes its **upwind** neighbour's moisture, rains + out more on windward upslopes (orographic), loses a multiplicative fraction per cell + (continentality) → leeward/interior drying. The raw field is near-binary, so it's + **diffused** (`climateMoistureSmooth` passes) into transition zones, then normalized to + `sMoist∈[0,1]` by anchoring the **median land precip → 0.5** (robust to orographic spikes). + +`classifyBiomes()` reads `sTemp` + `sMoist` (not a latitude hack) → rain-shadow/interior +deserts emerge; 13 biomes incl. polar Ice; wetlands require water adjacency. All biome & +climate thresholds are tunable `biome*` / `climate*` keys in `planet.cfg`. + +## Headless testing + +Engine is raylib-free, so logic is tested without a display. Build/run: +``` +g++ -std=c++17 -O2 -Isrc/sim test_logic.cpp src/sim/IcoSphere.cpp src/sim/Planet.cpp \ + src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp src/sim/PlanetErosion.cpp \ + src/sim/PlanetHydrology.cpp src/sim/PlanetBiomes.cpp src/sim/PlanetClimate.cpp \ + src/sim/PlanetIO.cpp -o /tmp/t && /tmp/t +``` +(add new `src/sim/*.cpp` to that list as stages are added). `Planet::step()` passes are +data-parallel + double-buffered → bit-identical for any OpenMP thread count (determinism). diff --git a/docs/fauna-flora-plan.md b/docs/fauna-flora-plan.md new file mode 100644 index 0000000..ffe2bb6 --- /dev/null +++ b/docs/fauna-flora-plan.md @@ -0,0 +1,82 @@ +# Next step — Fauna & Flora (derived carrying-capacity densities) + +> Status: **planned, not yet implemented.** This is the next World-Creation stage after +> climate. Approved design captured here so it can be picked up verbatim (e.g. after the +> project moves to another machine — the auto-memory does not travel, this doc does). + +## Decision summary + +- **Derived carrying-capacity map**, not a live simulation: vegetation + wildlife + **density scalars** recomputed each tick from the climate fields + terrain, exactly like + biomes/climate. It answers *"what could live here"*, not *"what is living/evolving here"*. +- The actual **time-evolving ecosystem** (population growth, migration, predator-prey, + extinction, seasons, species) is reserved for the future **Live World** real-time mode + (year/season timescales — not the My geological clock). +- Fauna uses a **herbivore + carnivore** food-chain split. +- Fields are **derived / not saved** (deterministic functions of climate) → **no + save-format change**. New config knobs ride the v6 self-describing save with no break. + +## Part A — Biosphere fields (`src/sim/`, raylib-free) + +- **`Planet.hpp`**: declare `void computeBiosphere();` + accessors + `const std::vector& vegetation()/herbivores()/carnivores() const`. Add scratch + members `std::vector sVeg, sHerb, sCarn;` (not saved). +- **New `src/sim/PlanetBiosphere.cpp`** — `Planet::computeBiosphere()` (run *after* + `computeClimate()`; reads `sTemp` + `sMoist`): + - **Vegetation (flora)** `sVeg[i] ∈ [0,1]`: `0` on water (`elev ≤ seaLevel`). On land, an + NPP-style **Liebig minimum** of a temperature factor and a moisture factor: + - `tF = clamp((sTemp[i] − bioVegTempMin)/(bioVegTempOpt − bioVegTempMin), 0, 1)` + - `mF = clamp(sMoist[i]/bioVegMoistRef, 0, 1)` + - `sVeg = min(tF, mF)` → lush warm-wet, ~0 in ice/desert/alpine (cold- or water-limited), + matching the biomes. + - **Herbivores** `sHerb[i] = clamp(sVeg[i] · bioHerbProductivity, 0, 1)` (capacity scales + with plant productivity). + - **Carnivores** `sCarn[i]`: present only where prey is abundant, ramping with it: + `sCarn = sHerb > bioCarnPreyMin ? clamp((sHerb − bioCarnPreyMin)/(1 − bioCarnPreyMin) + · bioCarnScale, 0, 1) : 0` → concentrated in the richest regions (an energy-pyramid feel, + visually distinct from the broader herbivore pattern). +- **`src/sim/Planet.cpp` `generate()`**: call `computeBiosphere()` after `classifyBiomes()`. + +## Part B — Config knobs (`PlanetTypes.hpp` + `PlanetIO.cpp`) + +Add to `PlanetConfig` + `CONFIG_FIELDS` + `validateConfig`: + +| key | default | meaning | +|---|---|---| +| `bioVegTempMin` | −5 °C | below this, no plant growth | +| `bioVegTempOpt` | 15 °C | at/above this, temperature isn't limiting | +| `bioVegMoistRef` | 0.5 | normalized moisture at which water isn't limiting | +| `bioHerbProductivity` | 0.9 | herbivore capacity per unit vegetation | +| `bioCarnPreyMin` | 0.30 | min herbivore density to support carnivores | +| `bioCarnScale` | 1.0 | carnivore density ramp above the prey threshold | + +## Part C — Rendering (`src/render/`) + +- **`Colors.{hpp,cpp}`**: `ColorMode` += `Vegetation`, `Herbivores`, `Carnivores`; add + `vegColor` (barren tan → lush green), `herbColor` (pale → amber/orange), `carnColor` + (pale → red/violet); extend `colorModeName`. +- **`Viewer.cpp`**: `recolor()` adds the three cases (read the accessors, empty-guarded); + `refreshView()` calls `planet.computeBiosphere()` after `classifyBiomes()`. +- **`ViewerInput.cpp`**: `KEY_EIGHT` → Vegetation, `KEY_NINE` → Herbivores, `KEY_ZERO` → + Carnivores (each `recolor()`). The top-center view-name label already shows the mode. +- **`ViewerRender.cpp` `renderHUD`**: controls line gains `8 veg 9 herbiv 0 carniv`. +- **`Panels.cpp` `cellInfo()`**: add a life line — `vegetation X% herbivores Y% + carnivores Z%` (guarded on the arrays being sized). + +## Verification (headless — extend the `src/sim` test set with `PlanetBiosphere.cpp`) + +- **Vegetation:** 0 on ocean & Ice cells; near-0 on Desert; high (>0.6) on warm-wet Forest; + all in [0,1]; finite; deterministic. +- **Fauna food chain:** `sHerb` correlates with `sVeg`; `sCarn[i] > 0 ⇒ sHerb[i] > + bioCarnPreyMin` and carnivore-positive cells are a strict subset of herbivore-rich cells; + all in [0,1]; deterministic. +- `test_logic`, biomes+config, climate suites still pass; app builds clean; no save change. +- **GUI:** `8` lush green continents thinning to barren in deserts/poles/peaks; `9` grazer + density following productivity; `0` predators concentrated in the richest belts; cell info + shows vegetation/herbivore/carnivore %. + +## Out of scope (future / Live World) + +Time-evolving populations (growth, migration, predator-prey, extinction), seasons, and +species/typed organisms — all part of the future **Live World** real-time simulation that +runs the finished planet at hours/days/weeks/months with dynamic weather + life. diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..63cc9b8 --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,11 @@ +#include "Viewer.hpp" + +// Program entry point: build the viewer, initialize it (window, layout, config, +// first world; honours --seed / --config), and run the frame loop. All the work +// lives in the engine (src/sim) and the viewer (src/render). +int main(int argc, char** argv) { + Viewer v; + if (!v.init(argc, argv)) return 1; + v.run(); + return 0; +} diff --git a/src/render/Colors.cpp b/src/render/Colors.cpp new file mode 100644 index 0000000..907aebf --- /dev/null +++ b/src/render/Colors.cpp @@ -0,0 +1,123 @@ +#include "Colors.hpp" +#include +#include + +Color elevationColor(double e, double seaLevel) { + if (e < seaLevel) { + // Water: deep -> shallow blue. + double t = std::clamp((e + 11000.0) / (seaLevel + 11000.0), 0.0, 1.0); + return Color{ (unsigned char)(10 + 20 * t), + (unsigned char)(30 + 90 * t), + (unsigned char)(80 + 120 * t), 255 }; + } + // Land: green -> brown -> white by height. + double t = std::clamp(e / 9000.0, 0.0, 1.0); + if (t < 0.4) { double u = t / 0.4; + return Color{ (unsigned char)(60 + 80 * u), (unsigned char)(140 - 30 * u), + (unsigned char)(50), 255 }; } + if (t < 0.75) { double u = (t - 0.4) / 0.35; + return Color{ (unsigned char)(140 - 30 * u), (unsigned char)(110 - 40 * u), + (unsigned char)(50 + 10 * u), 255 }; } + double u = (t - 0.75) / 0.25; + return Color{ (unsigned char)(110 + 145 * u), (unsigned char)(70 + 185 * u), + (unsigned char)(60 + 195 * u), 255 }; +} + +Color plateColor(int id) { + float h = std::fmod(id * 0.61803398875f, 1.0f) * 360.0f; + return ColorFromHSV(h, 0.65f, 0.85f); +} + +Color ageColor(double age, double maxAge) { + double t = std::clamp(age / std::max(1.0, maxAge), 0.0, 1.0); + return Color{ (unsigned char)(40 + 200 * t), (unsigned char)(40), + (unsigned char)(200 - 160 * t), 255 }; +} + +Color crustColor(bool oceanic) { + return oceanic ? Color{ 30, 60, 160, 255 } // oceanic: deep blue + : Color{ 160, 130, 70, 255 }; // continental: warm brown +} + +Color lakeColor() { return Color{ 40, 200, 210, 255 }; } // bright turquoise (vs ocean blue) + +Color biomeColor(Biome b) { + switch (b) { + case Biome::Ocean: return Color{ 20, 60, 120, 255 }; // deep blue + case Biome::Ice: return Color{ 235, 240, 250, 255 }; // white (polar caps / snow) + case Biome::Lake: return Color{ 40, 200, 210, 255 }; // bright turquoise + case Biome::Beach: return Color{ 222, 210, 150, 255 }; // pale sand + case Biome::Wetland: return Color{ 70, 115, 95, 255 }; // dark teal-green (swamp/bayou) + case Biome::Grassland: return Color{ 130, 185, 80, 255 }; // light green + case Biome::Savanna: return Color{ 185, 180, 85, 255 }; // yellow-green + case Biome::Desert: return Color{ 214, 184, 120, 255 }; // tan + case Biome::Forest: return Color{ 40, 110, 50, 255 }; // dark green + case Biome::Taiga: return Color{ 55, 105, 85, 255 }; // blue-green (boreal) + case Biome::Tundra: return Color{ 155, 165, 150, 255 }; // pale grey-green + case Biome::Hills: return Color{ 120, 135, 70, 255 }; // olive + case Biome::Mountains: return Color{ 135, 125, 115, 255 }; // grey-brown + } + return Color{ 255, 0, 255, 255 }; // unreachable; flags an unmapped biome +} + +const char* biomeName(Biome b) { + switch (b) { + case Biome::Ocean: return "Ocean"; + case Biome::Ice: return "Ice cap"; + case Biome::Lake: return "Lake"; + case Biome::Beach: return "Beach"; + case Biome::Wetland: return "Wetland"; + case Biome::Grassland: return "Grassland"; + case Biome::Savanna: return "Savanna"; + case Biome::Desert: return "Desert"; + case Biome::Forest: return "Forest"; + case Biome::Taiga: return "Taiga"; + case Biome::Tundra: return "Tundra"; + case Biome::Hills: return "Hills"; + case Biome::Mountains: return "Mountains"; + } + return "?"; +} + +const char* colorModeName(ColorMode m) { + switch (m) { + case ColorMode::Elevation: return "Elevation"; + case ColorMode::Plate: return "Plates"; + case ColorMode::Age: return "Crust age"; + case ColorMode::Crust: return "Crust type"; + case ColorMode::Biome: return "Biome"; + case ColorMode::Temperature: return "Temperature"; + case ColorMode::Precip: return "Precipitation"; + } + return "?"; +} + +// Temperature ramp over ~[-40, 40] C: deep blue -> cyan -> green -> yellow -> red. +Color tempColor(double celsius) { + double t = std::clamp((celsius + 40.0) / 80.0, 0.0, 1.0); // 0 cold .. 1 hot + // 4 segments between 5 control colors. + static const unsigned char key[5][3] = { + { 30, 40, 130 }, // -40 C deep blue + { 60, 160, 210 }, // -20 C cyan + { 90, 190, 90 }, // 0 C green + { 225, 200, 70 }, // +20 C yellow + { 210, 70, 50 }, // +40 C red + }; + double s = t * 4.0; int k = std::min(3, (int)s); double f = s - k; + auto L = [&](int c){ return (unsigned char)(key[k][c] + (key[k + 1][c] - key[k][c]) * f); }; + return Color{ L(0), L(1), L(2), 255 }; +} + +// Precipitation ramp over normalized [0,1]: tan (dry) -> green -> teal/blue (wet). +Color precipColor(double moist01) { + double t = std::clamp(moist01, 0.0, 1.0); + static const unsigned char key[4][3] = { + { 205, 180, 120 }, // 0.00 dry tan + { 170, 185, 90 }, // 0.33 scrub + { 70, 160, 90 }, // 0.66 green + { 40, 120, 190 }, // 1.00 wet blue + }; + double s = t * 3.0; int k = std::min(2, (int)s); double f = s - k; + auto L = [&](int c){ return (unsigned char)(key[k][c] + (key[k + 1][c] - key[k][c]) * f); }; + return Color{ L(0), L(1), L(2), 255 }; +} diff --git a/src/render/Colors.hpp b/src/render/Colors.hpp new file mode 100644 index 0000000..0f942aa --- /dev/null +++ b/src/render/Colors.hpp @@ -0,0 +1,23 @@ +#pragma once +#include "raylib.h" +#include "PlanetTypes.hpp" // Biome + +// Cell color mapping for the viewer. Pure functions of cell properties. + +enum class ColorMode { Elevation, Plate, Age, Crust, Biome, Temperature, Precip }; + +Color elevationColor(double e, double seaLevel); +Color plateColor(int id); +Color ageColor(double age, double maxAge); +Color crustColor(bool oceanic); +// Phase-2.5 inland water (lakes): a bright turquoise, clearly distinct from ocean. +Color lakeColor(); +// Phase-3 biome palette + display name. +Color biomeColor(Biome b); +const char* biomeName(Biome b); +// Human-readable name of a color mode (for the active-view label). +const char* colorModeName(ColorMode m); +// Phase-3 climate: temperature in deg C (blue cold -> red hot); precip normalized 0..1 +// (tan dry -> green -> blue wet). +Color tempColor(double celsius); +Color precipColor(double moist01); diff --git a/src/render/Map2D.cpp b/src/render/Map2D.cpp new file mode 100644 index 0000000..1e8dd80 --- /dev/null +++ b/src/render/Map2D.cpp @@ -0,0 +1,69 @@ +#include "Map2D.hpp" +#include "rlgl.h" +#include "Projection.hpp" +#include +#include + +void buildMap2D(const Planet& p, Rectangle r, Map2D& m) { + double hw = EqualEarth::halfWidth(), hh = EqualEarth::halfHeight(); + size_t n = p.cells.size(); + m.pos.resize(n); m.lon.resize(n); m.lat.resize(n); + for (size_t i = 0; i < n; ++i) { + double lon, lat; dirToLonLat(p.cells[i].unit, lon, lat); + double x, y; EqualEarth::forward(lon, lat, x, y); + m.pos[i] = Vector2{ (float)(r.x + (x / hw * 0.5 + 0.5) * r.width), + (float)(r.y + (0.5 - y / hh * 0.5) * r.height) }; + m.lon[i] = lon; m.lat[i] = lat; + } +} + +double wrapPi(double l) { while (l > M_PI) l -= 2*M_PI; while (l < -M_PI) l += 2*M_PI; return l; } + +Vector2 projLonLat(double lon, double lat, double lonOffset, Rectangle r) { + double hw = EqualEarth::halfWidth(), hh = EqualEarth::halfHeight(); + double x, y; EqualEarth::forward(wrapPi(lon + lonOffset), lat, x, y); + return Vector2{ (float)(r.x + (x / hw * 0.5 + 0.5) * r.width), + (float)(r.y + (0.5 - y / hh * 0.5) * r.height) }; +} +Vector2 mapScreen(const Map2D& m, int idx, Rectangle r, double lonOffset) { + return projLonLat(m.lon[idx], m.lat[idx], lonOffset, r); +} + +void drawMap2D(const Planet& p, const std::vector& vc, + const Map2D& m, Rectangle r, double lonOffset) { + double hw = EqualEarth::halfWidth(); + auto px = [&](double lon, double lat) -> float { + double x, y; EqualEarth::forward(lon, lat, x, y); + return (float)(r.x + (x / hw * 0.5 + 0.5) * r.width); + }; + const std::vector& tri = p.triIndices(); + rlDisableBackfaceCulling(); + rlBegin(RL_TRIANGLES); + for (size_t k = 0; k + 2 < tri.size(); k += 3) { + int v[3] = { tri[k], tri[k + 1], tri[k + 2] }; + double lo[3] = { wrapPi(m.lon[v[0]] + lonOffset), wrapPi(m.lon[v[1]] + lonOffset), + wrapPi(m.lon[v[2]] + lonOffset) }; + double mn = std::min({lo[0], lo[1], lo[2]}); + double mx = std::max({lo[0], lo[1], lo[2]}); + if (mx - mn <= M_PI) { // fast path (no wrap) + for (int t = 0; t < 3; ++t) { + rlColor4ub(vc[v[t]].r, vc[v[t]].g, vc[v[t]].b, 255); + rlVertex2f(px(lo[t], m.lat[v[t]]), m.pos[v[t]].y); + } + } else { // antimeridian seam + double ul[3] = { lo[0], lo[1], lo[2] }; // unwrap around v0 + double ref = ul[0]; + for (int t = 0; t < 3; ++t) { + while (ul[t] - ref > M_PI) ul[t] -= 2 * M_PI; + while (ref - ul[t] > M_PI) ul[t] += 2 * M_PI; + } + const double shift[3] = { 0.0, 2 * M_PI, -2 * M_PI }; // both edges; scissor clips + for (double sh : shift) + for (int t = 0; t < 3; ++t) { + rlColor4ub(vc[v[t]].r, vc[v[t]].g, vc[v[t]].b, 255); + rlVertex2f(px(ul[t] + sh, m.lat[v[t]]), m.pos[v[t]].y); + } + } + } + rlEnd(); +} diff --git a/src/render/Map2D.hpp b/src/render/Map2D.hpp new file mode 100644 index 0000000..5ab508d --- /dev/null +++ b/src/render/Map2D.hpp @@ -0,0 +1,23 @@ +#pragma once +#include "raylib.h" +#include "Planet.hpp" +#include + +// Equal Earth 2D map: per-cell screen positions + the projection/draw helpers. + +struct Map2D { + std::vector pos; // screen position per cell (fast path) + std::vector lon; // longitude per cell (radians, for seam detection) + std::vector lat; // latitude per cell (radians) +}; + +void buildMap2D(const Planet& p, Rectangle r, Map2D& m); + +double wrapPi(double l); +// Project a (lon,lat) onto the map with the longitude pan applied. +Vector2 projLonLat(double lon, double lat, double lonOffset, Rectangle r); +Vector2 mapScreen(const Map2D& m, int idx, Rectangle r, double lonOffset); + +// lonOffset pans the map east/west (radians); y is unchanged by the pan. +void drawMap2D(const Planet& p, const std::vector& vc, + const Map2D& m, Rectangle r, double lonOffset); diff --git a/src/render/Overlays.cpp b/src/render/Overlays.cpp new file mode 100644 index 0000000..9e1be25 --- /dev/null +++ b/src/render/Overlays.cpp @@ -0,0 +1,197 @@ +#include "Overlays.hpp" +#include "Map2D.hpp" // projLonLat (2D projection of overlays) +#include "Colors.hpp" // elevationColor (subgrid patch) +#include "rlgl.h" +#include "Projection.hpp" // dirToLonLat / lonLatToDir +#include +#include +#include // std::abs(int) + +void buildBorders(const Planet& p, float radius, + std::vector& real, std::vector& ridge) { + real.clear(); ridge.clear(); + auto midV = [&](int i, int j) -> Vector3 { + Vec3 m = ((p.cells[i].unit + p.cells[j].unit) * 0.5).normalized() * radius; + return Vector3{ (float)m.x, (float)m.y, (float)m.z }; + }; + auto isBaby = [&](int pid){ return pid >= 0 && pid < (int)p.plates.size() && p.plates[pid].baby; }; + auto emit = [&](const Vector3& a, const Vector3& b, int pA, int pB) { + std::vector& out = (isBaby(pA) || isBaby(pB)) ? ridge : real; + out.push_back(a); out.push_back(b); + }; + const std::vector& tri = p.triIndices(); + for (size_t k = 0; k + 2 < tri.size(); k += 3) { + int ia = tri[k], ib = tri[k + 1], ic = tri[k + 2]; + int pa = p.cells[ia].plateId, pb = p.cells[ib].plateId, pc = p.cells[ic].plateId; + if (pa == pb && pb == pc) continue; + if (pa != pb && pb != pc && pa != pc) { + Vec3 c = ((p.cells[ia].unit + p.cells[ib].unit + p.cells[ic].unit) + * (1.0 / 3.0)).normalized() * radius; + Vector3 C{ (float)c.x, (float)c.y, (float)c.z }; + emit(C, midV(ia, ib), pa, pb); + emit(C, midV(ib, ic), pb, pc); + emit(C, midV(ic, ia), pc, pa); + } else { + int lone, o1, o2; + if (pa == pb) { lone = ic; o1 = ia; o2 = ib; } + else if (pb == pc) { lone = ia; o1 = ib; o2 = ic; } + else { lone = ib; o1 = ia; o2 = ic; } + emit(midV(lone, o1), midV(lone, o2), p.cells[lone].plateId, p.cells[o1].plateId); + } + } +} + +void buildDriftArrows(const Planet& p, float radius, + std::vector& out, std::vector& labels) { + out.clear(); labels.clear(); + int np = (int)p.plates.size(); + if (np == 0) return; + std::vector sum(np, Vec3{0, 0, 0}); + std::vector cnt(np, 0); + for (const auto& c : p.cells) + if (c.plateId >= 0 && c.plateId < np) { sum[c.plateId] = sum[c.plateId] + c.unit; ++cnt[c.plateId]; } + std::vector centroid(np), vel(np); + double vmax = 1e-30; + for (int i = 0; i < np; ++i) { + if (cnt[i] == 0) continue; + centroid[i] = sum[i].normalized(); + Vec3 omega = p.plates[i].driftAxis * p.plates[i].driftSpeed; + vel[i] = omega.cross(centroid[i]); + vmax = std::max(vmax, vel[i].length()); + } + const double shaftMin = 0.12, shaftMax = 0.45; + auto V = [](const Vec3& v) { return Vector3{ (float)v.x, (float)v.y, (float)v.z }; }; + for (int i = 0; i < np; ++i) { + if (cnt[i] == 0 || vel[i].length() < 1e-12) continue; + Vec3 dir = vel[i].normalized(); + double len = shaftMin + (shaftMax - shaftMin) * (vel[i].length() / vmax); + Vec3 base = centroid[i] * (double)radius; + Vec3 tip = base + dir * len; + Vec3 perp = dir.cross(centroid[i]).normalized(); + double hl = len * 0.30; + Vec3 back = dir * -1.0; + Vec3 h1 = tip + (back * 0.8 + perp * 0.6) * hl; + Vec3 h2 = tip + (back * 0.8 - perp * 0.6) * hl; + out.push_back(V(base)); out.push_back(V(tip)); + out.push_back(V(tip)); out.push_back(V(h1)); + out.push_back(V(tip)); out.push_back(V(h2)); + + Vec3 lbl = centroid[i] * (radius + 0.015); // just above the planet surface at the plate centroid + labels.push_back({i, V(lbl)}); + } +} + +void buildRivers(const Planet& p, float radius, + std::vector& rivers, std::vector& bigRivers) { + rivers.clear(); bigRivers.clear(); + const std::vector& disc = p.discharge(); + const std::vector& flow = p.flowTo(); + if (disc.empty() || flow.empty()) return; + const double thr = p.cfg.riverThreshold; + auto V = [&](int i) { Vec3 m = p.cells[i].unit * (double)radius; return Vector3{ (float)m.x, (float)m.y, (float)m.z }; }; + for (int i = 0; i < (int)p.cells.size(); ++i) { + int d = flow[i]; + if (d < 0 || disc[i] < thr) continue; // not a river / reached the sea + (disc[i] > thr * 6.0 ? bigRivers : rivers).push_back(V(i)); + (disc[i] > thr * 6.0 ? bigRivers : rivers).push_back(V(d)); + } +} + +std::vector> buildGraticule() { + std::vector> g; + const double D = M_PI / 180.0; + for (int la = -60; la <= 60; la += 30) { // parallels + std::vector pl; + for (int lo = -180; lo <= 180; lo += 5) pl.push_back({ (float)(lo * D), (float)(la * D) }); + g.push_back(pl); + } + for (int lo = -180; lo < 180; lo += 30) { // meridians + std::vector pl; + for (int la = -85; la <= 85; la += 5) pl.push_back({ (float)(lo * D), (float)(la * D) }); + g.push_back(pl); + } + return g; +} +void drawGraticule3D(const std::vector>& g, float radius) { + rlBegin(RL_LINES); rlColor4ub(110, 125, 150, 150); + for (const auto& pl : g) + for (size_t i = 0; i + 1 < pl.size(); ++i) { + Vec3 a = lonLatToDir(pl[i].x, pl[i].y) * (double)radius; + Vec3 b = lonLatToDir(pl[i + 1].x, pl[i + 1].y) * (double)radius; + rlVertex3f((float)a.x, (float)a.y, (float)a.z); + rlVertex3f((float)b.x, (float)b.y, (float)b.z); + } + rlEnd(); +} +void drawGraticule2D(const std::vector>& g, Rectangle r, double lonOffset) { + rlBegin(RL_LINES); rlColor4ub(110, 125, 150, 150); + for (const auto& pl : g) + for (size_t i = 0; i + 1 < pl.size(); ++i) { + Vector2 pa = projLonLat(pl[i].x, pl[i].y, lonOffset, r); + Vector2 pb = projLonLat(pl[i + 1].x, pl[i + 1].y, lonOffset, r); + if (fabsf(pa.x - pb.x) > r.width * 0.5f) continue; + rlVertex2f(pa.x, pa.y); rlVertex2f(pb.x, pb.y); + } + rlEnd(); +} + +void drawGraticuleLabels2D(Rectangle r, double lonOffset) { + const double D = M_PI / 180.0; + const Color col{ 170, 185, 205, 220 }; + // Latitudes at the left edge (Equal Earth y depends only on lat). + const int lats[] = { 60, 30, 0, -30, -60 }; + for (int la : lats) { + Vector2 p = projLonLat(0.0, la * D, lonOffset, r); + const char* t = (la == 0) ? "0" : TextFormat("%d%c", std::abs(la), la > 0 ? 'N' : 'S'); + DrawText(t, (int)r.x + 3, (int)p.y - 6, 11, col); + } + // Longitudes along the bottom edge (skip any panned off the map). + const int lons[] = { -180, -120, -60, 0, 60, 120, 180 }; + int by = (int)(r.y + r.height) - 14; + for (int lo : lons) { + Vector2 p = projLonLat(lo * D, 0.0, lonOffset, r); + if (p.x < r.x + 2 || p.x > r.x + r.width - 2) continue; + const char* t = (lo == 0) ? "0" : TextFormat("%d%c", std::abs(lo), lo > 0 ? 'E' : 'W'); + DrawText(t, (int)p.x - 8, by, 11, col); + } +} + +void drawSegments2D(const std::vector& segs, Color col, float width, + Rectangle r, double lonOffset) { + if (segs.empty()) return; + rlSetLineWidth(width); rlBegin(RL_LINES); rlColor4ub(col.r, col.g, col.b, 255); + for (size_t i = 0; i + 1 < segs.size(); i += 2) { + Vec3 a = Vec3{segs[i].x, segs[i].y, segs[i].z}.normalized(); + Vec3 b = Vec3{segs[i + 1].x, segs[i + 1].y, segs[i + 1].z}.normalized(); + double alo, ala, blo, bla; dirToLonLat(a, alo, ala); dirToLonLat(b, blo, bla); + Vector2 pa = projLonLat(alo, ala, lonOffset, r), pb = projLonLat(blo, bla, lonOffset, r); + if (fabsf(pa.x - pb.x) > r.width * 0.5f) continue; + rlVertex2f(pa.x, pa.y); rlVertex2f(pb.x, pb.y); + } + rlEnd(); rlSetLineWidth(1.0f); +} + +void drawSubgrids(const std::vector>& sgs, + float visBase, float elevExagg, double seaLevel, float eps) { + for (const auto& sg : sgs) { + if (!sg || sg->res < 2) continue; + int R = sg->res; + rlBegin(RL_TRIANGLES); + for (int j = 0; j < R - 1; ++j) + for (int i = 0; i < R - 1; ++i) { + const SubCell* q[4] = { + &sg->sub[(size_t)j * R + i], &sg->sub[(size_t)j * R + i + 1], + &sg->sub[(size_t)(j + 1) * R + i + 1], &sg->sub[(size_t)(j + 1) * R + i] + }; + const int order[6] = { 0, 1, 2, 0, 2, 3 }; + for (int o : order) { + const SubCell* s = q[o]; + Color col = elevationColor(s->elevation, seaLevel); + float rr = visBase + (float)s->elevation * elevExagg + eps; + rlColor4ub(col.r, col.g, col.b, 255); + rlVertex3f((float)(s->unit.x * rr), (float)(s->unit.y * rr), (float)(s->unit.z * rr)); + } + } + rlEnd(); + } +} diff --git a/src/render/Overlays.hpp b/src/render/Overlays.hpp new file mode 100644 index 0000000..ba074f5 --- /dev/null +++ b/src/render/Overlays.hpp @@ -0,0 +1,43 @@ +#pragma once +#include "raylib.h" +#include "Planet.hpp" +#include +#include + +// World overlays drawn over the globe (3D) and/or the Equal Earth map (2D): +// plate borders, per-plate drift arrows + labels, Phase-3 rivers, lat/lon +// graticule, projected line segments, and the selected-cell subgrid patch. + +// ---- Plate borders (dual contour through boundary triangles) ---------------- +// Segments separating two real plates go to `real` (drawn yellow); any segment +// touching a baby (young spreading-ridge) plate goes to `ridge` (drawn red). +void buildBorders(const Planet& p, float radius, + std::vector& real, std::vector& ridge); + +// ---- Per-plate drift arrows ------------------------------------------------- +struct PlateLabel { int id; Vector3 pos; }; +void buildDriftArrows(const Planet& p, float radius, + std::vector& out, std::vector& labels); + +// ---- Phase-3 rivers: drainage-network segments (cell -> its downstream cell) - +// Split into normal and "big" rivers by discharge so they can be drawn at two +// line widths. Needs Planet::computeHydrology() to have been called. +void buildRivers(const Planet& p, float radius, + std::vector& rivers, std::vector& bigRivers); + +// ---- Lat/lon graticule ------------------------------------------------------ +// Polylines of (lon,lat) radians (Vector2.x=lon, .y=lat). +std::vector> buildGraticule(); +void drawGraticule3D(const std::vector>& g, float radius); +void drawGraticule2D(const std::vector>& g, Rectangle r, double lonOffset); +// Degree numbers along the 2D map edges (latitudes at the left, longitudes at the +// bottom). lonOffset is the map's east/west pan. Draw inside the map scissor. +void drawGraticuleLabels2D(Rectangle r, double lonOffset); + +// Project 3D sphere line segments (pairs of points) onto the 2D map. +void drawSegments2D(const std::vector& segs, Color col, float width, + Rectangle r, double lonOffset); + +// ---- Subgrid overlay (high-res patch over a selected cell) ------------------ +void drawSubgrids(const std::vector>& sgs, + float visBase, float elevExagg, double seaLevel, float eps); diff --git a/src/render/Panels.cpp b/src/render/Panels.cpp new file mode 100644 index 0000000..9aed36c --- /dev/null +++ b/src/render/Panels.cpp @@ -0,0 +1,162 @@ +#include "Panels.hpp" +#include "Colors.hpp" // elevationColor (subtile grid) +#include "Projection.hpp" // dirToLonLat +#include +#include +#include +#include + +// elev/age come from the display snapshot so the readout matches what is drawn. +static std::vector cellInfo(const Planet& p, int i, double elev, double age) { + const Cell& c = p.cells[i]; + double lon, lat; dirToLonLat(c.unit, lon, lat); + const Plate& pl = p.plates[c.plateId]; + const int n = (int)p.cells.size(); + auto sized = [&](const std::vector& v) { return (int)v.size() == n; }; + std::vector L; + L.push_back(std::string(TextFormat("Cell #%d", i))); + L.push_back(std::string(TextFormat("lat %+6.1f lon %+6.1f", lat * 180.0 / M_PI, lon * 180.0 / M_PI))); + L.push_back(std::string(TextFormat("elev %.0f m (%s)", elev, + elev < p.cfg.seaLevel ? "ocean" : "land"))); + L.push_back(std::string(TextFormat("plate %d (%s) crust %s", c.plateId, + pl.type == PlateType::Oceanic ? "Oceanic" : "Continental", + c.oceanic ? "Oceanic" : "Continental"))); + L.push_back(std::string(TextFormat("biome: %s", biomeName(c.biome)))); + // Climate (derived; present once computeClimate() has run). + if (sized(p.temperature()) && sized(p.moisture())) + L.push_back(std::string(TextFormat("temp %.1f C precip %.0f%%", + p.temperature()[i], p.moisture()[i] * 100.0))); + L.push_back(std::string(TextFormat("geoAge %.0f My neighbors %d", age, (int)c.neighbors.size()))); + // Hydrology (derived; present once routeFlow()/hydrology() has run). + if (sized(p.discharge()) && p.discharge()[i] > p.cfg.riverThreshold) + L.push_back(std::string(TextFormat("river: discharge %.0f", p.discharge()[i]))); + if (sized(p.lakeDepth()) && p.lakeDepth()[i] > p.cfg.biomeLakeMinDepth && elev > p.cfg.seaLevel) + L.push_back(std::string(TextFormat("lake: depth %.0f m", p.lakeDepth()[i]))); + return L; +} + +void drawDetailPanel(const Planet& p, const std::shared_ptr& sg, + int macro, double macroElev, double macroAge, + Rectangle panel, Rectangle grid, int hoveredSub) { + DrawRectangleRec(panel, Color{12, 14, 22, 235}); + DrawRectangleLinesEx(panel, 1, Color{120, 120, 150, 255}); + int tx = (int)panel.x + 10, ty = (int)panel.y + 8; + DrawText(TextFormat("Tile #%d", macro), tx, ty, 20, RAYWHITE); + DrawText("C: close", (int)(panel.x + panel.width) - 78, ty + 4, 14, Color{170, 170, 185, 255}); + ty += 28; + // Cramped above the subtile grid -> stop before overlapping it (the full list is + // always shown in the top-right hover panel, which has room). + for (auto& s : cellInfo(p, macro, macroElev, macroAge)) { + if (ty + 18 > (int)grid.y) break; + DrawText(s.c_str(), tx, ty, 15, Color{210, 210, 220, 255}); ty += 18; + } + + if (!sg || sg->res < 2) return; + int R = sg->res; + float cw = grid.width / R, ch = grid.height / R; + DrawText(TextFormat("Subtiles %dx%d (elevation; dim = neighbor)", R, R), + (int)grid.x, (int)grid.y - 18, 14, Color{200, 200, 210, 255}); + for (int j = 0; j < R; ++j) + for (int i = 0; i < R; ++i) { + const SubCell& s = sg->sub[(size_t)j * R + i]; + Color col = elevationColor(s.elevation, p.cfg.seaLevel); + if (s.nearestMacro != macro) { // territory of a neighbor + col.r = (unsigned char)(col.r * 0.55); col.g = (unsigned char)(col.g * 0.55); + col.b = (unsigned char)(col.b * 0.55); + } + DrawRectangle((int)(grid.x + i * cw), (int)(grid.y + j * ch), + (int)std::ceil(cw), (int)std::ceil(ch), col); + } + DrawRectangleLinesEx(grid, 1, Color{90, 90, 110, 255}); + + int by = (int)(grid.y + grid.height) + 6; + if (hoveredSub >= 0) { + int i = hoveredSub % R, j = hoveredSub / R; + DrawRectangleLinesEx(Rectangle{grid.x + i * cw, grid.y + j * ch, cw, ch}, 2, WHITE); + const SubCell& s = sg->sub[hoveredSub]; + double lon, lat; dirToLonLat(s.unit, lon, lat); + DrawText(TextFormat("subtile [%d,%d] elev %.0f m", i, j, s.elevation), + (int)panel.x + 10, by, 15, RAYWHITE); + DrawText(TextFormat("under macro #%d lat %+.2f lon %+.2f", + s.nearestMacro, lat * 180.0 / M_PI, lon * 180.0 / M_PI), + (int)panel.x + 10, by + 18, 14, Color{200, 200, 210, 255}); + } else { + DrawText("hover a subtile for detail", (int)panel.x + 10, by, 14, Color{170, 170, 185, 255}); + } +} + +void drawHoverPanel(const Planet& p, Rectangle r, int hovered, int selected) { + DrawRectangleRec(r, Color{12, 14, 22, 235}); + DrawRectangleLinesEx(r, 1, Color{120, 120, 150, 255}); + int x = (int)r.x + 18, y = (int)r.y + 14; + DrawText("Cell info", x, y, 24, RAYWHITE); + y += 46; + int shown = (hovered >= 0) ? hovered : selected; + if (shown < 0) { + DrawText("hover the 3D globe or the 2D map", x, y, 20, Color{170, 170, 185, 255}); + return; + } + if (hovered < 0) { DrawText("(selected tile)", x, y, 18, Color{210, 180, 120, 255}); y += 30; } + for (auto& s : cellInfo(p, shown, p.cells[shown].elevation, p.cells[shown].geoAge)) { + DrawText(s.c_str(), x, y, 24, Color{215, 220, 230, 255}); y += 32; + } +} + +void drawStats(const Planet& p, Rectangle r, double elapsedMy, bool drifting) { + DrawRectangleRec(r, Color{12, 14, 22, 235}); + DrawRectangleLinesEx(r, 1, Color{120, 120, 150, 255}); + int x = (int)r.x + 16, y = (int)r.y + 12; + DrawText("World statistics", x, y, 22, RAYWHITE); y += 38; + + int N = (int)p.cells.size(), np = (int)p.plates.size(); + double Rkm = p.cfg.radius / 1000.0; + double surfKm2 = 4.0 * M_PI * Rkm * Rkm, cellKm2 = surfKm2 / std::max(1, N); + std::vector pc(np, 0), pl(np, 0); + double seaLvl = p.cfg.seaLevel; + int landGeo = 0, cont = 0; double mn = 1e30, mx = -1e30, sum = 0; + for (const auto& c : p.cells) { + if (c.plateId >= 0 && c.plateId < np) { pc[c.plateId]++; if (!c.oceanic) pl[c.plateId]++; } + if (c.elevation > seaLvl) ++landGeo; // geographic land (above sea level) + if (!c.oceanic) ++cont; // continental crust + mn = std::min(mn, c.elevation); mx = std::max(mx, c.elevation); sum += c.elevation; + } + // Real plates vs. baby (young spreading-ridge) plates are counted separately. + int contPlates = 0, used = 0, babyPlates = 0, babyCells = 0; + for (int q = 0; q < np; ++q) { + if (pc[q] == 0) continue; + if (p.plates[q].baby) { ++babyPlates; babyCells += pc[q]; continue; } + ++used; + if (pl[q] * 2 > pc[q]) ++contPlates; + } + + auto L = [&](const char* s) { DrawText(s, x, y, 17, Color{210, 215, 225, 255}); y += 23; }; + L(TextFormat("Cells: %d cell area %.1fk km2 R %.0f km", N, cellKm2 / 1000.0, Rkm)); + L(TextFormat("Surface area: %.0f M km2", surfKm2 / 1.0e6)); + L(TextFormat("Plates: %d active %d continental / %d oceanic", used, contPlates, used - contPlates)); + L(TextFormat("Young ridges: %d strips %d cells", babyPlates, babyCells)); + int water = N - landGeo; + double wlRatio = landGeo > 0 ? (double)water / landGeo : 0.0; + L(TextFormat("Land %.0f%% Ocean %.0f%% (water:land %.2f:1)", + 100.0 * landGeo / N, 100.0 * water / N, wlRatio)); + L(TextFormat("Sea level: %+.0f m", seaLvl)); + L(TextFormat("Crust: %.0f%% continental / %.0f%% oceanic", 100.0 * cont / N, 100.0 * (N - cont) / N)); + L(TextFormat("Elevation: %.0f .. %.0f m mean %.0f m", mn, mx, sum / N)); + if (drifting) L(TextFormat("Sim time: %.0f My", elapsedMy)); + y += 8; + DrawText("plate cells size area speed land", x, y, 15, Color{150, 155, 170, 255}); y += 21; + + std::vector idx(np); for (int q = 0; q < np; ++q) idx[q] = q; + std::sort(idx.begin(), idx.end(), [&](int a, int b){ return pc[a] > pc[b]; }); + int rows = 0, rowMax = 13; + for (int q : idx) { + if (pc[q] == 0 || p.plates[q].baby) continue; // baby ridges summarised above + if (rows++ >= rowMax) break; + const char* ty = (pl[q] * 2 > pc[q]) ? "cont" : "ocn "; + DrawText(TextFormat("P%-2d %s %5dc %4.1f%% %6.1fM %4.1fcm/y %3.0f%%", q, ty, pc[q], + 100.0 * pc[q] / N, pc[q] * cellKm2 / 1.0e6, p.plates[q].speedCmYr, 100.0 * pl[q] / pc[q]), + x, y, 16, Color{200, 205, 220, 255}); + y += 21; + } + DrawText("click a tile to inspect its subtiles", + x, (int)(r.y + r.height) - 24, 14, Color{150, 150, 165, 255}); +} diff --git a/src/render/Panels.hpp b/src/render/Panels.hpp new file mode 100644 index 0000000..188177d --- /dev/null +++ b/src/render/Panels.hpp @@ -0,0 +1,20 @@ +#pragma once +#include "raylib.h" +#include "Planet.hpp" +#include + +// Right-column UI panels: the tile detail panel (subtiles grid), the hover/ +// selection info box, and the world-statistics panel. + +// Tile detail panel: tile info header + the selected cell's subgrid drawn as a +// flat, hoverable grid of subtiles (neighbor-owned subtiles dimmed). hoveredSub +// is the grid index under the cursor (-1 if none). +void drawDetailPanel(const Planet& p, const std::shared_ptr& sg, + int macro, double macroElev, double macroAge, + Rectangle panel, Rectangle grid, int hoveredSub); + +// Top-right quadrant: info for the hovered cell (or the selected one if none). +void drawHoverPanel(const Planet& p, Rectangle r, int hovered, int selected); + +// World statistics panel (shown in the subareas quadrant when no tile selected). +void drawStats(const Planet& p, Rectangle r, double elapsedMy, bool drifting); diff --git a/src/render/Picking.cpp b/src/render/Picking.cpp new file mode 100644 index 0000000..0ab5186 --- /dev/null +++ b/src/render/Picking.cpp @@ -0,0 +1,47 @@ +#include "Picking.hpp" +#include +#include + +Vec3 rayDirFromMouse(Vector3 pos, Vector3 target, float fovyDeg, + float mx, float my, float vw, float vh) { + Vec3 P{pos.x, pos.y, pos.z}, T{target.x, target.y, target.z}; + Vec3 f = (T - P).normalized(); + Vec3 r = f.cross(Vec3{0, 1, 0}).normalized(); + Vec3 u = r.cross(f); + double nx = 2.0 * mx / vw - 1.0; + double ny = 1.0 - 2.0 * my / vh; + double tanV = std::tan(fovyDeg * 0.5 * M_PI / 180.0); + double tanH = tanV * (double)vw / vh; + return (f + r * (nx * tanH) + u * (ny * tanV)).normalized(); +} + +bool raySphere(Vec3 o, Vec3 d, double R, Vec3& hitUnit) { + double b = 2.0 * o.dot(d); + double c = o.dot(o) - R * R; + double disc = b * b - 4.0 * c; + if (disc < 0) return false; + double sq = std::sqrt(disc); + double t = (-b - sq) * 0.5; + if (t < 0) t = (-b + sq) * 0.5; + if (t < 0) return false; + hitUnit = (o + d * t).normalized(); + return true; +} + +int nearestCell(const Planet& p, const Vec3& dir) { + int best = -1; double bd = -2.0; + for (size_t i = 0; i < p.cells.size(); ++i) { + double dd = p.cells[i].unit.dot(dir); + if (dd > bd) { bd = dd; best = (int)i; } + } + return best; +} + +double angBetween(const Vec3& a, const Vec3& b) { + return std::acos(std::clamp(a.dot(b), -1.0, 1.0)); +} + +Vec3 rotateZ(const Vec3& v, double deg) { + double r = deg * M_PI / 180.0, c = std::cos(r), s = std::sin(r); + return Vec3{ v.x * c - v.y * s, v.x * s + v.y * c, v.z }; +} diff --git a/src/render/Picking.hpp b/src/render/Picking.hpp new file mode 100644 index 0000000..a092198 --- /dev/null +++ b/src/render/Picking.hpp @@ -0,0 +1,19 @@ +#pragma once +#include "raylib.h" +#include "Vec3.hpp" +#include "Planet.hpp" + +// Mouse-picking helpers: ray from the 3D viewport, sphere intersection, nearest +// cell to a direction, and angular distance between two unit directions. + +// Custom camera ray for the 3D viewport (GetScreenToWorldRay assumes the full +// screen; the globe occupies only the top-left quadrant). +Vec3 rayDirFromMouse(Vector3 pos, Vector3 target, float fovyDeg, + float mx, float my, float vw, float vh); +bool raySphere(Vec3 o, Vec3 d, double R, Vec3& hitUnit); +int nearestCell(const Planet& p, const Vec3& dir); +double angBetween(const Vec3& a, const Vec3& b); + +// Rotate v by `deg` degrees CCW about the world +Z axis (matches rlRotatef(deg,0,0,1)). +// Used for the axial-tilt transform: +tilt maps model->world, -tilt world->model. +Vec3 rotateZ(const Vec3& v, double deg); diff --git a/src/render/Viewer.cpp b/src/render/Viewer.cpp new file mode 100644 index 0000000..48c7ce5 --- /dev/null +++ b/src/render/Viewer.cpp @@ -0,0 +1,254 @@ +#include "Viewer.hpp" +#include "Picking.hpp" // angBetween (rebuildSub) +#include "Projection.hpp" // EqualEarth (layout) +#include +#include +#include +#include + +bool Viewer::init(int argc, char** argv) { + uint32_t cliSeed = 0; // 0 = no --seed given + for (int a = 1; a < argc; ++a) { + if (!std::strcmp(argv[a], "--seed") && a + 1 < argc) + cliSeed = (uint32_t)std::strtoul(argv[++a], nullptr, 10); + else if (!std::strcmp(argv[a], "--config") && a + 1 < argc) + configPath = argv[++a]; + } + + SetConfigFlags(FLAG_MSAA_4X_HINT); + InitWindow(screenW, screenH, "Planet Sim - Phase 1: Tectonics"); + SetTargetFPS(60); + + // Layout: left column 70% wide (3D globe 60% h on top, 2D map 40% h below); + // right column 30% wide (cell info 50% h on top, subareas 50% below). + leftW = (int)(screenW * 0.70f); // 1344 + rightX = leftW; rightW = screenW - leftW; // 576 + rightH = screenH / 2; // 540 + + // Top-left: 3D globe (render texture, its own aspect). + view3DW = leftW; view3DH = (int)(screenH * 0.60f); // 1344 x 648 + rt3d = LoadRenderTexture(view3DW, view3DH); + SetTextureFilter(rt3d.texture, TEXTURE_FILTER_BILINEAR); + + // Bottom-left: 2D Equal Earth map, fit (keep aspect) into the 40% strip. + const int mapAreaY = view3DH, mapAreaH = screenH - view3DH; // (0,648) 1344 x 432 + int mapH = mapAreaH - 30; + int mapW = (int)(mapH * (EqualEarth::halfWidth() / EqualEarth::halfHeight())); + if (mapW > leftW - 30) { mapW = leftW - 30; mapH = (int)(mapW / (EqualEarth::halfWidth() / EqualEarth::halfHeight())); } + mapRect = Rectangle{ (float)((leftW - mapW) / 2), + (float)(mapAreaY + (mapAreaH - mapH) / 2), + (float)mapW, (float)mapH }; + + // Right column. + hoverRect = Rectangle{ (float)rightX + 8, 8.0f, (float)rightW - 16, (float)rightH - 16 }; + panelRect = Rectangle{ (float)rightX + 8, (float)rightH + 8, (float)rightW - 16, (float)rightH - 16 }; + const float panelHeader = 120.0f; + const float gridSide = std::min(panelRect.width - 40.0f, panelRect.height - panelHeader - 56.0f); + gridRect = Rectangle{ panelRect.x + (panelRect.width - gridSide) / 2.0f, + panelRect.y + panelHeader, gridSide, gridSide }; + + // Buttons (pause + Phase-3 prompt, centered in the 3D viewport). + pauseBtn = Rectangle{ 16.0f, (float)view3DH - 44.0f, 160.0f, 32.0f }; + const float pbW = 220.0f, pbH = 42.0f, pbGap = 24.0f; + pbCx = view3DW * 0.5f; pbCy = view3DH * 0.5f; + p3ContinueBtn = Rectangle{ pbCx - pbW - pbGap * 0.5f, pbCy + 8.0f, pbW, pbH }; + p3StartBtn = Rectangle{ pbCx + pbGap * 0.5f, pbCy + 8.0f, pbW, pbH }; + + cfg.subdivisions = 5; + if (!loadConfig(configPath, cfg)) saveConfig(configPath, cfg); // load, or create a default + if (cliSeed != 0) cfg.seed = cliSeed; // CLI --seed overrides config + std::string cfgErr = validateConfig(cfg); + if (!cfgErr.empty()) cfg = PlanetConfig{}; // revert to safe defaults + planet.generate(cfg); + + cam.position = {0, 0, 6}; cam.target = {0, 0, 0}; cam.up = {0, 1, 0}; + cam.fovy = 45; cam.projection = CAMERA_PERSPECTIVE; + + buildBorders(planet, borderR, borders, ridgeBorders); + buildDriftArrows(planet, driftR, driftArrows, plateLabels); + graticule = buildGraticule(); + buildMap2D(planet, mapRect, map2D); + phase3PromptAt = planet.cfg.phase3AfterMy; + + refreshView(); // colour the freshly generated (flat) world + return true; +} + +void Viewer::rebuildSub() { + subgrids.clear(); + if (selectedCell < 0) return; + subgrids.push_back(planet.makeSubGrid(selectedCell, subRes)); + const Cell& c = planet.cells[selectedCell]; + for (int nb : c.neighbors) subgrids.push_back(planet.makeSubGrid(nb, subRes)); + double ma = 0.0; + for (int nb : c.neighbors) ma += angBetween(c.unit, planet.cells[nb].unit); + ma /= std::max(1, c.neighbors.size()); + selectedThresh = ma * 1.4; +} + +void Viewer::selectCell(int idx) { + if (idx < 0) return; + if (idx == selectedCell) { selectedCell = -1; subgrids.clear(); return; } + selectedCell = idx; rebuildSub(); +} + +// Recolor the mesh + refresh the elevation range, read straight from cells. +void Viewer::recolor() { + double maxAge = 1.0; for (const auto& c : planet.cells) maxAge = std::max(maxAge, c.geoAge); + const std::vector& temp = planet.temperature(); + const std::vector& moist = planet.moisture(); // 0..1, already robustly normalized + vcolors.resize(planet.cells.size()); + for (size_t i = 0; i < planet.cells.size(); ++i) { + switch (mode) { + case ColorMode::Plate: { + int pid = planet.cells[i].plateId; + vcolors[i] = (pid >= 0 && pid < (int)planet.plates.size() && planet.plates[pid].baby) + ? Color{70, 80, 95, 255} // young spreading-ridge crust + : plateColor(pid); + break; + } + case ColorMode::Age: vcolors[i] = ageColor(planet.cells[i].geoAge, maxAge); break; + case ColorMode::Crust: vcolors[i] = crustColor(planet.cells[i].oceanic); break; + case ColorMode::Biome: vcolors[i] = biomeColor(planet.cells[i].biome); break; + case ColorMode::Temperature: vcolors[i] = temp.empty() ? Color{90,90,90,255} : tempColor(temp[i]); break; + case ColorMode::Precip: vcolors[i] = moist.empty() ? Color{90,90,90,255} : precipColor(moist[i]); break; + default: vcolors[i] = elevationColor(planet.cells[i].elevation, planet.cfg.seaLevel); + } + } + // Phase 3: shade filled basins above sea level as inland water (lakes). + const std::vector& lk = planet.lakeDepth(); + if (phase3 && !lk.empty()) + for (size_t i = 0; i < planet.cells.size(); ++i) + if (lk[i] > 20.0 && planet.cells[i].elevation > planet.cfg.seaLevel) + vcolors[i] = lakeColor(); + minE = planet.minElevation(); maxE = planet.maxElevation(); +} + +void Viewer::refreshView() { + if (phase3) planet.computeHydrology(); // refresh lakes/rivers for the view + planet.computeClimate(); // temperature + precipitation fields + planet.classifyBiomes(); // keep cell.biome current (reads the climate) + recolor(); + if (settled) { // Phase 2: plates moved -> boundaries moved + buildBorders(planet, borderR, borders, ridgeBorders); + buildDriftArrows(planet, driftR, driftArrows, plateLabels); + } + if (phase3) buildRivers(planet, riverR, rivers, bigRivers); + if (selectedCell >= 0) rebuildSub(); +} + +void Viewer::regenWorld() { // after generate(): geometry changed + buildBorders(planet, borderR, borders, ridgeBorders); + buildDriftArrows(planet, driftR, driftArrows, plateLabels); + buildMap2D(planet, mapRect, map2D); + selectedCell = -1; subgrids.clear(); + settled = false; settleRun = 0; formAccum = 0.0; stepCount = 0; paused = false; + planet.drifting = false; // Phase 1: original forming behavior + phase3 = false; phase3Prompt = false; phase3PromptAt = planet.cfg.phase3AfterMy; + rivers.clear(); bigRivers.clear(); + elapsedMy = 0.0; dtMy = 0.0; driftAccum = 0.0; + refreshView(); +} + +void Viewer::regen() { planet.generate(cfg); regenWorld(); } + +void Viewer::stepOnce() { // one tick + settle bookkeeping + maxChange = planet.step(); ++stepCount; + if (maxChange < settleThresh) { if (++settleRun >= settleNeed) settled = true; } + else settleRun = 0; +} + +void Viewer::pauseAction() { paused = !paused; } // pause/resume forming or drift + +void Viewer::setStatus(const std::string& m) { statusMsg = m; statusUntil = GetTime() + 3.0; } + +// F5: write seed + config + full planet state. F9: read it back and resume. +void Viewer::saveGame(const char* path) { + std::ofstream os(path, std::ios::binary); + if (!os) { setStatus("Save failed"); return; } + uint32_t ver = SAVE_VERSION; uint8_t st = settled ? 1 : 0; uint8_t p3 = phase3 ? 1 : 0; + os.write("PLSV", 4); + os.write(reinterpret_cast(&ver), sizeof ver); + os.write(reinterpret_cast(&elapsedMy), sizeof elapsedMy); + os.write(reinterpret_cast(&st), sizeof st); + os.write(reinterpret_cast(&driftRate), sizeof driftRate); + os.write(reinterpret_cast(&p3), sizeof p3); // v3: Phase-3 flag + planet.writeState(os); + setStatus(os ? std::string("Saved ") + path : "Save failed"); +} + +void Viewer::loadGame(const char* path) { + std::ifstream is(path, std::ios::binary); + if (!is) { setStatus(std::string("No ") + path); return; } + char magic[4] = {0}; uint32_t ver = 0; double em = 0; uint8_t st = 0; double dr = 4.0; uint8_t p3 = 0; + is.read(magic, 4); + is.read(reinterpret_cast(&ver), sizeof ver); + is.read(reinterpret_cast(&em), sizeof em); + is.read(reinterpret_cast(&st), sizeof st); + if (ver >= 2) is.read(reinterpret_cast(&dr), sizeof dr); + if (ver >= 3) is.read(reinterpret_cast(&p3), sizeof p3); + if (!is || std::memcmp(magic, "PLSV", 4) != 0 || ver > SAVE_VERSION) { setStatus("Load failed: bad file"); return; } + if (!planet.readState(is, ver >= 4)) { setStatus("Load failed: corrupt/mismatch"); return; } // v4: per-cell biome + cfg = planet.cfg; // adopt the loaded config + elapsedMy = em; settled = (st != 0); + planet.drifting = settled; // resume drift boosts iff mid-drift + phase3 = (p3 != 0); phase3Prompt = false; + phase3PromptAt = phase3 ? elapsedMy : (elapsedMy + planet.cfg.phase3AfterMy); + driftRate = dr; + settleRun = settleNeed; // keep the settled latch consistent + dtMy = settled ? planet.cflDtMy() : 0.0; + driftAccum = 0.0; formAccum = 0.0; + paused = true; selectedCell = -1; subgrids.clear(); + buildBorders(planet, borderR, borders, ridgeBorders); + buildDriftArrows(planet, driftR, driftArrows, plateLabels); + buildMap2D(planet, mapRect, map2D); + refreshView(); + setStatus(std::string("Loaded ") + path); +} + +// Advance the simulation this frame: Phase-1 forming (paced ticks toward +// equilibrium), or Phase-2 drift / Phase-3 drift+hydrology at a finer dt. +void Viewer::stepSim() { + if (!paused && !settled) { + // --- Phase 1: forming, paced ticks toward equilibrium ------------- + formAccum += GetFrameTime() * formRate; + int budget = 0; + while (formAccum >= 1.0 && budget < 2000) { + stepOnce(); formAccum -= 1.0; ++budget; + if (settled) break; + } + if (settled) { dtMy = planet.cflDtMy(); planet.drifting = true; } // entering Phase 2 + refreshView(); // live update so you watch the terrain rise + } else if (!paused && settled) { + // --- Phase 2 drift (and Phase 3 = drift + hydrology at a finer dt) -- + // Drift never stops; Phase 3 just uses a smaller timestep so each step + // advances fewer My (more steps before plates visibly move) while + // rivers/lakes/fluvial erosion resolve. + double dt = planet.cflDtMy() * (phase3 ? planet.cfg.phase3DtScale : 1.0); + dtMy = dt; + driftAccum += driftRate * GetFrameTime(); // accumulate across frames + const int guardMax = phase3 ? 60 : 500; // Phase-3 ticks are heavier + int guard = 0; bool advanced = false; + while (driftAccum >= dt && guard < guardMax) { + planet.advect(dt); planet.step(); planet.erode(dt); + if (phase3) planet.hydrology(dt); + elapsedMy += dt; driftAccum -= dt; ++guard; advanced = true; + // Timed Phase-3 invitation: pause + prompt once we cross the mark. + if (!phase3 && elapsedMy >= phase3PromptAt) { phase3Prompt = true; paused = true; break; } + } + if (driftAccum > 2.0 * dt) driftAccum = 2.0 * dt; // drop backlog (don't runaway) + if (advanced) refreshView(); // live: watch the world evolve + } +} + +void Viewer::run() { + while (!WindowShouldClose()) { + handleInput(); + stepSim(); + // A regenerate this frame may have shrunk the planet; keep indices valid. + if (hovered >= (int)planet.cells.size()) hovered = -1; + renderFrame(); + } + UnloadRenderTexture(rt3d); + CloseWindow(); +} diff --git a/src/render/Viewer.hpp b/src/render/Viewer.hpp new file mode 100644 index 0000000..1035fa6 --- /dev/null +++ b/src/render/Viewer.hpp @@ -0,0 +1,124 @@ +#pragma once +#include "raylib.h" +#include "Planet.hpp" +#include "Colors.hpp" // ColorMode +#include "Overlays.hpp" // PlateLabel +#include "Map2D.hpp" // Map2D +#include +#include +#include + +// The interactive viewer: owns all window/sim/view state and runs the frame +// loop. The old free-standing main() lived as one giant function with capturing +// lambdas; those lambdas are now methods and their captured locals are members, +// so the body splits cleanly across Viewer.cpp (setup + sim orchestration), +// ViewerInput.cpp (input/picking/keys) and ViewerRender.cpp (drawing). +struct Viewer { + // ---- Files / save format ------------------------------------------------ + static constexpr uint32_t SAVE_VERSION = 6; // v6: self-describing (text) config block; v4: +biome; v3: +phase3 + const char* CONFIG_PATH = "planet.cfg"; + const char* SAVE_PATH = "planet.save"; + std::string configPath = "planet.cfg"; // initial config (--config overrides) + + // ---- Window / layout (set in init) -------------------------------------- + int screenW = 1920, screenH = 1080; + int leftW = 0, rightX = 0, rightW = 0, rightH = 0; + int view3DW = 0, view3DH = 0; + RenderTexture2D rt3d{}; + Rectangle mapRect{}, hoverRect{}, panelRect{}, gridRect{}; + Rectangle pauseBtn{}, p3ContinueBtn{}, p3StartBtn{}; + float pbCx = 0.0f, pbCy = 0.0f; + + const float visBase = 2.0f; + const float elevExagg = 0.00000004f; + const float borderR = visBase + 0.004f; + const float driftR = visBase + 0.006f; + const float riverR = visBase + 0.005f; + const float gratR = visBase + 0.003f; + const int subRes = 16; + + // ---- Sim / world -------------------------------------------------------- + Planet planet; + PlanetConfig cfg; + + Camera3D cam{}; + float camYaw = 0.4f, camPitch = 0.3f, camDist = 6.0f; + + ColorMode mode = ColorMode::Elevation; + std::vector vcolors; + + std::vector borders, ridgeBorders; bool showBorders = true; + std::vector driftArrows; bool showDrift = true; + std::vector plateLabels; + std::vector rivers, bigRivers; bool showRivers = true; + std::vector> graticule; bool showGrat = false; + Map2D map2D; + + // Phase-1 forming model. + bool paused = false, settled = false; + long long stepCount = 0; + double maxChange = 0.0; + int settleRun = 0; + const double settleThresh = 2.0; // m/tick at or below which it's "settled" + const int settleNeed = 3; // consecutive settled ticks before pausing + const double formRate = 55.0; // forming ticks per second (watchable) + double formAccum = 0.0; + double minE = 0.0, maxE = 0.0; + + // Phase 2 drift. + double elapsedMy = 0.0, dtMy = 0.0, driftRate = 4.0, driftAccum = 0.0; + + // Phase 3 hydrology. + bool phase3 = false, phase3Prompt = false; + double phase3PromptAt = 0.0; + + // Selection + subgrid (phase 4/5 preview). + int selectedCell = -1; + double selectedThresh = 0.06; + std::vector> subgrids; + + // Transient on-screen status line. + std::string statusMsg; double statusUntil = 0.0; + + // Input state. + float dragDist = 0.0f; + double mapLon = 0.0; // 2D map longitude pan (radians) + bool pressInMap = false; // a drag that started on the map pans it + + // Per-frame picking state (set by handleInput, read by render). + Vector2 mp{}; + bool onPause = false; + int hovered = -1; + bool hasHoverSub = false; + SubCell hoverSub; + int hoveredSubIdx = -1; + + // ---- Lifecycle ---------------------------------------------------------- + bool init(int argc, char** argv); // window, layout, config, first world + void run(); // the frame loop (until window closes) + + // ---- Sim orchestration (Viewer.cpp) ------------------------------------- + void rebuildSub(); + void selectCell(int idx); + void recolor(); + void refreshView(); + void regenWorld(); // after generate(): geometry changed + void regen(); // generate(cfg) + regenWorld() + void stepOnce(); // one tick + settle bookkeeping + void pauseAction(); + void setStatus(const std::string& m); + void saveGame(const char* path); + void loadGame(const char* path); + void stepSim(); // advance forming / drift+hydrology this frame + + // ---- Input (ViewerInput.cpp) -------------------------------------------- + void handleInput(); + + // ---- Render (ViewerRender.cpp) ------------------------------------------ + void renderFrame(); + void renderGlobe3D(); + void renderMap2D(); + void renderPanels(); + void renderHUD(); + void renderPrompt(); +}; diff --git a/src/render/ViewerInput.cpp b/src/render/ViewerInput.cpp new file mode 100644 index 0000000..c6d7406 --- /dev/null +++ b/src/render/ViewerInput.cpp @@ -0,0 +1,143 @@ +#include "Viewer.hpp" +#include "Picking.hpp" +#include "Map2D.hpp" // wrapPi +#include "Projection.hpp" // EqualEarth, lonLatToDir +#include +#include + +// One frame of input: camera orbit/zoom, hover picking (panel subtile -> 3D ray +// -> 2D map), click-to-select, and key handling. Writes the per-frame picking +// members (mp/onPause/hovered/hasHoverSub/hoverSub/hoveredSubIdx) for the render +// pass, and may step / regenerate / save / load the world. +void Viewer::handleInput() { + mp = GetMousePosition(); + bool in3D = (mp.x < view3DW && mp.y < view3DH); // top-left quadrant + bool inMap = CheckCollisionPointRec(mp, mapRect); + bool inPanel = (selectedCell >= 0) && CheckCollisionPointRec(mp, panelRect); + onPause = CheckCollisionPointRec(mp, pauseBtn); + + // --- Camera input (LMB drag orbits; tracks drag distance for clicks) -- + if (IsMouseButtonPressed(MOUSE_BUTTON_LEFT)) { + if (phase3Prompt) { // modal: only the two buttons act + if (CheckCollisionPointRec(mp, p3StartBtn)) { + phase3 = true; phase3Prompt = false; paused = false; + phase3PromptAt = elapsedMy; setStatus("Hydrology started"); + refreshView(); + } else if (CheckCollisionPointRec(mp, p3ContinueBtn)) { + phase3Prompt = false; paused = false; + phase3PromptAt = elapsedMy + planet.cfg.phase3AfterMy; + setStatus("Continuing world-building"); + } + } else { + dragDist = 0.0f; + pressInMap = inMap; // drag started on the map -> pan it + if (onPause) pauseAction(); // clickable pause / re-evolve button + } + } + if (IsMouseButtonDown(MOUSE_BUTTON_LEFT) && !phase3Prompt) { + Vector2 d = GetMouseDelta(); + dragDist += fabsf(d.x) + fabsf(d.y); + if (in3D && !onPause) { // only orbit from the 3D quadrant + camYaw += d.x * 0.005f; + camPitch += d.y * 0.005f; + camPitch = std::clamp(camPitch, -1.5f, 1.5f); + } + if (pressInMap) // drag the map east/west + mapLon = wrapPi(mapLon + d.x * (2.0 * M_PI / mapRect.width)); + } + camDist -= GetMouseWheelMove() * 0.4f; + camDist = std::clamp(camDist, 2.6f, 14.0f); + cam.position = { camDist * cosf(camPitch) * sinf(camYaw), + camDist * sinf(camPitch), + camDist * cosf(camPitch) * cosf(camYaw) }; + + // --- Hover picking: tile panel subtile, else 3D ray, else 2D map ------ + // The globe is rendered tilted by axialTilt about world Z; the picking sphere is + // rotation-invariant, so un-rotate the world-space hit direction by -tilt to get + // the model-space direction used to match cells/subcells (hitModel). + hovered = -1; hasHoverSub = false; hoveredSubIdx = -1; + bool have3DHit = false; Vec3 hitUnit, hitModel; + if (inPanel) { + if (!subgrids.empty() && CheckCollisionPointRec(mp, gridRect)) { + const auto& sg = subgrids[0]; int R = sg->res; + int i = std::clamp((int)((mp.x - gridRect.x) / (gridRect.width / R)), 0, R - 1); + int j = std::clamp((int)((mp.y - gridRect.y) / (gridRect.height / R)), 0, R - 1); + hoveredSubIdx = j * R + i; + hoverSub = sg->sub[hoveredSubIdx]; hasHoverSub = true; // marks it on the globe + } + } else if (in3D) { + Vec3 d = rayDirFromMouse(cam.position, cam.target, cam.fovy, + mp.x, mp.y, (float)view3DW, (float)view3DH); + Vec3 o{cam.position.x, cam.position.y, cam.position.z}; + if (raySphere(o, d, visBase, hitUnit)) { + have3DHit = true; + hitModel = rotateZ(hitUnit, -planet.cfg.axialTilt); // world -> model (undo tilt) + hovered = nearestCell(planet, hitModel); + } + } else if (inMap) { + double nx = (mp.x - mapRect.x) / mapRect.width, ny = (mp.y - mapRect.y) / mapRect.height; + double X = (nx * 2.0 - 1.0) * EqualEarth::halfWidth(); + double Y = (1.0 - 2.0 * ny) * EqualEarth::halfHeight(); + double lon, lat; + if (EqualEarth::inverse(X, Y, lon, lat)) + hovered = nearestCell(planet, lonLatToDir(wrapPi(lon - mapLon), lat)); + } + if (have3DHit && !subgrids.empty()) { // prefer subcell under cursor + double best = -2.0; const SubCell* bs = nullptr; + for (auto& sg : subgrids) { + if (!sg) continue; + for (auto& s : sg->sub) { double dd = s.unit.dot(hitModel); if (dd > best) { best = dd; bs = &s; } } + } + if (bs && angBetween(bs->unit, hitModel) < selectedThresh) { hoverSub = *bs; hasHoverSub = true; } + } + + // --- Click = select a tile (ignored over panel/button / while dragging) - + if (IsMouseButtonReleased(MOUSE_BUTTON_LEFT) && dragDist < 6.0f && !inPanel && !onPause && !phase3Prompt && hovered >= 0) + selectCell(hovered); + + // --- Keys ------------------------------------------------------------- + if (IsKeyPressed(KEY_SPACE) && !phase3Prompt) pauseAction(); + if (IsKeyPressed(KEY_ONE)) { mode = ColorMode::Elevation; recolor(); } + if (IsKeyPressed(KEY_TWO)) { mode = ColorMode::Plate; recolor(); } + if (IsKeyPressed(KEY_THREE)) { mode = ColorMode::Age; recolor(); } + if (IsKeyPressed(KEY_FOUR)) { mode = ColorMode::Crust; recolor(); } + if (IsKeyPressed(KEY_FIVE)) { mode = ColorMode::Biome; recolor(); } + if (IsKeyPressed(KEY_SIX)) { mode = ColorMode::Temperature; recolor(); } + if (IsKeyPressed(KEY_SEVEN)) { mode = ColorMode::Precip; recolor(); } + if (IsKeyPressed(KEY_B)) showBorders = !showBorders; + if (IsKeyPressed(KEY_D)) showDrift = !showDrift; + if (IsKeyPressed(KEY_G)) showGrat = !showGrat; + if (IsKeyPressed(KEY_J)) showRivers = !showRivers; + if (IsKeyPressed(KEY_H) && settled) { // toggle Phase 3 (hydrology) + bool wasPrompt = phase3Prompt; + phase3 = !phase3; phase3Prompt = false; + if (wasPrompt) paused = false; // taking the choice resumes the sim + if (phase3) { phase3PromptAt = elapsedMy; setStatus("Hydrology ON"); } + else { phase3PromptAt = elapsedMy + planet.cfg.phase3AfterMy; rivers.clear(); bigRivers.clear(); setStatus("Hydrology OFF"); } + refreshView(); + } + if (IsKeyPressed(KEY_C)) { selectedCell = -1; subgrids.clear(); } + if (IsKeyPressed(KEY_R)) { cfg.seed = (uint32_t)(GetTime() * 100000) | 1; regen(); } + if (IsKeyPressed(KEY_S)) { stepOnce(); refreshView(); } // one tick (handy while paused/settled) + if (IsKeyPressed(KEY_F)) { // fast-forward to settled + if (!settled) { + while (!settled) stepOnce(); + dtMy = planet.cflDtMy(); planet.drifting = true; refreshView(); + } + } + if (IsKeyPressed(KEY_EQUAL) && cfg.subdivisions < 7) { cfg.subdivisions++; regen(); } + if (IsKeyPressed(KEY_MINUS) && cfg.subdivisions > 1) { cfg.subdivisions--; regen(); } + if (IsKeyPressed(KEY_F2)) { + if (loadConfig(CONFIG_PATH, cfg)) { + std::string cerr = validateConfig(cfg); + if (!cerr.empty()) { cfg = PlanetConfig{}; setStatus("Bad planet.cfg — using defaults"); } + regen(); + setStatus(cerr.empty() ? "Reloaded planet.cfg" : "Bad planet.cfg — using defaults"); + } else setStatus("No planet.cfg"); } + if (IsKeyPressed(KEY_F5)) saveGame(SAVE_PATH); + if (IsKeyPressed(KEY_F9)) loadGame(SAVE_PATH); + if (IsKeyPressed(KEY_F12)) { TakeScreenshot("screenshot.png"); setStatus("Screenshot saved to screenshot.png"); } + // Drift speed (Phase 2): My simulated per real second. + if (IsKeyPressed(KEY_RIGHT_BRACKET)) driftRate = std::min(driftRate * 1.5, 80.0); + if (IsKeyPressed(KEY_LEFT_BRACKET)) driftRate = std::max(driftRate / 1.5, 0.5); +} diff --git a/src/render/ViewerRender.cpp b/src/render/ViewerRender.cpp new file mode 100644 index 0000000..bce9a9f --- /dev/null +++ b/src/render/ViewerRender.cpp @@ -0,0 +1,262 @@ +#include "Viewer.hpp" +#include "Overlays.hpp" +#include "Map2D.hpp" +#include "Panels.hpp" +#include "Picking.hpp" // rotateZ (axial-tilt transform for labels) +#include "rlgl.h" +#include "Projection.hpp" // dirToLonLat (plate labels) +#include +#include +#include + +// Render the 3D globe into its own RenderTexture (its viewport != the screen). +void Viewer::renderGlobe3D() { + BeginTextureMode(rt3d); + ClearBackground(Color{8, 10, 16, 255}); + BeginMode3D(cam); + // Axial tilt: lean the whole globe (and everything drawn over it) by the + // obliquity about the world Z axis. Picking + plate labels rotate to match + // (see ViewerInput / renderFrame). The picking sphere is rotation-invariant. + rlPushMatrix(); + rlRotatef((float)planet.cfg.axialTilt, 0.0f, 0.0f, 1.0f); + const std::vector& tri = planet.triIndices(); + rlBegin(RL_TRIANGLES); + for (size_t k = 0; k + 2 < tri.size(); k += 3) { + int idx[3] = { tri[k], tri[k + 1], tri[k + 2] }; + for (int j = 0; j < 3; ++j) { + const Cell& cc = planet.cells[idx[j]]; + const Vec3& u = cc.unit; + float r = visBase + (float)cc.elevation * elevExagg; + const Color& col = vcolors[idx[j]]; + rlColor4ub(col.r, col.g, col.b, 255); + rlVertex3f((float)(u.x * r), (float)(u.y * r), (float)(u.z * r)); + } + } + rlEnd(); + + if (selectedCell >= 0) // detail patch + drawSubgrids(subgrids, visBase, elevExagg, planet.cfg.seaLevel, 0.0016f); + + if (showBorders && (!borders.empty() || !ridgeBorders.empty())) { + rlSetLineWidth(2.0f); rlBegin(RL_LINES); + rlColor4ub(255, 235, 90, 255); // real plate borders: yellow + for (size_t i = 0; i + 1 < borders.size(); i += 2) { + rlVertex3f(borders[i].x, borders[i].y, borders[i].z); + rlVertex3f(borders[i + 1].x, borders[i + 1].y, borders[i + 1].z); + } + rlColor4ub(220, 70, 60, 255); // young spreading ridges: red + for (size_t i = 0; i + 1 < ridgeBorders.size(); i += 2) { + rlVertex3f(ridgeBorders[i].x, ridgeBorders[i].y, ridgeBorders[i].z); + rlVertex3f(ridgeBorders[i + 1].x, ridgeBorders[i + 1].y, ridgeBorders[i + 1].z); + } + rlEnd(); rlSetLineWidth(1.0f); + } + if (showDrift && !driftArrows.empty()) { + rlSetLineWidth(2.5f); rlBegin(RL_LINES); rlColor4ub(90, 230, 255, 255); + for (size_t i = 0; i + 1 < driftArrows.size(); i += 2) { + rlVertex3f(driftArrows[i].x, driftArrows[i].y, driftArrows[i].z); + rlVertex3f(driftArrows[i + 1].x, driftArrows[i + 1].y, driftArrows[i + 1].z); + } + rlEnd(); rlSetLineWidth(1.0f); + } + if (phase3 && showRivers) { // Phase-3 river network + auto drawRiv = [&](const std::vector& segs, float w) { + if (segs.empty()) return; + rlSetLineWidth(w); rlBegin(RL_LINES); rlColor4ub(80, 170, 235, 255); + for (size_t i = 0; i + 1 < segs.size(); i += 2) { + rlVertex3f(segs[i].x, segs[i].y, segs[i].z); + rlVertex3f(segs[i + 1].x, segs[i + 1].y, segs[i + 1].z); + } + rlEnd(); rlSetLineWidth(1.0f); + }; + drawRiv(rivers, 1.5f); drawRiv(bigRivers, 3.0f); + } + if (showGrat) drawGraticule3D(graticule, gratR); + // Markers: selected (orange), hovered cell (yellow), hovered subcell (white). + if (selectedCell >= 0) { + Vec3 u = planet.cells[selectedCell].unit * (double)(visBase + 0.012f); + DrawSphere(Vector3{(float)u.x, (float)u.y, (float)u.z}, 0.03f, ORANGE); + } + if (hovered >= 0 && !hasHoverSub) { + Vec3 u = planet.cells[hovered].unit * (double)(visBase + 0.012f); + DrawSphere(Vector3{(float)u.x, (float)u.y, (float)u.z}, 0.022f, YELLOW); + } + if (hasHoverSub) { + Vec3 u = hoverSub.unit * (double)(visBase + 0.02f); + DrawSphere(Vector3{(float)u.x, (float)u.y, (float)u.z}, 0.012f, WHITE); + } + // Spin axis: a rod through the poles, extended beyond the surface (tilts with + // the globe since it's inside the rotated matrix). Pole caps mark N (red)/S (blue). + { + float ax = visBase + 0.6f; + rlSetLineWidth(2.5f); + rlBegin(RL_LINES); rlColor4ub(210, 220, 235, 255); + rlVertex3f(0.0f, -ax, 0.0f); rlVertex3f(0.0f, ax, 0.0f); + rlEnd(); rlSetLineWidth(1.0f); + DrawSphere(Vector3{0.0f, ax, 0.0f}, 0.05f, Color{230, 90, 80, 255}); // north + DrawSphere(Vector3{0.0f, -ax, 0.0f}, 0.05f, Color{80, 140, 230, 255}); // south + } + rlPopMatrix(); + EndMode3D(); + EndTextureMode(); +} + +// 2D Equal Earth map + its overlays (borders/drift/rivers/labels/markers). +void Viewer::renderMap2D() { + DrawRectangleRec(mapRect, Color{6, 8, 14, 255}); + BeginScissorMode((int)mapRect.x, (int)mapRect.y, (int)mapRect.width, (int)mapRect.height); + drawMap2D(planet, vcolors, map2D, mapRect, mapLon); + if (showGrat) { drawGraticule2D(graticule, mapRect, mapLon); drawGraticuleLabels2D(mapRect, mapLon); } + if (showBorders && !borders.empty()) drawSegments2D(borders, Color{255, 235, 90, 255}, 2.0f, mapRect, mapLon); + if (showBorders && !ridgeBorders.empty()) drawSegments2D(ridgeBorders, Color{220, 70, 60, 255}, 2.0f, mapRect, mapLon); + if (showDrift && !driftArrows.empty()) drawSegments2D(driftArrows, Color{90, 230, 255, 255}, 2.0f, mapRect, mapLon); + if (phase3 && showRivers) { + drawSegments2D(rivers, Color{80, 170, 235, 255}, 1.5f, mapRect, mapLon); + drawSegments2D(bigRivers, Color{80, 170, 235, 255}, 3.0f, mapRect, mapLon); + } + if (showDrift && !plateLabels.empty()) { + for (const auto& lbl : plateLabels) { + Vec3 u = Vec3{lbl.pos.x, lbl.pos.y, lbl.pos.z}.normalized(); + double lon, lat; dirToLonLat(u, lon, lat); + Vector2 lp = projLonLat(lon, lat, mapLon, mapRect); + const char* txt = TextFormat("P%d", lbl.id); + DrawText(txt, (int)lp.x + 4, (int)lp.y - 8, 12, RAYWHITE); + } + } + if (selectedCell >= 0) DrawCircleV(mapScreen(map2D, selectedCell, mapRect, mapLon), 5, ORANGE); + if (hovered >= 0) DrawCircleV(mapScreen(map2D, hovered, mapRect, mapLon), 4, YELLOW); + EndScissorMode(); + DrawRectangleLinesEx(mapRect, 1, Color{90, 90, 110, 255}); + DrawText("2D Equal Earth (hover, drag to pan)", (int)mapRect.x + 6, (int)mapRect.y + 4, 14, Color{200, 200, 210, 255}); +} + +// Right column: hover/selection info (top) + detail panel or world stats (bottom). +void Viewer::renderPanels() { + drawHoverPanel(planet, hoverRect, hovered, selectedCell); + if (selectedCell >= 0 && !subgrids.empty()) + drawDetailPanel(planet, subgrids[0], selectedCell, + planet.cells[selectedCell].elevation, planet.cells[selectedCell].geoAge, + panelRect, gridRect, hoveredSubIdx); + else + drawStats(planet, panelRect, elapsedMy, settled); +} + +// Top-left HUD text + the clickable pause button. +void Viewer::renderHUD() { + // Active view-mode label, centered at the top of the 3D viewport. + { + const char* vm = TextFormat("%s view", colorModeName(mode)); + int vw = MeasureText(vm, 22); + DrawText(vm, view3DW / 2 - vw / 2, 10, 22, Color{235, 225, 140, 255}); + } + + int y = 10; + auto line = [&](const std::string& s){ DrawText(s.c_str(), 12, y, 18, RAYWHITE); y += 22; }; + double fastest = 0.0; for (const auto& pl : planet.plates) fastest = std::max(fastest, pl.speedCmYr); + line(!settled ? "Planet Sim - World Creation: forming" + : phase3 ? "Planet Sim - World Creation: hydrology" + : "Planet Sim - World Creation: drift & erosion"); + line(TextFormat("Cells: %d Subdiv: %d CellWidth: %.0f km", + (int)planet.cells.size(), cfg.subdivisions, planet.cellWidthMeters() / 1000.0)); + line(TextFormat("Elevation: %.0f .. %.0f m", minE, maxE)); + if (!settled) + line(TextFormat("Forming terrain tick %lld max change %.1f m/tick%s", + stepCount, maxChange, paused ? " [PAUSED]" : "")); + else { + line(TextFormat("%s %.1f My elapsed %.1f My/s%s", + phase3 ? "Hydrology - drift, rivers & erosion" : "Drift & erosion", + elapsedMy, driftRate, paused ? " [PAUSED]" : "")); + line(TextFormat("dt %.2f My/step fastest plate %.1f cm/yr [ / ] speed", dtMy, fastest)); + if (phase3) { + int riverCells = 0, lakeCells = 0; const auto& dq = planet.discharge(); const auto& lk = planet.lakeDepth(); + for (size_t i = 0; i < planet.cells.size(); ++i) { + if (!dq.empty() && dq[i] > planet.cfg.riverThreshold) ++riverCells; + if (!lk.empty() && lk[i] > 20.0 && planet.cells[i].elevation > planet.cfg.seaLevel) ++lakeCells; + } + line(TextFormat("rivers: %d cells lakes: %d cells", riverCells, lakeCells)); + } + } + y += 8; + line("hover: cell info | click tile: open detail panel | C close"); + line("1 elev 2 plates 3 age 4 crust 5 biome 6 temp 7 precip"); + line(TextFormat("B borders [%s] | D vectors [%s] | G grid [%s] | J rivers [%s]", + showBorders ? "on" : "off", showDrift ? "on" : "off", showGrat ? "on" : "off", showRivers ? "on" : "off")); + line(TextFormat("SPACE pause | [ / ] speed | S step | F fast-fwd | H hydrology [%s] | R reseed | +/-", + phase3 ? "on" : "off")); + line("F5 save | F9 load | F12 screenshot | F2 reload planet.cfg"); + if (!statusMsg.empty() && GetTime() < statusUntil) { + y += 4; DrawText(statusMsg.c_str(), 12, y, 18, Color{120, 230, 140, 255}); y += 22; + } + + // Clickable pause button (bottom-left of the 3D quadrant). + DrawRectangleRec(pauseBtn, onPause ? Color{60, 70, 92, 255} : Color{28, 34, 46, 235}); + DrawRectangleLinesEx(pauseBtn, 1, Color{120, 120, 150, 255}); + const char* plbl = paused ? "> RESUME" : "|| PAUSE"; + Color plcol = paused ? Color{120, 230, 140, 255} : RAYWHITE; + int plw = MeasureText(plbl, 18); + DrawText(plbl, (int)(pauseBtn.x + (pauseBtn.width - plw) / 2), (int)pauseBtn.y + 7, 18, plcol); +} + +// Phase-3 transition prompt (modal overlay over the 3D viewport). +void Viewer::renderPrompt() { + if (!phase3Prompt) return; + DrawRectangle(0, 0, (int)view3DW, (int)view3DH, Color{0, 0, 0, 150}); + const char* q = TextFormat("Reached %.0f My of drift. Begin hydrology (rivers, lakes & erosion)?", elapsedMy); + int qw = MeasureText(q, 22); + DrawText(q, (int)(pbCx - qw / 2.0f), (int)(pbCy - 40.0f), 22, RAYWHITE); + auto drawBtn = [&](Rectangle b, const char* lbl, Color fill) { + bool hot = CheckCollisionPointRec(mp, b); + DrawRectangleRec(b, hot ? Color{70, 90, 120, 255} : fill); + DrawRectangleLinesEx(b, 1, Color{150, 150, 180, 255}); + int w = MeasureText(lbl, 18); + DrawText(lbl, (int)(b.x + (b.width - w) / 2.0f), (int)(b.y + 11.0f), 18, RAYWHITE); + }; + drawBtn(p3ContinueBtn, "Keep building", Color{40, 46, 60, 255}); + drawBtn(p3StartBtn, "Start hydrology", Color{30, 72, 60, 255}); +} + +// One full frame: globe texture, then composite + 3D labels + map + panels + +// HUD + prompt onto the screen. +void Viewer::renderFrame() { + renderGlobe3D(); + + BeginDrawing(); + ClearBackground(Color{8, 10, 16, 255}); + DrawTextureRec(rt3d.texture, Rectangle{0, 0, (float)view3DW, -(float)view3DH}, + Vector2{0, 0}, WHITE); + + // 3D plate labels (manually projected to match BeginMode3D's viewport exactly). + if (showDrift && !plateLabels.empty()) { + Vec3 camPos{cam.position.x, cam.position.y, cam.position.z}; + Vec3 camTgt{cam.target.x, cam.target.y, cam.target.z}; + Vec3 camUp {cam.up.x, cam.up.y, cam.up.z}; + Vec3 forward = (camTgt - camPos).normalized(); + Vec3 right = forward.cross(camUp).normalized(); + Vec3 up = right.cross(forward); + + double fovRad = cam.fovy * M_PI / 180.0; + double aspect = (double)view3DW / view3DH; + double projH = std::tan(fovRad * 0.5); // half-height of the view frustum (NDC) + double projW = projH * aspect; + + for (const auto& lbl : plateLabels) { + Vec3 lp = rotateZ(Vec3{lbl.pos.x, lbl.pos.y, lbl.pos.z}, planet.cfg.axialTilt); // tilt to match globe + if (lp.dot(camPos) <= 0.0) continue; // far hemisphere -> hidden by globe + Vec3 rel = lp - camPos; + double z = rel.dot(forward); + if (z <= 0.0) continue; + double xndc = rel.dot(right) / (projW * z); + double yndc = rel.dot(up) / (projH * z); + float sx = (float)((xndc * 0.5 + 0.5) * view3DW); + float sy = (float)((0.5 - yndc * 0.5) * view3DH); + DrawText(TextFormat("P%d", lbl.id), (int)sx + 6, (int)sy - 6, 16, RAYWHITE); + } + } + + renderMap2D(); + renderPanels(); + renderHUD(); + renderPrompt(); + + EndDrawing(); +} diff --git a/src/sim/IcoSphere.cpp b/src/sim/IcoSphere.cpp new file mode 100644 index 0000000..3fc4c64 --- /dev/null +++ b/src/sim/IcoSphere.cpp @@ -0,0 +1,70 @@ +#include "IcoSphere.hpp" +#include + +int IcoSphere::addMidpoint(int a, int b) { + int64_t lo = a < b ? a : b; + int64_t hi = a < b ? b : a; + int64_t key = (lo << 32) | hi; + auto it = midpointCache.find(key); + if (it != midpointCache.end()) return it->second; + + Vec3 mid = ((positions[a] + positions[b]) * 0.5).normalized(); + int idx = (int)positions.size(); + positions.push_back(mid); + midpointCache.emplace(key, idx); + return idx; +} + +void IcoSphere::build(int subdivisions) { + positions.clear(); + neighbors.clear(); + triIndices.clear(); + midpointCache.clear(); + + // Base icosahedron (golden-ratio rectangle vertices). + const double t = (1.0 + std::sqrt(5.0)) / 2.0; + std::vector base = { + {-1, t, 0}, { 1, t, 0}, {-1, -t, 0}, { 1, -t, 0}, + { 0, -1, t}, { 0, 1, t}, { 0, -1, -t}, { 0, 1, -t}, + { t, 0, -1}, { t, 0, 1}, {-t, 0, -1}, {-t, 0, 1} + }; + for (auto& v : base) positions.push_back(v.normalized()); + + std::vector faces = { + 0,11,5, 0,5,1, 0,1,7, 0,7,10, 0,10,11, + 1,5,9, 5,11,4, 11,10,2, 10,7,6, 7,1,8, + 3,9,4, 3,4,2, 3,2,6, 3,6,8, 3,8,9, + 4,9,5, 2,4,11, 6,2,10, 8,6,7, 9,8,1 + }; + + // Subdivide. + for (int s = 0; s < subdivisions; ++s) { + std::vector next; + next.reserve(faces.size() * 4); + for (size_t i = 0; i < faces.size(); i += 3) { + int a = faces[i], b = faces[i + 1], c = faces[i + 2]; + int ab = addMidpoint(a, b); + int bc = addMidpoint(b, c); + int ca = addMidpoint(c, a); + next.insert(next.end(), {a, ab, ca, + b, bc, ab, + c, ca, bc, + ab, bc, ca}); + } + faces.swap(next); + } + + triIndices = faces; + + // Build neighbor adjacency from triangle edges. + std::vector> adj(positions.size()); + for (size_t i = 0; i < faces.size(); i += 3) { + int a = faces[i], b = faces[i + 1], c = faces[i + 2]; + adj[a].insert(b); adj[a].insert(c); + adj[b].insert(a); adj[b].insert(c); + adj[c].insert(a); adj[c].insert(b); + } + neighbors.resize(positions.size()); + for (size_t i = 0; i < positions.size(); ++i) + neighbors[i].assign(adj[i].begin(), adj[i].end()); +} diff --git a/src/sim/IcoSphere.hpp b/src/sim/IcoSphere.hpp new file mode 100644 index 0000000..4f76c7c --- /dev/null +++ b/src/sim/IcoSphere.hpp @@ -0,0 +1,25 @@ +#pragma once +#include "Vec3.hpp" +#include +#include +#include + +// Generates a geodesic icosphere: subdivided icosahedron projected onto a +// unit sphere. Produces fixed vertices (cell centers live on faces in dual, +// but here we use the triangle vertices as cells for simplicity) and a +// per-vertex neighbor list. Geometry is fixed; properties flow over it. +struct IcoSphere { + std::vector positions; // unit-sphere positions per vertex + std::vector> neighbors; // adjacency per vertex + std::vector triIndices; // flat triangle index list (3 per tri) + + // Build an icosphere with the given subdivision level (0 = base icosahedron). + void build(int subdivisions); + +private: + int addMidpoint(int a, int b); + // Cache of shared edge midpoints, keyed by the packed (lo<<32 | hi) vertex + // pair, so each edge is split exactly once. Hash map -> O(1) lookup; a + // linear scan here made build() O(n^2) and froze high subdivisions. + std::unordered_map midpointCache; +}; diff --git a/src/sim/Planet.cpp b/src/sim/Planet.cpp new file mode 100644 index 0000000..b71112e --- /dev/null +++ b/src/sim/Planet.cpp @@ -0,0 +1,233 @@ +#include "Planet.hpp" +#include +#include + +// --- Core: generation, geometry, plate seeding, shared helpers, subgrid ------ +// The Planet class is implemented across several translation units (all sharing +// this one header): tectonics in PlanetTectonics.cpp, drift/plate-lifecycle in +// PlanetDrift.cpp, erosion in PlanetErosion.cpp, hydrology in PlanetHydrology.cpp, +// config/save I/O in PlanetIO.cpp. This file holds world generation plus the +// small helpers (RNG, drift velocity, plate speed) the others call. + +// xorshift32 -- deterministic, seedable. +uint32_t Planet::rnd() { + uint32_t x = rngState; + x ^= x << 13; x ^= x >> 17; x ^= x << 5; + rngState = x; + return x; +} +double Planet::rndf() { return (rnd() & 0xFFFFFF) / double(0x1000000); } + +void Planet::generate(const PlanetConfig& c) { + cfg = c; + rngState = c.seed ? c.seed : 1; + targetLand = -1; // recomputed at first advect (drift start) + driftIter = 0; + erodeIter = 0; + sPrevCount.clear(); + sStaleStreak.clear(); + sFreePlateIds.clear(); + + buildGeometry(); + for (auto& c : cells) { c.elevation = 0.0; c.plateId = -1; c.geoAge = 0.0; } + + assignPlates(); + seedInitialRelief(); + computeClimate(); // temperature + precipitation fields (biomes read these) + classifyBiomes(); // give the fresh world an initial biome per cell +} + +// Build the icosphere and copy fixed geometry (unit direction + neighbor +// adjacency) onto the cells. Shared by generate() and readState() (load). +void Planet::buildGeometry() { + sphere.build(cfg.subdivisions); + cells.clear(); + cells.resize(sphere.positions.size()); + for (size_t i = 0; i < cells.size(); ++i) { + cells[i].unit = sphere.positions[i]; + cells[i].neighbors = sphere.neighbors[i]; + } +} + +void Planet::assignPlates() { + plates.clear(); + plates.resize(cfg.plateCount); + + // Pick random seed cells, flood-fill plate ownership over neighbors. + std::vector frontier; + for (int p = 0; p < cfg.plateCount; ++p) { + int seed = rnd() % cells.size(); + cells[seed].plateId = p; + frontier.push_back(seed); + + plates[p].id = p; + plates[p].type = (rndf() < 0.6) ? PlateType::Oceanic + : PlateType::Continental; + plates[p].baby = false; + // Random rotation axis + real surface speed 1..maxDriftSpeed cm/yr. + randomizePlateDrift(plates[p]); + } + + // Multi-source BFS so plates grow at equal rate. + size_t head = 0; + while (head < frontier.size()) { + int cur = frontier[head++]; + int pid = cells[cur].plateId; + for (int nb : cells[cur].neighbors) { + if (cells[nb].plateId < 0) { + cells[nb].plateId = pid; + frontier.push_back(nb); + } + } + } +} + +void Planet::seedInitialRelief() { + // Continental plates sit higher; oceanic lower. Add mild noise. The bases + // double as the isostatic equilibrium each cell relaxes toward in step(). + for (auto& cell : cells) { + const Plate& pl = plates[cell.plateId]; + cell.oceanic = (pl.type == PlateType::Oceanic); // crust type now lives on the cell + // Seed an age spread on oceanic crust so the starting seafloor already has + // ridge->abyss variety (continental crust has no cooling-age depth). + cell.geoAge = cell.oceanic ? rndf() * cfg.seafloorSeedAge : 0.0; + double base = cell.oceanic ? oceanicBase(cell.geoAge) : cfg.continentBase; + double noise = (rndf() * 2 - 1) * 200.0; + cell.elevation = base + noise; + } +} + +// Velocity of a plate's material at position pos (tangential to sphere). +// v = omega x r, with omega = axis * speed. +Vec3 Planet::driftVelocity(int plateId, const Vec3& pos) const { + const Plate& pl = plates[plateId]; + Vec3 omega = pl.driftAxis * pl.driftSpeed; + return omega.cross(pos); +} + +// Seafloor subsidence (half-space cooling): oceanic crust deepens with age from +// the ridge toward a deep abyssal floor. Used as the relax target in step(). +double Planet::oceanicBase(double age) const { + return std::max(cfg.oceanBase, + cfg.ridgeDepth - cfg.seafloorSubsidence * std::sqrt(std::max(0.0, age))); +} + +// cm/yr -> the small driftSpeed the Phase-1 uplift stress uses + angSpeed (rad/My). +void Planet::setPlateSpeed(Plate& p, double cmYr) { + p.speedCmYr = cmYr; + p.driftSpeed = (cmYr / cfg.maxDriftSpeed) * 1e-3; + p.angSpeed = cmYr * 1.0e4 / cfg.radius; // cm/yr -> m/My -> rad/My +} + +// Random rotation axis + random surface speed 1..maxDriftSpeed cm/yr. +void Planet::randomizePlateDrift(Plate& p) { + Vec3 axis{rndf() * 2 - 1, rndf() * 2 - 1, rndf() * 2 - 1}; + p.driftAxis = axis.normalized(); + setPlateSpeed(p, 1.0 + rndf() * (cfg.maxDriftSpeed - 1.0)); +} + +// Get a plate slot: reuse a dead one (0 cells) if available, else append. Keeps +// the `plates` vector from growing without bound as rifts spawn baby plates. +int Planet::acquirePlate() { + if (!sFreePlateIds.empty()) { + int id = sFreePlateIds.back(); sFreePlateIds.pop_back(); + plates[id] = Plate{}; plates[id].id = id; + if (id < (int)sStaleStreak.size()) sStaleStreak[id] = 0; // reused slot: fresh streak + return id; + } + Plate np{}; np.id = (int)plates.size(); plates.push_back(np); // value-init: axis/speed zeroed + return np.id; +} + +// --- Subgrid generation ----------------------------------------------------- +namespace { + uint32_t hashU(uint32_t a) { a ^= a << 13; a ^= a >> 17; a ^= a << 5; return a; } + double latticeVal(int cell, int gx, int gy) { + uint32_t h = hashU((uint32_t)cell * 2654435761u ^ + hashU((uint32_t)(gx * 73856093) ^ (uint32_t)(gy * 19349663))); + return (h & 0xFFFFFF) / double(0x1000000) * 2.0 - 1.0; // [-1,1] + } + double smoothstep(double t) { return t * t * (3.0 - 2.0 * t); } + // Bilinear value noise on an integer lattice, smooth-interpolated. + double valueNoise(int cell, double fx, double fy) { + int x0 = (int)std::floor(fx), y0 = (int)std::floor(fy); + double tx = smoothstep(fx - x0), ty = smoothstep(fy - y0); + double v00 = latticeVal(cell, x0, y0), v10 = latticeVal(cell, x0 + 1, y0); + double v01 = latticeVal(cell, x0, y0 + 1), v11 = latticeVal(cell, x0 + 1, y0 + 1); + double a = v00 + (v10 - v00) * tx, b = v01 + (v11 - v01) * tx; + return a + (b - a) * ty; + } + double angDist(const Vec3& a, const Vec3& b) { + return std::acos(std::clamp(a.dot(b), -1.0, 1.0)); + } +} + +std::shared_ptr Planet::makeSubGrid(int cellIndex, int res) const { + auto sg = std::make_shared(); + sg->macroCell = cellIndex; + sg->res = res; + sg->sub.resize((size_t)res * res); + if (res < 2) return sg; + + const Cell& c = cells[cellIndex]; + const Vec3& n = c.unit; + + // Local tangent frame at the cell center. + Vec3 ref = (std::fabs(n.y) < 0.99) ? Vec3{0, 1, 0} : Vec3{1, 0, 0}; + Vec3 t = ref.cross(n).normalized(); + Vec3 b = n.cross(t).normalized(); + + // Patch reaches out to roughly the neighbor-cell centers. + double meanAng = 0.0; + for (int nb : c.neighbors) meanAng += angDist(n, cells[nb].unit); + double half = (c.neighbors.empty() ? 0.1 : meanAng / c.neighbors.size()); + + // Macro set whose elevations the patch blends: this cell + its neighbors. + std::vector macro; macro.reserve(c.neighbors.size() + 1); + macro.push_back(cellIndex); + for (int nb : c.neighbors) macro.push_back(nb); + + const double eps = (0.15 * half) * (0.15 * half) + 1e-9; // IDW smoothing + for (int j = 0; j < res; ++j) { + for (int i = 0; i < res; ++i) { + double u = ((double)i / (res - 1) * 2.0 - 1.0) * half; + double v = ((double)j / (res - 1) * 2.0 - 1.0) * half; + double r = std::sqrt(u * u + v * v); + Vec3 dir = (r < 1e-12) ? n + : (n * std::cos(r) + (t * (u / r) + b * (v / r)) * std::sin(r)).normalized(); + + // Inverse-distance-weighted blend of macro elevations. + double wsum = 0.0, esum = 0.0; int nearest = macro[0]; double best = 1e9; + for (int m : macro) { + double d = angDist(dir, cells[m].unit); + double w = 1.0 / (d * d + eps); + wsum += w; esum += w * cells[m].elevation; + if (d < best) { best = d; nearest = m; } + } + double elev = esum / wsum; + + // Fine sub-cell detail (two octaves of value noise, +/-~250 m). + double fx = (double)i / (res - 1), fy = (double)j / (res - 1); + double nz = valueNoise(cellIndex, fx * 4.0, fy * 4.0) + + valueNoise(cellIndex, fx * 8.0 + 11.3, fy * 8.0 + 7.7) * 0.5; + elev += (nz / 1.5) * 250.0; + + SubCell& s = sg->sub[(size_t)j * res + i]; + s.unit = dir; + s.elevation = elev; + s.nearestMacro = nearest; + } + } + return sg; +} + +double Planet::cellWidthMeters() const { + double area = 4.0 * M_PI * cfg.radius * cfg.radius; + return std::sqrt(area / std::max(1, cells.size())); +} +double Planet::minElevation() const { + double m = 1e30; for (auto& c : cells) m = std::min(m, c.elevation); return m; +} +double Planet::maxElevation() const { + double m = -1e30; for (auto& c : cells) m = std::max(m, c.elevation); return m; +} diff --git a/src/sim/Planet.hpp b/src/sim/Planet.hpp new file mode 100644 index 0000000..2d45be1 --- /dev/null +++ b/src/sim/Planet.hpp @@ -0,0 +1,139 @@ +#pragma once +#include "Vec3.hpp" +#include "IcoSphere.hpp" +#include "PlanetTypes.hpp" // Cell, Plate, SubGrid/SubCell, PlanetConfig +#include +#include +#include +#include +#include + +class Planet { +public: + PlanetConfig cfg; + std::vector cells; + std::vector plates; + + // Phase flag: false during Phase-1 forming (modest, original tectonics that + // settle), true during Phase-2 drift. Gates the increment-4 orogeny boosts + // (collision/arc uplift + isostatic persistence) so Phase 1 stays unchanged + // and tall persistent mountains only grow during drift (where erode() limits + // them). Not serialized -- the orchestrator sets it from the phase. + bool drifting = false; + + void generate(const PlanetConfig& c); + double step(); // one tectonic tick; returns max |elevation change| (m) this tick + + // Phase 2: stable timestep (My) from the fastest plate (CFL ~ half a cell). + double cflDtMy() const; + // Phase 2: advect plate membership + carried crust (plateId, elevation, + // oceanic, geoAge) over the fixed grid by dt My. Opening gaps become new + // young oceanic crust (spreading); overlaps subduct. + void advect(double dtMy); + + // Phase 2: erode the elevation field by dt My (highs wear down, sediment + // deposits downhill in basins / below sea level) and, every seaLevelEvery + // calls, nudge seaLevel toward landFractionTarget. Crust type is untouched. + void erode(double dtMy); + + // Phase 3: one hydrology tick over the fixed grid -- recompute the drainage + // network (depression-fill -> lakes, steepest-descent routing -> rivers, + // flow accumulation -> discharge) and apply mass-conserving fluvial erosion + // (stream-power incision + downstream sediment transport/deposition) to the + // elevation field by dt My. computeHydrology() does the routing only (no + // erosion) so the viewer can show rivers/lakes when paused / after load. + void hydrology(double dtMy); + void computeHydrology(); + + // Phase 3 (climate): compute per-cell temperature + precipitation fields from + // elevation, latitude and prevailing-wind orographic moisture transport (windward + // rain, leeward rain shadow, dry continental interiors). Derived (not saved); + // call before classifyBiomes(), which consumes these fields. + void computeClimate(); + const std::vector& temperature() const { return sTemp; } // deg C + const std::vector& precipitation() const { return sPrecip; } // relative units + const std::vector& moisture() const { return sMoist; } // 0..1 (median land -> 0.5) + + // Phase 3 (biomes): classify every cell into a Biome from elevation + the climate + // fields (temperature + normalized precipitation). Derived + written back into + // cell.biome (saved). Assumes computeClimate() ran this tick. Re-run as terrain evolves. + void classifyBiomes(); + // Derived hydrology fields (recomputed each route; not saved). Empty until + // the first computeHydrology()/hydrology() call. + const std::vector& lakeDepth() const { return sLakeDepth; } + const std::vector& discharge() const { return sDischarge; } + const std::vector& flowTo() const { return sFlowTo; } + + // Build a fine-resolution subgrid patch for one macro cell (phase 4/5 hook). + std::shared_ptr makeSubGrid(int cellIndex, int res) const; + + // Save/load the full simulation state (binary). readState rebuilds geometry + // from the saved cfg.subdivisions, so only dynamic per-cell fields are stored. + // Reloading resumes the simulation exactly (deterministic continuation). + void writeState(std::ostream& os) const; + // hasBiome: whether the stream carries the per-cell biome byte (save v4+). For + // older saves (v3) pass false -- biomes are reclassified after the cells load. + bool readState(std::istream& is, bool hasBiome = true); + + // Helpers for rendering / info. + double cellWidthMeters() const; // approx lateral cell spacing + double minElevation() const; + double maxElevation() const; + const std::vector& triIndices() const { return sphere.triIndices; } + +private: + IcoSphere sphere; + uint32_t rngState = 1; + int targetLand = -1; // land-cell count to conserve during drift + + uint32_t rnd(); + double rndf(); // [0,1) + void buildGeometry(); // build icosphere + per-cell unit/neighbors + void assignPlates(); + void seedInitialRelief(); + Vec3 driftVelocity(int plateId, const Vec3& pos) const; + double oceanicBase(double age) const; // age-dependent seafloor depth + + // Phase-2 plate-dynamics helpers (see PlanetConfig above). + void setPlateSpeed(Plate& p, double cmYr); // cm/yr -> driftSpeed + angSpeed + void randomizePlateDrift(Plate& p); // random axis + random speed + int acquirePlate(); // reuse a dead slot or append one + void splitPlate(int pid); // fission: cut a plate roughly in two + void maybeSplitPlates(const std::vector& cnt); + void kickStalemates(const std::vector& cnt); + void deleteEnclosedPlates(); // absorb plates ringed by one other + void fuseMiniPlates(); // cluster of mini plates steals + merges + void coalesceBabyPlates(); // merge connected baby cells into one id + void promoteBabyPlates(const std::vector& cnt); + void adjustSeaLevel(); // nudge seaLevel toward land target + + // Phase-3 hydrology helpers (see hydrology()). + void routeFlow(); // depression-fill -> lakes, flow, discharge + + int driftIter = 0; // counts advect() calls (gates periodic checks) + int erodeIter = 0; // counts erode() calls (gates sea-level control) + std::vector sPrevCount; // per-plate cell count at the previous check + std::vector sStaleStreak; // consecutive stuck windows per plate + std::vector sFreePlateIds; // dead plate slots free for reuse + + // Reusable scratch buffers for step()/erode() so they allocate nothing per tick. + std::vector sStress, sBelt, sBeltNext, sDelta, sSmoothed, sOldElev, sErode; + std::vector sSub, sOver, sColl; + + // Phase-3 hydrology scratch (derived from elevation each routeFlow(); not saved). + std::vector sFill, sLakeDepth, sDischarge; + std::vector sFlowTo, sHydroOrder; + + // Phase-3 climate scratch (derived each computeClimate(); not saved). sMoist is the + // 0..1-normalized precipitation the biome classifier reads. + std::vector sTemp, sPrecip, sMoist; + std::vector sWind; + std::vector sUpwind; +}; + +// Human-editable config file (key = value text). All PlanetConfig input +// parameters are written/read via one shared field table. Unknown keys ignored. +// validateConfig returns an empty string if the config is reasonable. +bool loadConfig(const std::string& path, PlanetConfig& cfg); +bool saveConfig(const std::string& path, const PlanetConfig& cfg); +std::string validateConfig(const PlanetConfig& cfg); diff --git a/src/sim/PlanetBiomes.cpp b/src/sim/PlanetBiomes.cpp new file mode 100644 index 0000000..a55e2f2 --- /dev/null +++ b/src/sim/PlanetBiomes.cpp @@ -0,0 +1,55 @@ +#include "Planet.hpp" +#include +#include + +// --- Phase 3 (biomes): per-cell biome classification ------------------------- +// Rule-based, reading the climate fields (computeClimate(): temperature sTemp + +// normalized precipitation sMoist) plus elevation and hydrology lakes. All +// thresholds are tunable in planet.cfg (the PlanetConfig `biome*` fields). Writes +// cell.biome (saved). Assumes computeClimate() ran this tick. + +void Planet::classifyBiomes() { + const int n = (int)cells.size(); + const double sea = cfg.seaLevel; + if ((int)sTemp.size() != n || (int)sMoist.size() != n) computeClimate(); // safety: ensure fields + + // Tunable thresholds (planet.cfg). Local aliases keep the classifier readable. + const double ICE_TEMP = cfg.biomeIceTemp, TUNDRA_TEMP = cfg.biomeTundraTemp; + const double TAIGA_TEMP = cfg.biomeTaigaTemp, SAVANNA_TEMP = cfg.biomeSavannaTemp; + const double MOUNTAIN_ELEV = cfg.biomeMountainElev, HILLS_ELEV = cfg.biomeHillsElev; + const double BEACH_BAND = cfg.biomeBeachBand, LOWLAND_ELEV = cfg.biomeLowlandElev; + const double WETLAND_MOIST = cfg.biomeWetlandMoist, DESERT_MOIST = cfg.biomeDesertMoist; + const double GRASS_MOIST = cfg.biomeGrassMoist, TAIGA_MOIST = cfg.biomeTaigaMoist; + const double LAKE_MIN_DEPTH = cfg.biomeLakeMinDepth; + const bool haveLake = !sLakeDepth.empty(); + + for (int i = 0; i < n; ++i) { + const Cell& c = cells[i]; + double elevAbove = std::max(0.0, c.elevation - sea); + double temp = sTemp[i]; // climate temperature (deg C) + double moist = sMoist[i]; // climate precipitation, normalized 0..1 + double lakeD = haveLake ? sLakeDepth[i] : 0.0; + bool adjOcean = false, adjWater = false; + for (int nb : c.neighbors) { + if (cells[nb].elevation <= sea) { adjOcean = true; adjWater = true; } + else if (haveLake && sLakeDepth[nb] > LAKE_MIN_DEPTH) adjWater = true; // lakeside + } + + Biome b; + if (temp < ICE_TEMP) b = Biome::Ice; // polar caps + snowcaps + else if (c.elevation <= sea) b = Biome::Ocean; + else if (lakeD > LAKE_MIN_DEPTH) b = Biome::Lake; // inland water above sea + else if (elevAbove <= BEACH_BAND && adjOcean) b = Biome::Beach; + else if (elevAbove > MOUNTAIN_ELEV) b = Biome::Mountains; + else if (elevAbove > HILLS_ELEV) b = Biome::Hills; + else { // lowland / plains + if (temp < TUNDRA_TEMP) b = Biome::Tundra; + else if (elevAbove < LOWLAND_ELEV && moist > WETLAND_MOIST && adjWater) b = Biome::Wetland; // swamps hug water + else if (temp < TAIGA_TEMP) b = (moist > TAIGA_MOIST) ? Biome::Taiga : Biome::Tundra; + else if (moist < DESERT_MOIST) b = Biome::Desert; + else if (moist < GRASS_MOIST) b = (temp > SAVANNA_TEMP) ? Biome::Savanna : Biome::Grassland; + else b = Biome::Forest; + } + cells[i].biome = b; + } +} diff --git a/src/sim/PlanetClimate.cpp b/src/sim/PlanetClimate.cpp new file mode 100644 index 0000000..4db9a6f --- /dev/null +++ b/src/sim/PlanetClimate.cpp @@ -0,0 +1,134 @@ +#include "Planet.hpp" +#include +#include +#include + +// --- Phase 3: climate (temperature + orographic precipitation) --------------- +// Continuous per-cell fields derived from elevation, latitude and prevailing +// (zonal) winds. Raylib-free, deterministic, not saved (recomputed each tick like +// hydrology). classifyBiomes() consumes sTemp + the normalized sMoist. + +namespace { + // Latitudinal rainfall belts: a wet equatorial ITCZ, dry subtropics (~25-30 deg + // -> deserts), a moist mid-latitude storm track, dry poles. (Moved here from + // PlanetBiomes; modulates how readily moving air rains out.) + double rainfallByLat(double latRad) { + double d = std::fabs(latRad) * 180.0 / M_PI; + double eq = std::exp(-((d - 0.0) * (d - 0.0)) / (14.0 * 14.0)); + double mid = std::exp(-((d - 52.0) * (d - 52.0)) / (18.0 * 18.0)); + return 0.15 + 0.55 * eq + 0.35 * mid; + } +} + +void Planet::computeClimate() { + const int n = (int)cells.size(); + const double sea = cfg.seaLevel; + sTemp.assign(n, 0.0); + sPrecip.assign(n, 0.0); + sMoist.assign(n, 0.0); + sWind.assign(n, Vec3{0, 0, 0}); + sUpwind.assign(n, -1); + + const double EQ = cfg.biomeEquatorTemp, POLE = cfg.biomePoleDrop; + const double LATEXP = cfg.biomeLatExp, LAPSE = cfg.biomeElevLapse; + auto isWater = [&](int i) { return cells[i].elevation <= sea; }; + + // 1. Temperature, prevailing wind (zonal bands), and the upwind neighbour. + std::vector belt(n); + const Vec3 up{0, 1, 0}; + for (int i = 0; i < n; ++i) { + double lat = std::asin(std::clamp(cells[i].unit.y, -1.0, 1.0)); + double elevAbove = std::max(0.0, cells[i].elevation - sea); + sTemp[i] = EQ - POLE * std::pow(std::fabs(lat) / (M_PI / 2.0), LATEXP) - LAPSE * elevAbove; + belt[i] = rainfallByLat(lat); + + Vec3 east = up.cross(cells[i].unit); // tangent toward +longitude + double el = east.length(); + if (el < 1e-6) continue; // near a pole: no zonal wind + east = east * (1.0 / el); + double absLatDeg = std::fabs(lat) * 180.0 / M_PI; + double sign = (absLatDeg < 30.0 || absLatDeg > 60.0) ? -1.0 : 1.0; // easterlies blow west + sWind[i] = east * sign; + + // Upwind neighbour: the one the wind comes from (step nb->i aligns with wind). + int best = -1; double bestA = -1e30; + for (int nb : cells[i].neighbors) { + Vec3 dv = cells[i].unit - cells[nb].unit; + double dl = dv.length(); if (dl < 1e-12) continue; + double a = dv.dot(sWind[i]) / dl; + if (a > bestA) { bestA = a; best = nb; } + } + sUpwind[i] = best; + } + + // 2. Steady-state moisture advection along the wind (iterative upwind differencing, + // double-buffered -> deterministic). Ocean cells are a moisture source; land + // cells take their upwind moisture, rain part of it out (more on windward upslopes) + // and lose a little per cell, so leeward/interior cells dry out. + const double oceanM = cfg.climateOceanMoisture, eff = cfg.climateRainEfficiency; + const double oro = cfg.climateOrographic, refH = std::max(1.0, cfg.climateOroRefHeight); + const double cont = cfg.climateContinentality; + std::vector M(n), Mn(n); + for (int i = 0; i < n; ++i) M[i] = isWater(i) ? oceanM : 0.0; + int passes = std::max(1, cfg.climateWindPasses); + for (int p = 0; p < passes; ++p) { + for (int i = 0; i < n; ++i) { + if (isWater(i)) { Mn[i] = oceanM; continue; } + int u = sUpwind[i]; + if (u < 0) { Mn[i] = M[i]; continue; } + double Mup = M[u]; + double upslope = std::max(0.0, cells[i].elevation - cells[u].elevation) / refH; + double rain = std::min(Mup, Mup * belt[i] * (1.0 + oro * upslope) * eff); + // Continentality is multiplicative (a fractional loss per cell) so deep + // interiors get progressively, smoothly drier instead of flooring at 0 -- + // which preserves a usable moisture gradient inland. + Mn[i] = std::max(0.0, Mup - rain) * (1.0 - cont); + } + M.swap(Mn); + } + + // 3. Precipitation from the converged incoming moisture (same rain formula). + for (int i = 0; i < n; ++i) { + if (isWater(i)) { sPrecip[i] = belt[i] * oceanM; } + else { + int u = sUpwind[i]; + if (u < 0) { sPrecip[i] = belt[i] * 0.1; } + else { + double Mup = M[u]; + double upslope = std::max(0.0, cells[i].elevation - cells[u].elevation) / refH; + sPrecip[i] = Mup * belt[i] * (1.0 + oro * upslope) * eff; + } + } + } + // Diffuse the precipitation field a few passes (atmospheric mixing). Pure upwind + // advection from a constant ocean source gives a near-binary field (saturated where + // the wind reaches the sea, ~0 elsewhere); diffusion turns the sharp coast/interior + // boundary into a gradient, so grassland/forest transition zones appear. + int smooth = std::max(0, cfg.climateMoistureSmooth); + if (smooth > 0) { + std::vector tmp(n); + const double keep = 0.45; // self-weight; the rest is the neighbour mean + for (int p = 0; p < smooth; ++p) { + for (int i = 0; i < n; ++i) { + double s = 0.0; int c = 0; + for (int nb : cells[i].neighbors) { s += sPrecip[nb]; ++c; } + tmp[i] = keep * sPrecip[i] + (1.0 - keep) * (c ? s / c : sPrecip[i]); + } + sPrecip.swap(tmp); + } + } + + // Normalize to 0..1 for the biome classifier by anchoring the MEDIAN land precip + // to ~0.5 (robust to the extreme orographic peaks, which would otherwise crush all + // lowland to "dry" if we divided by the max). Land = elevation above sea level. + std::vector landP; + landP.reserve(n); + for (int i = 0; i < n; ++i) if (!isWater(i)) landP.push_back(sPrecip[i]); + double ref = 1e-6; + if (!landP.empty()) { + size_t mid = landP.size() / 2; + std::nth_element(landP.begin(), landP.begin() + mid, landP.end()); + ref = std::max(1e-6, landP[mid] / 0.5); // median -> 0.5 + } + for (int i = 0; i < n; ++i) sMoist[i] = std::clamp(sPrecip[i] / ref, 0.0, 1.0); +} diff --git a/src/sim/PlanetDrift.cpp b/src/sim/PlanetDrift.cpp new file mode 100644 index 0000000..881c487 --- /dev/null +++ b/src/sim/PlanetDrift.cpp @@ -0,0 +1,347 @@ +#include "Planet.hpp" +#include +#include +#include + +// --- Phase 2: plate drift + plate lifecycle (Wilson cycle) ------------------- +// advect() moves plate membership + carried crust over the fixed grid, and every +// splitCheckEvery iterations runs the lifecycle (fission, stalemate kick, +// spreading "baby" plates, merging) so the plate count stays lively. + +double Planet::cflDtMy() const { + double maxAng = 1e-30; + for (const auto& p : plates) maxAng = std::max(maxAng, p.angSpeed); + double cellAngle = std::sqrt(4.0 * M_PI / std::max(1, cells.size())); // rad + return 0.5 * cellAngle / maxAng; // fastest plate advances ~half a cell / step +} + +// Accumulation advection (coherent) + an explicit land-conservation clamp. +// A boundary cell builds up `drift` = signed convergence distance with its +// dominant other-plate neighbor (+ closing, - rifting). +1 cell -> overrun +// (take that crust; oceanic subducts under a buoyant continent), -1 cell -> new +// ridge crust (spreading). A coherence cleanup keeps plates contiguous, and a +// conservation pass holds the total land-cell count at its drift-start value by +// rifting excess continental margins / accreting margins where land is short -- +// a pragmatic stand-in for true mass-conserving advection. +void Planet::advect(double dtMy) { + const int n = (int)cells.size(); + const double R = cfg.radius; + const double cellSize = cellWidthMeters(); + auto landCount = [&]{ int l = 0; for (auto& c : cells) if (!c.oceanic) ++l; return l; }; + if (targetLand < 0) targetLand = landCount(); + + std::vector oP(n); std::vector oE(n), oA(n); std::vector oO(n); + std::vector vel(n); + for (int i = 0; i < n; ++i) { + oP[i] = cells[i].plateId; oE[i] = cells[i].elevation; oA[i] = cells[i].geoAge; oO[i] = cells[i].oceanic; + const Plate& p = plates[cells[i].plateId]; + vel[i] = (p.driftAxis * p.angSpeed).cross(cells[i].unit) * R; + } + + for (int c = 0; c < n; ++c) { + cells[c].geoAge += dtMy; + double maxClose = -1e30; int invNb = -1; + for (int j : cells[c].neighbors) { + if (oP[j] == oP[c]) continue; + Vec3 dir = (cells[j].unit - cells[c].unit).normalized(); + double closing = (vel[c] - vel[j]).dot(dir); + if (closing > maxClose) { maxClose = closing; invNb = j; } + } + if (invNb < 0) continue; + if (cells[c].invader != oP[invNb]) { cells[c].drift = 0.0; cells[c].invader = oP[invNb]; } + cells[c].drift += maxClose * dtMy; + + if (cells[c].drift >= cellSize) { + bool invBaby = plates[oP[invNb]].baby; // young ridge crust can't override + if (invBaby || (!cells[c].oceanic && oO[invNb])) { // ocean subducts under continent + cells[c].drift -= cellSize; + } else { // overrun: take encroaching crust + cells[c].plateId = oP[invNb]; + cells[c].elevation = oE[invNb]; + cells[c].oceanic = oO[invNb]; + cells[c].geoAge = oA[invNb] + dtMy; + cells[c].drift -= cellSize; + } + cells[c].invader = -1; + } else if (cells[c].drift <= -cellSize) { // rift -> young crust on a baby plate + // The gap opened by divergence becomes young ocean crust belonging to a + // "baby" proto-plate: join an adjacent baby if one exists, else start one. + int babyId = -1; + for (int j : cells[c].neighbors) + if (cells[j].plateId >= 0 && plates[cells[j].plateId].baby) { babyId = cells[j].plateId; break; } + if (babyId < 0) { + babyId = acquirePlate(); + plates[babyId].baby = true; + plates[babyId].type = PlateType::Oceanic; + setPlateSpeed(plates[babyId], 0.0); // sits at the ridge, just grows + } + cells[c].plateId = babyId; + cells[c].oceanic = true; cells[c].elevation = cfg.ridgeDepth; cells[c].geoAge = 0.0; + cells[c].drift += cellSize; cells[c].invader = -1; + } + } + + // Coherence cleanup: a cell mostly surrounded by one other plate joins it. + std::vector p2(n); std::vector e2(n), a2(n); std::vector o2(n); + for (int i = 0; i < n; ++i) { p2[i] = cells[i].plateId; e2[i] = cells[i].elevation; a2[i] = cells[i].geoAge; o2[i] = cells[i].oceanic; } + for (int i = 0; i < n; ++i) { + if (cells[i].plateId >= 0 && plates[cells[i].plateId].baby) continue; // leave baby crust alone + int bestP = -1, bestC = 0, src = -1; + for (int j : cells[i].neighbors) { + int pj = p2[j], cnt = 0, s = -1; + for (int k : cells[i].neighbors) if (p2[k] == pj) { ++cnt; if (s < 0) s = k; } + if (cnt > bestC) { bestC = cnt; bestP = pj; src = s; } + } + if (bestP >= 0 && !plates[bestP].baby && bestP != p2[i] && bestC * 2 > (int)cells[i].neighbors.size()) { + cells[i].plateId = bestP; cells[i].elevation = e2[src]; cells[i].oceanic = o2[src]; + cells[i].geoAge = a2[src]; cells[i].drift = 0.0; cells[i].invader = -1; + } + } + + // Soft land band: land is free to vary within +/- landBand of targetLand (so + // volcanic islands and natural subduction can change it), and is nudged back + // only past the band edges -- with hard guards so it can never reach 0/100%. + // The neighbour-count lambdas use a pre-nudge snapshot of oceanic flags: reading + // the live state would chain-flip cells along index-order adjacency (snaking). + std::vector oSnap(n); + for (int i = 0; i < n; ++i) oSnap[i] = cells[i].oceanic; + auto oceanNbrs = [&](int i){ int c = 0; for (int j : cells[i].neighbors) if (oSnap[j]) ++c; return c; }; + auto landNbrs = [&](int i){ int c = 0, s = -1; for (int j : cells[i].neighbors) if (!oSnap[j]) { ++c; if (s < 0) s = j; } return std::pair{c, s}; }; + int lc = landCount(); + int lo = std::max((int)(targetLand * (1.0 - cfg.landBand)), std::max(1, n / 100)); + int hi = std::min((int)(targetLand * (1.0 + cfg.landBand)), n - std::max(1, n / 100)); + int diff = (lc > hi) ? (lc - hi) : (lc < lo ? (lc - lo) : 0); + for (int want = 4; want >= 2 && diff > 0; --want) // too much land: rift thinnest margins + for (int i = 0; i < n && diff > 0; ++i) + if (!cells[i].oceanic && !plates[cells[i].plateId].baby && oceanNbrs(i) >= want) { + cells[i].oceanic = true; cells[i].elevation = cfg.ridgeDepth; cells[i].geoAge = 0.0; --diff; + } + for (int want = 4; want >= 2 && diff < 0; --want) // too little land: accrete ocean margins + for (int i = 0; i < n && diff < 0; ++i) { + auto [ln, ls] = landNbrs(i); + if (cells[i].oceanic && !plates[cells[i].plateId].baby && ln >= want) { + cells[i].oceanic = false; cells[i].elevation = cells[ls].elevation; cells[i].geoAge = cells[ls].geoAge; ++diff; + } + } + + // Periodic plate dynamics: every splitCheckEvery iterations, kick stalemated + // plates, promote grown-up baby plates, and split oversized plates -- keeping + // the plate count lively instead of collapsing to one giant plate. + if (++driftIter % cfg.splitCheckEvery == 0) { + coalesceBabyPlates(); // merge scattered rift cells into ridge strips + deleteEnclosedPlates(); // absorb plates ringed by a single other + fuseMiniPlates(); // clusters of mini plates steal + merge + std::vector cnt(plates.size(), 0); + for (auto& c : cells) if (c.plateId >= 0 && c.plateId < (int)cnt.size()) ++cnt[c.plateId]; + kickStalemates(cnt); + promoteBabyPlates(cnt); + maybeSplitPlates(cnt); + // Refresh the per-plate snapshot + free-list (dead slots) for next window. + sPrevCount.assign(plates.size(), 0); + for (auto& c : cells) if (c.plateId >= 0) ++sPrevCount[c.plateId]; + sFreePlateIds.clear(); + for (int p = 0; p < (int)plates.size(); ++p) if (sPrevCount[p] == 0) sFreePlateIds.push_back(p); + } +} + +// Fission: a plate over splitFraction of all cells rifts in two with a +// probability that ramps from splitProbBase at the threshold (20%->5%, 21%->10%, +// ... 39%->100%). Checked once per window. +void Planet::maybeSplitPlates(const std::vector& cnt) { + const int n = (int)cells.size(); + const int existing = (int)cnt.size(); // don't re-split plates born this pass + for (int pid = 0; pid < existing; ++pid) { + if (plates[pid].baby || cnt[pid] < 8) continue; + double frac = (double)cnt[pid] / n; + if (frac < cfg.splitFraction) continue; + double overPct = (frac - cfg.splitFraction) * 100.0; + double p = std::clamp(cfg.splitProbBase + cfg.splitProbSlope * overPct, 0.0, 1.0); + if (rndf() < p) splitPlate(pid); + } +} + +// Cut plate `pid` roughly in half along a random great circle through its +// centroid; one half keeps `pid`, the other becomes a new plate with a fresh +// random drift vector. Crust type/elevation/age are preserved (land conserved). +void Planet::splitPlate(int pid) { + const int n = (int)cells.size(); + std::vector members; Vec3 centroid{0, 0, 0}; + for (int i = 0; i < n; ++i) + if (cells[i].plateId == pid) { members.push_back(i); centroid = centroid + cells[i].unit; } + if ((int)members.size() < 8) return; + centroid = centroid.normalized(); + // Random tangent at the centroid = the cut plane's normal (so the great circle + // passes through the centroid and bisects the blob). + Vec3 ref = (std::fabs(centroid.y) < 0.99) ? Vec3{0, 1, 0} : Vec3{1, 0, 0}; + Vec3 t = ref.cross(centroid).normalized(); + Vec3 b = centroid.cross(t).normalized(); + double ang = rndf() * 2.0 * M_PI; + Vec3 cutN = (t * std::cos(ang) + b * std::sin(ang)).normalized(); + + int newId = acquirePlate(); + plates[newId].type = plates[pid].type; + plates[newId].baby = false; + randomizePlateDrift(plates[newId]); + for (int i : members) if (cells[i].unit.dot(cutN) > 0.0) cells[i].plateId = newId; +} + +// A plate whose cell count barely changed over the last window (and that borders +// another plate) is stalemated: give it a new random direction and add momentum +// (a speed boost toward the cap) to break the deadlock. +void Planet::kickStalemates(const std::vector& cnt) { + sStaleStreak.resize(plates.size(), 0); + if (sPrevCount.empty()) return; // need a baseline window first + const int lim = (int)std::min(cnt.size(), sPrevCount.size()); + std::vector hasBoundary(plates.size(), 0); + for (int i = 0; i < (int)cells.size(); ++i) + for (int j : cells[i].neighbors) + if (cells[j].plateId != cells[i].plateId) { if (cells[i].plateId >= 0) hasBoundary[cells[i].plateId] = 1; break; } + for (int pid = 0; pid < lim; ++pid) { + if (plates[pid].baby || cnt[pid] < 8 || !hasBoundary[pid]) { sStaleStreak[pid] = 0; continue; } + int delta = std::abs(cnt[pid] - sPrevCount[pid]); + if (delta < (int)(cfg.stalemateEps * cnt[pid] + 0.5)) ++sStaleStreak[pid]; + else sStaleStreak[pid] = 0; + if (sStaleStreak[pid] >= cfg.stalemateWindows) { // stuck for several windows: kick + Vec3 axis{rndf() * 2 - 1, rndf() * 2 - 1, rndf() * 2 - 1}; + plates[pid].driftAxis = axis.normalized(); + setPlateSpeed(plates[pid], std::min(cfg.maxDriftSpeed, plates[pid].speedCmYr * cfg.stalemateBoost)); + sStaleStreak[pid] = 0; + } + } +} + +// Absorb any plate whose entire border touches a single other plate (an inclusion +// inside that plate): its cells join the surrounder, keeping their own crust. On a +// mutual pair (two plates each ringed only by the other) only the smaller is +// absorbed, so they can't just swap labels. Emptied ids are reclaimed. +void Planet::deleteEnclosedPlates() { + const int n = (int)cells.size(); + const int P = (int)plates.size(); + std::vector surr(P, -2); // -2 none seen, -1 multiple, >=0 the single other + std::vector cnt(P, 0); + for (int i = 0; i < n; ++i) { + int pi = cells[i].plateId; if (pi < 0) continue; + ++cnt[pi]; + for (int j : cells[i].neighbors) { + int pj = cells[j].plateId; + if (pj == pi || pj < 0) continue; + if (surr[pi] == -2) surr[pi] = pj; + else if (surr[pi] != pj) surr[pi] = -1; + } + } + std::vector target(P); + for (int pid = 0; pid < P; ++pid) { + target[pid] = pid; + int s = surr[pid]; + if (cnt[pid] == 0 || cnt[pid] == n || s < 0) continue; + if (cnt[pid] < cnt[s] || (cnt[pid] == cnt[s] && pid > s)) target[pid] = s; // absorb smaller + } + for (int i = 0; i < n; ++i) { + int pi = cells[i].plateId; + if (pi >= 0 && target[pi] != pi) { cells[i].plateId = target[pi]; cells[i].drift = 0.0; cells[i].invader = -1; } + } +} + +// When at least fuseMinPlates mini (small, non-baby) plates cluster together they +// fuse into one and steal a ring of cells from their largest big neighbour -- +// terrane amalgamation, the inverse of fission. One ring per check (bounded). +void Planet::fuseMiniPlates() { + const int n = (int)cells.size(); + std::vector cnt(plates.size(), 0); + for (auto& c : cells) if (c.plateId >= 0 && c.plateId < (int)cnt.size()) ++cnt[c.plateId]; + auto isMini = [&](int pid){ return pid >= 0 && pid < (int)plates.size() && !plates[pid].baby + && cnt[pid] > 0 && cnt[pid] < cfg.miniPlateCells; }; + std::vector seen(n, 0); + std::queue q; + for (int i = 0; i < n; ++i) { + if (seen[i] || !isMini(cells[i].plateId)) continue; + std::vector comp, ids; + seen[i] = 1; q.push(i); + while (!q.empty()) { + int u = q.front(); q.pop(); comp.push_back(u); + int pu = cells[u].plateId; + if (std::find(ids.begin(), ids.end(), pu) == ids.end()) ids.push_back(pu); + for (int v : cells[u].neighbors) + if (!seen[v] && isMini(cells[v].plateId)) { seen[v] = 1; q.push(v); } + } + if ((int)ids.size() < cfg.fuseMinPlates) continue; + // Fuse into the largest member (ties -> smallest id). + int target = ids[0]; + for (int id : ids) if (cnt[id] > cnt[target] || (cnt[id] == cnt[target] && id < target)) target = id; + for (int u : comp) { cells[u].plateId = target; cells[u].drift = 0.0; cells[u].invader = -1; } + // Steal one ring from the largest adjacent big (non-mini, non-baby) plate. + int bigP = -1; + for (int u : comp) + for (int v : cells[u].neighbors) { + int pv = cells[v].plateId; + if (pv == target || pv < 0 || plates[pv].baby || isMini(pv)) continue; + if (bigP < 0 || cnt[pv] > cnt[bigP]) bigP = pv; + } + if (bigP >= 0) + for (int u : comp) + for (int v : cells[u].neighbors) + if (cells[v].plateId == bigP) { cells[v].plateId = target; cells[v].drift = 0.0; cells[v].invader = -1; } + } +} + +// Group baby (young rift) cells into connected blobs. A blob of at least +// babyMinCells coalesces into one baby plate id (the smallest in the blob) -- a +// coherent mid-ocean-ridge strip that can grow to promotion size. Smaller blobs +// are advection noise (isolated single-cell rifts) and dissolve back into the +// neighbouring plate. Emptied ids become dead slots reclaimed by the free-list. +void Planet::coalesceBabyPlates() { + const int n = (int)cells.size(); + std::vector seen(n, 0); + auto isBaby = [&](int i){ return cells[i].plateId >= 0 && plates[cells[i].plateId].baby; }; + std::queue q; + for (int i = 0; i < n; ++i) { + if (seen[i] || !isBaby(i)) continue; + std::vector comp; int minId = cells[i].plateId; + seen[i] = 1; q.push(i); + while (!q.empty()) { + int u = q.front(); q.pop(); comp.push_back(u); + minId = std::min(minId, cells[u].plateId); + for (int v : cells[u].neighbors) + if (!seen[v] && isBaby(v)) { seen[v] = 1; q.push(v); } + } + if ((int)comp.size() >= cfg.babyMinCells) { + for (int u : comp) cells[u].plateId = minId; // keep: a real ridge strip + plates[minId].baby = true; + } else { + for (int u : comp) { // dissolve: noise -> adopt a neighbour + int src = -1; + for (int v : cells[u].neighbors) + if (cells[v].plateId >= 0 && !plates[cells[v].plateId].baby) { src = v; break; } + if (src < 0) continue; // surrounded only by baby cells: leave it + cells[u].plateId = cells[src].plateId; + cells[u].oceanic = cells[src].oceanic; + cells[u].elevation = cells[src].elevation; + cells[u].geoAge = cells[src].geoAge; + cells[u].drift = 0.0; cells[u].invader = -1; + } + } + } +} + +// A baby (spreading) plate that has grown past babyPromoteFrac of all cells +// becomes a real plate: it gets a random drift vector and grows volcanic +// landmass on its most interior cells (simulated volcanic activity). +void Planet::promoteBabyPlates(const std::vector& cnt) { + const int n = (int)cells.size(); + const int threshold = std::max(4, (int)(cfg.babyPromoteFrac * n)); + const int lim = (int)cnt.size(); + for (int pid = 0; pid < lim; ++pid) { + if (!plates[pid].baby || cnt[pid] < threshold) continue; + plates[pid].baby = false; + randomizePlateDrift(plates[pid]); + std::vector members; + for (int i = 0; i < n; ++i) if (cells[i].plateId == pid) members.push_back(i); + auto interior = [&](int i){ int s = 0; for (int j : cells[i].neighbors) if (cells[j].plateId == pid) ++s; return s; }; + std::sort(members.begin(), members.end(), [&](int a, int b){ return interior(a) > interior(b); }); + int makeLand = (int)(cfg.volcanicLandFrac * members.size()); + for (int k = 0; k < makeLand; ++k) { + cells[members[k]].oceanic = false; + cells[members[k]].elevation = cfg.volcanicElev; // young volcanic island + } + } +} diff --git a/src/sim/PlanetErosion.cpp b/src/sim/PlanetErosion.cpp new file mode 100644 index 0000000..b8dd6de --- /dev/null +++ b/src/sim/PlanetErosion.cpp @@ -0,0 +1,58 @@ +#include "Planet.hpp" +#include +#include + +// --- Phase 2: erosion + sea level ------------------------------------------- + +// Mass-conserving downhill sediment transport: along each edge the higher cell +// gives material to the lower one, faster above sea level (subaerial weathering) +// than below (slow submarine). Highs wear down toward an uplift<->erosion +// equilibrium; the sediment piles in basins and below sea level, building shelves +// and deltas and slowly filling trenches. One double-buffered gather pass (each +// cell writes only its own index -> bit-identical for any thread count). +void Planet::erode(double dtMy) { + const int n = (int)cells.size(); + sErode.resize(n); + const double sea = cfg.seaLevel; + + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + double ei = cells[i].elevation, delta = 0.0; + for (int j : cells[i].neighbors) { + double ej = cells[j].elevation; + if (ei > ej) { // i is higher: i -> j outflow + double rate = (ei > sea) ? cfg.erosionLandRate : cfg.erosionSeaRate; + delta -= std::min(rate * dtMy, 0.5) * (ei - ej); + } else if (ej > ei) { // j is higher: j -> i inflow + double rate = (ej > sea) ? cfg.erosionLandRate : cfg.erosionSeaRate; + delta += std::min(rate * dtMy, 0.5) * (ej - ei); + } + } + sErode[i] = delta; + } + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) + cells[i].elevation = std::clamp(cells[i].elevation + sErode[i], -11000.0, 9000.0); + + if (++erodeIter % cfg.seaLevelEvery == 0) adjustSeaLevel(); +} + +// Gradual eustatic controller: if the geographic land fraction is outside a +// deadband around the target, nudge sea level by a small fixed step (raising it +// floods land, lowering exposes it). Called infrequently (seaLevelEvery) so the +// coastline drifts slowly, not in a jump. A nudge is only taken if it actually +// reduces the error -- otherwise (e.g. a big mass of cells at one elevation, where +// a full step would overshoot) it rests at the closest a fixed step allows instead +// of oscillating back and forth across that "cliff". +void Planet::adjustSeaLevel() { + const int n = (int)cells.size(); + if (n == 0) return; + auto landFracAt = [&](double sl) { + int a = 0; for (const auto& c : cells) if (c.elevation > sl) ++a; return (double)a / n; + }; + double err = landFracAt(cfg.seaLevel) - cfg.landFractionTarget; + if (std::fabs(err) <= cfg.seaLevelTol) return; // within deadband: rest + double cand = cfg.seaLevel + (err > 0 ? cfg.seaLevelStep : -cfg.seaLevelStep); + double candErr = landFracAt(cand) - cfg.landFractionTarget; + if (std::fabs(candErr) < std::fabs(err)) cfg.seaLevel = cand; // nudge only if it helps +} diff --git a/src/sim/PlanetHydrology.cpp b/src/sim/PlanetHydrology.cpp new file mode 100644 index 0000000..5ef4487 --- /dev/null +++ b/src/sim/PlanetHydrology.cpp @@ -0,0 +1,104 @@ +#include "Planet.hpp" +#include +#include +#include + +// --- Phase 3: hydrology (rivers, lakes, fluvial erosion) -------------------- + +namespace { + struct PQItem { double fill; int idx; }; + // min-heap on fill, index tie-break -> deterministic ordering. + struct PQGreater { + bool operator()(const PQItem& a, const PQItem& b) const { + return a.fill > b.fill || (a.fill == b.fill && a.idx > b.idx); + } + }; +} + +// Build the drainage network from the current elevation field (deterministic, +// not saved): depression-fill (priority-flood + epsilon tilt) -> sFill/lakeDepth, +// steepest-descent over the filled surface -> sFlowTo, flow accumulation -> +// sDischarge. Ocean cells (elev <= seaLevel) are outlets/sinks. +void Planet::routeFlow() { + const int n = (int)cells.size(); + const double sea = cfg.seaLevel; + const double INF = 1e18, EPS = 1e-3; // eps tilt (m) so flats still drain + + sFill.assign(n, INF); + sLakeDepth.assign(n, 0.0); + sDischarge.assign(n, 0.0); + sFlowTo.assign(n, -1); + + std::priority_queue, PQGreater> pq; + std::vector done(n, 0); + int outlets = 0; + for (int i = 0; i < n; ++i) + if (cells[i].elevation <= sea) { sFill[i] = cells[i].elevation; done[i] = 1; pq.push({sFill[i], i}); ++outlets; } + if (outlets == 0) { // no ocean: seed the single lowest cell + int lo = 0; for (int i = 1; i < n; ++i) if (cells[i].elevation < cells[lo].elevation) lo = i; + sFill[lo] = cells[lo].elevation; done[lo] = 1; pq.push({sFill[lo], lo}); + } + while (!pq.empty()) { + PQItem c = pq.top(); pq.pop(); + for (int nb : cells[c.idx].neighbors) { + if (done[nb]) continue; + double f = std::max(cells[nb].elevation, c.fill + EPS); + sFill[nb] = f; done[nb] = 1; pq.push({f, nb}); + } + } + // Lake depth + flow direction (lowest-fill neighbour). + for (int i = 0; i < n; ++i) { + sLakeDepth[i] = std::max(0.0, sFill[i] - cells[i].elevation); + if (cells[i].elevation <= sea) { sFlowTo[i] = -1; continue; } // ocean = sink + int best = -1; double bestF = sFill[i]; + for (int nb : cells[i].neighbors) + if (sFill[nb] < bestF) { bestF = sFill[nb]; best = nb; } + sFlowTo[i] = best; + } + // Accumulation order: descending fill (so upstream is processed first). + sHydroOrder.resize(n); + for (int i = 0; i < n; ++i) sHydroOrder[i] = i; + std::sort(sHydroOrder.begin(), sHydroOrder.end(), [&](int a, int b) { + return sFill[a] > sFill[b] || (sFill[a] == sFill[b] && a < b); + }); + for (int idx : sHydroOrder) { + if (cells[idx].elevation > sea) sDischarge[idx] += cfg.rainfall; // rain on land + int d = sFlowTo[idx]; + if (d >= 0) sDischarge[d] += sDischarge[idx]; + } +} + +void Planet::computeHydrology() { routeFlow(); } + +// One hydrology tick: route, then mass-conserving stream-power incision with +// downstream sediment transport/deposition. Carving rivers, building deltas at +// mouths, slowly filling lakes. sum(elevation) is conserved (sink deposition). +void Planet::hydrology(double dtMy) { + routeFlow(); + const int n = (int)cells.size(); + const double m = cfg.riverDischargeExp, nexp = cfg.riverSlopeExp; + std::vector load(n, 0.0); // sediment arriving at each cell + for (int idx : sHydroOrder) { // upstream -> downstream + double carried = load[idx]; + int d = sFlowTo[idx]; + if (d < 0) { cells[idx].elevation += carried; continue; } // sink: deposit all + Vec3 ui = cells[idx].unit, ud = cells[d].unit; + double ang = std::acos(std::clamp(ui.dot(ud), -1.0, 1.0)); + double dist = std::max(1.0, cfg.radius * ang); + double drop = std::max(0.0, sFill[idx] - sFill[d]); + double S = drop / dist; + double cap = cfg.riverTransport * sDischarge[idx] * S; + if (carried > cap) { // over capacity -> deposit + double dep = (carried - cap) * cfg.depFrac; + cells[idx].elevation += dep; carried -= dep; + } else { // under capacity -> incise + double inc = cfg.riverIncision * std::pow(sDischarge[idx], m) * std::pow(S, nexp) * dtMy; + inc = std::min(inc, 0.5 * drop); // never invert the routed gradient + if (inc > 0.0) { cells[idx].elevation -= inc; carried += inc; } + } + load[d] += carried; // pass remaining sediment downstream + } + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) + cells[i].elevation = std::clamp(cells[i].elevation, -11000.0, 9000.0); +} diff --git a/src/sim/PlanetIO.cpp b/src/sim/PlanetIO.cpp new file mode 100644 index 0000000..ce9919c --- /dev/null +++ b/src/sim/PlanetIO.cpp @@ -0,0 +1,263 @@ +#include "Planet.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +// --- Config file (text) + save/load (binary) -------------------------------- + +// One shared field table so saveConfig/loadConfig can never drift apart. +// D = double field, I = int field, U = uint32 field. +#define CONFIG_FIELDS(D, I, U) \ + D(radius) D(seaLevel) D(axialTilt) D(continentBase) D(oceanBase) D(upliftGain) D(relax) \ + D(collisionFactor) D(arcFactor) D(isostaticPersist) D(rootScale) \ + D(seafloorSubsidence) D(seafloorSeedAge) \ + D(maxDriftSpeed) D(ridgeDepth) D(splitFraction) D(splitProbBase) D(splitProbSlope) \ + D(stalemateEps) D(stalemateBoost) D(babyPromoteFrac) D(volcanicLandFrac) \ + D(volcanicElev) D(landBand) D(erosionLandRate) D(erosionSeaRate) \ + D(landFractionTarget) D(seaLevelStep) D(seaLevelTol) \ + D(phase3AfterMy) D(phase3DtScale) D(rainfall) D(riverThreshold) D(riverIncision) \ + D(riverDischargeExp) D(riverSlopeExp) D(riverTransport) D(depFrac) \ + D(biomeEquatorTemp) D(biomePoleDrop) D(biomeLatExp) D(biomeElevLapse) \ + D(biomeIceTemp) D(biomeTundraTemp) D(biomeTaigaTemp) D(biomeSavannaTemp) \ + D(biomeMountainElev) D(biomeHillsElev) D(biomeBeachBand) D(biomeLowlandElev) \ + D(biomeWetlandMoist) D(biomeDesertMoist) D(biomeGrassMoist) D(biomeTaigaMoist) \ + D(biomeLakeMinDepth) \ + D(climateOceanMoisture) D(climateRainEfficiency) D(climateOrographic) \ + D(climateOroRefHeight) D(climateContinentality) \ + I(subdivisions) I(plateCount) I(beltWidth) I(splitCheckEvery) I(stalemateWindows) \ + I(miniPlateCells) I(fuseMinPlates) I(babyMinCells) I(seaLevelEvery) \ + I(climateWindPasses) I(climateMoistureSmooth) \ + U(seed) + +// Write all config fields as `key = value` lines (no header). Shared by the text +// config file (saveConfig) and the self-describing config block embedded in saves. +static void writeConfigFields(std::ostream& os, const PlanetConfig& cfg) { + #define WRITE(name) os << #name " = " << cfg.name << "\n"; + CONFIG_FIELDS(WRITE, WRITE, WRITE) + #undef WRITE +} + +// Parse `key = value` lines from any stream into cfg (unknown keys ignored, missing +// keys keep cfg's existing value). Shared by loadConfig + readState. This is why +// adding/removing config fields no longer breaks saves. +static void parseConfigStream(std::istream& is, PlanetConfig& cfg) { + auto trim = [](std::string& s) { + size_t a = s.find_first_not_of(" \t\r\n"), b = s.find_last_not_of(" \t\r\n"); + if (a == std::string::npos) s.clear(); else s = s.substr(a, b - a + 1); + }; + std::string line; + while (std::getline(is, line)) { + size_t hash = line.find('#'); if (hash != std::string::npos) line.resize(hash); + size_t eq = line.find('='); if (eq == std::string::npos) continue; + std::string key = line.substr(0, eq), val = line.substr(eq + 1); + trim(key); trim(val); + if (key.empty() || val.empty()) continue; + #define D(name) if (key == #name) { try { cfg.name = std::stod(val); } catch (...) {} continue; } + #define I(name) if (key == #name) { try { cfg.name = std::stoi(val); } catch (...) {} continue; } + #define U(name) if (key == #name) { try { cfg.name = (uint32_t)std::stoul(val); } catch (...) {} continue; } + CONFIG_FIELDS(D, I, U) + #undef D + #undef I + #undef U + } +} + +bool saveConfig(const std::string& path, const PlanetConfig& cfg) { + std::ofstream os(path); + if (!os) return false; + os.precision(15); // enough for the (nice, decimal) defaults; trailing zeros trimmed + os << "# Planet config -- edit values, then reload in-app (F2) or restart.\n"; + os << "# key = value; '#' starts a comment; unknown keys are ignored.\n\n"; + writeConfigFields(os, cfg); + return (bool)os; +} + +bool loadConfig(const std::string& path, PlanetConfig& cfg) { + std::ifstream is(path); + if (!is) return false; + parseConfigStream(is, cfg); + return true; +} + +std::string validateConfig(const PlanetConfig& cfg) { + auto rng = [](double v, double lo, double hi, const char* name) -> std::string { + if (v >= lo && v <= hi) return {}; + return std::string(name) + " = " + std::to_string(v) + " (expected " + + std::to_string(lo) + ".." + std::to_string(hi) + ")"; + }; + auto irng = [](int v, int lo, int hi, const char* name) -> std::string { + if (v >= lo && v <= hi) return {}; + return std::string(name) + " = " + std::to_string(v) + " (expected " + + std::to_string(lo) + ".." + std::to_string(hi) + ")"; + }; + std::vector bad; + auto E = [&](const std::string& s) { if (!s.empty()) bad.push_back(s); }; + + E(rng(cfg.radius, 1.0e3, 1.0e8, "radius")); + E(rng(cfg.seaLevel, -11000.0, 9000.0, "seaLevel")); + E(rng(cfg.axialTilt, 0.0, 180.0, "axialTilt")); + E(rng(cfg.continentBase, -2000.0, 6000.0, "continentBase")); + E(rng(cfg.oceanBase, -11000.0, 1000.0, "oceanBase")); + E(rng(cfg.upliftGain, 100.0, 1.0e7, "upliftGain")); + E(rng(cfg.relax, 0.001, 0.5, "relax")); + E(rng(cfg.collisionFactor, 0.0, 20.0, "collisionFactor")); + E(rng(cfg.arcFactor, 0.0, 20.0, "arcFactor")); + E(rng(cfg.isostaticPersist, 0.0, 0.95, "isostaticPersist")); + E(rng(cfg.rootScale, 100.0, 20000.0, "rootScale")); + E(rng(cfg.seafloorSubsidence, 0.0, 2000.0, "seafloorSubsidence")); + E(rng(cfg.seafloorSeedAge, 0.0, 1000.0, "seafloorSeedAge")); + E(rng(cfg.maxDriftSpeed, 0.1, 100.0, "maxDriftSpeed")); + E(rng(cfg.ridgeDepth, -8000.0, 0.0, "ridgeDepth")); + E(rng(cfg.splitFraction, 0.0, 1.0, "splitFraction")); + E(rng(cfg.splitProbBase, 0.0, 1.0, "splitProbBase")); + E(rng(cfg.splitProbSlope, 0.0, 1.0, "splitProbSlope")); + E(rng(cfg.stalemateEps, 0.0, 1.0, "stalemateEps")); + E(rng(cfg.stalemateBoost, 1.0, 5.0, "stalemateBoost")); + E(rng(cfg.babyPromoteFrac, 0.001, 0.5, "babyPromoteFrac")); + E(rng(cfg.volcanicLandFrac, 0.0, 1.0, "volcanicLandFrac")); + E(rng(cfg.volcanicElev, -1000.0, 5000.0, "volcanicElev")); + E(rng(cfg.landBand, 0.0, 1.0, "landBand")); + E(rng(cfg.erosionLandRate, 0.0, 1.0, "erosionLandRate")); + E(rng(cfg.erosionSeaRate, 0.0, 1.0, "erosionSeaRate")); + E(rng(cfg.landFractionTarget, 0.01, 0.99, "landFractionTarget")); + E(rng(cfg.seaLevelStep, 1.0, 2000.0, "seaLevelStep")); + E(rng(cfg.seaLevelTol, 0.001, 0.5, "seaLevelTol")); + E(rng(cfg.phase3AfterMy, 0.0, 1.0e6, "phase3AfterMy")); + E(rng(cfg.phase3DtScale, 0.001, 1.0, "phase3DtScale")); + E(rng(cfg.rainfall, 0.0, 1.0e6, "rainfall")); + E(rng(cfg.riverThreshold, 0.0, 1.0e9, "riverThreshold")); + E(rng(cfg.riverIncision, 0.0, 1.0e6, "riverIncision")); + E(rng(cfg.riverDischargeExp, 0.0, 5.0, "riverDischargeExp")); + E(rng(cfg.riverSlopeExp, 0.0, 5.0, "riverSlopeExp")); + E(rng(cfg.riverTransport, 0.0, 1.0e6, "riverTransport")); + E(rng(cfg.depFrac, 0.0, 1.0, "depFrac")); + E(rng(cfg.biomeEquatorTemp, -50.0, 80.0, "biomeEquatorTemp")); + E(rng(cfg.biomePoleDrop, 0.0, 150.0, "biomePoleDrop")); + E(rng(cfg.biomeLatExp, 0.1, 6.0, "biomeLatExp")); + E(rng(cfg.biomeElevLapse, 0.0, 0.05, "biomeElevLapse")); + E(rng(cfg.biomeIceTemp, -60.0, 20.0, "biomeIceTemp")); + E(rng(cfg.biomeTundraTemp, -60.0, 40.0, "biomeTundraTemp")); + E(rng(cfg.biomeTaigaTemp, -60.0, 40.0, "biomeTaigaTemp")); + E(rng(cfg.biomeSavannaTemp, -20.0, 60.0, "biomeSavannaTemp")); + E(rng(cfg.biomeMountainElev, 0.0, 11000.0, "biomeMountainElev")); + E(rng(cfg.biomeHillsElev, 0.0, 11000.0, "biomeHillsElev")); + E(rng(cfg.biomeBeachBand, 0.0, 2000.0, "biomeBeachBand")); + E(rng(cfg.biomeLowlandElev, 0.0, 11000.0, "biomeLowlandElev")); + E(rng(cfg.biomeWetlandMoist, 0.0, 1.0, "biomeWetlandMoist")); + E(rng(cfg.biomeDesertMoist, 0.0, 1.0, "biomeDesertMoist")); + E(rng(cfg.biomeGrassMoist, 0.0, 1.0, "biomeGrassMoist")); + E(rng(cfg.biomeTaigaMoist, 0.0, 1.0, "biomeTaigaMoist")); + E(rng(cfg.biomeLakeMinDepth, 0.0, 5000.0, "biomeLakeMinDepth")); + E(rng(cfg.climateOceanMoisture, 0.0, 1.0e3, "climateOceanMoisture")); + E(rng(cfg.climateRainEfficiency, 0.0, 1.0, "climateRainEfficiency")); + E(rng(cfg.climateOrographic, 0.0, 50.0, "climateOrographic")); + E(rng(cfg.climateOroRefHeight, 1.0, 1.0e5, "climateOroRefHeight")); + E(rng(cfg.climateContinentality, 0.0, 1.0, "climateContinentality")); + E(irng(cfg.subdivisions, 0, 7, "subdivisions")); + E(irng(cfg.plateCount, 1, 100, "plateCount")); + E(irng(cfg.beltWidth, 1, 12, "beltWidth")); + E(irng(cfg.splitCheckEvery, 1, 10000, "splitCheckEvery")); + E(irng(cfg.stalemateWindows, 1, 100, "stalemateWindows")); + E(irng(cfg.miniPlateCells, 1, 10000, "miniPlateCells")); + E(irng(cfg.fuseMinPlates, 2, 50, "fuseMinPlates")); + E(irng(cfg.babyMinCells, 1, 1000, "babyMinCells")); + E(irng(cfg.seaLevelEvery, 1, 100000, "seaLevelEvery")); + E(irng(cfg.climateWindPasses, 1, 1000, "climateWindPasses")); + E(irng(cfg.climateMoistureSmooth, 0, 100, "climateMoistureSmooth")); + + if (cfg.oceanBase >= cfg.continentBase) + bad.push_back("oceanBase >= continentBase (ocean floor must be below continents)"); + + if (bad.empty()) return {}; + std::string msg = "Bad config:"; + for (auto& s : bad) msg += "\n " + s; + return msg; +} + +namespace { + template void writePod(std::ostream& os, const T& v) { + static_assert(std::is_trivially_copyable::value, "writePod needs a POD type"); + os.write(reinterpret_cast(&v), sizeof(T)); + } + template void readPod(std::istream& is, T& v) { + static_assert(std::is_trivially_copyable::value, "readPod needs a POD type"); + is.read(reinterpret_cast(&v), sizeof(T)); + } + template void writeVec(std::ostream& os, const std::vector& v) { + static_assert(std::is_trivially_copyable::value, "writeVec needs POD elements"); + uint64_t n = v.size(); writePod(os, n); + if (n) os.write(reinterpret_cast(v.data()), (std::streamsize)(n * sizeof(T))); + } + template void readVec(std::istream& is, std::vector& v) { + static_assert(std::is_trivially_copyable::value, "readVec needs POD elements"); + uint64_t n = 0; readPod(is, n); v.resize((size_t)n); + if (n) is.read(reinterpret_cast(v.data()), (std::streamsize)(n * sizeof(T))); + } +} + +// Full simulation state. Geometry (unit/neighbors) is NOT stored -- it is rebuilt +// from cfg.subdivisions on load -- so only the dynamic per-cell fields are saved. +void Planet::writeState(std::ostream& os) const { + // Config is stored as a self-describing key=value text block (length-prefixed), + // not a raw POD dump, so adding/removing config fields never breaks old saves + // (unknown keys ignored, missing keys keep their defaults). precision(17) = + // max_digits10 for double, so values round-trip exactly (deterministic resume). + std::ostringstream cfgss; cfgss.precision(17); + writeConfigFields(cfgss, cfg); + std::string cfgText = cfgss.str(); + uint64_t clen = cfgText.size(); writePod(os, clen); + os.write(cfgText.data(), (std::streamsize)clen); + writePod(os, rngState); + writePod(os, driftIter); + writePod(os, erodeIter); + writePod(os, targetLand); + uint64_t nc = cells.size(); writePod(os, nc); + for (const Cell& c : cells) { + writePod(os, c.elevation); writePod(os, c.plateId); + uint8_t oc = c.oceanic ? 1 : 0; writePod(os, oc); + writePod(os, c.geoAge); writePod(os, c.drift); writePod(os, c.invader); + uint8_t bm = (uint8_t)c.biome; writePod(os, bm); // save v4: per-cell biome + } + writeVec(os, plates); + writeVec(os, sPrevCount); + writeVec(os, sStaleStreak); + writeVec(os, sFreePlateIds); +} + +bool Planet::readState(std::istream& is, bool hasBiome) { + // Read the length-prefixed key=value config block (see writeState). A default + // PlanetConfig is parsed over, so fields absent from an older save keep their + // current defaults. The length guard rejects pre-v6 (raw-POD-config) saves. + uint64_t clen = 0; readPod(is, clen); + if (!is || clen > 1000000) return false; + std::string cfgText(clen, '\0'); + is.read(&cfgText[0], (std::streamsize)clen); + if (!is) return false; + PlanetConfig c; + { std::istringstream cis(cfgText); parseConfigStream(cis, c); } + cfg = c; + buildGeometry(); // rebuild unit/neighbors from cfg.subdivisions + readPod(is, rngState); + readPod(is, driftIter); + readPod(is, erodeIter); + readPod(is, targetLand); + uint64_t nc = 0; readPod(is, nc); + if (!is || nc != cells.size()) return false; // subdivisions mismatch / corrupt file + for (Cell& cell : cells) { + readPod(is, cell.elevation); readPod(is, cell.plateId); + uint8_t oc = 0; readPod(is, oc); cell.oceanic = (oc != 0); + readPod(is, cell.geoAge); readPod(is, cell.drift); readPod(is, cell.invader); + if (hasBiome) { uint8_t bm = 0; readPod(is, bm); cell.biome = (Biome)bm; } // save v4 + } + readVec(is, plates); + readVec(is, sPrevCount); + readVec(is, sStaleStreak); + readVec(is, sFreePlateIds); + if (!hasBiome) classifyBiomes(); // old (v3) save: reclassify from loaded state + return (bool)is; +} diff --git a/src/sim/PlanetTectonics.cpp b/src/sim/PlanetTectonics.cpp new file mode 100644 index 0000000..86fc1e6 --- /dev/null +++ b/src/sim/PlanetTectonics.cpp @@ -0,0 +1,164 @@ +#include "Planet.hpp" +#include +#include + +// --- Phase 1/2 tectonics: one stress->uplift->relax tick --------------------- +// Data-parallel: every pass writes only its own cell index (double-buffered +// where it reads a field it writes), so the OpenMP loops are bit-identical for +// any thread count -- determinism preserved. + +double Planet::step() { + const int n = (int)cells.size(); + + // Remember elevation so we can report the largest change this tick (used by + // the viewer to detect when the world has settled into equilibrium). + sOldElev.resize(n); + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) sOldElev[i] = cells[i].elevation; + + // Per-boundary-cell shaping factors (multipliers on upliftGain). + const double collideFactor = 1.0; // continent interior baseline (no special boundary) + const double oceanArcFactor = 0.5; // oceanic island arc (mild) + const double trenchFactor = 2.4; // depth of the subduction trench + const double riftFactor = 0.5; // gentle deepening of divergent rifts + const double beltDecay = 0.6; // stress falloff per ring inland + + // Each pass below is data-parallel: every cell reads its own + neighbor + // data and writes ONLY its own index (double-buffered where it reads a field + // it also writes), so the OpenMP loops give bit-identical results for any + // thread count -- determinism is preserved. + + // 1. Boundary stress. For each cell touching another plate, the signed + // convergence (positive = plates closing, negative = pulling apart), + // plus flags marking subduction geometry from the cell's own crust type. + // Scratch buffers persist across calls (sized once); cleared per step. + sStress.assign(n, 0.0); + sSub.assign(n, 0); // oceanic cell facing continental -> trench + sOver.assign(n, 0); // continental cell facing oceanic -> arc + sColl.assign(n, 0); // continental cell facing continental -> collision + std::vector& stress = sStress; + std::vector& subducting = sSub; + std::vector& overriding = sOver; + std::vector& colliding = sColl; + + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + const Cell& a = cells[i]; + Vec3 va = driftVelocity(a.plateId, a.unit); + + double convergence = 0.0; + int boundaryNeighbors = 0; + bool facesOceanic = false, facesContinental = false; + + for (int nb : a.neighbors) { + const Cell& b = cells[nb]; + if (b.plateId == a.plateId) continue; + ++boundaryNeighbors; + Vec3 dir = (b.unit - a.unit).normalized(); + Vec3 vb = driftVelocity(b.plateId, b.unit); + convergence += (va - vb).dot(dir); + if (b.oceanic) facesOceanic = true; + else facesContinental = true; + } + if (boundaryNeighbors == 0) continue; + stress[i] = convergence / boundaryNeighbors; + + bool oceanic = a.oceanic; + if (oceanic && facesContinental) subducting[i] = 1; + if (!oceanic && facesOceanic) overriding[i] = 1; + if (!oceanic && facesContinental) colliding[i] = 1; // continent-continent orogeny + } + + // 2. Spread the convergent stress a few cell-rings into the plate interior + // so mountain belts have realistic width and flanks, with the peak left + // at the boundary (dilation with geometric decay, not blurring). Only + // the positive component spreads -- subduction trenches stay narrow. + // Boundary cells (sources) are re-anchored after each ring so they never + // accumulate stress from neighbouring boundaries, keeping peaks sharp. + std::vector& belt = sBelt; belt.resize(n); + std::vector& next = sBeltNext; next.resize(n); + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) belt[i] = std::max(0.0, stress[i]); + for (int it = 0; it < cfg.beltWidth; ++it) { + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + double v = belt[i]; + for (int nb : cells[i].neighbors) v = std::max(v, belt[nb] * beltDecay); + next[i] = v; + } + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) if (stress[i] > 0.0) next[i] = stress[i]; + belt.swap(next); + } + + // 3. Convert stress -> elevation change, shaped by plate geometry. + const double gain = cfg.upliftGain; + std::vector& delta = sDelta; delta.resize(n); + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + bool continental = !cells[i].oceanic; + double d = 0.0; + + if (subducting[i]) { + // Narrow, deep oceanic trench where this crust dives under a continent. + d -= std::max(0.0, stress[i]) * gain * trenchFactor; + } else if (belt[i] > 0.0) { + // Orogeny boosts (collision/arc) are drift-only; forming uses the + // original mild factors so Phase 1 settles as it always did. + double f; + if (colliding[i]) f = drifting ? cfg.collisionFactor : 1.0; // continent-continent + else if (overriding[i]) f = drifting ? cfg.arcFactor : 1.1; // continental arc + else f = continental ? collideFactor : oceanArcFactor; + d += belt[i] * gain * f; + } + if (stress[i] < 0.0) // divergent: rift valley / spreading + d += stress[i] * gain * riftFactor; + + delta[i] = d; + } + + // 4. Apply uplift, then isostatic relaxation toward the crust's base + // elevation. Equilibrium sits at base + delta/relax, which bounds growth + // well below the clamp instead of railing in a single tick. + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + double base, relaxEff; + if (cells[i].oceanic) { + base = oceanicBase(cells[i].geoAge); // deepens with crustal age (always on) + relaxEff = cfg.relax; + } else { + base = cfg.continentBase; + // Drift-only: thick (high) continental crust resists isostatic + // relaxation, so collision ranges stand and become erosion-limited + // instead of snapping back. Forming uses full relax (original). + if (drifting) { + double over = std::clamp((cells[i].elevation - cfg.continentBase) / cfg.rootScale, 0.0, 1.0); + relaxEff = cfg.relax * (1.0 - cfg.isostaticPersist * over); + } else { + relaxEff = cfg.relax; + } + } + cells[i].elevation += delta[i]; + cells[i].elevation += (base - cells[i].elevation) * relaxEff; + // (geoAge is crust age in My, advanced by advect() in Phase 2.) + } + + // 5. Gentle diffusion (erosion / sediment transport) so relief stays smooth. + std::vector& smoothed = sSmoothed; smoothed.resize(n); + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) { + double sum = cells[i].elevation; int cnt = 1; + for (int nb : cells[i].neighbors) { sum += cells[nb].elevation; ++cnt; } + smoothed[i] = cells[i].elevation * 0.96 + (sum / cnt) * 0.04; + } + #pragma omp parallel for schedule(static) if(n > 20000) + for (int i = 0; i < n; ++i) + cells[i].elevation = std::clamp(smoothed[i], -11000.0, 9000.0); + + // Largest elevation change this tick -> 0 as the world reaches equilibrium. + double maxChange = 0.0; + #pragma omp parallel for schedule(static) reduction(max:maxChange) if(n > 20000) + for (int i = 0; i < n; ++i) + maxChange = std::max(maxChange, std::fabs(cells[i].elevation - sOldElev[i])); + return maxChange; +} diff --git a/src/sim/PlanetTypes.hpp b/src/sim/PlanetTypes.hpp new file mode 100644 index 0000000..15cb07f --- /dev/null +++ b/src/sim/PlanetTypes.hpp @@ -0,0 +1,190 @@ +#pragma once +#include "Vec3.hpp" +#include +#include +#include + +// Data structures shared across the Planet engine (raylib-free). The Planet +// class itself lives in Planet.hpp; these are the per-cell / per-plate / config +// types it operates on. Geometry never moves -- properties flow over the fixed +// grid (see CLAUDE.md core principle). + +// Fine-resolution subgrid (phases 4/5 hook). Generated on demand for one macro +// cell: a high-res patch of the sphere around that cell, with elevation blended +// from the cell + its neighbors so it transitions smoothly across boundaries. +struct SubCell { + Vec3 unit; // direction on the unit sphere + double elevation = 0.0; // meters (blended + fine detail noise) + int nearestMacro = -1;// macro cell (this cell or a neighbor) it lies under +}; +struct SubGrid { + int macroCell = -1; + int res = 0; // grid is res x res, row-major + std::vector sub; +}; + +enum class PlateType { Oceanic, Continental }; + +// Phase-3 (climate & biomes) classification of a cell, derived from elevation, +// latitude (temperature) and hydrology/coast (moisture). Stored per cell (uint8, +// serialized) so Phase-4 civilization can read it. Keep Ocean == 0 so a default- +// constructed cell reads as ocean. Extend by appending new entries (don't reorder +// -- the numeric value is saved). +enum class Biome : uint8_t { + Ocean, Ice, Lake, Beach, Wetland, Grassland, Savanna, + Desert, Forest, Taiga, Tundra, Hills, Mountains +}; + +struct Plate { + int id = 0; + PlateType type = PlateType::Oceanic; // initial crust type seeded onto cells + Vec3 driftAxis; // rotation axis (unit) for tangential drift on sphere + double driftSpeed; // small value used by the Phase-1 uplift stress + double angSpeed = 0; // Phase-2 advection: rotation rate in radians / My + double speedCmYr = 0;// surface drift speed in cm/year (for display) + bool baby = false; // Phase-2: young spreading proto-plate, not yet promoted +}; + +// A fixed cell on the planet. Geometry never moves; properties flow. +struct Cell { + Vec3 unit; // unit-sphere direction (fixed) + double elevation = 0.0; // meters, relative to sea level (continuous, fine) + int plateId = -1; + double geoAge = 0.0; // My since this crust was (re)formed at a ridge + bool oceanic = true;// crust type travels WITH the cell (Phase-2 advection) + Biome biome = Biome::Ocean; // Phase-3 climate/biome classification (derived) + + // Phase-2 advection accumulator: signed convergence distance built up with + // the dominant other-plate neighbor (+ encroaching, - rifting), and which + // plate is encroaching. + double drift = 0.0; + int invader = -1; + + std::vector neighbors; + + // Phase 4/5 hook: optional fine-resolution subgrid. Null until needed. + std::shared_ptr subgrid; +}; + +struct PlanetConfig { + double radius = 6.371e6; // meters (Earth default) + int subdivisions = 5; // icosphere level + double seaLevel = 0.0; // meters + int plateCount = 12; + uint32_t seed = 1337; + double axialTilt = 23.44; // obliquity, degrees (Earth ~23.4). Visual tilt of + // the 3D globe + spin axis; groundwork for seasons. + + // --- Tectonic tuning (Phase 1) ------------------------------------------ + // Relief builds gradually toward an isostatic equilibrium instead of + // saturating: per-tick uplift competes with a relaxation pull toward the + // plate's base elevation, so peaks asymptote at base + uplift/relax. + double continentBase = 300.0; // m, resting elevation of continental crust + double oceanBase = -6000.0; // m, deep abyssal floor (oldest oceanic crust) + double upliftGain = 1.3e5; // m/tick per unit convergence stress + int beltWidth = 3; // cell-rings a mountain belt spreads inland + double relax = 0.02; // isostatic relaxation toward base, per tick + + // --- Orogeny: collision uplift + isostatic persistence (Phase 2 inc. 4) -- + // Continent-continent collisions build the tallest ranges; thick (high) + // continental crust resists isostatic relaxation, so ranges stand and are + // erosion-limited rather than snapping back to continentBase. + double collisionFactor = 1.8; // continent-continent uplift multiplier (Himalaya) + double arcFactor = 1.4; // continental subduction-arc uplift (Andes) + double isostaticPersist = 0.85; // how much high crust resists relax (0..<1) + double rootScale = 2500.0;// m above continentBase where persistence saturates + + // --- Plate drift (Phase 2) ---------------------------------------------- + double maxDriftSpeed = 20.0; // cm/year; fastest plate (Earth is 1-10) + double ridgeDepth = -2500.0; // m, elevation of brand-new crust at a ridge + + // --- Seafloor aging -> depth (Phase 2 inc. 4) --------------------------- + // Oceanic crust subsides as it ages (half-space cooling): depth = + // ridgeDepth - seafloorSubsidence * sqrt(geoAge), clamped at oceanBase. + double seafloorSubsidence = 280.0; // m per sqrt(My) of crustal age + double seafloorSeedAge = 80.0; // My, initial oceanic age spread at generation + + // --- Plate dynamics (Phase 2): fission, stalemate kick, spreading plates -- + // Periodic checks run every `splitCheckEvery` drift iterations. A plate over + // `splitFraction` of all cells may rift in two with probability + // splitProbBase + splitProbSlope * (percentOverThreshold) + // (20%->5%, 21%->10%, ... 39%->100%). A plate whose cell count barely changed + // over a window gets a stalemate "kick" (new direction + a speed boost). Young + // rift crust grows on a "baby" proto-plate; once it reaches `babyPromoteFrac` + // of all cells it becomes a real plate and grows volcanic landmass. + // (Related spreading/volcanic fields are at bottom for binary save compat.) + int splitCheckEvery = 10; // drift iterations between periodic checks + double splitFraction = 0.20; // plate share of cells that may rift apart + double splitProbBase = 0.05; // split probability at the threshold + double splitProbSlope = 0.05; // added per 1 percentage-point over threshold + double stalemateEps = 0.005; // |dCells|/cells below this over a window = stuck + int stalemateWindows = 4; // consecutive stuck windows required before a kick + double stalemateBoost = 1.5; // speed multiplier when kicking a stuck plate + int miniPlateCells = 50; // a non-baby plate smaller than this is "mini" + int fuseMinPlates = 3; // distinct mini plates in a cluster to fuse + steal + + // --- Erosion + sea level (Phase 2 increment 2) -------------------------- + // erode() moves sediment downhill (highs wear down, basins/seas fill); a + // proportional sea-level controller holds a target geographic land fraction. + double erosionLandRate = 0.08; // subaerial erosion fraction / My + double erosionSeaRate = 0.02; // submarine erosion fraction / My (slower) + double landFractionTarget = 0.30; // geographic land goal (cells above seaLevel) + double seaLevelStep = 100.0; // m, fixed nudge per adjustment when off target + double seaLevelTol = 0.02; // deadband (land-fraction) where sea level rests + int seaLevelEvery = 100; // erode calls between sea-level adjustments + + // --- Spreading & volcanic (Phase 2 plate dynamics, kept here for binary compat) + int babyMinCells = 4; // baby blobs smaller than this dissolve (noise) + double babyPromoteFrac = 0.007; // baby-patch size (x N cells) to become a plate + double volcanicLandFrac = 0.30; // fraction of a promoted patch turned into land + double volcanicElev = 400.0; // m, volcanic-island starting elevation + double landBand = 0.10; // soft land clamp: +/- around targetLand + + // --- Phase 3: hydrology (rivers, lakes, fluvial erosion) ---------------- + // Macro drainage network on the fixed grid: depression-fill -> lakes, + // steepest-descent routing -> rivers, stream-power incision + downstream + // sediment transport/deposition (mass-conserving). Drift keeps running but + // Phase 3 uses a finer timestep (cflDtMy * phase3DtScale). + double phase3AfterMy = 300.0; // My of drift before the Phase-3 prompt + double phase3DtScale = 0.2; // Phase-3 timestep = cflDtMy() * this (finer) + double rainfall = 1.0; // uniform precip per cell (drainage-area unit) + double riverThreshold = 50.0; // discharge above which a cell counts as a river + double riverIncision = 0.02; // K in stream-power incision K*Q^m*S^n*dt + double riverDischargeExp = 0.5; // m: discharge exponent in stream power + double riverSlopeExp = 1.0; // n: slope exponent in stream power + double riverTransport = 0.10; // transport-capacity coefficient (cap=this*Q*S) + double depFrac = 0.25; // fraction of excess load deposited per cell + + // --- Phase 3: biome classification thresholds (see PlanetBiomes.cpp) ----- + // Temperature model (deg C): equator-warm curve cooling super-linearly toward the + // poles minus an elevation lapse. Moisture comes from latitude belts + hydrology. + double biomeEquatorTemp = 30.0; // C at the equator, sea level + double biomePoleDrop = 58.0; // C drop from equator to pole + double biomeLatExp = 1.3; // >1 keeps mid-latitudes temperate (cold near poles) + double biomeElevLapse = 0.0060; // C lost per metre above sea level + double biomeIceTemp = -9.5; // below -> Ice (polar caps + glaciated peaks); raise = bigger caps + double biomeTundraTemp = 2.0; // below (and above ice) -> Tundra/Taiga + double biomeTaigaTemp = 10.0; // cool + wet -> boreal forest + double biomeSavannaTemp = 22.0; // warm + moderate moisture -> savanna + double biomeMountainElev= 3000.0; // m above sea level -> Mountains + double biomeHillsElev = 1200.0; // m above sea level -> Hills + double biomeBeachBand = 60.0; // m above sea level + adjacent ocean -> Beach + double biomeLowlandElev = 500.0; // wetlands only below this elevation + double biomeWetlandMoist= 0.72; // moisture above this (low lowland) -> Wetland + double biomeDesertMoist = 0.28; // moisture below this -> Desert + double biomeGrassMoist = 0.50; // moisture below this -> Grassland/Savanna, else Forest + double biomeTaigaMoist = 0.40; // cool + above this -> Taiga (else Tundra) + double biomeLakeMinDepth= 20.0; // filled-basin depth above sea level counting as a Lake + + // --- Phase 3: climate (orographic precipitation) -- see PlanetClimate.cpp -- + // Temperature reuses the biome* temperature fields above. Precipitation advects + // ocean moisture along prevailing (zonal) winds: it rains on windward upslopes and + // dries out leeward (rain shadow) and far inland (continentality). + double climateOceanMoisture = 1.0; // moisture air carries leaving the ocean (source) + double climateRainEfficiency = 0.5; // fraction of available moisture*belt that rains per cell + double climateOrographic = 3.0; // extra rain per unit normalized upslope (windward) + double climateOroRefHeight = 500.0; // m of upslope that counts as one orographic unit + double climateContinentality = 0.05; // moisture lost per land cell crossed (dries interiors) + int climateWindPasses = 50; // moisture-advection iterations (steady state) + int climateMoistureSmooth = 12; // precipitation diffusion passes (wet/dry transition zones) +}; diff --git a/src/sim/Projection.hpp b/src/sim/Projection.hpp new file mode 100644 index 0000000..667db7c --- /dev/null +++ b/src/sim/Projection.hpp @@ -0,0 +1,71 @@ +#pragma once +#include "Vec3.hpp" +#include +#include + +// Equal Earth projection (Savric, Patterson & Jenny, 2018): an equal-area +// pseudocylindrical projection with a natural, low-distortion "globe" look. +// Forward is closed-form; the inverse uses a few Newton steps (cheap, used per +// mouse point for hover). Kept raylib-free so it can be unit-tested headless. +// +// Convention: y is "up". lat = asin(unit.y), lon = atan2(unit.z, unit.x). +namespace EqualEarth { + +constexpr double A1 = 1.340264; +constexpr double A2 = -0.081106; +constexpr double A3 = 0.000893; +constexpr double A4 = 0.003796; + +inline double M_() { return std::sqrt(3.0) / 2.0; } // sin(theta) scale +inline double K_() { return 2.0 * std::sqrt(3.0) / 3.0; } // x scale factor + +// dy/dtheta == x denominator; shared by forward and inverse. +inline double denom(double th2, double th6) { + return A1 + 3.0 * A2 * th2 + th6 * (7.0 * A3 + 9.0 * A4 * th2); +} + +// (lon,lat) radians -> projection (x,y). +inline void forward(double lon, double lat, double& x, double& y) { + double th = std::asin(M_() * std::sin(lat)); + double th2 = th * th, th6 = th2 * th2 * th2; + x = K_() * lon * std::cos(th) / denom(th2, th6); + y = th * (A1 + A2 * th2 + th6 * (A3 + A4 * th2)); // A1 th + A2 th^3 + A3 th^7 + A4 th^9 +} + +// projection (x,y) -> (lon,lat) radians. Returns false if (x,y) is outside the +// projected globe (so callers can reject hovers off the map). +inline bool inverse(double x, double y, double& lon, double& lat) { + double th = y; // good initial guess (y ~ A1*theta near 0) + for (int it = 0; it < 16; ++it) { + double th2 = th * th, th6 = th2 * th2 * th2; + double fy = th * (A1 + A2 * th2 + th6 * (A3 + A4 * th2)) - y; + double d = fy / denom(th2, th6); + th -= d; + if (std::fabs(d) < 1e-12) break; + } + double s = std::sin(th) / M_(); + if (s < -1.0 - 1e-9 || s > 1.0 + 1e-9) return false; + lat = std::asin(std::clamp(s, -1.0, 1.0)); + double th2 = th * th, th6 = th2 * th2 * th2; + double c = std::cos(th); + if (std::fabs(c) < 1e-12) return false; + lon = x * denom(th2, th6) / (K_() * c); + if (lon < -M_PI - 1e-6 || lon > M_PI + 1e-6) return false; + return true; +} + +// Half-extents of the projected map (x in [-halfW,halfW], y in [-halfH,halfH]). +inline double halfWidth() { double x, y; forward(M_PI, 0.0, x, y); return x; } +inline double halfHeight() { double x, y; forward(0.0, M_PI / 2.0, x, y); return y; } + +} // namespace EqualEarth + +// Sphere direction <-> geographic coordinates (y up). +inline void dirToLonLat(const Vec3& u, double& lon, double& lat) { + lat = std::asin(std::clamp(u.y, -1.0, 1.0)); + lon = std::atan2(u.z, u.x); +} +inline Vec3 lonLatToDir(double lon, double lat) { + double c = std::cos(lat); + return Vec3{ c * std::cos(lon), std::sin(lat), c * std::sin(lon) }; +} diff --git a/src/sim/Vec3.hpp b/src/sim/Vec3.hpp new file mode 100644 index 0000000..3db5728 --- /dev/null +++ b/src/sim/Vec3.hpp @@ -0,0 +1,30 @@ +#pragma once +#include + +// Simple 3D vector for planet geometry (metric, double precision). +struct Vec3 { + double x = 0.0, y = 0.0, z = 0.0; + + Vec3() = default; + Vec3(double x_, double y_, double z_) : x(x_), y(y_), z(z_) {} + + Vec3 operator+(const Vec3& o) const { return {x + o.x, y + o.y, z + o.z}; } + Vec3 operator-(const Vec3& o) const { return {x - o.x, y - o.y, z - o.z}; } + Vec3 operator*(double s) const { return {x * s, y * s, z * s}; } + + double dot(const Vec3& o) const { return x * o.x + y * o.y + z * o.z; } + + Vec3 cross(const Vec3& o) const { + return {y * o.z - z * o.y, + z * o.x - x * o.z, + x * o.y - y * o.x}; + } + + double length() const { return std::sqrt(x * x + y * y + z * z); } + + Vec3 normalized() const { + double l = length(); + if (l <= 1e-300) return {0, 0, 0}; + return {x / l, y / l, z / l}; + } +}; diff --git a/test_logic.cpp b/test_logic.cpp new file mode 100644 index 0000000..fbdc98b --- /dev/null +++ b/test_logic.cpp @@ -0,0 +1,156 @@ +// Headless logic test for Phase 1 tectonics. No display / raylib needed. +// +// g++ -std=c++17 -O2 -Isrc/sim test_logic.cpp src/sim/IcoSphere.cpp \ +// src/sim/Planet.cpp src/sim/PlanetTectonics.cpp src/sim/PlanetDrift.cpp \ +// src/sim/PlanetErosion.cpp src/sim/PlanetHydrology.cpp src/sim/PlanetIO.cpp \ +// -o /tmp/t && /tmp/t +// +// Verifies the invariants documented in CLAUDE.md so Planet::step() and the +// icosphere can be changed with confidence without launching the window. + +#include "Planet.hpp" +#include "Projection.hpp" +#include +#include +#include +#include + +static int failures = 0; + +static double angBetweenDeg(const Vec3& a, const Vec3& b) { + return std::acos(std::clamp(a.dot(b), -1.0, 1.0)) * 180.0 / M_PI; +} + +static void check(bool cond, const char* what) { + std::printf(" [%s] %s\n", cond ? "PASS" : "FAIL", what); + if (!cond) ++failures; +} + +// Euler characteristic for the icosphere: V - E + F == 2 (sphere topology). +static bool eulerOk(const Planet& p, size_t vertCount) { + const std::vector& tri = p.triIndices(); + size_t F = tri.size() / 3; + // Each interior edge is shared by exactly 2 triangles; count unique edges. + std::map, int> edges; + for (size_t k = 0; k + 2 < tri.size(); k += 3) { + int a = tri[k], b = tri[k + 1], c = tri[k + 2]; + int e[3][2] = {{a, b}, {b, c}, {c, a}}; + for (auto& pr : e) { + int lo = std::min(pr[0], pr[1]), hi = std::max(pr[0], pr[1]); + edges[{lo, hi}]++; + } + } + size_t E = edges.size(); + long long euler = (long long)vertCount - (long long)E + (long long)F; + std::printf(" V=%zu E=%zu F=%zu V-E+F=%lld\n", vertCount, E, F, euler); + return euler == 2; +} + +int main() { + // ---- Geometry: icosphere level 5 ------------------------------------- + Planet planet; + PlanetConfig cfg; + cfg.subdivisions = 5; + cfg.seed = 1337; + planet.generate(cfg); + + std::printf("Geometry (level 5):\n"); + check(planet.cells.size() == 10242, "10242 cells at subdivision 5"); + check(eulerOk(planet, planet.cells.size()), "Euler characteristic V-E+F == 2"); + + // Vertex degrees: exactly 12 should have degree 5, the rest degree 6. + int deg5 = 0, deg6 = 0, other = 0; + for (auto& c : planet.cells) { + if (c.neighbors.size() == 5) ++deg5; + else if (c.neighbors.size() == 6) ++deg6; + else ++other; + } + std::printf(" degree5=%d degree6=%d other=%d\n", deg5, deg6, other); + check(deg5 == 12, "exactly 12 degree-5 vertices"); + check(other == 0, "all remaining vertices are degree 6"); + + double cw = planet.cellWidthMeters() / 1000.0; + std::printf(" cell width ~ %.0f km\n", cw); + check(cw > 200.0 && cw < 250.0, "cell width ~223 km"); + + // ---- Plates ---------------------------------------------------------- + bool allAssigned = true; + for (auto& c : planet.cells) + if (c.plateId < 0 || c.plateId >= cfg.plateCount) allAssigned = false; + check(allAssigned, "every cell assigned to a valid plate"); + check((int)planet.plates.size() == cfg.plateCount, "plateCount plates created"); + + // ---- Tectonics: run ~40 ticks ---------------------------------------- + for (int i = 0; i < 40; ++i) planet.step(); + + double lo = planet.minElevation(), hi = planet.maxElevation(); + std::printf("Tectonics after 40 ticks:\n"); + std::printf(" elevation range %.0f .. %.0f m\n", lo, hi); + check(hi > 2000.0, "clear mountains form (max > 2000 m)"); + check(lo < -6000.0, "deep trenches form (min < -6000 m)"); + + // Relief must be GRADED, not saturated to the clamp rails (the old bug: + // uplift so strong every boundary cell railed to +/- the clamp in 1 tick). + int pinned = 0, midLand = 0, midSea = 0; + for (auto& c : planet.cells) { + if (c.elevation >= 8999.0 || c.elevation <= -10999.0) ++pinned; + if (c.elevation > 800.0 && c.elevation < 2000.0) ++midLand; // belt flanks + if (c.elevation < -4500.0 && c.elevation > -6000.0) ++midSea; // trench flanks + } + double pinnedPct = 100.0 * pinned / planet.cells.size(); + std::printf(" pinned to clamp: %d (%.2f%%) flank cells: land=%d sea=%d\n", + pinned, pinnedPct, midLand, midSea); + check(pinnedPct < 2.0, "not saturated: <2% of cells pinned to clamp rails"); + check(midLand > 0 && midSea > 0, "graded relief: mountains/trenches have flanks"); + + bool finite = true; + for (auto& c : planet.cells) + if (!std::isfinite(c.elevation)) finite = false; + check(finite, "no NaN/Inf elevations (numerically stable)"); + + // ---- Determinism: same seed -> identical result ---------------------- + Planet p2; + p2.generate(cfg); + for (int i = 0; i < 40; ++i) p2.step(); + bool identical = (p2.cells.size() == planet.cells.size()); + for (size_t i = 0; identical && i < p2.cells.size(); ++i) + if (p2.cells[i].elevation != planet.cells[i].elevation) identical = false; + check(identical, "same seed reproduces identical world (deterministic)"); + + // ---- Equal Earth projection round-trip (used by the 2D map + hover) --- + { + double maxErr = 0.0; int rejected = 0; + for (const auto& c : planet.cells) { + double lon, lat; dirToLonLat(c.unit, lon, lat); + double x, y; EqualEarth::forward(lon, lat, x, y); + double lo2, la2; + if (!EqualEarth::inverse(x, y, lo2, la2)) { ++rejected; continue; } + Vec3 d2 = lonLatToDir(lo2, la2); + maxErr = std::max(maxErr, angBetweenDeg(c.unit, d2)); + } + std::printf(" Equal Earth round-trip: max err %.2e deg, rejected %d\n", maxErr, rejected); + check(maxErr < 1e-3 && rejected == 0, "Equal Earth forward/inverse round-trips"); + } + + // ---- Subgrid: continuity at center + neighbor coverage ---------------- + { + int cell = (int)planet.cells.size() / 2; + auto sg = planet.makeSubGrid(cell, 16); + bool sgFinite = true; int ownCell = 0, ownNbr = 0; + for (auto& s : sg->sub) { + if (!std::isfinite(s.elevation)) sgFinite = false; + if (s.nearestMacro == cell) ++ownCell; else ++ownNbr; + } + const SubCell& ctr = sg->sub[(16 / 2) * 16 + 16 / 2]; + double centerErr = std::fabs(ctr.elevation - planet.cells[cell].elevation); + std::printf(" subgrid: %zu subcells, center |diff| %.0f m, own %d / nbr %d\n", + sg->sub.size(), centerErr, ownCell, ownNbr); + check(sgFinite && sg->sub.size() == 256 && centerErr < 800.0 && ownCell > 0 && ownNbr > 0, + "subgrid is finite, continuous at center, and reaches neighbors"); + } + + std::printf("\n%s (%d failure%s)\n", + failures == 0 ? "ALL TESTS PASSED" : "TESTS FAILED", + failures, failures == 1 ? "" : "s"); + return failures == 0 ? 0 : 1; +}