Snowball Earth: why a dimmer sun freezes the planet and today's cannot thaw it

// game.js: dim the sun over an ice-capped Earth. A 90-band energy balance model runs live; turn the sun
// down or up and watch the ice line and the global mean temperature. Units: °C, years, % of today's sun.
// The notebook writes the parameters, the grid, today's state, the coastline, and the land mask into #data.
"use strict";

const DATA = JSON.parse(document.getElementById("data").textContent);
const {A, B, C, D, Q0, S2, a0, a2, b0, Ts} = DATA.params;
const {xe, x, dt_yr: DT, years_per_s: SPEED, view: VIEW} = DATA;
const N = x.length;
const Q_MIN = Number(document.getElementById("sun").min), Q_MAX = Number(document.getElementById("sun").max);   // the slider's range
const dx = x.map((_, i) => xe[i + 1] - xe[i]);
const P2 = x.map((v) => (3 * v * v - 1) / 2);
const S = P2.map((p) => 1 + S2 * p);
const K = x.slice(1).map((_, i) => D * (1 - xe[i + 1] ** 2) / (x[i + 1] - x[i]));   // conductance of each inner edge
const EDGE_LAT = xe.map((v) => Math.round(Math.asin(v) * 180 / Math.PI));

// ---- the state: everything that changes, in one object
const state = {T: Float64Array.from(DATA.T0), q: 100, year: 0, running: true, trail: [], spin: 0};

// ---- one tick: formula 1 on the bands, by forward Euler
const flux = new Float64Array(N);
function step(T, q, dt) {
  flux.fill(0);
  for (let i = 0; i < N - 1; i++) {                     // heat across the edge between band i and band i + 1
    const F = K[i] * (T[i + 1] - T[i]);
    flux[i] += F;
    flux[i + 1] -= F;
  }
  for (let i = 0; i < N; i++) {
    const a = T[i] < Ts ? b0 : a0 + a2 * P2[i];         // coalbedo: the share of sunlight absorbed, less on ice
    T[i] += dt / C * (flux[i] / dx[i] + q / 100 * Q0 * S[i] * a - (A + B * T[i]));
  }
}

// ---- the readouts
const globalMean = (T) => T.reduce((sum, t, i) => sum + t * dx[i], 0) / 2;   // bands weighted by area
function iceLine(T) {                                   // equatorward edge of the lowest ice band in the north
  for (let i = N / 2; i < N; i++) if (T[i] < Ts) return EDGE_LAT[i];
  return 90;
}
const minus = (s) => s.replace("-", "−");
const iceText = (e) => e === 0 ? "frozen over" : e === 90 ? "no ice" : `ice line ${e}°`;

// ---- the land mask: 720 × 360 cells of 0.5°, one bit each, rows from 90° N, columns from 180° W
const MASK = (() => {
  const bytes = atob(DATA.mask.bits), bits = new Uint8Array(DATA.mask.w * DATA.mask.h);
  for (let i = 0; i < bits.length; i++) bits[i] = (bytes.charCodeAt(i >> 3) >> (7 - (i & 7))) & 1;
  return bits;
})();

// ---- drawing
const $ = (id) => document.getElementById(id);
const views = {globe: $("globe"), strip: $("strip")};
const raster = document.createElement("canvas");
let dpr = 1, pixels = null;

function colors() {
  const css = getComputedStyle(document.documentElement);
  return Object.fromEntries(["bg", "fg", "ink", "accent", "second", "muted", "ice", "coast"].map((k) => [k, css.getPropertyValue(`--${k}`).trim()]));
}
const rgb = (hex) => [1, 3, 5].map((i) => parseInt(hex.slice(i, i + 2), 16));
const mix = (bg, fg, alpha) => rgb(bg).map((b, i) => Math.round((1 - alpha) * b + alpha * rgb(fg)[i]));

function layout() {
  dpr = window.devicePixelRatio || 1;
  for (const canvas of Object.values(views)) {
    const box = canvas.getBoundingClientRect();
    canvas.width = Math.round(box.width * dpr);
    canvas.height = Math.round(box.height * dpr);
  }
  // Every pixel of the disk, mapped back to latitude and to longitude east of the view's center, once per size.
  const box = views.globe.getBoundingClientRect();
  const n = Math.max(2, Math.min(720, Math.round((Math.min(box.width, box.height) - 8) * dpr)));
  raster.width = raster.height = n;
  const p = VIEW.lat0 * Math.PI / 180;
  pixels = {n, inside: new Uint8Array(n * n), band: new Uint8Array(n * n), row: new Uint16Array(n * n), lon: new Float32Array(n * n)};
  for (let j = 0; j < n; j++) for (let i = 0; i < n; i++) {
    const u = 2 * (i + 0.5) / n - 1, v = 1 - 2 * (j + 0.5) / n, z2 = 1 - u * u - v * v, k = j * n + i;
    if (z2 <= 0) continue;
    const z = Math.sqrt(z2);
    const lat = Math.asin(v * Math.cos(p) + z * Math.sin(p)) * 180 / Math.PI;
    pixels.inside[k] = 1;
    pixels.band[k] = Math.min(N - 1, Math.floor((lat + 90) / 2));
    pixels.row[k] = Math.min(DATA.mask.h - 1, Math.floor((90 - lat) * 2));
    pixels.lon[k] = Math.atan2(u, z * Math.cos(p) - v * Math.sin(p)) * 180 / Math.PI;
  }
  render();
}

