// Painted-map import — the uplift-painting workflow from the Salty terrain generator
// (Tools/Terrain, `terrain plan` / `terrain bake`) on World Orogen's sphere mesh.
//
// The author paints a flat equirectangular map where every colour is a *class*: a rate of rock
// uplift and an erodibility, never a height. A legend JSON beside the painting says what the
// colours mean. This module turns the two into a planet:
//
//   1. classify   — every pixel goes to its nearest legend colour (nearest, never "unmatched",
//                   so a JPEG halo or a stray pixel lands on something sensible; the report
//                   says how many were far from everything).
//   2. vote       — every Voronoi region takes the majority class of the pixels under it.
//   3. strokes    — the white outline an artist draws round every island dissolves into
//                   whichever real class is nearest; a white blob touching a pole is the ice cap.
//   4. coast      — the drawn shoreline is roughened by adding noise to the signed distance
//                   from it, because a drawn coast is a smooth curve and a real one is fractal.
//   5. uplift     — class rate × the planet's upland fabric (a massif is where one fabric,
//                   cut at a quantile of the *planet*, stands high) × a coastal-plain ramp ×
//                   a regional swell; erodibility = class k × the planet's rock field.
//   6. solve      — dh/dt = U − K·A^m·S, integrated implicitly up the drainage stack
//                   (Braun & Willett 2013) with the ocean fixed at sea level, until the land is
//                   in balance with its uplift. Rivers, divides and the valley hierarchy come out
//                   of the physics; the painting decides only where the land rises and how fast.
//
// Paint the uplift, never the height: a solve handed a painted surface erodes it into something
// else within a few hundred steps and throws the drainage network away.
//
// Everything here is pure computation with no DOM, so it runs inside planet-worker.js. The
// legend parser and the pixel classifier also run on the main thread, to show the match report
// before anything is solved.

import { SimplexNoise } from './simplex-noise.js';
import { elevToHeightKm } from './color-map.js';

export const DEFAULT_WARN_DISTANCE = 60;
export const DEFAULT_COASTAL_FLOOR_MM_YR = 0.02;
export const MAX_MASSIF_FRACTION = 0.6;
export const EARTH_CIRCUMFERENCE_KM = 40030;
export const EARTH_RADIUS_KM = 6371;

// Defaults for the planet block, matching RawContent/World/Planet.json in the Salty repo.
export const PLANET_DEFAULTS = {
    circumferenceKm: 100,
    massifWavelengthKm: 7,
    lithologyWavelengthKm: 8,
    lithology: [0.6, 1.0, 1.8],
    variation: 0.30,
    coastDetail: 0.35,
    steps: 200,
    peakKm: 4.5,
    oceanDepthKm: 4.0,
    seed: 7945,
    // The bake's geology grid - Planet.json's pipeline block - which the class table's angles are about
    // (painted-report.js). Not used by the solve here, whose relief is a scale.
    k: 5e-5, m: 0.5, cellM: 8, talusDeg: 35,
};

// ─── Legend ───────────────────────────────────────────────────────

/**
 * Parse a legend JSON object (the same schema Tools/Terrain reads) into a normalised legend.
 * Throws with a readable message when the legend cannot make a world.
 */
export function parseLegend(obj) {
    if (!obj || typeof obj !== 'object') throw new Error('The legend is not a JSON object');
    if (!Array.isArray(obj.classes) || obj.classes.length === 0) throw new Error('The legend has no "classes" array');
    if (obj.classes.length > 255) throw new Error('The legend has more than 255 classes');

    const classes = obj.classes.map((c, i) => {
        if (!c || typeof c !== 'object') throw new Error(`Class ${i} is not an object`);
        if (typeof c.name !== 'string' || !c.name) throw new Error(`Class ${i} has no name`);
        let rgb = null;
        if (Array.isArray(c.rgb) && c.rgb.length === 3) {
            rgb = c.rgb.map(v => Math.max(0, Math.min(255, Math.round(+v || 0))));
        }
        if (!rgb && !c.derived) throw new Error(`Class "${c.name}" has no rgb`);
        const sea = !!c.sea;
        let massif = null;
        if (c.massif && +c.massif.fraction > 0 && !sea) {
            massif = {
                floorMmYr: Math.max(0, +c.massif.floor_mm_yr || 0),
                fraction: Math.min(MAX_MASSIF_FRACTION, +c.massif.fraction),
            };
        }
        return {
            index: i,
            name: c.name,
            rgb: rgb || [0, 0, 0],
            sea,
            depthM: sea ? Math.max(0, +c.depth_m || 0) : 0,
            upliftMmYr: sea ? 0 : Math.max(0, +c.uplift_mm_yr || 0),
            kMult: (+c.k_mult > 0) ? +c.k_mult : 1,
            stroke: !!c.stroke,
            derived: !!c.derived,
            snow: !!c.snow,
            edgeClass: typeof c.edge_class === 'string' ? c.edge_class : '',
            edgeIndex: -1,
            massif,
            coastalPlainKm: Math.max(0, +c.coastal_plain_km || 0),
            coastalFloorMmYr: Math.max(0, +c.coastal_floor_mm_yr || 0),
            lithologyMix: (c.lithology_mix === undefined || c.lithology_mix === null) ? 1 : Math.max(0, Math.min(1, +c.lithology_mix)),
            raw: c,
        };
    });

    for (const c of classes) {
        if (!c.edgeClass) continue;
        const j = classes.findIndex(o => o.name === c.edgeClass);
        if (j < 0) throw new Error(`Class "${c.name}" names edge_class "${c.edgeClass}", which is not in the legend`);
        c.edgeIndex = j;
    }

    const paintable = classes.filter(c => !c.derived);
    if (!paintable.some(c => c.sea)) throw new Error('The legend has no sea class');
    if (!classes.some(c => !c.sea && !c.stroke)) throw new Error('The legend has no land class');

    const planet = readPlanetBlock(obj);

    return {
        classes,
        warnDistance: +obj.warn_distance > 0 ? +obj.warn_distance : DEFAULT_WARN_DISTANCE,
        planet,
        source: obj,
    };
}

