Stratigraphy: Drill a Core

// Stratigraphy Builder — drill-a-core
//
// Time on the y-axis (millions of years, Ma), space horizontal. As simulated
// time advances, sediment layers deposit downward from the *top* of the basin
// (newest at top, older at bottom — geological "deeper = older" convention).
// We store sediments column-by-column so faults and tilts that happened *after*
// a layer formed actually shear the layer at every x.
//
// Sea level controls facies: when sea level (set by mouseY) is HIGH the column
// at every x that sits below sea level deposits marine sediments (blue-greys,
// shales, limestones). When sea level is LOW the column deposits terrestrial
// sediments (sandy yellows / browns).
//
// Tectonic events: roughly every 15 Ma we pick either a gentle TILT (rotates
// future deposition's "deposition surface") or a FAULT (a vertical offset
// across a random x split — everything currently in the basin shifts on one
// side). Both produce clear unconformities and offsets you can see when you
// drill a core.
//
// Drilling: after ~25 Ma of basin time, the user can click any x to log a
// vertical core at that x. The drilled core highlights and a side-panel
// shows the stratigraphic column at that location (most recent at top).

const COLS = 220;          // horizontal resolution of the basin
const ROWS = 320;          // vertical resolution (deep stack of layers)
const MA_PER_FRAME = 0.06; // simulated time advance per frame
const LAYER_EVERY_MA = 0.5; // a fresh sediment band every ~0.5 Ma
const DRILL_UNLOCK_MA = 25; // can't drill until enough history exists
const TECTONIC_EVERY_MA = 15;

// Each cell stores the deposited layer band. 0 = empty (no sediment yet).
// Positive integer = bandId; we also keep a parallel facies/color table.
let grid;            // Int16Array COLS*ROWS, layer band id (0 = empty)
let facies;          // Uint8Array — 1=marine, 2=terrestrial, 3=transitional, per band id
let bandColor;       // Float32Array — per-band color jitter (0..1)
let bandMa;          // Float32Array — Ma at which each band formed
let nextBand;        // next band id to allocate
// per-column "fill height" — how many rows from the bottom are filled.
// We fill UP from row ROWS-1 toward row 0. Row 0 is the surface.
let columnFill;      // Int16Array COLS — number of filled rows at each column
let deformOffset;    // Float32Array COLS — fractional row offset applied to future deposits (tilts)
let ageMa;           // simulated time elapsed
let lastLayerMa;     // when we last formed a layer
let lastTectonicMa;  // when we last did a tectonic event
let cores;           // array of { x, atMa }
let lastTectonicLabel;
let lastTectonicLabelT;
let drillUnlockedFlash;

// Sea level state. seaRow is the basin row that water surface sits at
// (smaller = higher sea level). Driven smoothly by mouseY input.
let seaRow;
let seaRowTarget;

function init() {
  grid = new Int16Array(COLS * ROWS);
  facies = new Uint8Array(2048);
  bandColor = new Float32Array(2048);
  bandMa = new Float32Array(2048);
  nextBand = 1;
  columnFill = new Int16Array(COLS);
  deformOffset = new Float32Array(COLS);
  ageMa = 0;
  lastLayerMa = -LAYER_EVERY_MA;
  lastTectonicMa = 0;
  cores = [];
  lastTectonicLabel = "";
  lastTectonicLabelT = 0;
  drillUnlockedFlash = 0;
  seaRow = ROWS * 0.5;
  seaRowTarget = ROWS * 0.5;
}

function allocBand(faciesKind) {
  if (nextBand >= facies.length) {
    // grow tables
    const grow = facies.length * 2;
    const f = new Uint8Array(grow); f.set(facies); facies = f;
    const c = new Float32Array(grow); c.set(bandColor); bandColor = c;
    const m = new Float32Array(grow); m.set(bandMa); bandMa = m;
  }
  const id = nextBand++;
  facies[id] = faciesKind;
  bandColor[id] = Math.random();
  bandMa[id] = ageMa;
  return id;
}