function begin(canvas) {                                // a 2D context in CSS pixels
  const ctx = canvas.getContext("2d");
  ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
  return [ctx, canvas.width / dpr, canvas.height / dpr];
}

function drawGlobe(c, lon0, edge) {
  const [ctx, W, H] = begin(views.globe);
  ctx.fillStyle = c.bg;
  ctx.fillRect(0, 0, W, H);
  const R = Math.min(W, H) / 2 - 4, cx = W / 2, cy = H / 2;
  if (R < 10) return;
  // the raster: land from the mask, ice from the band's temperature
  const {n, inside, band, row, lon} = pixels, sea = mix(c.bg, c.second, 0.22), land = mix(c.bg, c.muted, 0.45), ice = rgb(c.ice);
  const rctx = raster.getContext("2d"), image = rctx.createImageData(n, n), px = image.data;
  for (let k = 0; k < n * n; k++) {
    if (!inside[k]) continue;
    const col = Math.floor((((lon[k] + lon0 + 180) % 360 + 360) % 360) * 2);
    const color = state.T[band[k]] < Ts ? ice : MASK[row[k] * DATA.mask.w + col] ? land : sea;
    px.set(color, 4 * k);
    px[4 * k + 3] = 255;
  }
  rctx.putImageData(image, 0, 0);
  ctx.drawImage(raster, cx - R, cy - R, 2 * R, 2 * R);

  // the orthographic projection forward: (lat, lon) to the disk, and whether the point faces the viewer
  const p = VIEW.lat0 * Math.PI / 180;
  const project = (lat, lng) => {
    const f = lat * Math.PI / 180, l = (lng - lon0) * Math.PI / 180;
    const X = Math.cos(f) * Math.cos(l), Y = Math.cos(f) * Math.sin(l), Z = Math.sin(f);
    return [cx + R * Y, cy - R * (Z * Math.cos(p) - X * Math.sin(p)), X * Math.cos(p) + Z * Math.sin(p) > 0];
  };
  const trace = (points) => {                           // a polyline, lifted wherever it goes behind the globe
    let down = false;
    for (const [lat, lng] of points) {
      const [u, v, front] = project(lat, lng);
      if (front && down) ctx.lineTo(u, v); else if (front) ctx.moveTo(u, v);
      down = front;
    }
  };
  ctx.strokeStyle = c.coast;
  ctx.lineWidth = 0.6;
  ctx.beginPath();
  for (const line of DATA.land) trace(Array.from({length: line.length / 2}, (_, i) => [line[2 * i + 1], line[2 * i]]));
  ctx.stroke();
  if (edge > 0 && edge < 90) {                          // the ice edge, in both hemispheres
    ctx.strokeStyle = c.accent;
    ctx.lineWidth = 2;
    ctx.beginPath();
    for (const lat of [edge, -edge]) trace(Array.from({length: 181}, (_, i) => [lat, 2 * i]));
    ctx.stroke();
  }
  ctx.strokeStyle = c.ink;
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.arc(cx, cy, R, 0, 2 * Math.PI);
  ctx.stroke();
}