/**
 * The "planet" block: the numbers Planet.json carries in the Go tool. Read from a legend that carries an
 * optional copy of it, or from Planet.json itself, which has the same keys plus the pipeline block the
 * class table's angles need. Missing keys keep `base`, which is PLANET_DEFAULTS for a legend.
 */
function readPlanetBlock(obj, base = PLANET_DEFAULTS) {
    const p = obj.planet || {};
    const pipe = obj.pipeline || {};
    const lith = pipe.lithology || p.lithology || null;
    const kmults = lith && Array.isArray(lith.k_multipliers) && lith.k_multipliers.length >= 2
        ? lith.k_multipliers.map(v => Math.max(0.05, +v || 1)) : (base.lithology || PLANET_DEFAULTS.lithology).slice();
    const num = (v, d) => (v !== undefined && v !== null && isFinite(+v)) ? +v : d;
    const fluvial = pipe.fluvial || {};
    const thermal = pipe.thermal || {};
    // The geology cell is the detail quad times the geology factor, as manifest.go derives it.
    const quadM = num(obj.quad_cm, 0) / 100;
    const geologyFactor = num(pipe.geology_factor, 0);
    return {
        circumferenceKm: Math.max(1, num(p.circumference_km, base.circumferenceKm)),
        massifWavelengthKm: Math.max(0, num(p.massif_wavelength_km, base.massifWavelengthKm)),
        lithologyWavelengthKm: Math.max(0, num(p.lithology_wavelength_km, base.lithologyWavelengthKm)),
        lithology: kmults,
        variation: Math.max(0, Math.min(0.6, num(p.uplift_variation, base.variation))),
        seed: num(obj.source && obj.source.seed, num(p.seed, base.seed)),
        k: Math.max(0, num(fluvial.k, base.k)),
        m: num(fluvial.m, base.m),
        cellM: quadM > 0 && geologyFactor > 0 ? quadM * geologyFactor : num(p.cell_m, base.cellM),
        talusDeg: num(thermal.talus_deg, base.talusDeg),
        overlayLegend: typeof p.overlay_legend === 'string' ? p.overlay_legend : (base.overlayLegend || ''),
        overlay: typeof p.overlay === 'string' ? p.overlay : (base.overlay || ''),
    };
}

/**
 * Read Planet.json - the Go tool's own manifest - onto a parsed legend, so its planet block, its lithology
 * multipliers, its seed and the geology grid's constants travel without being retyped. Keys it lacks keep
 * what the legend had. Returns the merged block.
 */
export function applyPlanetManifest(legend, manifest) {
    if (!manifest || typeof manifest !== 'object') throw new Error('Planet.json is not a JSON object');
    if (!manifest.planet || typeof manifest.planet !== 'object') throw new Error('Planet.json has no "planet" block; the painted path reads a painted planet');
    legend.planet = readPlanetBlock(manifest, legend.planet);
    return legend.planet;
}

/** The class's plain-at-the-waterline rate, never above the class rate. */
export function plainFloorMmYr(c) {
    let floor = c.coastalFloorMmYr > 0 ? c.coastalFloorMmYr : DEFAULT_COASTAL_FLOOR_MM_YR;
    if (floor > c.upliftMmYr) floor = c.upliftMmYr;
    return floor;
}

/**
 * Write a legend back out as JSON text, keeping every key of the loaded file and only replacing
 * the numbers the table edits, so the file still reads in Tools/Terrain with its commentary intact.
 */
export function serializeLegend(legend) {
    const out = JSON.parse(JSON.stringify(legend.source));
    for (let i = 0; i < legend.classes.length; i++) {
        const c = legend.classes[i];
        const raw = out.classes[i];
        if (!raw) continue;
        if (c.sea) raw.depth_m = c.depthM;
        else {
            raw.uplift_mm_yr = c.upliftMmYr;
            raw.k_mult = c.kMult;
            if (c.massif && raw.massif) {
                raw.massif.floor_mm_yr = c.massif.floorMmYr;
                raw.massif.fraction = c.massif.fraction;
            }
        }
    }
    return JSON.stringify(out, null, 2);
}

