/*{
  "DESCRIPTION": "A tank of clear water with ink dropped into it, run as an actual fluid simulation rather than a picture of one. Six passes: velocity is advected through itself, forced by the emitters and by the stirrer, given back the curl that advection steals (vorticity confinement), then made incompressible by a damped Jacobi pressure solve before it is allowed to carry any dye. The plumes, the mushroom caps, the shed vortex pairs and the long thin filaments are all things the solver does -- not noise faded in and out to look like smoke. Turn the pressure solve off and the ink collapses into a flat structureless smear, which is the honest test of whether one is really running. Nothing here waits to be touched: the emitters move on their own and the stirrer is a vortex that roams the tank on a figure of eight, so it keeps churning by itself. STIR sets the middle of that beat and WANDER sets how far it strays, down to nought if you would rather it stayed put. EMIT chooses where the ink comes from: fountains along the floor, rain from the ceiling, orbit, collide, ring, wander, or still if you want nothing in the tank but the stirrer. STYLE chooses how the tank is photographed -- ink in water, glow, schlieren, vorticity, speed, thermal, chrome, oil film, neon, marbling, aurora, x-ray, nebula, wet paper, contour and streamlines -- and the last one, divergence, shows the solver grading itself: red where the flow is compressing, blue where it is expanding, and near black everywhere the pressure solve has done its job.",
  "CREDIT": "CC0",
  "ISFVSN": "2",
  "CATEGORIES": ["fluid", "simulation", "generative", "feedback"],
  "INPUTS": [
    { "NAME": "uStir", "TYPE": "point2D", "DEFAULT": [0.5, 0.3] },
    { "NAME": "uStirF", "TYPE": "float", "DEFAULT": 0.35, "MIN": -3.0, "MAX": 3.0 },
    { "NAME": "uStirR", "TYPE": "float", "DEFAULT": 0.14, "MIN": 0.02, "MAX": 0.6 },
    { "NAME": "uStirInk", "TYPE": "float", "DEFAULT": 0.0, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uWander", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0, "MAX": 1.0 },

    { "NAME": "uSpeed", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.1, "MAX": 3.0 },
    { "NAME": "uForce", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.0, "MAX": 3.0 },
    { "NAME": "uVort", "TYPE": "float", "DEFAULT": 1.1, "MIN": 0.0, "MAX": 4.0 },
    { "NAME": "uDrag", "TYPE": "float", "DEFAULT": 0.12, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uBuoy", "TYPE": "float", "DEFAULT": 0.8, "MIN": -3.0, "MAX": 3.0 },
    { "NAME": "uSwirl", "TYPE": "float", "DEFAULT": 0.0, "MIN": -2.0, "MAX": 2.0 },
    { "NAME": "uCoarse", "TYPE": "float", "DEFAULT": 2.5, "MIN": 1.0, "MAX": 8.0 },
    { "NAME": "uSharp", "TYPE": "bool", "DEFAULT": true },
    { "NAME": "uWalls", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1],
      "LABELS": ["tank", "wrap"] },

    { "NAME": "uInk", "TYPE": "float", "DEFAULT": 1.3, "MIN": 0.0, "MAX": 3.0 },
    { "NAME": "uFade", "TYPE": "float", "DEFAULT": 0.10, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uRadius", "TYPE": "float", "DEFAULT": 0.030, "MIN": 0.004, "MAX": 0.15 },
    { "NAME": "uSources", "TYPE": "long", "DEFAULT": 4, "VALUES": [1, 2, 3, 4, 5, 6, 7, 8],
      "LABELS": ["1", "2", "3", "4", "5", "6", "7", "8"] },
    { "NAME": "uEmit", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1, 2, 3, 4, 5, 6],
      "LABELS": ["fountains", "rain", "orbit", "collide", "ring", "wander", "still"] },
    { "NAME": "uPalette", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1, 2, 3, 4, 5, 6, 7],
      "LABELS": ["spectrum", "inks", "sunset", "ocean", "neon", "mono", "ember", "candy"] },
    { "NAME": "uSat", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "uSeed", "TYPE": "float", "DEFAULT": 7.0, "MIN": 0.0, "MAX": 64.0 },
    { "NAME": "uReset", "TYPE": "event" },

    { "NAME": "uStyle", "TYPE": "long",
      "DEFAULT": 0, "VALUES": [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],
      "LABELS": ["ink in water", "glow", "schlieren", "vorticity", "speed", "thermal",
                 "chrome", "oil film", "neon", "marbling", "aurora", "x-ray",
                 "nebula", "wet paper", "contour", "streamlines", "divergence"] },
    { "NAME": "uGain", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.0, "MAX": 4.0 },
    { "NAME": "uGamma", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.4, "MAX": 2.4 },
    { "NAME": "uEdge", "TYPE": "float", "DEFAULT": 0.8, "MIN": 0.0, "MAX": 3.0 },
    { "NAME": "uGlow", "TYPE": "float", "DEFAULT": 0.5, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "uGrain", "TYPE": "float", "DEFAULT": 0.10, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uVign", "TYPE": "float", "DEFAULT": 0.35, "MIN": 0.0, "MAX": 1.0 }
  ],
  "PASSES": [
    { "TARGET": "velA", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "prsA", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "prsB", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "velB", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "dye",  "PERSISTENT": true, "FLOAT": true },
    { }
  ]
}*/

