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 <noreply@anthropic.com>
This commit is contained in:
Jonas Reith 2026-06-23 15:08:25 +02:00
commit acc0e5eec9
36 changed files with 4498 additions and 0 deletions

21
.gitignore vendored Normal file
View File

@ -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

168
BUILD.md Normal file
View File

@ -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<id> 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.

520
CLAUDE.md Normal file
View File

@ -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<id>` 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> 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<id>` 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.

54
CMakeLists.txt Normal file
View File

@ -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()

94
docs/design-notes.md Normal file
View File

@ -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).

82
docs/fauna-flora-plan.md Normal file
View File

@ -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<double>& vegetation()/herbivores()/carnivores() const`. Add scratch
members `std::vector<double> 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.

11
src/main.cpp Normal file
View File

@ -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;
}

123
src/render/Colors.cpp Normal file
View File

@ -0,0 +1,123 @@
#include "Colors.hpp"
#include <algorithm>
#include <cmath>
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 };
}

23
src/render/Colors.hpp Normal file
View File

@ -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);

69
src/render/Map2D.cpp Normal file
View File

@ -0,0 +1,69 @@
#include "Map2D.hpp"
#include "rlgl.h"
#include "Projection.hpp"
#include <algorithm>
#include <cmath>
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<Color>& 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<int>& 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();
}

23
src/render/Map2D.hpp Normal file
View File

@ -0,0 +1,23 @@
#pragma once
#include "raylib.h"
#include "Planet.hpp"
#include <vector>
// Equal Earth 2D map: per-cell screen positions + the projection/draw helpers.
struct Map2D {
std::vector<Vector2> pos; // screen position per cell (fast path)
std::vector<double> lon; // longitude per cell (radians, for seam detection)
std::vector<double> 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<Color>& vc,
const Map2D& m, Rectangle r, double lonOffset);

197
src/render/Overlays.cpp Normal file
View File

@ -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 <algorithm>
#include <cmath>
#include <cstdlib> // std::abs(int)
void buildBorders(const Planet& p, float radius,
std::vector<Vector3>& real, std::vector<Vector3>& 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<Vector3>& out = (isBaby(pA) || isBaby(pB)) ? ridge : real;
out.push_back(a); out.push_back(b);
};
const std::vector<int>& 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<Vector3>& out, std::vector<PlateLabel>& labels) {
out.clear(); labels.clear();
int np = (int)p.plates.size();
if (np == 0) return;
std::vector<Vec3> sum(np, Vec3{0, 0, 0});
std::vector<int> 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<Vec3> 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<Vector3>& rivers, std::vector<Vector3>& bigRivers) {
rivers.clear(); bigRivers.clear();
const std::vector<double>& disc = p.discharge();
const std::vector<int>& 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<std::vector<Vector2>> buildGraticule() {
std::vector<std::vector<Vector2>> g;
const double D = M_PI / 180.0;
for (int la = -60; la <= 60; la += 30) { // parallels
std::vector<Vector2> 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<Vector2> 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<std::vector<Vector2>>& 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<std::vector<Vector2>>& 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<Vector3>& 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<std::shared_ptr<SubGrid>>& 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();
}
}

43
src/render/Overlays.hpp Normal file
View File

@ -0,0 +1,43 @@
#pragma once
#include "raylib.h"
#include "Planet.hpp"
#include <vector>
#include <memory>
// 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<Vector3>& real, std::vector<Vector3>& ridge);
// ---- Per-plate drift arrows -------------------------------------------------
struct PlateLabel { int id; Vector3 pos; };
void buildDriftArrows(const Planet& p, float radius,
std::vector<Vector3>& out, std::vector<PlateLabel>& 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<Vector3>& rivers, std::vector<Vector3>& bigRivers);
// ---- Lat/lon graticule ------------------------------------------------------
// Polylines of (lon,lat) radians (Vector2.x=lon, .y=lat).
std::vector<std::vector<Vector2>> buildGraticule();
void drawGraticule3D(const std::vector<std::vector<Vector2>>& g, float radius);
void drawGraticule2D(const std::vector<std::vector<Vector2>>& 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<Vector3>& segs, Color col, float width,
Rectangle r, double lonOffset);
// ---- Subgrid overlay (high-res patch over a selected cell) ------------------
void drawSubgrids(const std::vector<std::shared_ptr<SubGrid>>& sgs,
float visBase, float elevExagg, double seaLevel, float eps);

162
src/render/Panels.cpp Normal file
View File

@ -0,0 +1,162 @@
#include "Panels.hpp"
#include "Colors.hpp" // elevationColor (subtile grid)
#include "Projection.hpp" // dirToLonLat
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>
// elev/age come from the display snapshot so the readout matches what is drawn.
static std::vector<std::string> 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<double>& v) { return (int)v.size() == n; };
std::vector<std::string> 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<SubGrid>& 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<int> 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<int> 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});
}

20
src/render/Panels.hpp Normal file
View File

@ -0,0 +1,20 @@
#pragma once
#include "raylib.h"
#include "Planet.hpp"
#include <memory>
// 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<SubGrid>& 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);