// ─── Pixel classification ─────────────────────────────────────────

/**
 * Assign every pixel of an RGBA image to its nearest paintable class.
 * Returns the class raster plus a match report.
 */
export function classifyImage(rgba, w, h, legend) {
    const cls = legend.classes;
    const n = cls.length;
    const paint = [];
    for (let i = 0; i < n; i++) if (!cls[i].derived) paint.push(i);
    const np = paint.length;
    const pr = new Int32Array(np), pg = new Int32Array(np), pb = new Int32Array(np);
    for (let k = 0; k < np; k++) {
        const c = cls[paint[k]];
        pr[k] = c.rgb[0]; pg[k] = c.rgb[1]; pb[k] = c.rgb[2];
    }

    // Colours are quantised to six bits a channel and classified once per cell. A legend's
    // classes sit tens of units apart (warn distance 60), so the ±2 of the cell is nothing.
    const cacheCls = new Uint8Array(1 << 18).fill(255);
    const cacheD2 = new Float32Array(1 << 18);

    const total = w * h;
    const out = new Uint8Array(total);
    const counts = new Int32Array(n);
    const warn2 = legend.warnDistance * legend.warnDistance;
    let far = 0, maxD2 = 0, maxAt = 0;

    for (let p = 0, o = 0; p < total; p++, o += 4) {
        const r = rgba[o], g = rgba[o + 1], b = rgba[o + 2];
        const key = ((r >> 2) << 12) | ((g >> 2) << 6) | (b >> 2);
        let c = cacheCls[key];
        if (c === 255) {
            const cr = (r & ~3) + 2, cg = (g & ~3) + 2, cb = (b & ~3) + 2;
            let best = 0, bestD = Infinity;
            for (let k = 0; k < np; k++) {
                const dr = cr - pr[k], dg = cg - pg[k], db = cb - pb[k];
                const d = dr * dr + dg * dg + db * db;
                if (d < bestD) { bestD = d; best = k; }
            }
            c = paint[best];
            cacheCls[key] = c;
            cacheD2[key] = bestD;
        }
        out[p] = c;
        counts[c]++;
        const d2 = cacheD2[key];
        if (d2 > warn2) far++;
        if (d2 > maxD2) { maxD2 = d2; maxAt = p; }
    }

    // The wrap: the left and right columns are the same meridian.
    let wrapDiffer = 0, wrapLandSea = 0;
    for (let y = 0; y < h; y++) {
        const a = out[y * w], b = out[y * w + w - 1];
        if (a !== b) {
            wrapDiffer++;
            if (cls[a].sea !== cls[b].sea) wrapLandSea++;
        }
    }

    return {
        classes: out,
        counts,
        total,
        far,
        maxDist: Math.sqrt(maxD2),
        maxAt: [maxAt % w, (maxAt / w) | 0],
        wrapRows: h,
        wrapDiffer,
        wrapLandSea,
    };
}

// ─── Region sampling ──────────────────────────────────────────────

const clamp1 = v => Math.max(-1, Math.min(1, v));

/**
 * Vote the class raster onto the mesh: every region takes the majority class of a box of
 * pixels the size of its own footprint, so a thin stroke or a JPEG halo never decides a cell.
 * Returns the class per region and the latitude per region (reused by the stroke pass).
 */
export function sampleClassesToMesh(mesh, r_xyz, classRaster, w, h, numClasses) {
    const N = mesh.numRegions;
    const r_class = new Uint8Array(N);
    const r_lat = new Float32Array(N);
    const spacing = Math.sqrt(4 * Math.PI / Math.max(1, N - 1)); // radians between neighbours
    const pxPerRadX = w / (2 * Math.PI);
    const pxPerRadY = h / Math.PI;
    const counts = new Int32Array(numClasses);
    const MAXS = 7;

    for (let r = 0; r < N; r++) {
        const x = r_xyz[3 * r], y = r_xyz[3 * r + 1], z = r_xyz[3 * r + 2];
        const lat = Math.asin(clamp1(y));
        const lon = Math.atan2(x, z);
        r_lat[r] = lat;
        const cx = (lon / Math.PI + 1) * 0.5 * w;
        const cy = (0.5 - lat / Math.PI) * h;
        const cosLat = Math.max(Math.cos(lat), 1e-3);
        let hx = 0.5 * spacing * pxPerRadX / cosLat;
        if (hx > w / 2) hx = w / 2;
        const hy = 0.5 * spacing * pxPerRadY;
        const nx = Math.min(MAXS, Math.max(1, Math.round(2 * hx)));
        const ny = Math.min(MAXS, Math.max(1, Math.round(2 * hy)));
        counts.fill(0);
        for (let j = 0; j < ny; j++) {
            const sy = ny === 1 ? cy : cy - hy + (2 * hy) * (j + 0.5) / ny;
            let py = Math.floor(sy);
            if (py < 0) py = 0; else if (py >= h) py = h - 1;
            const row = py * w;
            for (let i = 0; i < nx; i++) {
                const sx = nx === 1 ? cx : cx - hx + (2 * hx) * (i + 0.5) / nx;
                let px = Math.floor(sx);
                px = ((px % w) + w) % w;
                counts[classRaster[row + px]]++;
            }
        }
        let best = 0;
        for (let c = 1; c < numClasses; c++) if (counts[c] > counts[best]) best = c;
        r_class[r] = best;
    }
    return { r_class, r_lat };
}