// =============================================================================
//  ink tank -- Stam stable fluids in six ISF passes
// =============================================================================
//
//  velA  .xy velocity, advected and forced but still compressible
//        .z  1.0 once seeded    .w  forcing phase, 0..1
//  prsA  .x  pressure, sweep 1  .y  divergence  .z  curl  (divergence and curl
//  prsB  .x  pressure, sweep 2  .y  divergence  .z  curl   are carried, not
//                                                          recomputed -- they
//                                                          come free with sweep 1)
//  velB  .xy velocity with the pressure gradient taken out: what moves the ink
//        .z  curl               .w  the phase again, so the render pass reads
//                                   one texture instead of two
//  dye   .rgb ink               .a  1.0 once seeded (a fresh target reads zero,
//                                   which is how the first frame knows it is first)
//
//  Pass order is the whole trick. velA advects velB from LAST frame, which was
//  divergence free; prsA and prsB relax against THIS frame's velA, each sweep reading
//  the one before it, with prsA warm-starting from last frame's prsB; velB and
//  dye both run after the solve, so everything downstream sees a legal
//  incompressible field. Targets are double buffered, so a pass reading its own
//  target at a neighbouring texel gets last frame's value everywhere -- which is
//  what lets the dye advection sample itself.
//
//  Three properties of this renderer shaped the code more than anything else,
//  and all three were measured before a line of the solver was written:
//
//  1. FLOAT pass targets are HALF float, not 32 bit. That is the constraint that
//     decides how the pressure solve has to be written. Relaxation converges by
//     accumulating corrections; the moment a correction is smaller than the
//     storage quantum the sweep stalls and pressure stops existing. The textbook
//     form stores p with a factor of the cell size squared in it, which at this
//     grid is around 3e-5 -- underneath half float's smallest normal number, so
//     it flushes to zero and the solver silently does nothing. Everything below
//     is therefore written in CELL units, where the cell size cancels out of
//     both the sweep and the projection (see passPrs) and every stored quantity
//     is order 1. A coarser lattice also helps directly: the relative size of a
//     Jacobi correction goes as 1/L^2 in cells, so uCoarse is not only a
//     performance dial, it is what keeps the solve above the noise floor.
//
//  2. A pass target is always the size of the window; a PASS WIDTH/HEIGHT is
//     ignored. So the lattice would be four times finer at 1080p than on the 512
//     poster and the same sliders would give a different fluid on every canvas.
//     gSP scales the stencil with the window instead, so the tank is always
//     about 512/uCoarse cells tall wherever it runs.
//
//  3. FRAMEINDEX is round(TIME*60) rather than a count of frames drawn, and TIME
//     does not advance at all during the 120 frame poster warm up. Anything that
//     had to move per frame would freeze. So the sim keeps its own phase in
//     velA.w and every moving thing reads that. It is stored as a wrapped phase
//     rather than a clock in seconds for reason 1 again: a half float seconds
//     counter stops advancing entirely at about 35 seconds, because by then one
//     frame's DT is smaller than the gap between representable numbers.
//
// =============================================================================

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#else
precision mediump float;
#endif

#define DT_BASE (1.0 / 60.0)
#define PERIOD  12.0
// The dye luminance that reads as "fully inked". Not a taste constant: the ink
// spreads as it is carried, and measured over the emitter modes the 95th
// percentile of luminance lands at 0.13 to 0.22 with a peak near 0.5 -- so the
// styles below, whose curves are all written for a 0..1 density, were being fed
// the bottom tenth of their range and came out nearly black. Normalising here
// once is better than seventeen magic numbers.
#define DREF    0.20
#define RELAX   0.5   // damped Jacobi; see passPrs for why this is not 1.0

// ---- per-fragment globals, filled by setup() --------------------------------
vec2  gR;     // render size, pixels
float gAsp;   // width / height
vec2  gTex;   // one pixel, in uv
vec2  gUV;    // this fragment, in uv
float gSP;    // stencil spacing, in pixels (see note 2 above)
vec2  gDX;    // one cell to the right, in uv
vec2  gDY;    // one cell up, in uv
float gDT;    // seconds per frame
float gPh;    // forcing phase, 0..1, one turn per PERIOD seconds

// World space is uv with x stretched by the aspect ratio, so y runs 0..1 and x
// runs 0..gAsp. One world unit is gR.y pixels on BOTH axes, which is what keeps
// every vortex round on a wide canvas.
vec2 worldP(vec2 uv){ return vec2(uv.x * gAsp, uv.y); }

void setup(){
  gR   = RENDERSIZE.xy;
  gAsp = gR.x / gR.y;
  gTex = 1.0 / gR;
  gUV  = gl_FragCoord.xy * gTex;
  gSP  = max(1.0, uCoarse * gR.y / 512.0);
  gDX  = vec2(gTex.x * gSP, 0.0);
  gDY  = vec2(0.0, gTex.y * gSP);
  gDT  = DT_BASE * uSpeed;
}

// ---- small maths ------------------------------------------------------------

float lum(vec3 c){ return dot(c, vec3(0.2126, 0.7152, 0.0722)); }

