planetsim/test_live.cpp

// Headless test for the Live World stage (insolation + live seasonal temperature).
// No display / raylib needed.
//
//   g++ -std=c++17 -O2 -Isrc/sim test_live.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/PlanetLive.cpp \
//       src/sim/PlanetBiota.cpp src/sim/PlanetFloraGen.cpp src/sim/PlanetFaunaGen.cpp \
//       src/sim/PlanetFungiGen.cpp src/sim/PlanetIO.cpp -o /tmp/tl && /tmp/tl
//
// Verifies: insolation range; the sub-solar hemisphere is lit and the night side dark;
// the solar declination tracks axialTilt (poles lit/dark at solstice, neutral at equinox);
// live temperature stays within the summer/winter band and is anti-phased across the
// hemispheres; the snow line advances in the winter hemisphere; and determinism.

#include "Planet.hpp"
#include "Projection.hpp"
#include <cstdio>
#include <cmath>
#include <algorithm>

static int failures = 0;
static void check(bool cond, const char* what) {
    std::printf("  [%s] %s\n", cond ? "PASS" : "FAIL", what);
    if (!cond) ++failures;
}

// Index of the cell whose unit direction is closest to dir.
static int nearestCell(const Planet& p, const Vec3& dir) {
    int best = 0; double bd = -2.0;
    for (int i = 0; i < (int)p.cells.size(); ++i) {
        double d = p.cells[i].unit.dot(dir);
        if (d > bd) { bd = d; best = i; }
    }
    return best;
}

int main() {
    Planet planet;
    PlanetConfig cfg;
    cfg.subdivisions = 5;
    cfg.seed = 1337;
    planet.generate(cfg);
    planet.computeClimate();                 // live season needs the climate fields
    const int n = (int)planet.cells.size();

    // North/south pole cells (extreme latitude) for declination checks.
    int north = 0, south = 0;
    for (int i = 0; i < n; ++i) {
        if (planet.cells[i].unit.y > planet.cells[north].unit.y) north = i;
        if (planet.cells[i].unit.y < planet.cells[south].unit.y) south = i;
    }
    const double tilt = cfg.axialTilt * M_PI / 180.0;

    std::printf("Live World: insolation\n");
    // --- Range + lit/dark hemispheres (arbitrary time, summer solstice) -------
    planet.computeInsolation(0.25, 0.30);
    const std::vector<double>& sun = planet.insolation();
    bool inRange = true; int lit = 0; double mx = 0.0;
    for (int i = 0; i < n; ++i) {
        if (sun[i] < 0.0 || sun[i] > 1.0) inRange = false;
        if (sun[i] > 0.0) ++lit;
        mx = std::max(mx, sun[i]);
    }
    check(inRange, "insolation in [0,1]");
    check(lit > (int)(0.40 * n) && lit < (int)(0.60 * n), "about half the planet is in daylight");
    // Sub-solar cell is fully lit; its antipode is dark.
    double decl = tilt * std::sin(2.0 * M_PI * 0.25);
    double lon  = M_PI * (1.0 - 2.0 * 0.30);
    Vec3 sd = lonLatToDir(lon, decl);
    int subsolar = nearestCell(planet, sd);
    int antipode = nearestCell(planet, sd * -1.0);
    check(sun[subsolar] > 0.99, "sub-solar cell is fully lit");
    check(sun[antipode] == 0.0, "antipodal (midnight) cell is dark");

    std::printf("Live World: declination tracks axial tilt\n");
    // Equinox (decl=0): both poles near the terminator (~0). Solstice: summer pole lit,
    // winter pole in polar night.
    planet.computeInsolation(0.0, 0.0);                     // equinox
    double npEq = planet.insolation()[north], spEq = planet.insolation()[south];
    planet.computeInsolation(0.25, 0.0);                    // northern summer solstice
    double npSol = planet.insolation()[north], spSol = planet.insolation()[south];
    check(npEq < 0.1 && spEq < 0.1, "equinox: both poles near the terminator");
    check(npSol > 0.3, "summer solstice: summer pole sees the midnight sun");
    check(spSol == 0.0, "summer solstice: winter pole is in polar night");
    check(std::fabs(npSol - std::sin(tilt)) < 0.05, "polar insolation ~ sin(axialTilt)");

    std::printf("Live World: seasonal temperature\n");
    // At northern summer (doy=0.25) live temp = summer in the north, winter in the south;
    // everywhere it stays within the [winter, summer] band.
    planet.computeLiveSeason(0.25);
    const std::vector<double>& lt   = planet.liveTemp();
    const std::vector<double>& summ = planet.summerTemp();
    const std::vector<double>& wint = planet.winterTemp();
    bool inBand = true; double nLiveSummer = 0.0; int nN = 0;
    for (int i = 0; i < n; ++i) {
        if (lt[i] < wint[i] - 1e-6 || lt[i] > summ[i] + 1e-6) inBand = false;
        if (planet.cells[i].unit.y > 0.3) { nLiveSummer += lt[i]; ++nN; }
    }
    nLiveSummer /= std::max(1, nN);
    check(inBand, "live temp within [winter, summer] for every cell");
    check(std::fabs(lt[north] - summ[north]) < 0.5, "northern summer solstice -> north at its summer temp");
    check(std::fabs(lt[south] - wint[south]) < 0.5, "northern summer solstice -> south at its winter temp");
    // Half a year later the northern hemisphere is colder (anti-phase).
    planet.computeLiveSeason(0.75);
    double nLiveWinter = 0.0; nN = 0;
    for (int i = 0; i < n; ++i)
        if (planet.cells[i].unit.y > 0.3) { nLiveWinter += planet.liveTemp()[i]; ++nN; }
    nLiveWinter /= std::max(1, nN);
    check(nLiveWinter < nLiveSummer - 1.0, "northern hemisphere colder in its winter than its summer");

    std::printf("Live World: snow line advances in winter\n");
    auto snowCountNorth = [&](double doy) {
        planet.computeLiveSeason(doy); const std::vector<double>& t = planet.liveTemp();
        int c = 0;
        for (int i = 0; i < n; ++i)
            if (planet.cells[i].unit.y > 0.0 && planet.cells[i].elevation > cfg.seaLevel
                && t[i] < cfg.snowTemp) ++c;
        return c;
    };
    int snowSummer = snowCountNorth(0.25), snowWinter = snowCountNorth(0.75);
    std::printf("    north land snow cells: summer %d, winter %d\n", snowSummer, snowWinter);
    check(snowWinter > snowSummer, "more northern land under snow in winter than summer");

    std::printf("Live World: determinism\n");
    Planet p2; p2.generate(cfg); p2.computeClimate();
    p2.computeInsolation(0.37, 0.61); p2.computeLiveSeason(0.37);
    planet.computeInsolation(0.37, 0.61); planet.computeLiveSeason(0.37);
    bool same = true;
    for (int i = 0; i < n; ++i)
        if (p2.insolation()[i] != planet.insolation()[i] || p2.liveTemp()[i] != planet.liveTemp()[i])
            same = false;
    check(same, "same seed + time -> identical insolation & live temp");

    std::printf(failures ? "\nSOME LIVE CHECKS FAILED (%d)\n" : "\nALL LIVE CHECKS PASSED\n", failures);
    return failures ? 1 : 0;
}