/**
 * Resolve stroke classes on the mesh. A stroke component touching a pole becomes its
 * edge_class (the ice cap painted in the same white as the outlines); everything else
 * dissolves into the nearest non-stroke class by BFS.
 */
export function resolveStrokes(mesh, r_class, r_lat, legend) {
    const cls = legend.classes;
    const N = mesh.numRegions;
    const isStroke = new Uint8Array(cls.length);
    let any = false;
    for (let i = 0; i < cls.length; i++) if (cls[i].stroke) { isStroke[i] = 1; any = true; }
    if (!any) return { edgeAssigned: 0, dissolved: 0 };

    const { adjOffset, adjList } = mesh;
    const spacing = Math.sqrt(4 * Math.PI / Math.max(1, N - 1));
    const poleLat = Math.PI / 2 - 2.5 * spacing;
    const queue = new Int32Array(N);
    const visited = new Uint8Array(N);
    let edgeAssigned = 0, dissolved = 0;

    // Connected components of each stroke class; the ones at a pole become the edge class.
    for (let r = 0; r < N; r++) {
        if (visited[r] || !isStroke[r_class[r]]) continue;
        const ci = r_class[r];
        const edge = cls[ci].edgeIndex;
        let head = 0, tail = 0, touches = false;
        queue[tail++] = r; visited[r] = 1;
        while (head < tail) {
            const c = queue[head++];
            if (Math.abs(r_lat[c]) > poleLat) touches = true;
            for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
                const nb = adjList[j];
                if (!visited[nb] && r_class[nb] === ci) { visited[nb] = 1; queue[tail++] = nb; }
            }
        }
        if (touches && edge >= 0) {
            for (let k = 0; k < tail; k++) r_class[queue[k]] = edge;
            edgeAssigned += tail;
        }
    }

    // Dissolve what is left into the nearest real class.
    const assigned = new Uint8Array(N);
    let head = 0, tail = 0;
    for (let r = 0; r < N; r++) {
        if (!isStroke[r_class[r]]) { assigned[r] = 1; queue[tail++] = r; }
    }
    while (head < tail) {
        const c = queue[head++];
        for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
            const nb = adjList[j];
            if (!assigned[nb]) { assigned[nb] = 1; r_class[nb] = r_class[c]; queue[tail++] = nb; dissolved++; }
        }
    }
    return { edgeAssigned, dissolved };
}

// ─── Coast ────────────────────────────────────────────────────────

/**
 * Hop distance from the coast: for land, to the nearest sea region; for sea, to the nearest
 * land region. 1 means adjacent. 0 means the planet has no coast at all.
 */
export function coastDistance(mesh, r_land) {
    const N = mesh.numRegions;
    const { adjOffset, adjList } = mesh;
    const hop = new Int32Array(N);
    const queue = new Int32Array(N);
    let head = 0, tail = 0;
    for (let r = 0; r < N; r++) {
        const land = r_land[r];
        for (let j = adjOffset[r], jEnd = adjOffset[r + 1]; j < jEnd; j++) {
            if (r_land[adjList[j]] !== land) { hop[r] = 1; queue[tail++] = r; break; }
        }
    }
    while (head < tail) {
        const c = queue[head++];
        const land = r_land[c];
        const d = hop[c] + 1;
        for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
            const nb = adjList[j];
            if (hop[nb] === 0 && r_land[nb] === land) { hop[nb] = d; queue[tail++] = nb; }
        }
    }
    return hop;
}

/** Connected components of land regions; returns per-region component id and each component's max coast hop. */
function landComponents(mesh, r_land, hop) {
    const N = mesh.numRegions;
    const { adjOffset, adjList } = mesh;
    const comp = new Int32Array(N).fill(-1);
    const compMax = [];
    const queue = new Int32Array(N);
    for (let r = 0; r < N; r++) {
        if (!r_land[r] || comp[r] >= 0) continue;
        const id = compMax.length;
        let head = 0, tail = 0, mx = 0;
        queue[tail++] = r; comp[r] = id;
        while (head < tail) {
            const c = queue[head++];
            if (hop[c] > mx) mx = hop[c];
            for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
                const nb = adjList[j];
                if (r_land[nb] && comp[nb] < 0) { comp[nb] = id; queue[tail++] = nb; }
            }
        }
        compMax.push(mx);
    }
    return { comp, compMax };
}

/**
 * Roughen the drawn coastline by adding fractal noise to the signed hop distance from it.
 * `amount` 0..1 is up to three cells of shift; an islet may lose at most two thirds of its
 * width so an archipelago does not vanish. Returns a new land mask.
 */
