/*{
  "DESCRIPTION": "Serpent Field: three raymarched snakes flocking in real 3D. A persistent boid solver flies the heads with a bounded turn rate and a lateral swim stroke, a following chain drags the bodies along the exact path each head travelled, and the render pass marches a tapered swept-sphere distance field with parallel-transported scales, mutual coloured bounce light, and volumetric bloom.",
  "CATEGORIES": ["generative", "3d", "raymarch"],
  "INPUTS": [
    { "NAME": "speed", "TYPE": "float", "DEFAULT": 0.90, "MIN": 0.0, "MAX": 1.8 },
    { "NAME": "girth", "TYPE": "float", "DEFAULT": 0.062, "MIN": 0.028, "MAX": 0.110 },
    { "NAME": "bodyLength", "TYPE": "float", "DEFAULT": 1.55, "MIN": 0.60, "MAX": 2.60 },
    { "NAME": "flock", "TYPE": "float", "DEFAULT": 0.72, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "chase", "TYPE": "float", "DEFAULT": 0.46, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "avoid", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "undulate", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 1.4 },
    { "NAME": "swim", "TYPE": "float", "DEFAULT": 1.20, "MIN": 0.2, "MAX": 3.0 },
    { "NAME": "mess", "TYPE": "float", "DEFAULT": 0.50, "MIN": 0.0, "MAX": 1.2 },
    { "NAME": "scaleDepth", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "gloss", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "glow", "TYPE": "float", "DEFAULT": 0.60, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "seed", "TYPE": "float", "DEFAULT": 11.0, "MIN": 0.0, "MAX": 30.0 }
  ],
  "PASSES": [
    { "TARGET": "heads", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "spine", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "bounds", "FLOAT": true },
    {}
  ]
}*/

// GLSL ES 1.00 gives fragment stages no default float precision and makes highp
// optional there, so a conforming host declares `precision mediump float;` -
// about ten bits of mantissa. This shader integrates head positions across
// thousands of frames and then raymarches against them, where mediump would
// both quantize the flight path and band the surface normals. Ask for highp
// where the hardware advertises it; the guard is required because naming highp
// on a part that lacks it is a compile error.
#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI 3.14159265359
#define TAU 6.28318530718

#define SNAKES 3
#define NODES 33          // spine samples per snake, node 0 is the snout
#define SEGS 32           // NODES - 1 capsules
#define GROUPS 4          // bounding spheres per snake
#define GROUPSPAN 8       // capsules covered by one bounding sphere (GROUPS*GROUPSPAN == SEGS)

#define SCALE_RINGS 15.0  // scales around the circumference; integer keeps the wrap seamless
#define SCALE_PITCH 46.0  // scale rows per unit of arc length

// How far past the snout the tongue can reach. Declared up here because the
// bounds pass has to inflate the head's group sphere by it, and that pass is
// written before the render pass that draws the tongue.
#define TONGUE_REACH 0.17

// How far ahead the head starts steering around a body. It has to exceed the
// turn radius - cruise divided by the turn rate, about 0.37 here - or the head
// gets its warning too late to act on it and only the contact solver saves it.
#define AVOID_REACH 0.20

// Clear air left between two surfaces by the contact solver. Small enough to
// read as touching, large enough that neither the march's hit threshold nor
// half-float position quantisation can close it.
#define CONTACT_MARGIN 0.014

// Nodes this close along the same body are neighbours, not obstacles.
#define SELF_SKIP 4.5

// ---------------------------------------------------------------------------
// Small shared helpers
// ---------------------------------------------------------------------------

float hash11(float p){ return fract(sin(p*127.1 + seed*19.73)*43758.5453); }

// Bounding a 3D quantity has to bound its magnitude, not its components: clamp()
// on a vec3 bounds it to a cube, so anything saturating gets pushed into the
// corners and its direction snaps onto the box diagonals. Limiting the length
// leaves direction untouched.
vec3 limitLength(vec3 v, float limit){
  float m = length(v);
  return m > limit ? v*(limit/max(m, 1e-5)) : v;
}

vec3 safeNormalize(vec3 v, vec3 fallback){
  float m = length(v);
  return m > 1e-5 ? v/m : fallback;
}

// A divergence-free wander field, written as the exact curl of an analytic sine
// potential rather than as finite differences of noise. Being divergence free
// is what stops the heads from piling into the same attracting sheet: the field
// stirs without ever having a sink for them to fall into.
vec3 curlField(vec3 p, float t){
  vec3 a = vec3(p.y*1.7 + t*0.61, p.z*1.3 - t*0.44, p.x*1.9 + t*0.53);
  vec3 b = vec3(p.z*2.6 - t*0.37, p.x*2.9 + t*0.71, p.y*2.3 - t*0.66);
  vec3 f = vec3(-1.3*cos(a.y), -1.9*cos(a.z), -1.7*cos(a.x));
  vec3 g = vec3(-2.9*cos(b.y), -2.3*cos(b.z), -2.6*cos(b.x));
  return f + g*0.38;
}

// Body radius along the snake. u = 0 at the snout, 1 at the tip of the tail.
//
// The shape of this curve is most of what separates "snake" from "worm" at a
// glance. A snake holds a near-constant diameter for the front half and only
// then tapers, so the taper is deferred rather than applied over the whole
// body; a plain pow() falloff from the head gives the carrot silhouette that
// the first render produced. The tail still has to reach an actual point, or
// it reads as a cut tube.
// Two failure modes this curve has to avoid, both found by looking at renders:
//
//   - If the neck thins faster than the snout bump decays, the profile dips
//     between node 0 and node 2 and the head detaches into a lollipop ball. The
//     two terms have to be tuned against each other, not independently.
//   - If the taper reaches zero exactly at u = 1, the last capsule is a needle
//     a whole segment long and the tail renders as a stray hair. Ending the
//     taper past the end of the body leaves a small rounded tip instead.
//   - A head built as a single bump at u = 0 is one segment long and wider than
//     it is long, so it renders as a blunt tube end with a cap. A head needs to
//     span about three segments and have three distinct parts: a narrow snout
//     tip, a broad skull behind the eyes, and a slight neck.
float radiusAt(float u){
  u = clamp(u, 0.0, 1.0);
  float taper = pow(1.0 - smoothstep(0.38, 1.02, u), 0.52);
  float neck  = 0.86 + 0.14*smoothstep(0.0, 0.12, u);
  float skull = 0.62*exp(-((u - 0.028)*(u - 0.028))/0.00090);
  float snoutTaper = 0.42 + 0.58*smoothstep(0.0, 0.045, u);
  return girth*snoutTaper*(taper*neck + skull*0.45);
}

float segmentLength(){ return bodyLength/float(SEGS); }

// Cruise speed is a constant rather than an input because the body solver needs
// to know exactly how far the head advanced this frame in order to slide itself
// forward by the same arc.
#define CRUISE 0.88

vec3 catmull(vec3 p0, vec3 p1, vec3 p2, vec3 p3, float t){
  vec3 a = 2.0*p1;
  vec3 b = p2 - p0;
  vec3 c = 2.0*p0 - 5.0*p1 + 4.0*p2 - p3;
  vec3 d = -p0 + 3.0*p1 - 3.0*p2 + p3;
  return 0.5*(a + b*t + c*t*t + d*t*t*t);
}

