planetsim/src/sim/Planet.cpp

#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
    computeBiotaDensity();   // derived flora/fauna/funga density (population is on-demand)
}

// 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];
    }
    sBiota.assign(cells.size(), {});   // empty biota population until generateBiota()
    sHasBiota = false;
}

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