export function roughenCoast(mesh, r_xyz, r_land, amount, seed, r_jitter = null) {
    const N = mesh.numRegions;
    const out = new Uint8Array(r_land);
    if (amount <= 0) return { r_land: out, flipped: 0 };
    const hop = coastDistance(mesh, r_land);
    const amp = amount * 3;
    // The overlay's coast_jitter marks scale the reach per region (painted-overlay.js): 0 pins the shore as
    // painted, above 1 chews it harder. The early-out below has to use the largest reach any of them asks.
    let ampMax = amp;
    if (r_jitter) {
        let jm = 1;
        for (let r = 0; r < N; r++) if (r_jitter[r] > jm) jm = r_jitter[r];
        ampMax = amp * jm;
    }
    const { comp, compMax } = landComponents(mesh, r_land, hop);
    const noise = new SimplexNoise(seed * 31 + 17);
    // Bays about eight cells wide, with four octaves below that.
    const spacing = Math.sqrt(4 * Math.PI / Math.max(1, N - 1));
    const F = 1 / (8 * spacing);
    let flipped = 0;
    for (let r = 0; r < N; r++) {
        if (hop[r] === 0) continue;
        const d = hop[r] - 0.5;
        if (d > ampMax + 1) continue;
        let a = r_jitter ? amp * r_jitter[r] : amp;
        if (a <= 0 || d > a + 1) continue;
        if (r_land[r]) {
            const cap = 0.66 * compMax[comp[r]];
            if (cap < a) a = cap;
            if (a < 0.5) continue;
        }
        const n = noise.fbm(r_xyz[3 * r] * F, r_xyz[3 * r + 1] * F, r_xyz[3 * r + 2] * F, 4, 0.5);
        const s = (r_land[r] ? d : -d) + 2 * n * a;
        const land = s > 0 ? 1 : 0;
        if (land !== r_land[r]) { out[r] = land; flipped++; }
    }
    return { r_land: out, flipped };
}

/**
 * After the coast moves, a region can be land wearing a sea class or the other way round.
 * Give each such region the class of the nearest region that agrees with its new type.
 */
export function reconcileClasses(mesh, r_class, r_land, legend) {
    const cls = legend.classes;
    const N = mesh.numRegions;
    const { adjOffset, adjList } = mesh;
    const queue = new Int32Array(N);
    const done = new Uint8Array(N);
    let changed = 0;
    for (const wantLand of [1, 0]) {
        done.fill(0);
        let head = 0, tail = 0;
        for (let r = 0; r < N; r++) {
            const consistent = (cls[r_class[r]].sea ? 0 : 1) === r_land[r];
            if (consistent) { done[r] = 1; if (r_land[r] === wantLand) queue[tail++] = r; }
            else if (r_land[r] !== wantLand) done[r] = 1; // the other pass's problem
        }
        while (head < tail) {
            const c = queue[head++];
            for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
                const nb = adjList[j];
                if (!done[nb]) { done[nb] = 1; r_class[nb] = r_class[c]; queue[tail++] = nb; changed++; }
            }
        }
        // A region nothing reached (an all-sea planet turned to land somewhere) takes the first
        // class of the wanted type.
        for (let r = 0; r < N; r++) {
            if (done[r]) continue;
            const idx = cls.findIndex(c => !c.derived && !c.stroke && (c.sea ? 0 : 1) === wantLand);
            if (idx >= 0) { r_class[r] = idx; changed++; }
            done[r] = 1;
        }
    }
    return changed;
}

// ─── Uplift field ─────────────────────────────────────────────────

/** Exact rank of every value in 0..1 over the whole array: the share of the planet standing below it. */
export function rankField(values) {
    const N = values.length;
    const idx = new Uint32Array(N);
    for (let i = 0; i < N; i++) idx[i] = i;
    idx.sort((a, b) => values[a] - values[b]);
    const rank = new Float32Array(N);
    const denom = Math.max(1, N - 1);
    for (let i = 0; i < N; i++) rank[idx[i]] = i / denom;
    return rank;
}

function massifShape(rank, fraction) {
    const lo = 1 - 1.5 * fraction;
    const hi = 1 - 0.5 * fraction;
    const t = (rank - lo) / (hi - lo);
    if (t <= 0) return 0;
    if (t >= 1) return 1;
    return t * t * (3 - 2 * t);
}

/** The class rate where the fabric is high, the floor where it is low, cut in rank so `fraction` means what it says. */
export function massifRate(floorMmYr, rateMmYr, rank, fraction) {
    return floorMmYr + (rateMmYr - floorMmYr) * massifShape(rank, fraction);
}

function sampleFbm(noise, r_xyz, N, F, octaves, gain) {
    const out = new Float32Array(N);
    for (let r = 0; r < N; r++) {
        out[r] = noise.fbm(r_xyz[3 * r] * F, r_xyz[3 * r + 1] * F, r_xyz[3 * r + 2] * F, octaves, gain);
    }
    return out;
}

/**
 * Build the uplift-rate field (mm/yr per region) and the erodibility multiplier from the
 * classes and the legend. Distances in the legend are in the painted planet's kilometres and
 * are scaled to Orogen's Earth-sized globe by 40030 / circumferenceKm.
 */