// These are linear-space albedos and they have to be chosen against the output
// transform, not by eye in isolation. The final gamma raises small channel
// values hard - a linear 0.115 leaves as 0.38 - so an albedo that looks like a
// deep red as numbers arrives on screen as salmon. Everything read pale pink,
// mint and periwinkle until the off-channels came down by roughly five times.
vec3 snakeColor(int s){
  if(s == 0) return vec3(1.00, 0.022, 0.035);
  if(s == 1) return vec3(0.030, 0.95, 0.105);
  return vec3(0.032, 0.105, 1.00);
}

// ---------------------------------------------------------------------------
// State access. Every target is RENDERSIZE-sized because the renderer ignores
// per-pass WIDTH/HEIGHT, so the simulation lives in the first few texels and
// every other fragment early-outs.
//   heads : column = snake. row 0 = position + init flag, row 1 = velocity,
//           row 2 = roll reference (the snake's own "up").
//   spine : column = node, row = snake. xyz = position, w = init flag.
//   bounds: row 0, column = snake*GROUPS + group. xyz = centre, w = radius.
// ---------------------------------------------------------------------------

vec4 headCell(int s, int row){
  return IMG_NORM_PIXEL(heads, (vec2(float(s), float(row)) + 0.5)/RENDERSIZE);
}

vec4 spineCell(int j, int s){
  return IMG_NORM_PIXEL(spine, (vec2(float(j), float(s)) + 0.5)/RENDERSIZE);
}

vec3 spineAtF(float col, int s){
  return IMG_NORM_PIXEL(spine, (vec2(col, float(s)) + 0.5)/RENDERSIZE).xyz;
}

// The body as a curve that can be sampled at a fractional node index. Catmull-
// Rom rather than a straight lerp: the body is resampled from itself every
// single frame, so any interpolation error is applied dozens of times to the
// same piece of shape as it travels to the tail. A lerp's chord error alone
// flattened the swim wave by roughly forty percent on the way down.
vec3 pathAt(float fi, int s){
  float m = float(NODES - 1);
  fi = clamp(fi, 0.0, m);
  float f0 = floor(fi);
  float fr = fi - f0;
  return catmull(spineAtF(max(f0 - 1.0, 0.0), s),
                 spineAtF(f0, s),
                 spineAtF(min(f0 + 1.0, m), s),
                 spineAtF(min(f0 + 2.0, m), s), fr);
}

vec4 boundCell(int s, int k){
  return IMG_NORM_PIXEL(bounds, (vec2(float(s*GROUPS + k), 0.0) + 0.5)/RENDERSIZE);
}

// ---------------------------------------------------------------------------
// Cold-start pose
//
// glslop derives FRAMEINDEX from the requested poster time, so a still rendered
// at t=8 hands this shader a first-ever frame numbered near 480 and every
// "is this frame zero" test fails. State is therefore seeded from the buffer's
// own contents (the w flag), and the seed pose is a function of TIME so that a
// single cold frame at any moment is already a composed picture of three snakes
// mid-flight rather than three straight rods.
// ---------------------------------------------------------------------------

void seedPose(int s, float t, out vec3 pos, out vec3 dir){
  float fs = float(s);
  float ph = fs*(TAU/float(SNAKES)) + hash11(fs + 3.0)*1.7;
  float orbit = t*0.21 + ph;
  float rad = 1.12 + 0.20*sin(t*0.23 + ph*2.3);
  pos = vec3(cos(orbit)*rad,
             0.46*sin(orbit*0.77 + ph) + 0.22*sin(t*0.17 + ph*1.9),
             sin(orbit)*rad);
  dir = normalize(vec3(-sin(orbit)*rad, 0.30*cos(t*0.31 + ph), cos(orbit)*rad));
}

// A seeded body is laid back along the reverse of the heading with a swim wave
// already in it, so a cold frame shows serpents rather than sticks.
vec3 seedNode(int s, int j, float t){
  vec3 pos, dir;
  seedPose(s, t, pos, dir);
  vec3 side = safeNormalize(cross(dir, vec3(0.0, 1.0, 0.0)), vec3(1.0, 0.0, 0.0));
  vec3 up = cross(side, dir);
  float L = segmentLength();
  float fj = float(j);
  float a = fj*0.62 - t*2.1 + float(s)*1.3;
  return pos - dir*(L*fj) + side*(sin(a)*L*0.62) + up*(cos(a*0.63)*L*0.26);
}

// ===========================================================================
// PASS 0 - the boids
// ===========================================================================

// Camera state, filtered and stored rather than derived fresh in the render
// pass. Two reasons it belongs here:
//
//   - `extent` is a max over the whole flock, and a max switches which point it
//     is following. That is continuous in value but not in slope, so the framing
//     picked up a small kink every time the outermost node changed - invisible
//     in a still, and a persistent low-level unsteadiness in motion.
//   - Deriving it per pixel cost the render pass two dozen texture reads per
//     fragment to compute a number that is identical for the whole frame.
//
// It runs off last frame's body, since this pass precedes the body pass. A frame
// of camera lag is not observable; the filter's own time constants are hundreds
// of times longer.
void mainCamera(){
  float dt = clamp(TIMEDELTA, 1.0/240.0, 1.0/30.0)*speed;

  vec3 centre = vec3(0.0);
  float count = 0.0;
  for(int o = 0; o < SNAKES; o++){
    for(int n = 0; n < NODES; n++){
      vec4 c = spineCell(n, o);
      if(c.w < 0.5) continue;
      centre += c.xyz;
      count += 1.0;
    }
  }
  if(count < 0.5){ gl_FragColor = vec4(0.0, 0.0, 0.0, 4.2); return; }
  centre /= count;

  float extent = 0.0;
  for(int o = 0; o < SNAKES; o++){
    for(int n = 0; n < NODES; n++){
      vec4 c = spineCell(n, o);
      if(c.w < 0.5) continue;
      extent = max(extent, length(c.xyz - centre) + radiusAt(float(n)/float(SEGS)));
    }
  }

  // At focal f the half-height visible at distance d is d/(2f), so fitting a
  // group of radius `extent` needs d >= 2*f*extent - about 3.1 here.
  float rawDolly = clamp(extent*3.20, 2.9, 7.6);
  vec3 rawTarget = centre*0.92;

  vec4 prev = headCell(0, 3);
  // A live dolly is never below its own 2.9 floor, so w is its own cold flag.
  if(prev.w < 0.5){ gl_FragColor = vec4(rawTarget, rawDolly); return; }

  float kTarget = 1.0 - exp(-dt/0.40);
  float kDolly = 1.0 - exp(-dt/0.90);
  gl_FragColor = vec4(mix(prev.xyz, rawTarget, kTarget), mix(prev.w, rawDolly, kDolly));
}

