planetsim/src/sim/Planet.hpp

#pragma once
#include "Vec3.hpp"
#include "IcoSphere.hpp"
#include "PlanetTypes.hpp"   // Cell, Plate, SubGrid/SubCell, PlanetConfig
#include "PlanetBiota.hpp"   // BiotaKind, Organism, CellBiota
#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, annual mean
    const std::vector<double>& precipitation() const { return sPrecip; } // relative units
    const std::vector<double>& moisture()      const { return sMoist; }  // 0..1 (median land -> 0.5)
    const std::vector<double>& summerTemp()    const { return sTempSummer; } // deg C, warmest month
    const std::vector<double>& winterTemp()    const { return sTempWinter; } // deg C, coldest month

    // 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();

    // Biota stage (flora/fauna/funga). computeBiotaDensity() builds the derived
    // per-cell density scalars (0..1) each tick (like climate; not saved); call it
    // after classifyBiomes(). generateBiota() does the on-demand slot/point fill of
    // the discrete population into sBiota (saved) -- NOT called per tick. See
    // PlanetBiota.cpp + PlanetFloraGen/FaunaGen/FungiGen.cpp.
    void   computeBiotaDensity();
    void   generateBiota();
    bool   biotaPopulated() const;
    const std::vector<double>& floraDensity() const { return sFloraDensity; }
    const std::vector<double>& faunaDensity() const { return sFaunaDensity; }
    const std::vector<double>& fungaDensity() const { return sFungaDensity; }
    const std::vector<CellBiota>& biota() const { return sBiota; }

    // 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.
    // hasBiota: whether the stream carries the biota population block (save v7+).
    bool readState(std::istream& is, bool hasBiome = true, bool hasBiota = 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

    // Biota helpers (PlanetFloraGen/FaunaGen/FungiGen.cpp). compute*Density write the
    // derived scalars; fill* draw the per-cell population (nbr = already-filled,
    // same-biome neighbours, for regional consistency).
    void  computeFloraDensity();
    void  computeFaunaDensity();
    void  computeFungaDensity();
    double neighbourhoodPrey(int i) const;       // mean fauna density over i + neighbours
    std::vector<Organism> fillFlora(int i, const std::vector<int>& nbr, uint32_t& rng);
    std::vector<Organism> fillFauna(int i, const std::vector<int>& nbr, uint32_t& rng);
    std::vector<Organism> fillFunga(int i, const std::vector<int>& nbr, uint32_t& rng);

    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. sTempSummer/sTempWinter are
    // the obliquity-driven seasonal extremes around the annual mean sTemp (see Seasons).
    std::vector<double>  sTemp, sPrecip, sMoist, sTempSummer, sTempWinter;
    std::vector<Vec3>    sWind;
    std::vector<int>     sUpwind;

    // Biota: derived density scalars (0..1; recomputed each tick, not saved) and the
    // on-demand discrete population (saved). sHasBiota latches once generated/loaded.
    std::vector<double>     sFloraDensity, sFaunaDensity, sFungaDensity;
    std::vector<CellBiota>  sBiota;
    bool                    sHasBiota = false;
};

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