// A persistent buffer never forgets, so one NaN is permanent: it advects into
// its neighbours and the tank is dead until the page reloads. Every value that
// goes into a target passes through here first. Written as a comparison rather
// than clamp() because clamp(NaN) is undefined and usually just keeps the NaN.
float safe1(float x, float lim){ return (x > -lim && x < lim) ? x : 0.0; }
vec2  safe2(vec2 v, float lim){ return vec2(safe1(v.x, lim), safe1(v.y, lim)); }
vec3  safe3(vec3 v, float lim){ return vec3(safe1(v.x, lim), safe1(v.y, lim), safe1(v.z, lim)); }

// Only ever called with small arguments (source index, cell index), where the
// sin hash is stable across drivers. Anything that decided the geometry of the
// sim would need an exact integer hash instead; here it picks seed blobs, paper
// fibre and film grain.
float hash1(float n){ return fract(sin(n * 12.9898 + uSeed * 7.13) * 43758.5453); }
vec2  hash2(float n){ return vec2(hash1(n), hash1(n + 17.31)); }

float vnoise(vec2 x){
  vec2 i = floor(x), f = fract(x);
  f = f * f * (3.0 - 2.0 * f);
  float a = hash1(dot(i, vec2(1.0, 57.0)));
  float b = hash1(dot(i + vec2(1.0, 0.0), vec2(1.0, 57.0)));
  float c = hash1(dot(i + vec2(0.0, 1.0), vec2(1.0, 57.0)));
  float e = hash1(dot(i + vec2(1.0, 1.0), vec2(1.0, 57.0)));
  return mix(mix(a, b, f.x), mix(c, e, f.x), f.y);
}

// A soft signed squash, standing in for tanh() (which ISF does not allow).
float squash(float x){ return x / (1.0 + abs(x)); }

vec3 cosPal(vec3 a, vec3 b, vec3 c, vec3 d, float t){
  return a + b * cos(6.283185 * (c * t + d));
}

vec3 pal(float t){
  t = fract(t);
  vec3 c;
  if      (uPalette == 0) c = cosPal(vec3(0.50), vec3(0.50), vec3(1.0), vec3(0.00, 0.33, 0.67), t);
  else if (uPalette == 1) c = cosPal(vec3(0.55, 0.50, 0.50), vec3(0.45, 0.50, 0.50), vec3(1.0), vec3(0.10, 0.45, 0.80), t);
  else if (uPalette == 2) c = cosPal(vec3(0.62, 0.42, 0.36), vec3(0.42, 0.36, 0.30), vec3(1.0, 1.0, 0.55), vec3(0.02, 0.15, 0.30), t);
  else if (uPalette == 3) c = cosPal(vec3(0.28, 0.44, 0.55), vec3(0.26, 0.32, 0.38), vec3(1.0), vec3(0.55, 0.62, 0.72), t);
  else if (uPalette == 4) c = cosPal(vec3(0.55, 0.42, 0.60), vec3(0.50, 0.48, 0.52), vec3(2.0, 1.0, 0.0), vec3(0.50, 0.20, 0.25), t);
  else if (uPalette == 5) c = vec3(0.86 + 0.14 * cos(6.283185 * t));
  else if (uPalette == 6) c = cosPal(vec3(0.52, 0.28, 0.14), vec3(0.48, 0.32, 0.18), vec3(1.0), vec3(0.00, 0.12, 0.24), t);
  else                    c = cosPal(vec3(0.70, 0.52, 0.62), vec3(0.30, 0.42, 0.32), vec3(1.0), vec3(0.00, 0.25, 0.50), t);
  c = max(c, 0.0);
  return mix(vec3(lum(c)), c, uSat);
}

// Black body-ish ramp shared by the thermal and speed views.
vec3 ironbow(float t){
  t = clamp(t, 0.0, 1.0);
  vec3 c = vec3(1.9 * t, 1.6 * t * t, 3.4 * t * t * t * t);
  return c / (1.0 + c * 0.55);
}

// ---- sampling ---------------------------------------------------------------
// One place decides what happens off the edge of the tank, so the walls behave
// the same for advection, for the pressure stencil and for the render taps.
// Clamping is not just a safety net: it makes the pressure solve Neumann at the
// wall -- the ghost cell equals the wall cell, so the normal gradient is zero,
// which is exactly the boundary condition a closed tank wants.
vec2 edgeUV(vec2 uv){
  if (uWalls == 1) return fract(uv);
  return clamp(uv, 0.5 * gTex, 1.0 - 0.5 * gTex);
}
vec4 tap(sampler2D s, vec2 uv){ return IMG_NORM_PIXEL(s, edgeUV(uv)); }

// Step backwards along a velocity, in uv. x is divided by the aspect because a
// world unit sideways covers less uv than a world unit vertically.
vec2 trace(vec2 uv, vec2 v, float dt){
  return edgeUV(uv - vec2(v.x / gAsp, v.y) * dt);
}