void mainHeads(){
  float col = floor(gl_FragCoord.x);
  float row = floor(gl_FragCoord.y);
  if(row > 2.5){
    if(row > 3.5 || col > 0.5){ gl_FragColor = vec4(0.0); return; }
    mainCamera();
    return;
  }
  if(col > float(SNAKES) - 0.5){ gl_FragColor = vec4(0.0); return; }

  int s = int(col);
  float t = TIME*speed;
  float dt = clamp(TIMEDELTA, 1.0/240.0, 1.0/30.0)*speed;

  vec4 posCell = headCell(s, 0);
  vec4 velCell = headCell(s, 1);
  vec4 rollCell = headCell(s, 2);

  vec3 pos = posCell.xyz;
  vec3 vel = velCell.xyz;
  vec3 roll = rollCell.xyz;

  // The swim phase is integrated with the simulation's own clamped dt and kept
  // as a wrapped fraction of a cycle, rather than being read straight off TIME.
  //
  // Driving it from TIME couples the stroke to wall-clock time while the head
  // position is integrated with a dt that is deliberately clamped for
  // stability. The two agree at normal frame rates and diverge badly anywhere
  // else - glslop's poster schedule steps at t/120, which for an eight second
  // poster is 0.067 per step against a clamp of 0.033. The head then covers
  // half the distance per stroke cycle that it should, the wavelength halves,
  // the amplitude halves with it, and the animals arrive at the canonical
  // renderer visibly straighter than they are in the studio. Accumulating the
  // phase alongside the position keeps the stroke locked to the motion it is
  // supposed to describe, at any frame rate and under any host's schedule.
  //
  // It is wrapped rather than accumulated without bound because these targets
  // are RGBA16F: an unwrapped phase would reach a magnitude where the mantissa
  // can no longer resolve one frame's increment, and the stroke would quietly
  // freeze after a few minutes of playback.
  float phase = fract(velCell.w);

  // Contents decide whether this is a cold buffer - never the frame number.
  if(posCell.w < 0.5){
    vec3 sp, sd;
    seedPose(s, t, sp, sd);
    pos = sp;
    vel = sd*0.9;
    roll = safeNormalize(cross(cross(sd, vec3(0.0, 1.0, 0.0)), sd), vec3(0.0, 1.0, 0.0));
    phase = hash11(float(s)*7.31 + 2.0);
  }

  phase = fract(phase + swim*dt);

  vec3 vdir = safeNormalize(vel, vec3(0.0, 0.0, 1.0));

  // Roll is parallel transported through time instead of rebuilt per frame from
  // a world axis. A frame rebuilt from a fixed reference flips whenever the
  // heading passes that reference, and the flip would snap every scale on the
  // body a quarter turn in one frame. Transport has no such singularity.
  roll = safeNormalize(roll - vdir*dot(roll, vdir), vec3(0.0, 1.0, 0.0));

  // Pure transport is free to leave a snake flying belly-up for long stretches,
  // which reads as an error even though it is physically consistent. A weak
  // righting pull toward world up fixes it. The pull is faded out as the heading
  // approaches vertical, where "up" stops being a meaningful reference - which
  // is also the only place the pull could snap rather than drift.
  vec3 worldUp = vec3(0.0, 1.0, 0.0);
  float upright = 1.0 - abs(dot(vdir, worldUp));
  vec3 rightTo = safeNormalize(worldUp - vdir*dot(worldUp, vdir), roll);
  roll = safeNormalize(mix(roll, rightTo, clamp(upright*2.6*dt, 0.0, 0.10)), roll);
  vec3 side = safeNormalize(cross(vdir, roll), vec3(1.0, 0.0, 0.0));

  vec3 force = vec3(0.0);

  // Orbit. A tangential push around a slowly tipping axis is what turns three
  // free boids into a shared circulation instead of three unrelated wanderers.
  vec3 axis = normalize(vec3(0.22*sin(t*0.11 + float(s)), 1.0, 0.18*cos(t*0.09)));
  force += cross(axis, pos)*1.15;

  // Containment, plus a weak outward push near the middle so they never all
  // collapse through the origin at once.
  float r = length(pos);
  force -= pos*(smoothstep(0.70, 1.30, r)*4.8);
  force += safeNormalize(pos, vec3(1.0, 0.0, 0.0))*((1.0 - smoothstep(0.20, 0.62, r))*1.5);

  // Flocking against the other heads.
  for(int o = 0; o < SNAKES; o++){
    if(o == s) continue;
    vec4 oPos = headCell(o, 0);
    vec4 oVel = headCell(o, 1);
    if(oPos.w < 0.5) continue;
    vec3 d = pos - oPos.xyz;
    float dl = max(length(d), 1e-4);
    vec3 dn = d/dl;
    // The equilibrium spacing is where these two balance. Pushed out to a full
    // body length the three snakes simply orbit as separate objects; bringing it
    // in to roughly half that is what lets them cross, coil past one another and
    // light each other up, which is the whole point of flying them as a flock.
    force += dn*(flock*1.85/(0.20 + dl*dl));                       // separation
    force -= dn*(flock*1.05*smoothstep(0.35, 1.20, dl));           // cohesion
    force += oVel.xyz*(flock*0.30);                                // alignment
  }

  // Cyclic pursuit: red chases green, green chases blue, blue chases red, each
  // aiming at a point trailing the one ahead rather than at its head.
  //
  // Symmetric cohesion alone gives three animals loosely orbiting a common
  // centre, which reads as three unrelated swimmers that happen to share a
  // frame. Making the attraction a cycle gives the group an intent: the classic
  // three-body pursuit spirals into a rotating braid, and because each snake is
  // aiming behind its target it falls into the wake rather than colliding with
  // it. Separation still has the final say when they close.
  int chased = s + 1;
  if(chased >= SNAKES) chased = 0;
  vec4 chasedPos = headCell(chased, 0);
  vec4 chasedVel = headCell(chased, 1);
  if(chasedPos.w > 0.5){
    vec3 lead = safeNormalize(chasedVel.xyz, vec3(0.0, 0.0, 1.0));
    vec3 aim = (chasedPos.xyz - lead*0.55) - pos;
    force += safeNormalize(aim, vec3(0.0))*(chase*2.30);
  }

  // Body avoidance - the thing head-to-head separation cannot do.
  //
  // A body here is a passive trail: only the head steers, and every node behind
  // it just slides along the path the head already laid down. So two bodies will
  // drift straight through one another no matter how well the heads avoid each
  // other, because nothing in the body is able to move out of the way. The head
  // has to steer around whole *bodies*, its own included, and it has to start
  // early enough that a rate-limited turn can still clear them.
  //
  // The falloff is quadratic and reaches exactly zero at the reach. A linear
  // cutoff would put a kink in the steering force every time a node crossed the
  // boundary, and with ninety-nine nodes streaming past that reads as a
  // permanent tremor rather than as avoidance.
  for(int o = 0; o < SNAKES; o++){
    for(int n = 0; n < NODES; n++){
      // A snake's own neck is legitimately right behind its own head.
      if(o == s && float(n) < 7.5) continue;
      vec4 cell = spineCell(n, o);
      if(cell.w < 0.5) continue;
      vec3 d = pos - cell.xyz;
      float dl = length(d);
      float reach = AVOID_REACH + radiusAt(float(n)/float(SEGS));
      if(dl > reach || dl < 1e-4) continue;
      float k = 1.0 - dl/reach;
      force += (d/dl)*(k*k*avoid*6.0);
    }
  }

  // Wander.
  force += curlField(pos*0.85, t)*(mess*1.35);

  // Bounded turn rate on the *travel* heading only. Steering by a capped angle
  // rather than by a capped force is what keeps the flight readable: the head
  // can always commit to a heading, it just cannot snap to one.
  vec3 desired = safeNormalize(force, vdir);
  float cosang = clamp(dot(vdir, desired), -1.0, 1.0);
  float ang = acos(cosang);
  float maxTurn = (2.35 + flock*0.6)*dt;
  float k = ang > 1e-4 ? min(maxTurn/ang, 1.0) : 1.0;
  vec3 travel = safeNormalize(mix(vdir, desired, k), vdir);

  // The swim stroke, applied after the steer instead of as another force.
  //
  // Steering the head from side to side and letting the body inherit the S is
  // how a real snake works, and it means the undulation needs no per-node frame
  // - the wave is simply the shape of the path the head left behind. But fed in
  // as a force it does not survive: a stroke tight enough to fit two waves in a
  // body needs a turn radius of a few centimetres, which is an angular rate far
  // above anything a flock should be allowed to steer at. Feeding it in as a
  // force therefore either gets clipped away by the turn cap, or forces the cap
  // so high that the flocking itself becomes twitchy. Rotating the already
  // steered heading decouples the two: the flock decides where the snake is
  // going, the stroke decides how it gets there.
  // The lateral weight is a tangent, so it sets the maximum angle the heading
  // ever departs from travel. That angle has to be large - around fifty degrees
  // at the default - because the path amplitude it produces is
  // (angle * wavelength / 2pi), and anything below about twice the body radius
  // is simply swallowed by the thickness of the animal. At the first attempt's
  // 24 degrees the wave amplitude came out at 0.043 against a radius of 0.062
  // and all three snakes rendered as straight rods.
  float stroke = phase*TAU + float(s)*2.09;
  vec3 strokeSide = safeNormalize(cross(travel, roll), side);
  vec3 strokeUp = cross(strokeSide, travel);
  vec3 heading = safeNormalize(travel + (strokeSide*sin(stroke)
                                        + strokeUp*(cos(stroke)*0.30))*(undulate*2.00), travel);

  // Constant cruise speed. A varying speed makes the body bunch and stretch; a
  // constant one lets it hold its length and just flow.
  pos += heading*(CRUISE*dt);

  // What gets stored is `travel`, the smooth flight direction - NOT the stroked
  // heading the snake actually moved along.
  //
  // Storing the stroked heading makes the stroke an increment on top of itself:
  // next frame's steer starts from an already-deviated direction and adds
  // another full deviation to it. Because the stroke phase moves only about
  // eleven degrees per frame, consecutive frames add near-identical rotations,
  // which compound into a corkscrew instead of cancelling into a weave. It came
  // out as a helix of radius 0.022 against a body radius of 0.062 - the wave was
  // real, and entirely hidden inside the animal. Keeping travel separate makes
  // the stroke an absolute offset from a smooth path, which is what it is.
  vel = travel*CRUISE;

  // Contact resolution for the head itself.
  //
  // The body solver pushes every node out of every surface, but node 0 is not
  // its to move - it is written straight from this pass - so the head was the
  // one part of the animal with steering and no contact response. It is also the
  // leading edge, the part actually driving into things, which is why it
  // accounted for essentially all of the residual interpenetration once the
  // bodies were solved. Correcting it here keeps `heads` and `spine` agreeing.
  float headRadius = radiusAt(0.0);
  vec3 headPush = vec3(0.0);
  for(int o = 0; o < SNAKES; o++){
    if(spineCell(0, o).w < 0.5) continue;
    vec3 a = spineCell(0, o).xyz;
    for(int n = 0; n < SEGS; n++){
      vec3 b = spineCell(n + 1, o).xyz;
      if(o == s && float(n) < SELF_SKIP + 3.5){ a = b; continue; }
      vec3 ba = b - a;
      float h = clamp(dot(pos - a, ba)/max(dot(ba, ba), 1e-8), 0.0, 1.0);
      vec3 d = pos - (a + ba*h);
      float dl = length(d);
      float need = headRadius + CONTACT_MARGIN
                 + mix(radiusAt(float(n)/float(SEGS)), radiusAt(float(n + 1)/float(SEGS)), h);
      if(dl < need && dl > 1e-5) headPush += (d/dl)*(need - dl);
      a = b;
    }
  }
  pos += headPush*0.42;

  // Hard backstop so a bad frame time can never fling a head out of the scene.
  if(length(pos) > 1.95) pos = normalize(pos)*1.95;

  if(row < 0.5)      gl_FragColor = vec4(pos, 1.0);
  else if(row < 1.5) gl_FragColor = vec4(vel, phase);
  else               gl_FragColor = vec4(roll, 1.0);
}