function depositLayer() {
  // For each column, decide facies based on whether that column's *current
  // deposition surface* sits below the sea level. seaRow uses basin-row
  // coordinates where 0 is "high in the air" and ROWS is the basin floor.
  // A column's deposition surface = (ROWS - columnFill[x]) + deformOffset[x].
  // If surface > seaRow it is submerged (marine); else terrestrial.
  // We allocate up to 3 bands for this Ma slice (marine / terrestrial /
  // transitional) and pick one per column so the layer reads as a single
  // time-band but with facies variation across space.
  const marineId = allocBand(1);
  const terrId   = allocBand(2);
  const transId  = allocBand(3);
  // randomize thickness 1..3 rows
  const thickness = 1 + ((Math.random() * 2.7) | 0);
  for (let x = 0; x < COLS; x++) {
    const surfaceRow = (ROWS - columnFill[x]) + deformOffset[x];
    let id;
    if (surfaceRow > seaRow + 6) id = marineId;
    else if (surfaceRow > seaRow - 6) id = transId;
    else id = terrId;
    // Fill `thickness` cells upward from current fill top.
    for (let t = 0; t < thickness; t++) {
      const fill = columnFill[x];
      if (fill >= ROWS) break;
      const row = ROWS - 1 - fill;
      grid[row * COLS + x] = id;
      columnFill[x] = fill + 1;
    }
  }
}

function applyTilt() {
  // Rotate future deposition surface across the basin. We just bias
  // deformOffset linearly; positive on one side means "this side is being
  // pushed down so accumulates more thickness".
  const slope = (Math.random() - 0.5) * 0.06; // rows per column
  const mid = COLS * 0.5;
  for (let x = 0; x < COLS; x++) deformOffset[x] += (x - mid) * slope;
  lastTectonicLabel = `tilt ${slope > 0 ? "→ down-right" : "→ down-left"}`;
  lastTectonicLabelT = 240;
}

function applyFault() {
  // Vertical offset across a random x. Shift all deposited cells on the
  // hangingwall side up or down by `offset` rows. We also shift columnFill.
  const fx = (COLS * (0.2 + Math.random() * 0.6)) | 0;
  const offset = (1 + ((Math.random() * 6) | 0)) * (Math.random() < 0.5 ? -1 : 1);
  const side = Math.random() < 0.5 ? "L" : "R";
  for (let x = 0; x < COLS; x++) {
    const onSide = side === "L" ? x < fx : x > fx;
    if (!onSide) continue;
    // shift this column's filled cells by `offset` rows (down = +offset)
    if (offset > 0) {
      // shift downward: rows near the top move toward bottom — but we only
      // have ROWS rows, so cells falling off the bottom are lost.
      for (let r = ROWS - 1; r >= 0; r--) {
        const src = r - offset;
        grid[r * COLS + x] = src >= 0 ? grid[src * COLS + x] : 0;
      }
      columnFill[x] = Math.min(ROWS, columnFill[x] + offset);
    } else {
      const off = -offset;
      // shift upward: rows shift toward 0; cells exposed at the bottom become 0
      for (let r = 0; r < ROWS; r++) {
        const src = r + off;
        grid[r * COLS + x] = src < ROWS ? grid[src * COLS + x] : 0;
      }
      columnFill[x] = Math.max(0, columnFill[x] - off);
    }
  }
  lastTectonicLabel = `fault @ x≈${(fx / COLS * 100) | 0}%  (${side}-side ${offset > 0 ? "down" : "up"} ${Math.abs(offset)})`;
  lastTectonicLabelT = 240;
}

function maybeTectonic() {
  if (ageMa - lastTectonicMa < TECTONIC_EVERY_MA) return;
  lastTectonicMa = ageMa;
  if (Math.random() < 0.5) applyTilt(); else applyFault();
}

// Color a band cell. Returns an "rgb(...)" string. We mix the band's facies
// base palette with the per-band jitter so each layer is recognizable.
function colorFor(bandId) {
  if (!bandId) return "#0a0e14";
  const f = facies[bandId];
  const j = bandColor[bandId];
  if (f === 1) {
    // marine — blue-grey to dark slate
    const r = 50 + (j * 30) | 0;
    const g = 70 + (j * 35) | 0;
    const b = 100 + (j * 55) | 0;
    return `rgb(${r},${g},${b})`;
  } else if (f === 2) {
    // terrestrial — sandy yellow-brown
    const r = 150 + (j * 70) | 0;
    const g = 110 + (j * 55) | 0;
    const b = 65 + (j * 40) | 0;
    return `rgb(${r},${g},${b})`;
  } else {
    // transitional — greenish-grey (coastal, shale-mudstone)
    const r = 90 + (j * 30) | 0;
    const g = 105 + (j * 35) | 0;
    const b = 85 + (j * 25) | 0;
    return `rgb(${r},${g},${b})`;
  }
}

