tectonics.js
138 lines · in natura-src-1.0.0.zip as natura/tectonics.js · a module · its entry in the module reference · its test: tectonics-test.js
// tectonics.js -- Natura TECTONICS (SYSTEMS.md #10, ALGORITHMS.md): the deepest, slowest rung -- the
// ground itself is not fixed. The crust is broken into PLATES that drift, and where two plates CONVERGE
// their edges crumple upward into mountain ranges (the Himalaya, the Andes); where they pull APART the
// crust thins and RIFTS into valleys and new sea; and away from the edges the old highlands erode slowly
// back down (isostasy). Run it over millions of years and a flat craton grows a cordillera along its
// convergent margins. This is the boundary condition UNDER the terrain -- a continent that worldgen draws
// in an instant, tectonics makes EVOLVE. The bottom of the cross-scale stack: plate motion over Myr sets
// the elevation that weather rains on, that rivers carve, that life lives on.
//
// A Natura system spec (init/step/fields/events). Reads only the clock (deep time) + a seeded rng.
// Deterministic + node-testable. UMD (node test / window.Tectonics).
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Tectonics = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
// bilinear sample of a grid field at fractional cell (fi,fj), edge-clamped (for the drift backtrace)
function sampleE(F, N, fi, fj) {
fi = clamp(fi, 0, N - 1.001); fj = clamp(fj, 0, N - 1.001);
const i0 = fi | 0, j0 = fj | 0, tx = fi - i0, ty = fj - j0;
const a = F[j0 * N + i0], b = F[j0 * N + i0 + 1], c = F[(j0 + 1) * N + i0], d = F[(j0 + 1) * N + i0 + 1];
return a * (1 - tx) * (1 - ty) + b * tx * (1 - ty) + c * (1 - tx) * ty + d * tx * ty;
}
// value-noise fBm, for DOMAIN-WARPING the plate Voronoi so the boundaries are organic + irregular
// (real plate margins wiggle; a bare Voronoi gives straight seams) rather than clean geometric lines.
function nz(ix, iy, s) { let h = (Math.imul(ix | 0, 374761393) + Math.imul(iy | 0, 668265263) + Math.imul(s | 0, 2147483647)) | 0; h = Math.imul(h ^ (h >>> 13), 1274126177); return ((h ^ (h >>> 16)) >>> 0) / 4294967296; }
function vnz(x, y, s) { const ix = Math.floor(x), iy = Math.floor(y), fx = x - ix, fy = y - iy, u = fx * fx * (3 - 2 * fx), v = fy * fy * (3 - 2 * fy);
const a = nz(ix, iy, s), b = nz(ix + 1, iy, s), c = nz(ix, iy + 1, s), d = nz(ix + 1, iy + 1, s); return (a * (1 - u) + b * u) * (1 - v) + (c * (1 - u) + d * u) * v; }
function fbm2(x, y, s) { return 0.5 * vnz(x, y, s) + 0.32 * vnz(x * 2.03, y * 2.03, s + 7) + 0.18 * vnz(x * 4.1, y * 4.1, s + 13); }
// make a Tectonics system. config: { N, plates, upliftK, erodeK, elevation(seed field), driftK, name }.
function makeTectonics(config) {
config = config || {};
const N = config.N || 64;
const nPlates = config.plates != null ? config.plates : 5;
const upliftK = config.upliftK != null ? config.upliftK : 0.03; // uplift per convergence per Myr
const erodeK = config.erodeK != null ? config.erodeK : 0.004; // isostatic decay of high ground per Myr
const driftK = config.driftK != null ? config.driftK : 0.00013; // plate DRIFT speed (world-units per Myr): the continents move -- moderate, so the crust slides coherently rather than smearing
const upliftSpread = config.upliftSpread != null ? config.upliftSpread : 2; // radius over which convergence uplift is spread -> BROAD belts, not knife-edge lines
const smoothK = config.smoothK != null ? config.smoothK : 0.1; // light per-step elevation smoothing so ranges + coasts read soft, not as harsh scored lines
const warpK = config.warpK != null ? config.warpK : 0.28; // domain-warp strength for organic (non-straight) plate boundaries
const outSmooth = config.outSmooth != null ? config.outSmooth : 2; // passes of smoothing on the OUTPUT elevation (see fields): the internal crust carries cell-scale advection noise; consumers should see a coherent landform
const elevFn = config.elevation || ((x, y) => 0.05); // starting heightfield (a low craton)
const gx = (i) => (i / (N - 1)) * 2 - 1;
return {
name: config.name || 'tectonics',
domain: 'physical',
scale: 'planetary',
params: config,
init(rng) {
// seed the plates (Voronoi), give each a drift velocity, assign every cell to its nearest plate
const seeds = [];
for (let p = 0; p < nPlates; p++) seeds.push({ x: rng() * 2 - 1, y: rng() * 2 - 1, vx: (rng() - 0.5), vy: (rng() - 0.5) });
const plate = new Int8Array(N * N), elev = new Float32Array(N * N);
// a fixed organic DOMAIN WARP for the plate boundaries -- the seams wiggle through it as they migrate
const ws = (rng() * 1e6) | 0, warpX = new Float32Array(N * N), warpY = new Float32Array(N * N);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const x = gx(i), y = gx(j), k = j * N + i;
warpX[k] = warpK * (fbm2(x * 1.6 + 3, y * 1.6, ws) - 0.5); warpY[k] = warpK * (fbm2(x * 1.6, y * 1.6 + 3, ws + 41) - 0.5); }
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i, wx = gx(i) + warpX[k], wy = gx(j) + warpY[k]; let best = 0, bd = 1e9;
for (let p = 0; p < nPlates; p++) { const dx = wx - seeds[p].x, dy = wy - seeds[p].y, d = dx * dx + dy * dy; if (d < bd) { bd = d; best = p; } }
plate[k] = best; elev[k] = elevFn(gx(i), gx(j));
}
return { N, rng, seeds, plate, elev, baseline: Float32Array.from(elev), elevOut: Float32Array.from(elev), warpX, warpY, active: 0 };
},
step(st, ctx, dt) {
const dMyr = clamp(dt / (1e6 * 365 * 86400), 0, 50); // millions of years this step
if (dMyr <= 0) return st;
const { plate, elev, seeds, baseline, warpX, warpY } = st;
// 1. the PLATES DRIFT AS BODIES: move each seed (bounce at the edges) and RE-VORONOI, so the
// boundaries themselves MIGRATE -- a mountain belt or a rift is born, travels, and dies, and the
// whole tectonic pattern keeps reshuffling instead of pinning fault lines to fixed spots forever.
if (driftK > 0) {
const d = driftK * dMyr, gs = (N - 1) / 2;
for (const s of seeds) { s.x += s.vx * d; s.y += s.vy * d;
if (s.x < -1.15 || s.x > 1.15) s.vx *= -1; if (s.y < -1.15 || s.y > 1.15) s.vy *= -1; }
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i, x = gx(i) + warpX[k], y = gx(j) + warpY[k]; let best = 0, bd = 1e9;
for (let p = 0; p < seeds.length; p++) { const dx = x - seeds[p].x, dy = y - seeds[p].y, dd = dx * dx + dy * dy; if (dd < bd) { bd = dd; best = p; } }
plate[k] = best;
}
const nE = new Float32Array(N * N), nB = new Float32Array(N * N); // advect the crust WITH its (moving) plate
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const k = j * N + i, v = seeds[plate[k]];
nE[k] = sampleE(elev, N, i - v.vx * d * gs, j - v.vy * d * gs); nB[k] = sampleE(baseline, N, i - v.vx * d * gs, j - v.vy * d * gs); }
elev.set(nE); baseline.set(nB);
}
// 2. CONVERGENCE at the (migrated) boundaries, then SPREAD it over a few cells -> broad mountain
// BELTS instead of knife-edge lines.
const raw = new Float32Array(N * N); let active = 0;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i, pA = plate[k], vA = seeds[pA]; let conv = 0;
for (let dy = -1; dy <= 1; dy++) for (let dx = -1; dx <= 1; dx++) {
if (!dx && !dy) continue; const ni = i + dx, nj = j + dy; if (ni < 0 || ni >= N || nj < 0 || nj >= N) continue;
const pB = plate[nj * N + ni]; if (pB === pA) continue; const vB = seeds[pB], nl = Math.hypot(dx, dy);
conv += ((vA.vx - vB.vx) * dx + (vA.vy - vB.vy) * dy) / nl;
}
raw[k] = conv; if (conv > 0.05) active++;
}
const up = new Float32Array(N * N), R = upliftSpread;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { let s = 0, n = 0;
for (let dy = -R; dy <= R; dy++) for (let dx = -R; dx <= R; dx++) { const ni = i + dx, nj = j + dy; if (ni < 0 || ni >= N || nj < 0 || nj >= N) continue; s += raw[nj * N + ni]; n++; }
up[j * N + i] = s / n;
}
// 3. uplift the belts / rift the divergences (rifts SOFTENED), erode toward baseline
for (let k = 0; k < N * N; k++) { const u = up[k] > 0 ? up[k] * 1.8 : up[k] * 0.4; // x1.8 restores height lost to the spread; rifts x0.4 (gentler trenches)
elev[k] += upliftK * u * dMyr; elev[k] -= erodeK * (elev[k] - baseline[k]) * dMyr; elev[k] = clamp(elev[k], -0.4, 1.2); }
// 4. a light smoothing so ranges + coasts read soft, not as harsh scored lines
if (smoothK > 0) { const nb = new Float32Array(N * N);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const k = j * N + i;
const il = Math.max(0, i - 1), ir = Math.min(N - 1, i + 1), jd = Math.max(0, j - 1), ju = Math.min(N - 1, j + 1);
nb[k] = elev[k] * (1 - smoothK) + smoothK * 0.25 * (elev[j * N + il] + elev[j * N + ir] + elev[jd * N + i] + elev[ju * N + i]); }
elev.set(nb); }
st.active = active;
// OUTPUT elevation: a smoothed copy the consumers read. The per-plate semi-Lagrangian advection injects
// cell-scale noise at every (migrating) plate boundary every step -- faster than the internal smoothing
// can remove it -- which would dither the biome map and speckle the hillshade. The landforms (belts,
// basins, coasts) are low-frequency and survive; the noise does not. The rough field stays internal.
{ const out = st.elevOut, t2 = new Float32Array(N * N); out.set(elev);
for (let it = 0; it < outSmooth; it++) {
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const il = Math.max(0, i - 1), ir = Math.min(N - 1, i + 1), jd = Math.max(0, j - 1), ju = Math.min(N - 1, j + 1);
t2[j * N + i] = 0.2 * (out[j * N + i] + out[j * N + il] + out[j * N + ir] + out[jd * N + i] + out[ju * N + i]); }
out.set(t2); } }
return st;
},
// OUT: elevation (the evolving heightfield, smoothed for consumers -- terrain/worldgen/biome build on it), plate
fields(st) { return { N: st.N, elevation: st.elevOut, plate: st.plate }; },
events(st) { return st.active > 8 ? [{ kind: 'orogeny', cells: st.active }] : []; }
};
}
return { makeTectonics };
});