// ===========================================================================
// PASS 1 - the following body
//
// The body is not a chain and not a spring. Node j is defined as the point
// exactly j*segmentLength of true arc back along the path the head has actually
// travelled, found by walking last frame's curve and accumulating real
// distance. Everything else was tried first and none of it survives contact
// with a long body:
//
//   - Follow-the-leader (each node pulled to a fixed distance behind the next)
//     is a pursuit curve, and a pursuit curve low-passes the path it chases.
//     Compounded over thirty-two links it erased a thirteen-segment swim wave
//     entirely and rendered all three snakes as rigid rods.
//   - Resampling the old curve at a fixed index offset convects the shape
//     correctly but has no control over spacing at all: the correction term
//     comes out identical for every node, so it slides the whole body without
//     ever changing the gaps, and the bodies crumpled into clusters of spheres.
//
// Measuring real arc length has no such feedback loop to get wrong. The spacing
// is exact every frame by construction, and it stays exact through frame-rate
// changes, speed changes and live edits to the body length.
// ===========================================================================

// Walk back from this frame's head along last frame's curve until `target` of
// arc length has been covered. The first segment - new head to old head - is
// the piece of path that was created this frame, and it is straight, so a lerp
// is exact there. Further back the spline is used instead of the chord, because
// this resampling happens to the same piece of shape on every frame of its trip
// to the tail and chord error would quietly eat the wave.
vec3 bodyPoint(int s, float target){
  vec3 newHead = headCell(s, 0).xyz;
  vec3 prev = newHead;
  float acc = 0.0;
  for(int i = 0; i < NODES; i++){
    vec3 p = spineAtF(float(i), s);
    float seg = length(p - prev);
    if(acc + seg >= target){
      float f = (target - acc)/max(seg, 1e-6);
      if(i == 0) return mix(newHead, p, f);
      return pathAt(float(i - 1) + f, s);
    }
    acc += seg;
    prev = p;
  }
  return prev;
}