// Semi-Lagrangian advection, optionally sharpened.
//
// Plain backtracing is unconditionally stable but eats detail: every step is a
// bilinear blend, so a plume turns to fog in a few seconds. The uSharp branch is
// MacCormack -- advect the result forwards again, and half the round trip error
// is the diffusion, so subtract it. Unlimited, that overshoots into new extrema
// and detonates, so the corrected value is clamped into the range of the cells
// it was actually interpolated from. That clamp is the difference between a
// sharp fluid and a screen full of NaN.
vec4 advect(sampler2D src, sampler2D vf, vec2 uv, float dt){
  vec2 v0 = tap(vf, uv).xy;
  vec2 u1 = trace(uv, v0, dt);
  vec4 q1 = tap(src, u1);
  if (!uSharp) return q1;
  vec2 v1 = tap(vf, u1).xy;
  vec4 q2 = tap(src, trace(u1, -v1, dt));
  vec4 q  = q1 + 0.5 * (tap(src, uv) - q2);
  vec4 a = tap(src, u1 + gDX), b = tap(src, u1 - gDX);
  vec4 c = tap(src, u1 + gDY), d = tap(src, u1 - gDY);
  return clamp(q, min(min(a, b), min(min(c, d), q1)),
                  max(max(a, b), max(max(c, d), q1)));
}

// ---- the ink sources --------------------------------------------------------
// Position, direction and colour of source i. Everything moves on harmonics of
// gPh, so the forcing is exactly periodic and the phase can wrap forever without
// the sources jumping. The flow never repeats regardless: the fluid carries its
// own state, and two identical forcings a period apart meet completely different
// water.
void source(int i, out vec2 sp, out vec2 sd, out vec3 sc){
  float n  = float(uSources);
  float fi = float(i);
  float a  = (fi + 0.5) / n;
  float w  = 6.283185 * gPh;
  float j  = 6.283185 * hash1(fi * 3.71 + 0.5);

  sc = pal(a + 0.13 * sin(w + j));

  if (uEmit == 0) {                                   // fountains on the floor
    sp = vec2(a * gAsp + 0.05 * sin(w * 2.0 + j), 0.05);
    sd = normalize(vec2(0.45 * sin(w * 3.0 + j * 1.7), 1.0));
  } else if (uEmit == 1) {                            // rain from the ceiling
    sp = vec2(a * gAsp + 0.05 * sin(w * 2.0 + j), 0.95);
    sd = normalize(vec2(0.45 * sin(w * 3.0 + j * 1.7), -1.0));
  } else if (uEmit == 2) {                            // orbiting, blowing sideways
    float th = 6.283185 * a + w;
    sp = vec2(gAsp, 1.0) * 0.5 + 0.28 * vec2(cos(th), sin(th));
    sd = vec2(-sin(th), cos(th));
  } else if (uEmit == 3) {                            // two banks facing off
    float side = mod(fi, 2.0);
    float row  = (floor(fi * 0.5) + 0.5) / max(1.0, ceil(n * 0.5));
    sp = vec2(mix(0.05, gAsp - 0.05, side), 0.15 + 0.7 * fract(row + 0.12 * sin(w + j)));
    sd = vec2(mix(1.0, -1.0, side), 0.0);
  } else if (uEmit == 4) {                            // a ring firing inward
    float th = 6.283185 * a + 0.25 * sin(w + j);
    sp = vec2(gAsp, 1.0) * 0.5 + 0.40 * vec2(cos(th), sin(th));
    sd = -vec2(cos(th), sin(th));
  } else {                                            // wandering
    sp = vec2(gAsp, 1.0) * (0.5 + 0.34 * vec2(sin(w + j), cos(w * 1.31 + j * 2.1)));
    sd = vec2(cos(w * 0.7 + j), sin(w * 0.9 + j));
  }
}

// Gaussian falloff around a source; r is the world-space radius.
float blob(vec2 p, vec2 c, float r){
  vec2 d = p - c;
  return exp(-dot(d, d) / (r * r));
}

vec2 emitForce(vec2 p){
  vec2 f = vec2(0.0);
  if (uEmit != 6) {
    for (int i = 0; i < 8; i++){
      if (i >= uSources) break;
      vec2 sp, sd; vec3 sc;
      source(i, sp, sd, sc);
      f += sd * blob(p, sp, uRadius);
    }
  }
  return f * uForce * 5.0;
}

// Where the spoon is this frame.
//
// There is no mouse over the canvas here -- point2D is a control widget, not a
// pointer -- so a stirrer that only sat where it was placed would be a static
// lump in a shader that is supposed to run itself. uStir sets the middle of its
// beat and uWander sets how far it roams from there, on a figure of eight.
// Both harmonics are integer multiples of the phase, so the path closes exactly
// when the phase wraps and the spoon never jumps; the water it finds on the
// next lap is different every time, which is what keeps the loop from reading
// as a loop.
vec2 stirP(){
  float w = 6.283185 * gPh;
  vec2 c = worldP(uStir) + uWander * 0.30 * vec2(sin(w), sin(2.0 * w + 1.1));
  return clamp(c, vec2(0.08), vec2(gAsp - 0.08, 0.92));
}

vec3 inkAt(vec2 p){
  vec3 c = vec3(0.0);
  if (uEmit != 6) {
    for (int i = 0; i < 8; i++){
      if (i >= uSources) break;
      vec2 sp, sd; vec3 sc;
      source(i, sp, sd, sc);
      c += sc * blob(p, sp, uRadius);
    }
  }
  if (uStirInk > 0.0) {
    c += pal(gPh * 0.7) * uStirInk * blob(p, stirP(), uStirR * 0.55);
  }
  return c * uInk * 2.4;
}

