planetsim/src/sim/PlanetWeather.cpp

#include "Planet.hpp"
#include <algorithm>
#include <cmath>
#include <vector>

// Live World weather: a dynamic per-cell humidity / cloud / rain cycle advanced on the live
// clock (geometry fixed -- these are fields flowed over the grid, like climate, but time-varying).
// One step: evaporate over warm sunlit seas -> advect humidity & cloud along the prevailing wind
// -> condense the supersaturated air into cloud (extra on windward upslopes) -> rain out the
// thick cloud -> dissipate. Reads the static climate scaffolding (sTemp/sWind/sUpwind/sMoist set
// by computeClimate) and the live sInsolation (computeInsolation). Deterministic; saved (v10).

void Planet::initWeather() {
    const int n = (int)cells.size();
    sHumidity.assign(n, 0.0);
    sCloud.assign(n, 0.0);
    sRain.assign(n, 0.0);
    const double sea = cfg.seaLevel;
    const bool haveM = ((int)sMoist.size() == n);
    for (int i = 0; i < n; ++i) {
        if (cells[i].elevation <= sea) sHumidity[i] = 0.9;                       // saturated marine air
        else sHumidity[i] = haveM ? (0.2 + 0.5 * sMoist[i]) : 0.3;              // land: from climatology
    }
    sStorms.clear(); sStormNextId = 1;
    sWeatherRng = cfg.seed ? (cfg.seed ^ 0x5701A123u) : 0x5701A123u;            // separate RNG
    sHasWeather = true;
}

WeatherSnapshot Planet::captureWeather() const {
    WeatherSnapshot s;
    s.humidity = sHumidity; s.cloud = sCloud; s.rain = sRain;
    s.storms = sStorms; s.rng = sWeatherRng; s.nextId = sStormNextId;
    s.volcanoes = volcanoes; s.volRng = sVolRng;
    s.settlementPop.reserve(settlements.size());          // civ: only population is mutable
    for (const Settlement& st : settlements) s.settlementPop.push_back(st.population);
    s.settlementAllegiance = sSettleAllegiance;           // civ Step 5: conquest state
    s.wars = wars; s.warRng = sWarRng; s.warNextId = sWarNextId;
    s.diplomacy = diplomacy;                              // civ Step 6: standing realm relations
    return s;
}

void Planet::restoreWeather(const WeatherSnapshot& s) {
    sHumidity = s.humidity; sCloud = s.cloud; sRain = s.rain;
    sStorms = s.storms; sWeatherRng = s.rng; sStormNextId = s.nextId;
    volcanoes = s.volcanoes; sVolRng = s.volRng;
    if (s.settlementPop.size() == settlements.size())     // restore populations (set is fixed)
        for (size_t k = 0; k < settlements.size(); ++k) settlements[k].population = s.settlementPop[k];
    sSettleAllegiance = s.settlementAllegiance;            // civ Step 5: restore conquest state
    if (sSettleAllegiance.size() != settlements.size()) sSettleAllegiance.assign(settlements.size(), -1);
    wars = s.wars; sWarRng = s.warRng ? s.warRng : sWarRng; sWarNextId = s.warNextId ? s.warNextId : sWarNextId;
    diplomacy = s.diplomacy;                              // civ Step 6: restore realm relations
    sHasWeather = !sHumidity.empty();
}