void mainSpine(){
  float col = floor(gl_FragCoord.x);
  float row = floor(gl_FragCoord.y);
  if(col > float(NODES) - 0.5 || row > float(SNAKES) - 0.5){ gl_FragColor = vec4(0.0); return; }

  int j = int(col);
  int s = int(row);
  float t = TIME*speed;
  float L = segmentLength();

  // Node 0 is the head itself, taken from this frame's boid pass.
  if(j == 0){
    vec4 h = headCell(s, 0);
    vec3 p = h.w < 0.5 ? seedNode(s, 0, t) : h.xyz;
    gl_FragColor = vec4(p, 1.0);
    return;
  }

  vec4 selfCell = spineCell(j, s);
  if(selfCell.w < 0.5){ gl_FragColor = vec4(seedNode(s, j, t), 1.0); return; }

  vec3 self = bodyPoint(s, float(j)*L);

  // Contact resolution. Steering gets the heads out of each other's way, but two
  // bodies that are already overlapping have no steering of their own, and a
  // tail cannot decide to move. So after the arc-length walk each node is pushed
  // out of any surface it is inside - other snakes' and, further along, its own.
  //
  // This is safe to do *after* the walk precisely because the walk re-measures
  // true arc length from the stored curve every frame: the displacement becomes
  // part of the recorded path, next frame's walk re-parameterises over it, and
  // node spacing stays exact. The deformation then convects toward the tail on
  // its own, which is what a real snake sliding over another one looks like.
  // Tested against the other body's *surface*, not against its nodes. Node-to-
  // node is not enough: two segments can cross between their endpoints with
  // every node pair still comfortably apart, and that is exactly the case that
  // renders as one snake passing through another. Measuring to the swept capsule
  // closes the gap and also gives the correct push direction - away from the
  // nearest point on the surface rather than away from an arbitrary sample.
  float ur = float(j)/float(SEGS);
  float selfRadius = radiusAt(ur);
  vec3 push = vec3(0.0);
  for(int o = 0; o < SNAKES; o++){
    if(spineCell(0, o).w < 0.5) continue;
    vec3 a = spineCell(0, o).xyz;
    for(int n = 0; n < SEGS; n++){
      vec3 b = spineCell(n + 1, o).xyz;
      // Segments along this same body near this node are neighbours, not
      // obstacles. Compare as floats: ES 1.00 has no integer overload of abs().
      if(o == s && abs(float(n) + 0.5 - float(j)) < SELF_SKIP){ a = b; continue; }
      vec3 ba = b - a;
      float h = clamp(dot(self - a, ba)/max(dot(ba, ba), 1e-8), 0.0, 1.0);
      vec3 d = self - (a + ba*h);
      float dl = length(d);
      float need = selfRadius + CONTACT_MARGIN
                 + mix(radiusAt(float(n)/float(SEGS)), radiusAt(float(n + 1)/float(SEGS)), h);
      if(dl < need && dl > 1e-5) push += (d/dl)*(need - dl);
      a = b;
    }
  }

  // Relaxed rather than applied whole: the node this one is being pushed off is
  // being pushed off it in the same frame, so resolving the full overlap on both
  // sides at once overshoots and the two buzz against each other.
  self += push*0.50;

  gl_FragColor = vec4(self, 1.0);
}

// ===========================================================================
// PASS 2 - bounding spheres
//
// A single sphere around a whole snake is a useless bound: the body is long and
// thin, so the sphere is mostly empty and the march pays full price the moment
// it enters. Four spheres per snake keep the empty-space stepping honest and
// cost twelve texels a frame.
// ===========================================================================

void mainBounds(){
  float col = floor(gl_FragCoord.x);
  float row = floor(gl_FragCoord.y);
  if(col > float(SNAKES*GROUPS) - 0.5 || row > 0.5){ gl_FragColor = vec4(0.0); return; }

  float fs = floor(col/float(GROUPS));
  float fk = col - fs*float(GROUPS);
  int s = int(fs);
  int k = int(fk);

  vec3 centre = vec3(0.0);
  for(int i = 0; i <= GROUPSPAN; i++){
    centre += spineCell(k*GROUPSPAN + i, s).xyz;
  }
  centre /= float(GROUPSPAN + 1);

  float rad = 0.0;
  for(int i = 0; i <= GROUPSPAN; i++){
    int node = k*GROUPSPAN + i;
    float u = float(node)/float(SEGS);
    rad = max(rad, length(spineCell(node, s).xyz - centre) + radiusAt(u));
  }

  // The head's group has to cover the tongue as well, or the coarse gate in the
  // march will skip past it.
  if(k == 0) rad += TONGUE_REACH;

  gl_FragColor = vec4(centre, rad + 0.015);
}

// ===========================================================================
// PASS 3 - the render
// ===========================================================================

// Distance to one tapered capsule. h comes back so the caller can recover where
// along the segment the closest point sat.
float capsuleDistance(vec3 p, vec3 a, vec3 b, float ra, float rb, out float h){
  vec3 pa = p - a;
  vec3 ba = b - a;
  h = clamp(dot(pa, ba)/max(dot(ba, ba), 1e-8), 0.0, 1.0);
  return length(pa - ba*h) - mix(ra, rb, h);
}

// Tongue flick. Timed from TIME rather than from the integrated phase because
// it is pure decoration - nothing downstream depends on it staying in step with
// the body's motion, and a poster is entitled to catch it at any point.
float tongueExtend(int s){
  float fp = fract(TIME*0.31 + hash11(float(s)*3.71 + 5.0));
  return smoothstep(0.0, 0.13, fp)*smoothstep(0.34, 0.18, fp);
}

// Exact distance to one snake: the union of capsules that share endpoints is
// the sphere swept along the polyline, so the surface stays smooth across
// joints instead of creasing.
float snakeDistance(vec3 p, int s, out float isTongue){
  float best = 1e9;
  isTongue = 0.0;
  vec3 head = spineCell(0, s).xyz;
  vec3 a = head;
  vec3 neck = a;
  for(int j = 0; j < SEGS; j++){
    vec3 b = spineCell(j + 1, s).xyz;
    if(j == 0) neck = b;
    float h;
    float d = capsuleDistance(p, a, b,
                              radiusAt(float(j)/float(SEGS)),
                              radiusAt(float(j + 1)/float(SEGS)), h);
    best = min(best, d);
    a = b;
  }

  // Gate on a sphere around the snout so the tongue costs nothing for the vast
  // majority of march steps, which are nowhere near a head.
  float ext = tongueExtend(s);
  if(ext > 0.01 && length(p - head) < TONGUE_REACH + 0.02){
    vec3 dir = safeNormalize(head - neck, vec3(0.0, 0.0, 1.0));
    vec3 side = safeNormalize(cross(dir, vec3(0.0, 1.0, 0.0)), vec3(1.0, 0.0, 0.0));
    float len = TONGUE_REACH*ext;
    vec3 root = head + dir*(radiusAt(0.0)*0.55);
    vec3 split = root + dir*(len*0.58);
    float h;
    float dt0 = capsuleDistance(p, root, split, 0.0125, 0.0075, h);
    vec3 tipA = split + dir*(len*0.42) + side*(len*0.26);
    vec3 tipB = split + dir*(len*0.42) - side*(len*0.26);
    float dt1 = min(capsuleDistance(p, split, tipA, 0.0072, 0.0026, h),
                    capsuleDistance(p, split, tipB, 0.0072, 0.0026, h));
    float dtongue = min(dt0, dt1);
    if(dtongue < best){ best = dtongue; isTongue = 1.0; }
  }
  return best;
}

// Coarse-to-fine scene distance. The group spheres give a valid lower bound, so
// returning one of them when it is still far is conservative - the march simply
// takes a shorter step than it could have - while skipping twenty capsules.
#define FINE_GATE 0.075