47
src/render/Picking.cpp Normal file
View File

@ -0,0 +1,47 @@
#include "Picking.hpp"
#include <algorithm>
#include <cmath>
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 };
}

19
src/render/Picking.hpp Normal file
View File

@ -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);

254
src/render/Viewer.cpp Normal file
View File

@ -0,0 +1,254 @@
#include "Viewer.hpp"
#include "Picking.hpp" // angBetween (rebuildSub)
#include "Projection.hpp" // EqualEarth (layout)
#include <algorithm>
#include <cmath>
#include <cstring>
#include <fstream>
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<size_t>(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<double>& temp = planet.temperature();
const std::vector<double>& 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<double>& 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<const char*>(&ver), sizeof ver);
os.write(reinterpret_cast<const char*>(&elapsedMy), sizeof elapsedMy);
os.write(reinterpret_cast<const char*>(&st), sizeof st);
os.write(reinterpret_cast<const char*>(&driftRate), sizeof driftRate);
os.write(reinterpret_cast<const char*>(&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<char*>(&ver), sizeof ver);
is.read(reinterpret_cast<char*>(&em), sizeof em);
is.read(reinterpret_cast<char*>(&st), sizeof st);
if (ver >= 2) is.read(reinterpret_cast<char*>(&dr), sizeof dr);
if (ver >= 3) is.read(reinterpret_cast<char*>(&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();
}

124
src/render/Viewer.hpp Normal file
View File

@ -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 <vector>
#include <string>
#include <memory>
// 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<Color> vcolors;
std::vector<Vector3> borders, ridgeBorders; bool showBorders = true;
std::vector<Vector3> driftArrows; bool showDrift = true;
std::vector<PlateLabel> plateLabels;
std::vector<Vector3> rivers, bigRivers; bool showRivers = true;
std::vector<std::vector<Vector2>> 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<std::shared_ptr<SubGrid>> 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();
};

143
src/render/ViewerInput.cpp Normal file
View File

@ -0,0 +1,143 @@
#include "Viewer.hpp"
#include "Picking.hpp"
#include "Map2D.hpp" // wrapPi
#include "Projection.hpp" // EqualEarth, lonLatToDir
#include <algorithm>
#include <cmath>
// 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);
}

262
src/render/ViewerRender.cpp Normal file
View File

@ -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 <algorithm>
#include <cmath>
#include <vector>
// 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<int>& 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<Vector3>& 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();
}

70
src/sim/IcoSphere.cpp Normal file
View File

@ -0,0 +1,70 @@
#include "IcoSphere.hpp"
#include <set>
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<Vec3> 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<int> 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<int> 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<std::set<int>> 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());
}

25
src/sim/IcoSphere.hpp Normal file
View File

@ -0,0 +1,25 @@
#pragma once
#include "Vec3.hpp"
#include <vector>
#include <cstdint>
#include <unordered_map>
// 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<Vec3> positions; // unit-sphere positions per vertex
std::vector<std::vector<int>> neighbors; // adjacency per vertex
std::vector<int> 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<int64_t, int> midpointCache;
};

233
src/sim/Planet.cpp Normal file
View File

@ -0,0 +1,233 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
// --- 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<int> 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<SubGrid> Planet::makeSubGrid(int cellIndex, int res) const {
auto sg = std::make_shared<SubGrid>();
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<int> 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<size_t>(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;
}

139
src/sim/Planet.hpp Normal file
View File

@ -0,0 +1,139 @@
#pragma once
#include "Vec3.hpp"
#include "IcoSphere.hpp"
#include "PlanetTypes.hpp" // Cell, Plate, SubGrid/SubCell, PlanetConfig
#include <vector>
#include <memory>
#include <cstdint>
#include <string>
#include <iosfwd>
class Planet {
public:
PlanetConfig cfg;
std::vector<Cell> cells;
std::vector<Plate> 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<double>& temperature() const { return sTemp; } // deg C
const std::vector<double>& precipitation() const { return sPrecip; } // relative units
const std::vector<double>& 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<double>& lakeDepth() const { return sLakeDepth; }
const std::vector<double>& discharge() const { return sDischarge; }
const std::vector<int>& flowTo() const { return sFlowTo; }
// Build a fine-resolution subgrid patch for one macro cell (phase 4/5 hook).
std::shared_ptr<SubGrid> 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<int>& 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<int>& cnt);
void kickStalemates(const std::vector<int>& 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<int>& 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<int> sPrevCount; // per-plate cell count at the previous check
std::vector<int> sStaleStreak; // consecutive stuck windows per plate
std::vector<int> sFreePlateIds; // dead plate slots free for reuse
// Reusable scratch buffers for step()/erode() so they allocate nothing per tick.
std::vector<double> sStress, sBelt, sBeltNext, sDelta, sSmoothed, sOldElev, sErode;
std::vector<uint8_t> sSub, sOver, sColl;
// Phase-3 hydrology scratch (derived from elevation each routeFlow(); not saved).
std::vector<double> sFill, sLakeDepth, sDischarge;
std::vector<int> sFlowTo, sHydroOrder;
// Phase-3 climate scratch (derived each computeClimate(); not saved). sMoist is the
// 0..1-normalized precipitation the biome classifier reads.
std::vector<double> sTemp, sPrecip, sMoist;
std::vector<Vec3> sWind;
std::vector<int> 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);