function drawStrip(c) {                                 // the player's path: ice line against the sun
  const [ctx, W, H] = begin(views.strip);
  ctx.fillStyle = c.bg;
  ctx.fillRect(0, 0, W, H);
  if (W < 80 || H < 60) return;
  const fs = parseFloat(getComputedStyle(document.body).fontSize);   // the page's text size, 14 px or more
  const m = {l: 2.4 * fs, r: 0.8 * fs, t: 2.3 * fs, b: 2.7 * fs}, w = W - m.l - m.r, h = H - m.t - m.b;
  const X = (q) => m.l + (q - Q_MIN) / (Q_MAX - Q_MIN) * w, Y = (e) => m.t + h - e / 90 * h;
  ctx.font = `${fs}px system-ui, sans-serif`;
  ctx.fillStyle = c.fg;
  ctx.strokeStyle = c.muted;
  ctx.lineWidth = 1;
  ctx.textBaseline = "top";
  for (let q = Math.ceil(Q_MIN / 20) * 20; q <= Q_MAX; q += 20) {   // labels at the ends align inward, clear of the 0
    ctx.textAlign = q === Q_MIN ? "left" : q === Q_MAX ? "right" : "center";
    ctx.fillText(q, X(q), m.t + h + 0.3 * fs);
  }
  ctx.textAlign = "center";
  ctx.fillText("sun / % of today", m.l + w / 2, m.t + h + 1.45 * fs);
  ctx.textAlign = "right";
  ctx.textBaseline = "middle";
  for (const e of [0, 30, 60, 90]) ctx.fillText(`${e}`, m.l - 0.4 * fs, Y(e));
  ctx.textAlign = "left";
  ctx.textBaseline = "top";
  ctx.fillText("ice line / °", 2, 0.4 * fs);
  ctx.globalAlpha = 0.25;
  ctx.beginPath();
  for (const e of [30, 60, 90]) { ctx.moveTo(m.l, Y(e)); ctx.lineTo(m.l + w, Y(e)); }
  ctx.stroke();
  ctx.globalAlpha = 1;
  ctx.beginPath();
  ctx.moveTo(m.l, m.t);
  ctx.lineTo(m.l, m.t + h);
  ctx.lineTo(m.l + w, m.t + h);
  ctx.stroke();
  ctx.fillStyle = c.accent;
  for (const [q, e] of state.trail) {
    ctx.beginPath();
    ctx.arc(X(q), Y(e), 2.5, 0, 2 * Math.PI);
    ctx.fill();
  }
  ctx.strokeStyle = c.accent;                           // where the planet is now
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  ctx.arc(X(state.q), Y(iceLine(state.T)), 5, 0, 2 * Math.PI);
  ctx.stroke();
}

function render() {
  const c = colors(), edge = iceLine(state.T);
  drawGlobe(c, VIEW.lon0 - state.spin, edge);
  drawStrip(c);
  $("mean").textContent = `global mean ${minus(globalMean(state.T).toFixed(1))} °C`;
  $("ice").textContent = iceText(edge);
  $("year").textContent = `year ${Math.floor(state.year)}`;
  $("sun").value = state.q;
  $("sun-val").textContent = `${state.q} %`;
  $("play").textContent = state.running ? "Pause" : "Play";
}

// ---- controls, keys, and the loop
function set(partial) {
  if (partial.T) partial = {...partial, T: Float64Array.from(partial.T)};
  Object.assign(state, partial);
  render();
}
const clamp = (v) => Math.min(Q_MAX, Math.max(Q_MIN, v));
const nudge = (sign) => set({q: clamp(state.q + sign)});
const togglePlay = () => set({running: !state.running});
const reset = () => set({T: DATA.T0, q: 100, year: 0, trail: []});

$("sun").addEventListener("input", (e) => set({q: clamp(Number(e.target.value))}));
$("minus").addEventListener("click", () => nudge(-1));
$("plus").addEventListener("click", () => nudge(1));
$("play").addEventListener("click", togglePlay);
$("reset").addEventListener("click", reset);
document.addEventListener("keydown", (e) => {
  const actions = {" ": togglePlay, ArrowLeft: () => nudge(-1), ArrowRight: () => nudge(1), r: reset};
  const act = actions[e.key.length === 1 ? e.key.toLowerCase() : e.key];
  const native = (e.target instanceof HTMLInputElement && e.key.startsWith("Arrow"))
    || (e.target instanceof HTMLButtonElement && e.key === " ");   // the focused control handles these
  if (act && !native) {
    e.preventDefault();
    act();
  }
});

let last = null, owed = 0;
function frame(now) {                                   // called by the browser once per screen refresh
  const elapsed = last === null ? 0 : Math.min((now - last) / 1000, 0.05);   // capped: a background tab does not jump
  last = now;
  if (state.running) {
    state.spin = (state.spin + 360 * elapsed / VIEW.period_s) % 360;
    owed += SPEED * elapsed / DT;                       // model steps due: 25 years per second of play
    for (; owed >= 1; owed--) {
      step(state.T, state.q, DT);
      const before = state.year;
      state.year += DT;
      if (Math.floor(state.year / 2) > Math.floor(before / 2)) {   // a dot on the strip every 2 model years
        state.trail.push([state.q, iceLine(state.T)]);
        if (state.trail.length > 2000) state.trail.shift();
      }
    }
  }
  render();
  requestAnimationFrame(frame);
}

const resized = new ResizeObserver(layout);             // the boxes change with the frame and with the readouts
for (const id of ["globe-view", "strip-view"]) resized.observe($(id));
window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", render);
layout();
requestAnimationFrame(frame);
window.game = {state, set, render};