// The spoon: a rotating vortex, dragging the water round wherever it is.
vec2 stirForce(vec2 p){
  vec2 d = p - stirP();
  float r = length(d);
  if (r < 1e-4) return vec2(0.0);
  return (vec2(-d.y, d.x) / r) * uStirF * blob(p, stirP(), uStirR) * 2.2;
}

// Free-slip walls: kill only the component heading out of the tank, so ink
// slides along the glass instead of sticking to it.
vec2 walls(vec2 uv, vec2 v){
  if (uWalls == 1) return v;
  vec2 b = gTex * gSP * 1.5;
  if (uv.x <       b.x) v.x = max(v.x, 0.0);
  if (uv.x > 1.0 - b.x) v.x = min(v.x, 0.0);
  if (uv.y <       b.y) v.y = max(v.y, 0.0);
  if (uv.y > 1.0 - b.y) v.y = min(v.y, 0.0);
  return v;
}

// ---- initial conditions -----------------------------------------------------
// The tank does not start still. A stream function built from a few sines, then
// differentiated, gives a velocity field that is divergence free by construction
// -- the solver's own invariant, satisfied before the solver has run once. The
// ink therefore starts curling on frame one instead of sitting in a dead pool,
// which also means the poster is a developed fluid and not an empty rectangle.
float potential(vec2 p){
  float s = 0.0, f = 3.0, a = 1.0;
  for (int k = 0; k < 4; k++){
    vec2 o = hash2(float(k) * 5.13) * 6.283185;
    s += a * sin(p.x * f + o.x) * sin(p.y * f * 1.13 + o.y);
    f *= 2.03; a *= 0.5;
  }
  return s;
}

vec2 seedVelocity(vec2 p){
  float e = 0.01;
  float dx = potential(p + vec2(e, 0.0)) - potential(p - vec2(e, 0.0));
  float dy = potential(p + vec2(0.0, e)) - potential(p - vec2(0.0, e));
  return vec2(dy, -dx) * (0.05 / (2.0 * e));
}

vec3 seedDye(vec2 p){
  vec3 c = vec3(0.0);
  for (int i = 0; i < 7; i++){
    float fi = float(i);
    vec2  q  = vec2(hash1(fi * 1.77 + 0.3) * gAsp, 0.1 + 0.8 * hash1(fi * 2.91 + 5.1));
    float r  = mix(0.05, 0.13, hash1(fi * 4.43 + 9.0));
    c += pal(hash1(fi * 6.17 + 2.0)) * blob(p, q, r) * 1.3;
  }
  return c;
}

// =============================================================================
//  pass 0 -- advect and force the velocity
// =============================================================================
void passVel(){
  vec4 prev = IMG_THIS_PIXEL(velA);
  bool first = (prev.z < 0.5) || uReset;

  // Advance the phase FIRST, so every later pass this frame reads the same phase
  // the forces were evaluated at.
  gPh = fract(prev.w + gDT / PERIOD);

  vec2 p = worldP(gUV);
  vec2 v;

  if (first) {
    gPh = 0.0;
    v = seedVelocity(p);
  } else {
    // Plain backtracing here, not MacCormack. The sharpening is nine taps and
    // measurably does nothing for the picture on the VELOCITY field: velocity is
    // smooth, confinement puts back the curl advection loses, and what the eye
    // reads as detail is all in the dye, which does still get the sharp path.
    vec4 vb = tap(velB, gUV);
    v = tap(velB, trace(gUV, vb.xy, gDT)).xy;

    // Vorticity confinement. Advection quietly drains angular momentum out of
    // small eddies; this pushes it back by nudging velocity along the gradient
    // of |curl|, so filaments keep curling instead of dissolving into fog.
    float cL = tap(velB, gUV - gDX).z, cR = tap(velB, gUV + gDX).z;
    float cB = tap(velB, gUV - gDY).z, cT = tap(velB, gUV + gDY).z;
    vec2  gc = vec2(abs(cR) - abs(cL), abs(cT) - abs(cB));
    float gn = length(gc);
    if (gn > 1e-5) {
      // f = eps * h * (N x omega). The curl stored in velB.z is a raw
      // difference, which IS h*omega, so the cell size has already cancelled and
      // what is left is eps * curl * dt. Dropping that dt is not a scale tweak:
      // confinement injects energy in proportion to the curl it finds, so at 27x
      // strength it becomes a feedback loop that doubles the energy every few
      // seconds and tears the tank apart. It looks fine for the first hundred
      // frames, which is exactly long enough to ship it.
      vec2 N = gc / gn;
      v += uVort * vec2(N.y, -N.x) * vb.z * gDT * 6.0;
    }

    v += emitForce(p) * gDT;
    v += stirForce(p) * gDT;
    v.y += uBuoy * lum(tap(dye, gUV).rgb) * gDT * 1.6;

    vec2 q = p - vec2(gAsp, 1.0) * 0.5;
    v += vec2(-q.y, q.x) * uSwirl * gDT * 1.5;

    v *= exp(-uDrag * gDT * 4.0);
  }

  v = walls(gUV, safe2(v, 40.0));
  gl_FragColor = vec4(v, 1.0, gPh);
}