void Planet::stepWeather(double dtHours) {
    const int n = (int)cells.size();
    if (!sHasWeather || (int)sHumidity.size() != n || (int)sCloud.size() != n || (int)sRain.size() != n)
        initWeather();
    if (dtHours <= 0.0) return;                 // paused: hold the current sky
    if ((int)sTemp.size() != n) return;         // need the climate fields

    const double sea = cfg.seaLevel;
    auto isOcean = [&](int i) { return cells[i].elevation <= sea; };
    const double cw = std::max(1.0, cellWidthMeters());
    double advFrac = std::clamp(cfg.weatherWindKmh * 1000.0 * dtHours / cw, 0.0, 1.0);
    const bool haveSun = ((int)sInsolation.size() == n);
    const bool haveUp  = ((int)sUpwind.size() == n);

    // 1. Advect humidity downwind (upwind differencing) + evaporate over warm sunlit ocean.
    std::vector<double> nh(n);
    for (int i = 0; i < n; ++i) {
        double hUp = (haveUp && sUpwind[i] >= 0) ? sHumidity[sUpwind[i]] : sHumidity[i];
        double h = sHumidity[i] * (1.0 - advFrac) + hUp * advFrac;
        if (isOcean(i)) {
            double tf  = std::clamp((sTemp[i] + 2.0) / 30.0, 0.0, 1.0);          // warm seas evaporate more
            double sun = haveSun ? (0.5 + 0.5 * sInsolation[i]) : 0.7;           // daytime boost
            double target = 0.55 + 0.45 * tf;                                    // marine humidity target
            double rate = 1.0 - std::exp(-cfg.weatherEvapRate * sun * dtHours);
            if (target > h) h += (target - h) * rate;
        }
        nh[i] = h;
    }
    sHumidity.swap(nh);

    // 2. Advect cloud (it drifts with the wind too).
    std::vector<double> nc(n);
    for (int i = 0; i < n; ++i) {
        double cUp = (haveUp && sUpwind[i] >= 0) ? sCloud[sUpwind[i]] : sCloud[i];
        nc[i] = sCloud[i] * (1.0 - advFrac) + cUp * advFrac;
    }
    sCloud.swap(nc);

    // 3. Condense (saturation + orographic lift) -> rain -> dissipate, per cell.
    const double condR = 1.0 - std::exp(-cfg.weatherCondense * dtHours);
    const double rainR = 1.0 - std::exp(-cfg.weatherRainRate * dtHours);
    const double dissR = 1.0 - std::exp(-cfg.weatherCloudDissip * dtHours);
    const double invDt = 1.0 / dtHours;
    for (int i = 0; i < n; ++i) {
        double sat = std::max(0.05, cfg.weatherSatBase + cfg.weatherSatTempCoef * std::max(0.0, sTemp[i]));
        double cond = 0.0;
        double excess = sHumidity[i] - sat;
        if (excess > 0.0) cond += excess * condR;                               // convective/thermal
        if (haveUp && sUpwind[i] >= 0) {                                        // orographic (windward)
            double up = cells[i].elevation - cells[sUpwind[i]].elevation;
            if (up > 0.0) cond += sHumidity[i] * std::min(1.0, up * cfg.weatherOrographic) * condR;
        }
        cond = std::min(cond, sHumidity[i]);
        sHumidity[i] -= cond;
        sCloud[i]    += cond;

        double rain = 0.0;
        if (sCloud[i] > cfg.weatherRainThresh) {
            rain = (sCloud[i] - cfg.weatherRainThresh) * rainR;
            sCloud[i] -= rain;
        }
        double diss = sCloud[i] * dissR;
        sCloud[i]    -= diss;
        sHumidity[i] += diss * 0.5;                                             // half re-evaporates

        sRain[i] = rain * invDt;                                                // intensity (per hour)
        if (sHumidity[i] < 0.0) sHumidity[i] = 0.0;
        sCloud[i] = std::clamp(sCloud[i], 0.0, 1.5);
    }

    // 4. Moving weather systems (lows / hurricanes / typhoons). Drifting agents that travel with
    //    the steering wind and stamp cloud/rain onto the grid, so the sky visibly evolves.
    const bool haveWind = ((int)sWind.size() == n);
    auto wrnd = [&]() { uint32_t x = sWeatherRng; x ^= x << 13; x ^= x >> 17; x ^= x << 5; sWeatherRng = x; return x; };
    auto wrf  = [&]() { return (wrnd() & 0xFFFFFFu) / double(0x1000000); };
    const Vec3 worldUp{0, 1, 0};
    const double D2R = M_PI / 180.0;

    // 4a. Genesis: over warm tropical ocean (5..25 deg) or a mid-latitude (30..62 deg) ocean low.
    if ((int)sStorms.size() < cfg.weatherSystemMax) {
        double pSpawn = 1.0 - std::exp(-cfg.weatherSpawnRate * dtHours);
        if (wrf() < pSpawn) {
            int bestIdx = -1; double bestScore = 0.0; bool bestTrop = false;
            for (int t = 0; t < 8; ++t) {
                int ci = (int)(wrnd() % (uint32_t)n);
                if (cells[ci].elevation > sea) continue;
                double absdeg = std::fabs(std::asin(std::clamp(cells[ci].unit.y, -1.0, 1.0))) / D2R;
                double score = 0.0; bool trop = false;
                if (absdeg > 5.0 && absdeg < 25.0 && sTemp[ci] >= cfg.weatherTropicalSST) { score = 0.6 + 0.4 * wrf(); trop = true; }
                else if (absdeg >= 30.0 && absdeg <= 62.0) { score = 0.3 + 0.3 * wrf(); }
                if (score > bestScore) { bestScore = score; bestIdx = ci; bestTrop = trop; }
            }
            if (bestIdx >= 0) {
                WeatherSystem ws;
                ws.pos = cells[bestIdx].unit;
                ws.strength = 0.15;
                ws.radius = cfg.weatherSystemRadius * (bestTrop ? 0.8 : 1.25);
                ws.life = bestTrop ? (120.0 + 180.0 * wrf()) : (60.0 + 90.0 * wrf());
                ws.spin = (cells[bestIdx].unit.y >= 0.0) ? 1.0 : -1.0;
                ws.tropical = bestTrop;
                ws.id = sStormNextId++;
                sStorms.push_back(ws);
            }
        }
    }

    // 4b. Move, intensify and cull each system.
    for (size_t s = 0; s < sStorms.size(); ) {
        WeatherSystem& ws = sStorms[s];
        int nc = 0; double nd = -2.0;                          // nearest cell to the system
        for (int i = 0; i < n; ++i) { double d = cells[i].unit.dot(ws.pos); if (d > nd) { nd = d; nc = i; } }
        const Vec3& nrm = ws.pos;
        Vec3 steer = (haveWind && sWind[nc].length() > 1e-9) ? sWind[nc].normalized() : Vec3{0, 0, 0};
        Vec3 northT = worldUp - nrm * worldUp.dot(nrm); double nl = northT.length();
        if (nl > 1e-9) northT = northT * (1.0 / nl);
        double poleSign = (nrm.y >= 0.0) ? 1.0 : -1.0;
        Vec3 vel = steer + northT * (poleSign * 0.35);         // steering + poleward recurve
        vel = vel - nrm * vel.dot(nrm);                        // keep tangent
        double vl = vel.length();
        if (vl > 1e-9) {
            double dAng = cfg.weatherSystemSpeed * 1000.0 * dtHours / std::max(1.0, cfg.radius);
            Vec3 vdir = vel * (1.0 / vl);
            ws.pos = (nrm * std::cos(dAng) + vdir * std::sin(dAng)).normalized();
        }
        bool overWarmSea = (cells[nc].elevation <= sea) && (sTemp[nc] >= cfg.weatherTropicalSST - 4.0);
        if (ws.tropical) {
            if (overWarmSea) ws.strength += (1.0 - ws.strength) * (1.0 - std::exp(-0.05 * dtHours));
            else             ws.strength -= ws.strength * (1.0 - std::exp(-0.15 * dtHours));
        } else {
            double frac = std::min(1.0, ws.age / std::max(1.0, ws.life));
            ws.strength = 0.2 + 0.6 * std::sin(frac * M_PI);   // rise then fade
            if (cells[nc].elevation > sea) ws.strength *= 0.7; // weaker over land
        }
        ws.strength = std::clamp(ws.strength, 0.0, 1.0);
        ws.age += dtHours;
        if (ws.age > ws.life || (ws.tropical && ws.strength < 0.05 && cells[nc].elevation > sea)) {
            sStorms[s] = sStorms.back(); sStorms.pop_back();   // swap-remove dead system
        } else ++s;
    }

    // 4c. Stamp each system's cloud/rain shield onto the grid (Gaussian-ish core falloff).
    if (!sStorms.empty()) {
        std::vector<double> cosR(sStorms.size());
        for (size_t s = 0; s < sStorms.size(); ++s) cosR[s] = std::cos(std::min(M_PI, sStorms[s].radius));
        for (int i = 0; i < n; ++i) {
            double moist = 0.3 + 0.7 * std::clamp(sHumidity[i], 0.0, 1.0);
            for (size_t s = 0; s < sStorms.size(); ++s) {
                double dot = cells[i].unit.dot(sStorms[s].pos);
                if (dot < cosR[s]) continue;                   // outside the system radius
                double d = std::acos(std::clamp(dot, -1.0, 1.0));
                double fall = 1.0 - d / sStorms[s].radius; fall *= fall;
                double st = sStorms[s].strength;
                sCloud[i] += st * fall * cfg.weatherSystemCloud * dtHours * moist;
                sRain[i]  += st * fall * cfg.weatherSystemRain * moist;
            }
            sCloud[i] = std::clamp(sCloud[i], 0.0, 1.5);
        }
    }
}