export function buildUpliftField(mesh, r_xyz, r_class, r_land, legend, opts) {
    const cls = legend.classes;
    const N = mesh.numRegions;
    const circ = Math.max(1, opts.circumferenceKm);
    const scale = EARTH_CIRCUMFERENCE_KM / circ;
    const avgEdgeKm = Math.PI * EARTH_RADIUS_KM / Math.sqrt(N);
    const seed = opts.seed | 0;
    const hop = coastDistance(mesh, r_land);

    // A noise wavelength of λ painted kilometres is F = circ / (2π λ) noise units per radian.
    const freqFor = km => circ / (2 * Math.PI * Math.max(1e-6, km));

    let rank = null;
    if (opts.massifWavelengthKm > 0 && cls.some(c => c.massif)) {
        const fab = sampleFbm(new SimplexNoise(seed * 7 + 1), r_xyz, N, freqFor(opts.massifWavelengthKm), 5, 0.45);
        rank = rankField(fab);
    }

    let rock = null;
    const kmults = opts.lithology || [];
    if (opts.lithologyWavelengthKm > 0 && kmults.length >= 2 && cls.some(c => c.lithologyMix > 0)) {
        const fab = sampleFbm(new SimplexNoise(seed * 7 + 2), r_xyz, N, freqFor(opts.lithologyWavelengthKm), 4, 0.5);
        const rr = rankField(fab);
        rock = new Float32Array(N);
        const types = kmults.length;
        for (let r = 0; r < N; r++) rock[r] = kmults[Math.min(types - 1, Math.floor(rr[r] * types))];
    }

    let swell = null;
    const variation = Math.max(0, opts.variation || 0);
    if (variation > 0) swell = sampleFbm(new SimplexNoise(seed * 7 + 3), r_xyz, N, freqFor(25), 4, 0.5);

    const r_rate = new Float32Array(N);
    const r_k = new Float32Array(N);
    let rateMax = 0, landCount = 0;
    for (let r = 0; r < N; r++) {
        const c = cls[r_class[r]];
        let k = c.kMult;
        if (rock && c.lithologyMix > 0) k *= 1 + c.lithologyMix * (rock[r] - 1);
        r_k[r] = k;
        if (!r_land[r]) continue;
        landCount++;
        let rate = c.upliftMmYr;
        if (rank && c.massif) rate = massifRate(c.massif.floorMmYr, rate, rank[r], c.massif.fraction);
        if (c.coastalPlainKm > 0) {
            const plainKm = c.coastalPlainKm * scale;
            const shoreKm = Math.max(0, hop[r] - 0.5) * avgEdgeKm;
            let t = Math.min(1, shoreKm / plainKm);
            t = t * t * (3 - 2 * t);
            const floor = plainFloorMmYr(c);
            if (rate > floor) rate = floor + (rate - floor) * t;
        }
        if (swell) {
            let sw = 0.5 + swell[r] * 0.8;
            if (sw < 0) sw = 0; else if (sw > 1) sw = 1;
            rate *= 1 + variation * (2 * sw - 1);
        }
        r_rate[r] = rate;
        if (rate > rateMax) rateMax = rate;
    }
    return { r_rate, r_k, rateMax, avgEdgeKm, scale, hop, landCount, massifRank: rank };
}

// ─── The solve ────────────────────────────────────────────────────

function hash01(a, b) {
    let x = (Math.imul(a, 374761393) + Math.imul(b, 668265263)) | 0;
    x = Math.imul(x ^ (x >>> 13), 1274126177);
    x ^= x >>> 16;
    return (x >>> 0) / 4294967296;
}

/** Binary min-heap over region indices with a copied key, sized once. */
class RegionHeap {
    constructor(capacity) {
        this.keys = new Float64Array(capacity);
        this.items = new Int32Array(capacity);
        this.size = 0;
    }
    clear() { this.size = 0; }
    push(item, key) {
        let i = this.size++;
        const keys = this.keys, items = this.items;
        while (i > 0) {
            const p = (i - 1) >> 1;
            if (keys[p] <= key) break;
            keys[i] = keys[p]; items[i] = items[p];
            i = p;
        }
        keys[i] = key; items[i] = item;
    }
    pop() {
        const keys = this.keys, items = this.items;
        const top = items[0];
        const n = --this.size;
        if (n > 0) {
            const key = keys[n], item = items[n];
            let i = 0;
            while (true) {
                let l = 2 * i + 1;
                if (l >= n) break;
                const r = l + 1;
                if (r < n && keys[r] < keys[l]) l = r;
                if (keys[l] >= key) break;
                keys[i] = keys[l]; items[i] = items[l];
                i = l;
            }
            keys[i] = key; items[i] = item;
        }
        return top;
    }
}

/**
 * Solve the stream-power equation on the mesh from an uplift field.
 *
 * Units are dimensionless: U is the rate as a fraction of the largest class rate, K is the
 * erodibility multiplier, A is drainage area in cells and lengths are in mean edges, so a
 * divide one cell from the sea at full rate stands about one unit high. For n = 1 the steady
 * state is linear in U/K, so the relief is set afterwards by a single scale (see toElevation)
 * and the shape of the land — where the rivers run, how the valleys nest, how far a coast is
 * from its divide — is what the solve decides.
 *
 * Each step: priority-flood so every land cell has a downhill path to the ocean, steepest
 * receivers, a donor stack, drainage area down the stack, the implicit update up it, and a
 * touch of hillslope diffusion so the divides are rounded rather than needles.
 */