// =============================================================================
//  passes 1..3 -- three Jacobi sweeps of the pressure Poisson equation
// =============================================================================
//  We want v - grad(p) to be divergence free, which means laplacian(p) = div(v).
//  The textbook discretisation is p = (sum of neighbours - div * h^2) / 4 with
//  div measured as a difference over 2h. Both of those h factors are poison in a
//  half float target: they push p down to around 3e-5, below the smallest normal
//  half float, so it flushes to zero and no amount of iterating does anything.
//
//  Writing it in CELL units makes every h cancel, which is what keeps the stored
//  numbers order 1 and every sweep's correction representable:
//
//      P = (PL + PR + PB + PT - dvg) / 4
//
//  The other half of getting this right is WHICH differences dvg and the
//  gradient use, and it is not a detail. The obvious choice is a central
//  difference for both. Compose them and you get
//
//      div(grad(p))[i] = (p[i+2] - 2p[i] + p[i-2]) / 4
//
//  which is the Laplacian at spacing TWO, not the compact five point Laplacian
//  the sweep above actually solves. The gradient being subtracted is therefore
//  not the gradient of the field that was solved for, and the mismatch lives
//  entirely in the checkerboard mode -- which the solver cannot see and so can
//  never remove. Measured on a fixed field, that scheme reduces divergence by
//  1.5x and then STOPS: 300 sweeps and 3000 sweeps give the same answer. It
//  looks like an iteration count problem and no amount of iterating fixes it.
//
//  So dvg is a FORWARD difference and the projection is a BACKWARD one:
//
//      dvg    = (vx(i+1) - vx(i)) + (vy(j+1) - vy(j))
//      v     -= (P(i) - P(i-1), P(j) - P(j-1))
//
//  Now div(grad(p)) is exactly p[i+1] - 2p[i] + p[i-1]: the same operator the
//  sweep inverts, so at convergence the projection removes the divergence
//  completely rather than asymptotically missing most of it. This is the MAC
//  staggering every serious solver uses, and it costs one tap less than the
//  wrong version.
//
//  Two sweeps a frame, warm started from last frame's answer, is not a
//  converged solve -- but the field it is chasing barely moves between frames,
//  so the error that matters is the change, not the total, and the "divergence"
//  style shows what is left.
// =============================================================================
void passPrs(sampler2D src, bool measure){
  vec4 a = IMG_THIS_PIXEL(velA);
  gPh = a.w;

  vec4  s0 = IMG_THIS_PIXEL(src);
  float dvg, crl;
  if (measure) {
    // The two taps this needs are also the two the curl needs, so taking curl
    // here as a forward difference makes it free and saves the projection pass
    // four taps of its own.
    vec2 vR = tap(velA, gUV + gDX).xy;
    vec2 vT = tap(velA, gUV + gDY).xy;
    dvg = (vR.x - a.x) + (vT.y - a.y);
    crl = (vR.y - a.y) - (vT.x - a.x);
  } else {
    dvg = s0.y; crl = s0.z;             // carried, not recomputed
  }

  float pL = tap(src, gUV - gDX).x, pR = tap(src, gUV + gDX).x;
  float pB = tap(src, gUV - gDY).x, pT = tap(src, gUV + gDY).x;
  float pj = (pL + pR + pB + pT - dvg) * 0.25;
  // Damped Jacobi. Taking the full step is not just slower to settle, it is
  // UNSTABLE here: with only a few sweeps a frame the pressure never catches up,
  // and subtracting an unconverged gradient injects more divergence than it
  // removes, which feeds back through the advection and tears the tank apart in
  // about twenty frames. Modelled over the whole feedback loop, half steps leave
  // divergence at 0.086 of the flow's own shear where full steps blow up and the
  // old central scheme sat at 1.54.
  float p = mix(s0.x, pj, RELAX);
  if (uReset) p = 0.0;

  gl_FragColor = vec4(safe1(p, 4000.0), safe1(dvg, 40.0), safe1(crl, 40.0), 1.0);
}

// =============================================================================
//  pass 4 -- subtract the pressure gradient, and store the curl
// =============================================================================
void passProj(){
  vec4 a = IMG_THIS_PIXEL(velA);
  gPh = a.w;

  vec2 v = a.xy;
  // backward difference, the adjoint of the forward one dvg was measured with
  float p0 = IMG_THIS_PIXEL(prsB).x;
  v -= vec2(p0 - tap(prsB, gUV - gDX).x, p0 - tap(prsB, gUV - gDY).x);
  v = walls(gUV, safe2(v, 40.0));

  // Curl comes from velA rather than from the projected field, which is not an
  // approximation: projection subtracts a gradient, and the curl of a gradient
  // is zero, so both fields have the same curl. It was already computed in the
  // first pressure sweep and carried down the chain.
  gl_FragColor = vec4(v, IMG_THIS_PIXEL(prsB).z, a.w);
}

// =============================================================================
//  pass 5 -- carry the ink
// =============================================================================
void passDye(){
  vec4 a = IMG_THIS_PIXEL(velA);
  gPh = a.w;

  vec4 prev = IMG_THIS_PIXEL(dye);
  vec3 d;
  if (prev.a < 0.5 || uReset) {
    d = seedDye(worldP(gUV));
  } else {
    d = advect(dye, velB, gUV, gDT).rgb;
    d *= exp(-uFade * gDT * 3.0);
    d += inkAt(worldP(gUV)) * gDT;
  }
  gl_FragColor = vec4(clamp(safe3(d, 64.0), 0.0, 16.0), 1.0);
}