55
src/sim/PlanetBiomes.cpp Normal file
View File

@ -0,0 +1,55 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
// --- 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;
}
}

134
src/sim/PlanetClimate.cpp Normal file
View File

@ -0,0 +1,134 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
#include <vector>
// --- 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<double> 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<double> 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<double> 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<double> 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);
}

347
src/sim/PlanetDrift.cpp Normal file
View File

@ -0,0 +1,347 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
#include <queue>
// --- 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<size_t>(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<int> oP(n); std::vector<double> oE(n), oA(n); std::vector<uint8_t> oO(n);
std::vector<Vec3> 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<int> p2(n); std::vector<double> e2(n), a2(n); std::vector<uint8_t> 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<uint8_t> 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<int,int>{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<int> 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<int>& 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<int> 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<int>& 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<uint8_t> 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<int> surr(P, -2); // -2 none seen, -1 multiple, >=0 the single other
std::vector<int> 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<int> 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<int> 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<char> seen(n, 0);
std::queue<int> q;
for (int i = 0; i < n; ++i) {
if (seen[i] || !isMini(cells[i].plateId)) continue;
std::vector<int> 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<char> seen(n, 0);
auto isBaby = [&](int i){ return cells[i].plateId >= 0 && plates[cells[i].plateId].baby; };
std::queue<int> q;
for (int i = 0; i < n; ++i) {
if (seen[i] || !isBaby(i)) continue;
std::vector<int> 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<int>& 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<int> 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
}
}
}

58
src/sim/PlanetErosion.cpp Normal file
View File

@ -0,0 +1,58 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
// --- 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
}

104
src/sim/PlanetHydrology.cpp Normal file
View File

@ -0,0 +1,104 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
#include <queue>
// --- 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<PQItem, std::vector<PQItem>, PQGreater> pq;
std::vector<uint8_t> 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<double> 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);
}

263
src/sim/PlanetIO.cpp Normal file
View File

@ -0,0 +1,263 @@
#include "Planet.hpp"
#include <algorithm>
#include <string>
#include <vector>
#include <type_traits>
#include <istream>
#include <ostream>
#include <fstream>
#include <sstream>
// --- 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<std::string> 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 <class T> void writePod(std::ostream& os, const T& v) {
static_assert(std::is_trivially_copyable<T>::value, "writePod needs a POD type");
os.write(reinterpret_cast<const char*>(&v), sizeof(T));
}
template <class T> void readPod(std::istream& is, T& v) {
static_assert(std::is_trivially_copyable<T>::value, "readPod needs a POD type");
is.read(reinterpret_cast<char*>(&v), sizeof(T));
}
template <class T> void writeVec(std::ostream& os, const std::vector<T>& v) {
static_assert(std::is_trivially_copyable<T>::value, "writeVec needs POD elements");
uint64_t n = v.size(); writePod(os, n);
if (n) os.write(reinterpret_cast<const char*>(v.data()), (std::streamsize)(n * sizeof(T)));
}
template <class T> void readVec(std::istream& is, std::vector<T>& v) {
static_assert(std::is_trivially_copyable<T>::value, "readVec needs POD elements");
uint64_t n = 0; readPod(is, n); v.resize((size_t)n);
if (n) is.read(reinterpret_cast<char*>(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;
}

164
src/sim/PlanetTectonics.cpp Normal file
View File

@ -0,0 +1,164 @@
#include "Planet.hpp"
#include <algorithm>
#include <cmath>
// --- 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<double>& stress = sStress;
std::vector<uint8_t>& subducting = sSub;
std::vector<uint8_t>& overriding = sOver;
std::vector<uint8_t>& 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<double>& belt = sBelt; belt.resize(n);
std::vector<double>& 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<double>& 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<double>& 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;
}

190
src/sim/PlanetTypes.hpp Normal file
View File

@ -0,0 +1,190 @@
#pragma once
#include "Vec3.hpp"
#include <vector>
#include <memory>
#include <cstdint>
// 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<SubCell> 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<int> neighbors;
// Phase 4/5 hook: optional fine-resolution subgrid. Null until needed.
std::shared_ptr<SubGrid> 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)
};

71
src/sim/Projection.hpp Normal file
View File

@ -0,0 +1,71 @@
#pragma once
#include "Vec3.hpp"
#include <cmath>
#include <algorithm>
// 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) };
}

30
src/sim/Vec3.hpp Normal file
View File

@ -0,0 +1,30 @@
#pragma once
#include <cmath>
// 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};
}
};

156
test_logic.cpp Normal file
View File

@ -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 <cstdio>
#include <map>
#include <cmath>
#include <algorithm>
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<int>& tri = p.triIndices();
size_t F = tri.size() / 3;
// Each interior edge is shared by exactly 2 triangles; count unique edges.
std::map<std::pair<int,int>, 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;
}