export function solveUplift(mesh, neighborDist, r_land, r_rate, r_k, rateMax, params, onProgress) {
    const N = mesh.numRegions;
    const { adjOffset, adjList } = mesh;
    const steps = Math.max(1, params.steps | 0);
    const dt = 1;
    const m = params.m ?? 0.5;
    const alpha = params.diffusion ?? 0.04;
    const seed = params.seed | 0;
    const eps = 1e-4;

    let sumL = 0;
    for (let i = 0; i < adjList.length; i++) sumL += neighborDist[i];
    const meanEdge = sumL / Math.max(1, adjList.length);
    const invMean = 1 / meanEdge;

    const U = new Float32Array(N);
    const h = new Float64Array(N);
    for (let r = 0; r < N; r++) {
        if (!r_land[r]) continue;
        U[r] = rateMax > 0 ? r_rate[r] / rateMax : 0;
        // A little initial relief, scaled by the rate, only so the first routing has something
        // to bite on. The solve produces the relief; starting from ridges means tearing them down.
        h[r] = 0.05 * U[r] * (0.75 + 0.5 * hash01(r, seed)) + 1e-3 * hash01(r, seed + 1);
    }

    const receiver = new Int32Array(N);
    const recvLen = new Float32Array(N);
    const donorOff = new Int32Array(N + 1);
    const donorList = new Int32Array(N);
    const cursor = new Int32Array(N);
    const stack = new Int32Array(N);
    const area = new Float32Array(N);
    const closed = new Uint8Array(N);
    const tmp = new Float64Array(N);
    const heap = new RegionHeap(N);

    function flood(step) {
        closed.fill(0);
        heap.clear();
        for (let r = 0; r < N; r++) {
            if (!r_land[r]) { closed[r] = 1; heap.push(r, h[r]); }
        }
        while (heap.size > 0) {
            const c = heap.pop();
            const hc = h[c];
            for (let j = adjOffset[c], jEnd = adjOffset[c + 1]; j < jEnd; j++) {
                const nb = adjList[j];
                if (closed[nb]) continue;
                closed[nb] = 1;
                // The epsilon is scattered per cell and per step, so a filled flat has no
                // gradient the router could read as the flood's own traversal order.
                if (h[nb] <= hc) h[nb] = hc + eps * (0.5 + hash01(nb, step * 7919 + 13));
                heap.push(nb, h[nb]);
            }
        }
    }

    function receivers() {
        for (let r = 0; r < N; r++) {
            if (!r_land[r]) { receiver[r] = r; recvLen[r] = meanEdge; continue; }
            const hr = h[r];
            let best = -1, bestS = 0, bestJ = -1;
            for (let j = adjOffset[r], jEnd = adjOffset[r + 1]; j < jEnd; j++) {
                const nb = adjList[j];
                const drop = hr - h[nb];
                if (drop <= 0) continue;
                const s = drop / neighborDist[j];
                if (s > bestS) { bestS = s; best = nb; bestJ = j; }
            }
            if (best < 0) { receiver[r] = r; recvLen[r] = meanEdge; }
            else { receiver[r] = best; recvLen[r] = neighborDist[bestJ]; }
        }
    }

    function buildStack() {
        donorOff.fill(0);
        for (let i = 0; i < N; i++) { const r = receiver[i]; if (r !== i) donorOff[r + 1]++; }
        for (let i = 0; i < N; i++) donorOff[i + 1] += donorOff[i];
        cursor.set(donorOff.subarray(0, N));
        for (let i = 0; i < N; i++) { const r = receiver[i]; if (r !== i) donorList[cursor[r]++] = i; }
        let tail = 0;
        for (let i = 0; i < N; i++) if (receiver[i] === i) stack[tail++] = i;
        for (let read = 0; read < tail; read++) {
            const c = stack[read];
            for (let d = donorOff[c], dEnd = donorOff[c + 1]; d < dEnd; d++) stack[tail++] = donorList[d];
        }
        return tail;
    }

    function accumulate(len) {
        area.fill(1);
        for (let k = len - 1; k >= 0; k--) {
            const i = stack[k];
            const r = receiver[i];
            if (r !== i) area[r] += area[i];
        }
    }

    function update(len) {
        for (let k = 0; k < len; k++) {
            const i = stack[k];
            if (!r_land[i]) continue;
            const r = receiver[i];
            if (r === i) { h[i] += dt * U[i]; continue; }
            const L = recvLen[i] * invMean;
            const f = r_k[i] * dt * Math.pow(area[i], m) / L;
            const hr = h[r];
            let next = (h[i] + dt * U[i] + f * hr) / (1 + f);
            if (next < hr) next = hr;
            h[i] = next;
        }
    }

    function diffuse() {
        if (alpha <= 0) return;
        for (let r = 0; r < N; r++) {
            if (!r_land[r]) continue;
            let sum = 0, cnt = 0;
            for (let j = adjOffset[r], jEnd = adjOffset[r + 1]; j < jEnd; j++) { sum += h[adjList[j]]; cnt++; }
            tmp[r] = cnt > 0 ? h[r] + alpha * (sum / cnt - h[r]) : h[r];
        }
        for (let r = 0; r < N; r++) if (r_land[r]) h[r] = tmp[r];
    }

    const t0 = (typeof performance !== 'undefined') ? performance.now() : Date.now();
    let stackLen = 0;
    const report = Math.max(1, Math.floor(steps / 20));
    for (let step = 0; step < steps; step++) {
        flood(step);
        receivers();
        stackLen = buildStack();
        accumulate(stackLen);
        update(stackLen);
        diffuse();
        if (onProgress && (step % report === 0 || step === steps - 1)) onProgress(step + 1, steps);
    }
    // One last fill and routing, so area and receiver describe the surface that is returned.
    flood(steps);
    receivers();
    stackLen = buildStack();
    accumulate(stackLen);

    const basin = new Int32Array(N);
    for (let k = 0; k < stackLen; k++) {
        const i = stack[k];
        const r = receiver[i];
        basin[i] = r === i ? i : basin[r];
    }
    const t1 = (typeof performance !== 'undefined') ? performance.now() : Date.now();
    return { h, area, receiver, recvLen, basin, meanEdge, solveMs: t1 - t0 };
}