// =============================================================================
//  pass 6 -- photograph the tank
// =============================================================================
void passRender(){
  vec4 D = IMG_THIS_PIXEL(dye);
  vec4 V = IMG_THIS_PIXEL(velB);
  gPh = V.w;

  float ex = uGain / DREF;      // dye -> 0..1 working density
  vec3  d  = max(D.rgb, 0.0) * ex;
  float L  = lum(d);
  float sp = length(V.xy);
  float cu = V.z;

  // Density gradient: the schlieren image, and the surface normal for anything
  // that wants to look like a lit material. Kept in per-cell units (no divide by
  // the cell size) so the styles look the same at any window size.
  float lR = lum(tap(dye, gUV + gDX).rgb) * ex, lL = lum(tap(dye, gUV - gDX).rgb) * ex;
  float lT = lum(tap(dye, gUV + gDY).rgb) * ex, lB = lum(tap(dye, gUV - gDY).rgb) * ex;
  vec2  g  = vec2(lR - lL, lT - lB) * 0.5;
  vec3  n  = normalize(vec3(-g * uEdge * 12.0, 1.0));

  // A cheap four tap bloom, but only for the four styles that light up.
  vec3 blur = vec3(0.0);
  if (uStyle == 1 || uStyle == 8 || uStyle == 10 || uStyle == 12) {
    blur = (tap(dye, gUV + gDX * 7.0).rgb + tap(dye, gUV - gDX * 7.0).rgb
          + tap(dye, gUV + gDY * 7.0).rgb + tap(dye, gUV - gDY * 7.0).rgb) * 0.25 * ex;
  }

  vec3 col;

  if (uStyle == 0) {
    // Ink in water, backlit. Beer-Lambert through the dye, so thick ink goes
    // black and thin ink only tints -- the reason this reads as liquid and not
    // as paint.
    vec3 back = mix(vec3(0.97, 0.97, 0.95), vec3(0.55, 0.62, 0.72), gUV.y);
    col = back * exp(-d * 2.6);
    col *= 0.85 + 0.30 * dot(n, normalize(vec3(-0.5, 0.7, 0.6)));
  } else if (uStyle == 1) {
    col = d * 0.55 + d * d * 0.10 + blur * uGlow * 0.5;
  } else if (uStyle == 2) {
    // Schlieren: a lab photograph of density gradient, dark where ink shears.
    col = vec3(1.0 - squash(length(g) * 9.0 * uEdge)) * vec3(0.99, 0.99, 1.0);
  } else if (uStyle == 3) {
    float c = squash(cu * 22.0 * uEdge);
    col = mix(vec3(0.04, 0.05, 0.08),
              c > 0.0 ? vec3(1.00, 0.42, 0.20) : vec3(0.22, 0.55, 1.00), abs(c));
    col += d * 0.10;
  } else if (uStyle == 4) {
    col = ironbow(sp * 2.2 * uEdge) * 1.15;
  } else if (uStyle == 5) {
    col = ironbow(L * 0.92);
  } else if (uStyle == 6) {
    // Chrome: the dye is a height field, the tank a studio with two lights.
    vec3 r = reflect(vec3(0.0, 0.0, -1.0), n);
    vec3 env = mix(vec3(0.05, 0.06, 0.09), vec3(0.55, 0.62, 0.72), r.y * 0.5 + 0.5);
    env += vec3(1.00, 0.96, 0.90) * pow(max(dot(r, normalize(vec3(-0.55, 0.72, 0.42))), 0.0), 48.0);
    env += vec3(0.35, 0.45, 0.65) * pow(max(dot(r, normalize(vec3(0.62, -0.35, 0.70))), 0.0), 16.0);
    col = env * (0.35 + 0.65 * smoothstep(0.0, 0.35, L));
  } else if (uStyle == 7) {
    // Oil film: thickness beats against wavelength, so the colour cycles with
    // the depth of the ink rather than being painted on it.
    float th = L * 5.0 + 0.35;
    col = (0.5 + 0.5 * cos(6.283185 * (th + vec3(0.00, 0.33, 0.67))))
        * smoothstep(0.0, 0.10, L) * (0.55 + 0.45 * n.z);
    col += vec3(0.02, 0.03, 0.05);
  } else if (uStyle == 8) {
    float e = squash(length(g) * 7.0 * uEdge);
    col = pal(fract(L * 0.45 + gPh)) * e * 1.6 + d * 0.06 + blur * uGlow * 0.35;
  } else if (uStyle == 9) {
    // Suminagashi: the density field cut into bands, so the ink reads as the
    // combed rings of paper marbling.
    float band = fract(L * 5.0);
    float line = smoothstep(0.0, 0.09, band) * smoothstep(0.0, 0.09, 1.0 - band);
    col = mix(vec3(0.05, 0.05, 0.07), vec3(0.95, 0.93, 0.87), line);
    col *= 0.75 + 0.25 * n.z;
  } else if (uStyle == 10) {
    col = pal(fract(0.35 + L * 0.22 + gUV.y * 0.20)) * pow(L, 0.75) * 1.1;
    col += blur * uGlow * 0.45;
    col += vec3(0.01, 0.02, 0.04);
  } else if (uStyle == 11) {
    float t = 1.0 - exp(-L * 2.2);
    col = mix(vec3(0.02, 0.03, 0.05), vec3(0.80, 0.90, 1.00), t);
    col += vec3(0.35, 0.55, 0.85) * squash(length(g) * 4.0 * uEdge);
  } else if (uStyle == 12) {
    col = d * 0.75 + blur * uGlow * 0.6;
    float st = vnoise(gUV * gR.y * 0.30);
    col += vec3(0.9, 0.95, 1.0) * smoothstep(0.985, 1.0, st) * (1.0 - smoothstep(0.0, 0.25, L));
    col += vec3(0.015, 0.012, 0.030);
  } else if (uStyle == 13) {
    // Ink on wet paper: fibre in the sheet, and a darker rim where the wash
    // dried at its own edge.
    float fib = vnoise(gUV * vec2(gR.y * 0.8, gR.y * 0.22)) * 0.10
              + vnoise(gUV * vec2(gR.y * 0.20, gR.y * 0.9)) * 0.06;
    vec3 paper = vec3(0.94, 0.91, 0.84) - fib * 0.35;
    float rim = squash(length(g) * 6.0 * uEdge);
    col = paper * exp(-(L * 3.2 + rim * 0.9) * vec3(1.00, 1.06, 1.16));
  } else if (uStyle == 14) {
    // Contour map of the ink.
    float f  = L * 7.0;
    float ln = abs(fract(f) - 0.5) * 2.0;
    // fwidth() is banned in ISF, but g is already the change in L per cell, so
    // the screen-space width of a contour is that gradient rescaled. Same
    // result, no derivative instruction.
    float w  = clamp(28.0 * length(g) / gSP, 0.004, 0.45);
    col = pal(fract(floor(f) * 0.11 + 0.2)) * (0.28 + 0.5 * smoothstep(0.0, 1.0, L));
    col = mix(vec3(0.02, 0.02, 0.03), col, smoothstep(0.0, 0.02, L));
    col = mix(col, vec3(0.98), 1.0 - smoothstep(w, w * 2.4, ln));
  } else if (uStyle == 15) {
    // Streamlines: walk backwards along the flow, averaging a fixed noise field
    // as you go. Points on the same streamline average the same samples, so the
    // noise smears into brush strokes that ARE the velocity field.
    // Two things have to be true or this reads as television static rather than
    // brush strokes: the noise has to be FINER than the walk is long, and the
    // step has to follow the direction of the flow rather than its speed --
    // otherwise slow water gets a dot and fast water gets a streak, and the
    // texture reports speed twice instead of showing direction once.
    vec2  q = gUV;
    float acc = 0.0, wsum = 0.0;
    float stp = 3.0 / gR.y;                 // three pixels per step, either way
    for (int k = 0; k < 14; k++){
      float side = (k < 7) ? -1.0 : 1.0;
      if (k == 7) q = gUV;
      vec2 vv = tap(velB, q).xy;
      float sl = length(vv);
      vv = (sl > 1e-4) ? vv / sl : vec2(0.0, 1.0);
      q = edgeUV(q + side * vec2(vv.x / gAsp, vv.y) * stp);
      float w = 1.0 - abs(float(k) - 6.5) / 8.0;
      acc  += w * vnoise(q * gR.y * 0.42);
      wsum += w;
    }
    float lic = clamp((acc / wsum - 0.5) * (3.0 + 3.0 * uEdge) + 0.5, 0.0, 1.0);
    vec3 tint = mix(vec3(0.55, 0.60, 0.70), d * 1.4 + 0.15, smoothstep(0.0, 0.35, L));
    col = tint * (0.25 + 0.9 * lic);
  } else {
    // The solver, exposed. This is the residual the pressure sweeps have not
    // cleaned up yet: red where the flow is compressing, blue where it is
    // expanding. It should be near black except in a thin shell around whatever
    // is being forced hardest.
    float dv = squash(tap(prsB, gUV).y * 18.0 * uEdge);
    col = mix(vec3(0.02), dv > 0.0 ? vec3(1.0, 0.25, 0.15) : vec3(0.15, 0.45, 1.0), abs(dv));
  }

  // ---- finish -------------------------------------------------------------
  col = max(col, 0.0);
  if (uGamma != 1.0) col = pow(col, vec3(1.0 / uGamma));

  float grain = hash1(dot(floor(gl_FragCoord.xy), vec2(1.0, 311.7))
                    + floor(gPh * PERIOD * 60.0) * 0.017);
  col += (grain - 0.5) * uGrain * 0.13;

  vec2 vg = (gUV - 0.5) * vec2(gAsp, 1.0);
  col *= 1.0 - uVign * smoothstep(0.20, 0.95, dot(vg, vg) * 1.6);

  gl_FragColor = vec4(clamp(col, 0.0, 1.0), 1.0);
}

void main(){
  setup();
  if      (PASSINDEX == 0) passVel();
  else if (PASSINDEX == 1) passPrs(prsB, true);   // warm start, measure divergence
  else if (PASSINDEX == 2) passPrs(prsA, false);
  else if (PASSINDEX == 3) passProj();
  else if (PASSINDEX == 4) passDye();
  else                     passRender();
}