float mapScene(vec3 p, out float id){
  float best = 1e9;
  id = -1.0;
  for(int s = 0; s < SNAKES; s++){
    float sb = 1e9;
    for(int k = 0; k < GROUPS; k++){
      vec4 b = boundCell(s, k);
      sb = min(sb, length(p - b.xyz) - b.w);
    }
    if(sb > best) continue;
    if(sb > FINE_GATE){ best = sb; continue; }
    float isTongue;
    float d = snakeDistance(p, s, isTongue);
    // The fractional part carries "this is the tongue, not the body", so the
    // shading pass can give it flesh instead of scales without a second march.
    if(d < best){ best = d; id = float(s) + isTongue*0.5; }
  }
  return best;
}

float mapDistance(vec3 p){
  float id;
  return mapScene(p, id);
}

// Surface frame at the hit. Walked once, not per march step: the frame vector is
// parallel transported down the chain so the scale lattice cannot spin around
// the body wherever the heading happens to align with a world axis.
void surfaceFrame(vec3 p, int s, out float arc, out vec3 tangent, out vec3 axisU, out vec3 axisV){
  float bestD = 1e9;
  arc = 0.0;
  tangent = vec3(0.0, 0.0, 1.0);
  axisU = vec3(0.0, 1.0, 0.0);

  vec3 a = spineCell(0, s).xyz;
  vec3 b0 = spineCell(1, s).xyz;
  vec3 dir = safeNormalize(b0 - a, vec3(0.0, 0.0, 1.0));
  vec3 frame = headCell(s, 2).xyz;
  frame = safeNormalize(frame - dir*dot(frame, dir), vec3(0.0, 1.0, 0.0));

  float L = segmentLength();
  for(int j = 0; j < SEGS; j++){
    vec3 b = spineCell(j + 1, s).xyz;
    dir = safeNormalize(b - a, dir);
    frame = safeNormalize(frame - dir*dot(frame, dir), frame);
    float h;
    float d = capsuleDistance(p, a, b,
                              radiusAt(float(j)/float(SEGS)),
                              radiusAt(float(j + 1)/float(SEGS)), h);
    if(d < bestD){
      bestD = d;
      arc = (float(j) + h)*L;
      tangent = dir;
      axisU = frame;
    }
    a = b;
  }
  axisV = cross(tangent, axisU);
}

vec3 sceneNormal(vec3 p, float scale){
  vec2 e = vec2(1.0, -1.0)*scale;
  return normalize(e.xyy*mapDistance(p + e.xyy) +
                   e.yyx*mapDistance(p + e.yyx) +
                   e.yxy*mapDistance(p + e.yxy) +
                   e.xxx*mapDistance(p + e.xxx));
}

float ambientOcclusion(vec3 p, vec3 n){
  float occ = 0.0;
  float sca = 1.0;
  for(int i = 0; i < 5; i++){
    float h = 0.012 + 0.055*float(i);
    occ += (h - mapDistance(p + n*h))*sca;
    sca *= 0.72;
  }
  return clamp(1.0 - 2.4*occ, 0.0, 1.0);
}

float softShadow(vec3 ro, vec3 rd){
  float res = 1.0;
  float t = 0.035;
  for(int i = 0; i < 22; i++){
    float h = mapDistance(ro + rd*t);
    res = min(res, 9.0*h/t);
    t += clamp(h, 0.020, 0.20);
    if(res < 0.005 || t > 4.0) break;
  }
  return clamp(res, 0.0, 1.0);
}

// Span of the ray over the union of every group sphere. Pixels that miss all of
// them never enter the march at all, which is most of the frame.
vec2 sceneSpan(vec3 ro, vec3 rd){
  float tn = 1e9;
  float tf = -1e9;
  for(int s = 0; s < SNAKES; s++){
    for(int k = 0; k < GROUPS; k++){
      vec4 b = boundCell(s, k);
      vec3 oc = ro - b.xyz;
      float hb = dot(oc, rd);
      float c = dot(oc, oc) - b.w*b.w;
      float disc = hb*hb - c;
      if(disc < 0.0) continue;
      float sq = sqrt(disc);
      float t0 = -hb - sq;
      float t1 = -hb + sq;
      if(t1 < 0.0) continue;
      tn = min(tn, max(t0, 0.0));
      tf = max(tf, t1);
    }
  }
  return vec2(tn, tf);
}

// Volumetric bloom. Sampling the spine twice per segment and accumulating an
// inverse-square falloff against the ray is far cheaper than a second march and
// is what makes the bodies read as coloured light sources rather than plastic.
// Clamping the projection to the surface hit is the occlusion: anything behind
// the snake in front projects to the hit and lands far from the ray.
#define BLOOM_REACH 0.85

vec3 bloom(vec3 ro, vec3 rd, float tMax){
  vec3 sum = vec3(0.0);
  for(int s = 0; s < SNAKES; s++){
    float w = 0.0;
    for(int k = 0; k < GROUPS; k++){
      // Gate each group of eight segments on its bounding sphere. Without this
      // every pixel in the frame - overwhelmingly background - pays for all
      // ninety-six spine samples, and the glow ends up costing more than the
      // raymarch it decorates.
      vec4 bnd = boundCell(s, k);
      vec3 rel0 = bnd.xyz - ro;
      float tc = clamp(dot(rel0, rd), 0.0, tMax);
      if(length(rel0 - rd*tc) > bnd.w + BLOOM_REACH) continue;

      for(int i = 0; i < GROUPSPAN; i++){
        int j = k*GROUPSPAN + i;
        vec3 a = spineCell(j, s).xyz;
        vec3 b = spineCell(j + 1, s).xyz;
        for(int q = 0; q < 2; q++){
          vec3 pt = mix(a, b, 0.25 + 0.5*float(q));
          vec3 rel = pt - ro;
          float t = clamp(dot(rel, rd), 0.0, tMax);
          float d = length(rel - rd*t);
          w += radiusAt(float(j)/float(SEGS))/(0.005 + 6.5*d*d);
        }
      }
    }
    sum += snakeColor(s)*w;
  }
  return sum*0.0055;
}

// The scale lattice, in (angle around the body, arc length along it). The ring
// count is an integer so the pattern closes on itself at the atan seam, and the
// half-step row offset gives the staggered overlap real scales have.
float scaleField(float ang01, float arc, out vec2 grad, out float curve){
  vec2 q = vec2(ang01*SCALE_RINGS, arc*SCALE_PITCH);
  q.x += 0.5*floor(q.y);
  vec2 f = q - floor(q) - 0.5;
  vec2 w = vec2(f.x, f.y*1.30);
  float d = length(w);
  grad = w/max(d, 1e-4);
  curve = smoothstep(0.0, 0.46, d);
  return smoothstep(0.50, 0.20, d);
}