// ─── Scaling to Orogen's elevation ────────────────────────────────

// Inverse of elevToHeightKm on land: a table over 0..6 km, built once.
let _invTable = null;
const INV_BINS = 6000;
function invHeightKm(km) {
    if (!_invTable) {
        const table = new Float32Array(INV_BINS + 1);
        let t = 0;
        const dtStep = 1 / 65536;
        for (let b = 0; b <= INV_BINS; b++) {
            const target = 6 * b / INV_BINS;
            while (t < 1 && elevToHeightKm(t) < target) t += dtStep;
            table[b] = Math.min(1, t);
        }
        _invTable = table;
    }
    if (km <= 0) return 0;
    if (km >= 6) return 1;
    const x = km / 6 * INV_BINS;
    const b = Math.floor(x);
    const f = x - b;
    return _invTable[b] * (1 - f) + _invTable[Math.min(INV_BINS, b + 1)] * f;
}

/**
 * Turn the solved units into Orogen's elevation field, plus the derived layers.
 * Land is scaled so its 99.5th percentile stands at peakKm; sea takes the class depth,
 * scaled so the deepest class sits at oceanDepthKm, with a short ramp down from the shore.
 */
export function toElevation(mesh, solved, r_land, r_class, legend, hop, opts) {
    const cls = legend.classes;
    const N = mesh.numRegions;
    const { h, area, receiver, recvLen, basin, meanEdge } = solved;
    const avgEdgeKm = opts.avgEdgeKm;

    const landVals = [];
    for (let r = 0; r < N; r++) if (r_land[r]) landVals.push(h[r]);
    landVals.sort((a, b) => a - b);
    const p995 = landVals.length ? landVals[Math.min(landVals.length - 1, Math.floor(landVals.length * 0.995))] : 0;
    const kmScale = p995 > 0 ? opts.peakKm / p995 : 0;

    let deepest = 0;
    for (const c of cls) if (c.sea && c.depthM > deepest) deepest = c.depthM;
    const depthScale = deepest > 0 ? (opts.oceanDepthKm * 1000) / deepest : 0;

    const r_elevation = new Float32Array(N);
    const slopeDeg = new Float32Array(N);
    const flowLog = new Float32Array(N);
    const basinOut = new Int32Array(N);
    let maxKm = 0;
    for (let r = 0; r < N; r++) {
        if (r_land[r]) {
            let km = h[r] * kmScale;
            if (km > 6) km = 6;
            if (km > maxKm) maxKm = km;
            let e = invHeightKm(km);
            if (e < 0.002) e = 0.002;
            r_elevation[r] = e;
            const rec = receiver[r];
            if (rec !== r) {
                const dz = (h[r] - h[rec]) * kmScale;
                const dx = recvLen[r] / meanEdge * avgEdgeKm;
                slopeDeg[r] = Math.atan2(Math.max(0, dz), Math.max(1e-6, dx)) * 180 / Math.PI;
            }
            flowLog[r] = Math.log10(Math.max(1, area[r]));
            basinOut[r] = basin[r];
        } else {
            const c = cls[r_class[r]];
            let depthKm = c.depthM / 1000 * depthScale;
            const f = Math.min(1, Math.max(0, (hop[r] - 0.5) / 2));
            depthKm *= f;
            if (depthKm < 0.005) depthKm = 0.005;
            let e = -depthKm / 10;
            if (e < -0.5) e = -0.5;
            r_elevation[r] = e;
            slopeDeg[r] = -1;
            flowLog[r] = -1;
            basinOut[r] = -1;
        }
    }
    return { r_elevation, slopeDeg, flowLog, basin: basinOut, kmScale, p995, maxKm };
}