#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;
    }
    // Upper clamp tracks peakSoftCapEnd (matches step()); the soft peak cap, not a
    // fixed 9000 m wall here, governs how high mountains stand.
    #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, cfg.peakSoftCapEnd);

    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
}