function tick({ ctx, width: W, height: H, input }) {
  // ----- input -----
  // mouseY in the basin region drives target sea level. We map y across the
  // basin draw region (computed below). For now, latch on input.mouseY:
  // a low mouseY (near top) = high sea level (seaRow small), and vice versa.
  const my = input && typeof input.mouseY === "number" ? input.mouseY : -1;
  if (my >= 0 && my <= H) {
    const t = Math.max(0, Math.min(1, my / H));
    seaRowTarget = ROWS * (0.15 + t * 0.7); // mouse drives sea level
  }
  // smooth follow so the sea level glides
  seaRow += (seaRowTarget - seaRow) * 0.05;

  // ----- simulate -----
  ageMa += MA_PER_FRAME;
  // deposit a layer band whenever enough simulated time has passed
  while (ageMa - lastLayerMa >= LAYER_EVERY_MA) {
    lastLayerMa += LAYER_EVERY_MA;
    depositLayer();
  }
  maybeTectonic();

  // unlock drill flash
  if (ageMa >= DRILL_UNLOCK_MA && drillUnlockedFlash === 0 && cores.length === 0) {
    drillUnlockedFlash = 180;
  }
  if (drillUnlockedFlash > 0) drillUnlockedFlash--;

  // ----- layout -----
  const pad = 10;
  const headerH = 28;
  const panelW = Math.min(170, Math.max(110, (W * 0.22) | 0));
  const basinX = pad;
  const basinY = pad + headerH;
  const basinW = W - panelW - pad * 3;
  const basinH = H - basinY - pad;
  const cellW = basinW / COLS;
  const cellH = basinH / ROWS;

  // handle clicks → drill. Map click x to column.
  let clickedCore = -1;
  if (input && typeof input.consumeClicks === "function") {
    const clicks = input.consumeClicks();
    if (clicks.length > 0 && ageMa >= DRILL_UNLOCK_MA) {
      for (let ci = 0; ci < clicks.length; ci++) {
        const cx0 = clicks[ci].x;
        const cy0 = clicks[ci].y;
        if (cx0 >= basinX && cx0 <= basinX + basinW && cy0 >= basinY && cy0 <= basinY + basinH) {
          const col = Math.max(0, Math.min(COLS - 1, ((cx0 - basinX) / cellW) | 0));
          // dedupe: replace nearest if within a few columns
          let replaced = false;
          for (let i = 0; i < cores.length; i++) {
            if (Math.abs(cores[i].x - col) < 4) { cores[i] = { x: col, atMa: ageMa }; replaced = true; clickedCore = i; break; }
          }
          if (!replaced) { cores.push({ x: col, atMa: ageMa }); clickedCore = cores.length - 1; }
          // cap to last 6 cores
          if (cores.length > 6) cores.shift();
        }
      }
    }
  }

  // ----- draw background -----
  ctx.fillStyle = "#0a0e14";
  ctx.fillRect(0, 0, W, H);

  // basin area background (atmosphere / future deposition zone)
  ctx.fillStyle = "#0d1218";
  ctx.fillRect(basinX, basinY, basinW, basinH);

  // ----- draw deposited cells as horizontal runs per row -----
  // For each row, scan columns and coalesce same-bandId runs to reduce fillRect calls.
  for (let r = 0; r < ROWS; r++) {
    let runStart = -1;
    let runBand = 0;
    const rowOff = r * COLS;
    for (let x = 0; x < COLS; x++) {
      const b = grid[rowOff + x];
      if (b !== runBand) {
        if (runBand !== 0 && runStart >= 0) {
          ctx.fillStyle = colorFor(runBand);
          ctx.fillRect(basinX + runStart * cellW, basinY + r * cellH, (x - runStart) * cellW + 0.5, cellH + 0.5);
        }
        runBand = b;
        runStart = x;
      }
    }
    if (runBand !== 0 && runStart >= 0) {
      ctx.fillStyle = colorFor(runBand);
      ctx.fillRect(basinX + runStart * cellW, basinY + r * cellH, (COLS - runStart) * cellW + 0.5, cellH + 0.5);
    }
  }

  // ----- sea level line -----
  const seaY = basinY + (seaRow / ROWS) * basinH;
  // water tint above the current basin surface but below sea level
  ctx.fillStyle = "rgba(70,120,180,0.10)";
  ctx.fillRect(basinX, basinY, basinW, Math.max(0, seaY - basinY));
  // animated sea line
  ctx.strokeStyle = "rgba(120,200,255,0.85)";
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  for (let x = 0; x <= basinW; x += 6) {
    const wy = seaY + Math.sin((x * 0.05) + ageMa * 1.5) * 1.6;
    if (x === 0) ctx.moveTo(basinX + x, wy); else ctx.lineTo(basinX + x, wy);
  }
  ctx.stroke();
  ctx.fillStyle = "rgba(120,200,255,0.9)";
  ctx.font = "10px monospace";
  ctx.textAlign = "right";
  ctx.fillText("sea level", basinX + basinW - 6, seaY - 4);

  // ----- drill highlights -----
  for (let i = 0; i < cores.length; i++) {
    const c = cores[i];
    const cx = basinX + c.x * cellW + cellW * 0.5;
    // glow
    ctx.fillStyle = "rgba(255,220,140,0.13)";
    ctx.fillRect(cx - 3, basinY, 6, basinH);
    // bright center line
    ctx.strokeStyle = i === clickedCore ? "#fff7c0" : "rgba(255,220,140,0.85)";
    ctx.lineWidth = 1;
    ctx.beginPath(); ctx.moveTo(cx, basinY); ctx.lineTo(cx, basinY + basinH); ctx.stroke();
    // label dot at top
    ctx.fillStyle = "#ffd86b";
    ctx.beginPath(); ctx.arc(cx, basinY - 2, 3, 0, Math.PI * 2); ctx.fill();
    ctx.fillStyle = "#0a0e14";
    ctx.font = "9px monospace";
    ctx.textAlign = "center";
    ctx.fillText(String(i + 1), cx, basinY + 1);
  }

  // basin frame
  ctx.strokeStyle = "#333";
  ctx.lineWidth = 1;
  ctx.strokeRect(basinX - 0.5, basinY - 0.5, basinW + 1, basinH + 1);

  // ----- header / HUD -----
  ctx.textAlign = "left";
  ctx.fillStyle = "#dde";
  ctx.font = "13px monospace";
  ctx.fillText("Stratigraphy: Drill a Core", pad, pad + 14);

  ctx.font = "11px monospace";
  // Layout the header strip left-to-right with measured widths so it doesn't
  // collide on narrow canvases.
  let hx = pad + ctx.measureText("Stratigraphy: Drill a Core").width + 16;
  ctx.fillStyle = "#9af";
  const ageStr = `${ageMa.toFixed(1)} Ma`;
  ctx.fillText(ageStr, hx, pad + 14);
  hx += ctx.measureText(ageStr).width + 16;

  // current facies indicator (what is depositing RIGHT NOW at the basin center)
  const midX = (COLS / 2) | 0;
  const curSurface = (ROWS - columnFill[midX]) + deformOffset[midX];
  let curFaciesLbl, curFaciesCol;
  if (curSurface > seaRow + 6) { curFaciesLbl = "marine"; curFaciesCol = "#7cf"; }
  else if (curSurface > seaRow - 6) { curFaciesLbl = "coastal"; curFaciesCol = "#9c9"; }
  else { curFaciesLbl = "terrestrial"; curFaciesCol = "#dba"; }
  const headerRightLimit = W - panelW - pad * 2 - 4;
  if (hx + 110 < headerRightLimit) {
    ctx.fillStyle = "#888";
    ctx.fillText("depositing:", hx, pad + 14);
    hx += ctx.measureText("depositing:").width + 4;
    ctx.fillStyle = curFaciesCol;
    ctx.fillText(curFaciesLbl, hx, pad + 14);
    hx += ctx.measureText(curFaciesLbl).width + 16;
  }
  if (hx + 60 < headerRightLimit) {
    ctx.fillStyle = "#888";
    ctx.fillText(`cores: ${cores.length}`, hx, pad + 14);
  }

  // hint line just below the header
  ctx.font = "10px monospace";
  if (ageMa < DRILL_UNLOCK_MA) {
    ctx.fillStyle = "#665";
    ctx.fillText(`drill unlocks at ${DRILL_UNLOCK_MA} Ma · move mouse Y to set sea level`,
      pad, pad + 26);
  } else {
    ctx.fillStyle = drillUnlockedFlash > 0
      ? `rgba(255,220,140,${0.4 + (drillUnlockedFlash / 180) * 0.6})`
      : "#776";
    ctx.fillText("click in the basin to drill a core · mouse Y sets sea level",
      pad, pad + 26);
  }

  // tectonic event toast
  if (lastTectonicLabelT > 0) {
    lastTectonicLabelT--;
    const alpha = Math.min(1, lastTectonicLabelT / 60);
    ctx.fillStyle = `rgba(255,160,90,${alpha})`;
    ctx.font = "11px monospace";
    ctx.textAlign = "right";
    ctx.fillText("⚡ " + lastTectonicLabel, W - panelW - pad * 2 - 4, pad + 26);
    ctx.textAlign = "left";
  }

  // ----- side panel: core logs -----
  const px0 = W - panelW - pad;
  const py0 = basinY;
  const pw = panelW;
  const ph = basinH;
  ctx.fillStyle = "#0c1118";
  ctx.fillRect(px0, py0, pw, ph);
  ctx.strokeStyle = "#222";
  ctx.strokeRect(px0 + 0.5, py0 + 0.5, pw - 1, ph - 1);
  ctx.fillStyle = "#ccd";
  ctx.font = "11px monospace";
  ctx.fillText("core logs", px0 + 8, py0 + 14);

  if (cores.length === 0) {
    ctx.fillStyle = "#555";
    ctx.font = "10px monospace";
    ctx.fillText(ageMa < DRILL_UNLOCK_MA
      ? `(wait ${(DRILL_UNLOCK_MA - ageMa).toFixed(1)} Ma)`
      : "click basin to drill", px0 + 8, py0 + 30);
  } else {
    // Render up to 4 most-recent cores side by side as vertical strips
    const show = cores.slice(-4);
    const stripW = (pw - 16) / show.length;
    const stripTop = py0 + 24;
    const stripBot = py0 + ph - 10;
    const stripH = stripBot - stripTop;
    for (let i = 0; i < show.length; i++) {
      const c = show[i];
      const sx = px0 + 8 + i * stripW;
      // Render the column at c.x: row 0..ROWS-1 mapped to stripTop..stripBot
      // Use a single-pixel-wide vertical scan, coalescing runs of same band.
      let runStart = -1;
      let runBand = 0;
      for (let r = 0; r < ROWS; r++) {
        const b = grid[r * COLS + c.x];
        if (b !== runBand) {
          if (runBand !== 0 && runStart >= 0) {
            const y0 = stripTop + (runStart / ROWS) * stripH;
            const y1 = stripTop + (r / ROWS) * stripH;
            ctx.fillStyle = colorFor(runBand);
            ctx.fillRect(sx, y0, stripW - 4, y1 - y0 + 0.5);
          }
          runBand = b;
          runStart = r;
        }
      }
      if (runBand !== 0 && runStart >= 0) {
        const y0 = stripTop + (runStart / ROWS) * stripH;
        const y1 = stripTop + stripH;
        ctx.fillStyle = colorFor(runBand);
        ctx.fillRect(sx, y0, stripW - 4, y1 - y0 + 0.5);
      }
      // strip frame + label
      ctx.strokeStyle = "#333";
      ctx.strokeRect(sx - 0.5, stripTop - 0.5, stripW - 4 + 1, stripH + 1);
      ctx.fillStyle = "#bbb";
      ctx.font = "9px monospace";
      ctx.textAlign = "center";
      const idxInArr = cores.length - show.length + i;
      ctx.fillText(`#${idxInArr + 1}`, sx + (stripW - 4) / 2, stripTop - 2);
      ctx.font = "8px monospace";
      ctx.fillStyle = "#666";
      ctx.fillText(`${c.atMa.toFixed(0)}Ma`, sx + (stripW - 4) / 2, stripBot + 8);
      ctx.textAlign = "left";
    }
  }

  // ----- sea level indicator gutter on the left edge of basin -----
  // small triangle at sea level on the y-axis to make the control obvious
  ctx.fillStyle = "rgba(120,200,255,0.9)";
  ctx.beginPath();
  ctx.moveTo(basinX - 1, seaY);
  ctx.lineTo(basinX - 8, seaY - 4);
  ctx.lineTo(basinX - 8, seaY + 4);
  ctx.closePath();
  ctx.fill();
}