Triple Pendulum Chaos

let st, trail, cx, cy;
const L = [90, 90, 90];
const m = [1, 1, 1];
const g = 9.81;

function init({ width, height }) {
  cx = width / 2;
  cy = height * 0.35;
  randomize();
  trail = [];
}

function randomize() {
  st = [
    Math.PI / 2 + (Math.random() - 0.5) * 2,
    Math.PI / 2 + (Math.random() - 0.5) * 2,
    Math.PI / 2 + (Math.random() - 0.5) * 2,
    0, 0, 0,
  ];
}

function deriv(s) {
  const [t1, t2, t3, w1, w2, w3] = s;
  const [L1, L2, L3] = L;
  const [m1, m2, m3] = m;
  const c12 = Math.cos(t1 - t2), s12 = Math.sin(t1 - t2);
  const c13 = Math.cos(t1 - t3), s13 = Math.sin(t1 - t3);
  const c23 = Math.cos(t2 - t3), s23 = Math.sin(t2 - t3);
  const M = [
    [(m1 + m2 + m3) * L1, (m2 + m3) * L2 * c12, m3 * L3 * c13],
    [(m2 + m3) * L1 * c12, (m2 + m3) * L2, m3 * L3 * c23],
    [m3 * L1 * c13, m3 * L2 * c23, m3 * L3],
  ];
  const b = [
    -(m2 + m3) * L2 * w2 * w2 * s12 - m3 * L3 * w3 * w3 * s13 - (m1 + m2 + m3) * g * Math.sin(t1),
    (m2 + m3) * L1 * w1 * w1 * s12 - m3 * L3 * w3 * w3 * s23 - (m2 + m3) * g * Math.sin(t2),
    m3 * L1 * w1 * w1 * s13 + m3 * L2 * w2 * w2 * s23 - m3 * g * Math.sin(t3),
  ];
  const det3 = (A) =>
    A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) -
    A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) +
    A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
  const D = det3(M);
  const Mx = [[b[0], M[0][1], M[0][2]], [b[1], M[1][1], M[1][2]], [b[2], M[2][1], M[2][2]]];
  const My = [[M[0][0], b[0], M[0][2]], [M[1][0], b[1], M[1][2]], [M[2][0], b[2], M[2][2]]];
  const Mz = [[M[0][0], M[0][1], b[0]], [M[1][0], M[1][1], b[1]], [M[2][0], M[2][1], b[2]]];
  return [w1, w2, w3, det3(Mx) / D, det3(My) / D, det3(Mz) / D];
}

function vadd(a, b, s) { return a.map((v, i) => v + b[i] * s); }

function rk4(s, h) {
  const k1 = deriv(s);
  const k2 = deriv(vadd(s, k1, h / 2));
  const k3 = deriv(vadd(s, k2, h / 2));
  const k4 = deriv(vadd(s, k3, h));
  return s.map((v, i) => v + (h * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i])) / 6);
}

function tick({ ctx, dt, time, width, height, input }) {
  if (input.consumeClicks().length > 0) { randomize(); trail = []; }
  cx = width / 2;
  cy = height * 0.35;

  const steps = 16;
  const h = Math.min(dt, 0.033) / steps;
  for (let i = 0; i < steps; i++) st = rk4(st, h);

  const [t1, t2, t3] = st;
  const x1 = cx + L[0] * Math.sin(t1), y1 = cy + L[0] * Math.cos(t1);
  const x2 = x1 + L[1] * Math.sin(t2), y2 = y1 + L[1] * Math.cos(t2);
  const x3 = x2 + L[2] * Math.sin(t3), y3 = y2 + L[2] * Math.cos(t3);

  trail.push({ x: x3, y: y3, t: time });
  if (trail.length > 400) trail.shift();

  ctx.fillStyle = "rgba(10,10,20,0.18)";
  ctx.fillRect(0, 0, width, height);

  ctx.lineWidth = 2;
  for (let i = 1; i < trail.length; i++) {
    const p = trail[i - 1], q = trail[i];
    const age = (trail.length - i) / trail.length;
    const hue = (q.t * 60 + i * 2) % 360;
    ctx.strokeStyle = `hsla(${hue},100%,60%,${1 - age})`;
    ctx.beginPath();
    ctx.moveTo(p.x, p.y);
    ctx.lineTo(q.x, q.y);
    ctx.stroke();
  }

  ctx.strokeStyle = "#ddd";
  ctx.lineWidth = 2;
  ctx.beginPath();
  ctx.moveTo(cx, cy);
  ctx.lineTo(x1, y1);
  ctx.lineTo(x2, y2);
  ctx.lineTo(x3, y3);
  ctx.stroke();

  ctx.fillStyle = "#888";
  ctx.beginPath();
  ctx.arc(cx, cy, 4, 0, Math.PI * 2);
  ctx.fill();

  const pts = [[x1, y1, "#5af"], [x2, y2, "#f5a"], [x3, y3, "#fa5"]];
  for (const [x, y, c] of pts) {
    ctx.fillStyle = c;
    ctx.beginPath();
    ctx.arc(x, y, 9, 0, Math.PI * 2);
    ctx.fill();
    ctx.strokeStyle = "#fff";
    ctx.lineWidth = 1.5;
    ctx.stroke();
  }

  ctx.fillStyle = "#aaa";
  ctx.font = "12px monospace";
  ctx.fillText("click to randomize", 10, height - 12);
}