void mainImage(){
  vec2 uv = (gl_FragCoord.xy - 0.5*RENDERSIZE)/min(RENDERSIZE.x, RENDERSIZE.y);
  float t = TIME*speed;

  // Camera. A slow drift plus a slight rise keeps the flock readable from
  // changing angles without ever cutting away.
  // The camera frames the flock rather than the world origin: three free
  // swimmers spread and gather by a factor of two or more, so a fixed orbit
  // alternates between one snake filling the frame and all three as specks.
  // Both the look-at and the distance are filtered in the heads pass.
  vec4 cam = headCell(0, 3);
  float dolly = cam.w > 0.5 ? cam.w : 4.2;
  float ca = t*0.085 + 0.7;
  float cy = 0.62 + 0.42*sin(t*0.061);
  vec3 ta = cam.xyz + vec3(0.0, 0.05*sin(t*0.09), 0.0);
  vec3 ro = ta + vec3(sin(ca)*dolly, cy, cos(ca)*dolly);
  vec3 fw = normalize(ta - ro);
  vec3 rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
  vec3 up = cross(rt, fw);
  vec3 rd = normalize(fw*1.55 + rt*uv.x + up*uv.y);

  vec3 background = mix(vec3(0.012, 0.014, 0.026), vec3(0.030, 0.026, 0.052),
                        smoothstep(-0.6, 0.7, rd.y));
  // A dead-flat void reads as a missing background rather than as depth. Two
  // crossed low-frequency waves are enough to give it structure at a level the
  // eye reads but does not resolve.
  float haze = sin(rd.x*3.1 + rd.y*2.3)*sin(rd.z*2.7 - rd.y*1.9);
  background += vec3(0.012, 0.013, 0.024)*(0.5 + 0.5*haze);
  background *= 1.0 - 0.35*smoothstep(0.35, 1.05, length(uv));

  vec3 colour = background;
  float tHit = -1.0;
  float hitId = -1.0;

  vec2 span = sceneSpan(ro, rd);
  if(span.x < span.y){
    float tt = span.x;
    float tEnd = min(span.y, 12.0);
    for(int i = 0; i < 96; i++){
      vec3 p = ro + rd*tt;
      float id;
      float d = mapScene(p, id);
      // Widening the hit threshold with distance keeps the silhouettes from
      // shimmering without softening the near geometry.
      if(d < 0.0009 + 0.0006*tt){ tHit = tt; hitId = id; break; }
      tt += d;
      if(tt > tEnd) break;
    }
  }

  if(tHit > 0.0 && hitId >= 0.0){
    // The whole part is the snake, the half carries the tongue flag.
    int s = int(floor(hitId + 0.001));
    float tongueHit = fract(hitId + 0.001) > 0.25 ? 1.0 : 0.0;
    vec3 p = ro + rd*tHit;
    vec3 n = sceneNormal(p, 0.0016 + 0.0009*tHit);
    vec3 base = snakeColor(s);

    float arc;
    vec3 tangent, axisU, axisV;
    surfaceFrame(p, s, arc, tangent, axisU, axisV);

    // Radial direction in the transported frame gives both the scale angle and
    // which side of the animal we are looking at. On the snout and the tail tip
    // the normal turns to face along the body, the radial component collapses,
    // and the angle it yields is noise - which is exactly where the eye and the
    // belly were smearing into a dark blot. Keep the length as a confidence and
    // fade every angular feature out with it.
    vec3 radRaw = n - tangent*dot(n, tangent);
    float radial = length(radRaw);
    vec3 rad = radial > 1e-4 ? radRaw/radial : axisU;
    float frameConf = smoothstep(0.18, 0.62, radial);
    float cu = dot(rad, axisU);
    float cv = dot(rad, axisV);
    float ang01 = atan(cv, cu)/TAU + 0.5;

    vec2 grad;
    float scaleCurve;
    float scaleMask = scaleField(ang01, arc, grad, scaleCurve);

    // Two places the body lattice must not be drawn. On a cap the angle sweeps a
    // full turn across a few pixels and the rings converge into a starburst; and
    // on the head the lattice is simply wrong - a snake's skull is covered in
    // large smooth scutes, not in the fine rows that run down its body.
    // Detail level-of-detail. A scale is about 1/SCALE_PITCH across in world
    // units; once that projects to fewer than a couple of pixels the lattice is
    // past Nyquist and every frame samples it at a slightly different phase,
    // which crawls. Fading it out by projected size costs nothing at the sizes
    // where the detail is actually legible and removes the shimmer on distant
    // animals and in small previews entirely.
    float scalePixels = (min(RENDERSIZE.x, RENDERSIZE.y)*1.55)/(SCALE_PITCH*max(tHit, 0.001));
    float detailLod = smoothstep(1.4, 4.0, scalePixels);
    float scaleWeight = frameConf*smoothstep(0.050, 0.150, arc)*detailLod;
    scaleMask = mix(0.55, scaleMask, scaleWeight);
    scaleCurve = scaleCurve*scaleWeight;

    // Dome each scale, then crease the join. Perturbing only at the cell
    // boundary - which is the obvious thing to do - engraves a lattice of
    // grooves and the animal reads as wireframe rather than as overlapping
    // plates. A real scale is a shallow dome, so the normal has to tilt
    // progressively across the whole cell, proportional to the distance from
    // its centre, with the crease as a much smaller term on top.
    float rim = 4.0*scaleMask*(1.0 - scaleMask);
    vec3 bitan = normalize(cross(n, tangent));
    vec3 bump = bitan*grad.x + tangent*grad.y;
    n = normalize(n - bump*((scaleCurve*0.40 + rim*0.20)*scaleDepth));

    float u = clamp(arc/max(bodyLength, 1e-3), 0.0, 1.0);

    // Belly. cu is signed against the transported "up", so the underside picks
    // up pale plates and the dorsal keeps the saturated colour. The band has to
    // stay narrow: at a wide threshold the cream wraps around far enough to be
    // the dominant tone from most angles and the whole animal goes chalky,
    // which is what the first render did to all three of them.
    float belly = smoothstep(-0.42, -0.95, cu)*frameConf;
    vec3 dorsal = base*(0.55 + 0.45*scaleMask);
    vec3 vent = mix(vec3(0.52, 0.46, 0.34), base, 0.58)*(0.70 + 0.30*scaleMask);
    vec3 albedo = mix(dorsal, vent, belly*0.90);

    // Dorsal saddles. Two incommensurate frequencies rather than one, so the
    // blotches never settle into an obvious repeat, and gated to the back so the
    // flanks stay clean and the belly stays pale - which is how almost every
    // real snake is marked, and is most of what separates a patterned animal
    // from a painted tube.
    float pat = sin(arc*13.0 + float(s)*2.1) + 0.62*sin(arc*29.0 - float(s)*1.3);
    float saddle = smoothstep(0.10, 1.05, pat);
    float dorsalness = smoothstep(-0.25, 0.55, cu)*frameConf;
    albedo *= 1.0 - 0.50*saddle*dorsalness*(1.0 - belly*0.85);
    // A pale margin where each saddle meets the ground colour.
    float edge = 4.0*saddle*(1.0 - saddle);
    albedo += base*edge*dorsalness*0.16;
    albedo *= 0.90 + 0.17*hash11(floor(arc*SCALE_PITCH) + float(s)*37.0);
    // Darken the join between plates. The normal bump alone only shows where a
    // light happens to rake across it; a little occlusion in the groove keeps
    // the scales legible on the shadow side too.
    albedo *= 1.0 - 0.30*rim*scaleDepth;

    // Head: a darker cap over the skull and a wet dark eye. The eye sits on the
    // side of the head (cv near +/-1) and slightly above the midline, so it is
    // placed against the transported frame rather than against the world.
    float headMask = 1.0 - smoothstep(0.012, 0.115, arc);
    albedo *= 1.0 - 0.11*headMask*smoothstep(-0.20, 0.55, cu);
    float onHead = smoothstep(-0.15, 0.30, cu)*frameConf;
    float eye = smoothstep(1.0, 0.40, length(vec2((arc - 0.042)/0.017,
                                                  (abs(cv) - 0.80)/0.30)))*onHead;
    albedo = mix(albedo, vec3(0.015, 0.012, 0.015), eye*0.95);
    // A catchlight, placed rather than lit. Without one the dark ellipse reads
    // as a hole in the head instead of as a wet eye.
    float glint = smoothstep(1.0, 0.20, length(vec2((arc - 0.037)/0.0060,
                                                    (abs(cv) - 0.865)/0.115)))*onHead;

    // Tongue: wet dark flesh, no scales, no markings. It borrows the body's
    // surface frame, which is meaningless out there, so every angular feature
    // has to be switched off rather than merely recoloured.
    albedo = mix(albedo, vec3(0.30, 0.018, 0.055), tongueHit);
    glint *= 1.0 - tongueHit;
    scaleMask = mix(scaleMask, 0.85, tongueHit);

    vec3 key = normalize(vec3(0.55, 0.78, 0.32));
    float ndl = max(dot(n, key), 0.0);
    // Surfaces facing away from the key are already unlit, and the only two
    // terms `shadow` scales - diffuse and the Blinn specular - are both zero
    // there anyway, so the shadow march is pure waste on roughly half the
    // visible surface. Skipping it is exact, not an approximation.
    float shadow = ndl > 0.0015 ? softShadow(p + n*0.012, key) : 0.0;
    float occ = ambientOcclusion(p, n);

    vec3 lit = albedo*(ndl*shadow*vec3(1.00, 0.96, 0.90)*1.30);

    // Hemisphere ambient rather than a flat constant, so an unlit belly still
    // shows which way it is facing.
    float upness = n.y*0.5 + 0.5;
    lit += albedo*occ*mix(vec3(0.030, 0.034, 0.052), vec3(0.105, 0.125, 0.180), upness)*1.5;

    // Cool back rim. With a single key and a black void behind them the three
    // bodies lose their silhouettes wherever they overlap; a dim light from
    // behind and to the side is what separates them from the dark and from each
    // other without lifting the whole frame.
    vec3 backLight = normalize(vec3(-0.45, -0.30, -0.85));
    lit += albedo*max(dot(n, backLight), 0.0)*vec3(0.24, 0.44, 0.88)*0.55;

    // Mutual coloured bounce. Every snake is a light for the others, which is
    // the whole reason to make them red, green and blue: the interaction is
    // visible on the surfaces, not just in the palette.
    //
    // This has to be an average over a snake's four group lights, not a sum.
    // Summed, two neighbours contribute eight lights at up to 0.85 each, which
    // is a six-fold overbright - and since the neighbours of any red snake are
    // green and blue, the three primaries reconstruct white. Every body in the
    // frame turned pale pink. The bounce is meant to tint, not to illuminate.
    for(int o = 0; o < SNAKES; o++){
      if(o == s) continue;
      for(int k = 0; k < GROUPS; k++){
        vec4 b = boundCell(o, k);
        vec3 ld = b.xyz - p;
        float dist = length(ld);
        ld /= max(dist, 1e-4);
        float fall = 1.0/(0.20 + dist*dist*3.2);
        lit += albedo*snakeColor(o)*max(dot(n, ld), 0.0)*fall*(0.45/float(GROUPS));
      }
    }

    // Specular, sharpened by the scale mask so the plates read as wet. Kept
    // deliberately weak: every scale is a dome, so each one presents a range of
    // normals and some part of nearly every scale catches the highlight. A
    // strength that looks right on a smooth surface turns the whole animal into
    // white speckle once the domes are in - which is what actually bleached the
    // bodies, not the coloured bounce.
    vec3 hv = normalize(key - rd);
    float shine = mix(28.0, 150.0, gloss);
    float spec = pow(max(dot(n, hv), 0.0), shine)*(0.35 + 0.65*scaleMask);
    lit += vec3(1.0, 0.97, 0.92)*spec*shadow*gloss*0.42;

    // Fresnel rim in the animal's own colour, so silhouettes stay identifiable
    // where they overlap, with a cool sheen mixed in at the very edge - keratin
    // scales go silver-blue at grazing angles rather than simply brighter.
    // Grazing angles get a thin-film hue shift on top of the cool sheen. Kept
    // to a third of the mix so the animal's identity stays red, green or blue -
    // the iridescence is a property of the keratin, not a change of species.
    float fres = pow(1.0 - max(dot(n, -rd), 0.0), 3.4);
    vec3 sheenHue = 0.5 + 0.5*cos(TAU*(vec3(0.00, 0.33, 0.67) + fres*1.35 + arc*0.9));
    vec3 sheen = mix(mix(base, vec3(0.55, 0.72, 1.00), fres*0.45), sheenHue, 0.33*gloss);
    lit += sheen*fres*(0.55 + 0.85*gloss);
    lit += vec3(0.92, 0.96, 1.00)*glint*0.85;

    // Forward scattering through the thin end of the animal. A tail tapering to
    // a couple of millimetres of translucent keratin lights up when the key is
    // behind it, and that single term does more to sell these as bodies with
    // volume than any amount of extra specular.
    float thin = 1.0 - smoothstep(0.16, 0.72, radiusAt(u)/max(girth, 1e-4));
    float through = pow(max(dot(rd, key), 0.0), 3.0);
    lit += base*through*thin*(0.85 + 0.55*gloss);

    lit *= occ;
    colour = mix(background, lit, exp(-max(tHit - 2.4, 0.0)*0.16));
  }

  colour += bloom(ro, rd, tHit > 0.0 ? tHit : 12.0)*glow;

  // Filmic curve then gamma. Working roughly linear up to here is what lets the
  // three saturated lights sit on top of one another without clipping to white.
  colour *= 0.92;
  vec3 a = colour*(2.51*colour + 0.03);
  vec3 b = colour*(2.43*colour + 0.59) + 0.14;
  colour = clamp(a/b, 0.0, 1.0);
  colour = pow(colour, vec3(0.4545));

  // One least-significant bit of noise before quantisation. Most of this frame
  // is a very dark, very shallow gradient, and eight bits across that range puts
  // the contour lines wide enough apart to see as bands. Dithering trades them
  // for noise below the level the eye resolves.
  float grain = hash11(gl_FragCoord.x*3.7 + gl_FragCoord.y*811.3 + 0.5) - 0.5;
  colour += grain/255.0;

  gl_FragColor = vec4(colour, 1.0);
}

void main(){
  if(PASSINDEX == 0){ mainHeads(); return; }
  if(PASSINDEX == 1){ mainSpine(); return; }
  if(PASSINDEX == 2){ mainBounds(); return; }
  mainImage();
}