/*{
  "DESCRIPTION": "Vein Atlas: a dark mineral plain cracked open by living light. Three colonies drain their own territory through a flow-accumulation network computed on drifting fractal terrain, so a channel that carries load thickens and a channel that loses its nutrient withers and seals over into cold stone. Each colony carries a vigour that tracks its share of the nutrients currently on the plain, which widens or collapses its territory - so the three trade dominance and the plain is never twice the same. Rendered as a marched heightfield lit almost entirely by its own emission, with pulses running the fissures and the glow scattering into the haze above.",
  "CATEGORIES": ["generative", "simulation", "organic"],
  "INPUTS": [
    { "NAME": "tempo", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.15, "MAX": 3.0 },
    { "NAME": "density", "TYPE": "float", "DEFAULT": 0.50, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "wither", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.05, "MAX": 2.0 },
    { "NAME": "nutrition", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "meander", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 2.4 },
    { "NAME": "carve", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0, "MAX": 3.0 },
    { "NAME": "swirl", "TYPE": "float", "DEFAULT": 0.85, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "sealing", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "relief", "TYPE": "float", "DEFAULT": 0.125, "MIN": 0.02, "MAX": 0.30 },
    { "NAME": "emissive", "TYPE": "float", "DEFAULT": 0.45, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "pulse", "TYPE": "float", "DEFAULT": 0.75, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "haze", "TYPE": "float", "DEFAULT": 0.85, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "bloom", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "seed", "TYPE": "float", "DEFAULT": 7.0, "MIN": 0.0, "MAX": 30.0 }
  ],
  "PASSES": [
    { "TARGET": "acc1", "FLOAT": true },
    { "TARGET": "acc", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "net", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "mem", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "glow", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "scene", "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;`.
// This shader integrates conductance and crust across thousands of frames and
// then differentiates a stored heightfield for normals; mediump would band both.
#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI  3.14159265359
#define TAU 6.28318530718

#define COLONIES 3
#define NUTS 9            // nutrient slots, all analytic from the sim clock

// Seconds of simulated time per clock epoch. The clock is stored as
// (fraction, epoch) rather than as one number because persistent targets are
// RGBA16F: a single accumulating seconds value reaches a magnitude where one
// frame's increment falls below the mantissa and the whole piece freezes.
#define CLOCK_UNIT 8.0
#define CLOCK_WRAP 1024.0   // epochs before wrap; 1024 is exact in half float

// Routing is restricted to neighbours that strictly decrease distance to the
// colony's heart. That single constraint is what makes the whole solver safe:
// heart distance is a strictly decreasing quantity along every flow path, so
// cycles are impossible and every drop reaches a heart, no matter how violent
// the terrain or how strong the carving feedback.
//
// The first version instead subtracted the carve term from a potential and took
// steepest descent on it. It looked correct for two seconds and then degraded:
// a cell whose conductance had lowered its own potential below its downstream
// neighbour became a sink, trapped its flow, and starved everything behind it.
// The field filled with stalled bright cells and the trees broke into diagonal
// streaks. Terrain amplitude and carve strength are free artistic parameters
// here precisely because neither can violate the constraint.
// How strongly the routing prefers the neighbour that gets closest to the heart
// fastest. It only breaks ties: the admissibility test above already guarantees
// progress, so this is free to be small. Large values make every path a
// straight radial ray, and because there are only eight directions to choose
// from those rays quantise into parallel bands.
#define PROGRESS_W 0.28

// Static per-cell tiebreak. Without it, every cell in a region where the
// terrain is locally flat picks the same one of eight directions and the field
// fills with parallel diagonal streaks - the classic D8 artifact. Hashed on the
// cell, not on time, or the network flickers.
#define JITTER 0.55

#define NUT_R 0.085       // nutrient blob radius in world units
#define WARM_BOOST 16.0   // cold-start dt multiplier, see clockUpdate()
#define WARM_RATE 0.25

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

float hash11(float p){ return fract(sin(p*127.1 + seed*19.73)*43758.5453); }
float hash12(vec2 p){ return fract(sin(dot(p, vec2(127.1, 311.7)) + seed*3.17)*43758.5453); }

// The sim grid is a square block of texels in the bottom-left of every target,
// with one control row directly above it. Deriving the edge from RENDERSIZE
// instead of fixing it means a render smaller than the block degrades to a
// coarser simulation rather than reading texels that do not exist.
float simEdge(){ return min(160.0, floor(min(RENDERSIZE.x, RENDERSIZE.y)) - 2.0); }

vec4 fetchTexel(sampler2D tex, vec2 texel){
  return texture2D(tex, (texel + 0.5)/RENDERSIZE);
}

// Every target is allocated with NEAREST filtering, so upsampling the sim grid
// to render resolution has to be done by hand. Smoothstepped weights rather
// than linear ones because the heightfield is differentiated for normals and
// plain bilinear is only C0 - the facet edges show up as creases in the light.
vec4 simBilerp(sampler2D tex, vec2 grid, float N){
  vec2 g = grid - 0.5;
  vec2 i = floor(g);
  vec2 f = g - i;
  f = f*f*(3.0 - 2.0*f);
  vec2 lim = vec2(N - 1.0);
  vec4 a = fetchTexel(tex, clamp(i, vec2(0.0), lim));
  vec4 b = fetchTexel(tex, clamp(i + vec2(1.0, 0.0), vec2(0.0), lim));
  vec4 c = fetchTexel(tex, clamp(i + vec2(0.0, 1.0), vec2(0.0), lim));
  vec4 d = fetchTexel(tex, clamp(i + vec2(1.0, 1.0), vec2(0.0), lim));
  return mix(mix(a, b, f.x), mix(c, d, f.x), f.y);
}

// D8 neighbourhood. A function rather than a const array because indexing a
// const array with a non-constant expression is not guaranteed in ES 1.00.
vec2 dirVec(float d){
  if(d < 0.5) return vec2( 1.0,  0.0);
  if(d < 1.5) return vec2( 1.0,  1.0);
  if(d < 2.5) return vec2( 0.0,  1.0);
  if(d < 3.5) return vec2(-1.0,  1.0);
  if(d < 4.5) return vec2(-1.0,  0.0);
  if(d < 5.5) return vec2(-1.0, -1.0);
  if(d < 6.5) return vec2( 0.0, -1.0);
  return vec2( 1.0, -1.0);
}

// Three D8 directions dirBits into one channel, base 9 so that 8 can mean "no
// outflow" (a cell that is a local minimum). Max value 728, and half float
// represents every integer to 2048 exactly, so this survives storage.
float packDirs(float a, float b, float c){ return a + 9.0*b + 81.0*c; }
float unpackDir(float p, int c){
  if(c == 0) return mod(p, 9.0);
  if(c == 1) return mod(floor(p/9.0), 9.0);
  return floor(p/81.0);
}

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

vec2 cellToWorld(vec2 texel, float N){ return ((texel + 0.5)/N)*2.0 - 1.0; }
vec2 worldToGrid(vec2 w, float N){ return (w*0.5 + 0.5)*N; }

// D8 routing can only ever move along an axis or a diagonal, so every channel
// is built from straight runs meeting at ninety and forty-five degrees. That is
// invisible in a flat plot of the field and unmissable once the channel is
// given depth and lit: the plain turns into blocky terraces.
//
// Warping the lookup breaks the staircase without touching the simulation. It
// has to be shared between the height bake and the shading, or the emission
// slides off the fissure it is supposed to be sitting in.
// Amplitude has to stay well under half the wavelength or the map folds, and a
// folded lookup samples the same channel several times over and lays a chain of
// discrete blobs along every vein. Two gentle octaves reach the same break-up
// with a Jacobian around 0.6.
vec2 warpWorld(vec2 w){
  vec2 a = vec2(vnoise(w*6.5 + 3.1) - 0.5, vnoise(w*6.5 + 19.7) - 0.5);
  vec2 b = vec2(vnoise(w*15.0 + 7.7) - 0.5, vnoise(w*15.0 + 31.3) - 0.5);
  return w + a*0.019 + b*0.007;
}

// Flow must not leave the plain, and the plain's edge must not be visible as an
// edge. Both are the same mask: conductance and rain fade out before the
// boundary, so the rim is uniform dead crust that the fog swallows.
float borderMask(vec2 w){
  float e = min(1.0 - abs(w.x), 1.0 - abs(w.y));
  return smoothstep(0.02, 0.26, e);
}

// ---------------------------------------------------------------------------
// Clock
// ---------------------------------------------------------------------------

vec4 clockCell(){ return fetchTexel(net, vec2(0.0, simEdge())); }

// Per-colony vigour, the piece's only global state besides the clock. It is
// what makes the three colonies a conflict rather than three independent trees:
// vigour tracks each colony's *share* of the nutrients currently on the plain,
// so it is zero sum, and it scales both how far that colony's network reaches
// and how far its territory border sits from its heart.
vec4 vigourCell(){ return fetchTexel(net, vec2(1.0, simEdge())); }

vec3 colonyVigour(){
  vec4 v = vigourCell();
  if(v.a < 0.5) return vec3(1.0);          // cold, or not yet written
  return clamp(v.rgb, 0.35, 2.6);
}

// Simulated seconds. Everything physical is driven from this rather than from
// TIME: glslop renders a poster by ramping to the requested time in 120 steps,
// so TIMEDELTA there can be an order of magnitude larger than a live frame's
// and anything read off the clock desynchronises from anything integrated.
float simClock(){
  vec4 c = clockCell();
  return (c.g + c.r)*CLOCK_UNIT;
}

float simDelta(){
  vec4 c = clockCell();
  float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
  return dtr*(1.0 + WARM_BOOST*c.a)*tempo;
}

// ---------------------------------------------------------------------------
// Terrain, hearts, nutrients
// ---------------------------------------------------------------------------


// The terrain is what makes the drainage dendritic instead of radial, and its
// slow drift is what keeps the network reorganising after the nutrients have
// settled into a rhythm. Each octave drifts a different way so the field
// evolves rather than translating.
// Amplitude falls faster than frequency rises, so the gradient is dominated by
// the large scales. With the usual 0.5 falloff every octave contributes about
// the same slope and the routing decision is made by per-cell jitter, which
// reads as a noisy path rather than as a meander.
float terrain(vec2 w, float t){
  float s = 0.0, amp = 0.5, f = 2.4;
  vec2 dr = vec2(0.013, -0.009);
  for(int o = 0; o < 4; o++){
    s += amp*vnoise(w*f + dr*t*f);
    f *= 2.05;
    amp *= 0.44;
    dr = vec2(-dr.y, dr.x)*1.31;
  }
  return s;
}

// One shared random rotation, not one per colony: an independent offset per
// heart collapses two of them together often enough that a colony gets squeezed
// out entirely, which is what happened on the first run.
// Curl of a scalar noise, i.e. a divergence-free direction field. Steepest
// descent on a smooth slope points the same way for every cell in a
// neighbourhood, and with only eight directions available that prints long
// straight runs meeting at right angles - the network reads as circuit traces.
// Biasing the choice along a curl field instead makes the paths sinuous, and
// because it cannot override the heart-distance constraint it cannot stall the
// flow or close a loop.
vec2 curlDir(vec2 w, float t){
  float e = 0.030;
  vec2 dr = vec2(0.0, t*0.008);
  float p0 = vnoise(w*3.1 + dr) + 0.5*vnoise(w*7.4 + dr*1.7);
  float px = vnoise((w + vec2(e, 0.0))*3.1 + dr) + 0.5*vnoise((w + vec2(e, 0.0))*7.4 + dr*1.7);
  float py = vnoise((w + vec2(0.0, e))*3.1 + dr) + 0.5*vnoise((w + vec2(0.0, e))*7.4 + dr*1.7);
  vec2 g = vec2(px - p0, py - p0);
  vec2 c = vec2(g.y, -g.x);
  return c/max(length(c), 1e-5);
}

vec2 heartPos(int c, float t){
  float fc = float(c);
  float a = TAU*fc/3.0 + t*0.019 + hash11(3.0)*TAU;
  float r = 0.30 + 0.075*sin(t*0.031 + fc*2.1);
  vec2 p = vec2(cos(a), sin(a))*r;
  // An independent slow wander on top of the shared orbit. Without it the three
  // hearts turn in lockstep, the territory partition is rigid, and the only
  // thing that ever changes at large scale is which way the whole plain faces.
  p += vec2(vnoise(vec2(t*0.021, fc*7.3 + 2.0)) - 0.5,
            vnoise(vec2(fc*3.1 + 5.0, t*0.021 + 11.0)) - 0.5)*0.30;
  return p;
}

// Nutrient slot k as a pure function of the sim clock: position, and an
// amplitude envelope that ramps in, holds, and decays. No state, so there is no
// cold-start branch to go unrun, and no way for a poster to render an
// uninitialised nutrient table.
vec3 nutrient(int k, float t){
  float fk = float(k);
  float period = 24.0 + 14.0*hash11(fk*3.7 + 2.0);
  float u = (t + hash11(fk*9.1 + 5.0)*period)/period;
  float e = mod(floor(u), 64.0);   // bounded so the hash stays well conditioned
  float ph = fract(u);
  vec2 pos = vec2(hash11(e*13.1 + fk*57.3), hash11(e*29.7 + fk*17.9))*1.44 - 0.72;
  float amt = smoothstep(0.0, 0.13, ph)*(1.0 - smoothstep(0.52, 0.98, ph));
  return vec3(pos, amt);
}

// Soft nearest-heart partition. Territory is a smoothed Voronoi of the drifting
// hearts, so the visible competition is the interdigitation along the seams
// rather than a hard line.
// Returns the three shares in xyz and the winning *effective* distance in w.
//
// Distance is divided by vigour rather than the exponent being weighted by it:
// a multiplicatively weighted partition moves the border a long way for a
// modest strength difference, where weighting the softmax only nudges the seam.
// The same effective distance then sets how far the network reaches, so a
// colony that is winning visibly grows and a starving one visibly pulls back.
vec4 heartShare(vec2 w, float t, vec3 vig){
  vec3 s = pow(max(vig, vec3(0.05)), vec3(0.45));
  vec3 d;
  d.x = length(w - heartPos(0, t))/s.x;
  d.y = length(w - heartPos(1, t))/s.y;
  d.z = length(w - heartPos(2, t))/s.z;
  float m = min(d.x, min(d.y, d.z));
  // Sharp. A soft partition gives every colony a little rain everywhere, so all
  // three networks span the whole plain and overlay into mud; only the seam
  // should be shared.
  vec3 e = exp(-(d - vec3(m))*22.0);
  return vec4(e/max(e.x + e.y + e.z, 1e-5), m);
}

vec3 colonyColour(int c){
  if(c == 0) return vec3(1.00, 0.42, 0.075);   // molten gold
  if(c == 1) return vec3(0.13, 0.82, 1.00);    // cold cyan
  return vec3(0.72, 0.16, 1.00);               // violet
}

// ---------------------------------------------------------------------------
// PASS 0 - acc1: terrain, nutrient attraction, and D8 flow routing
// ---------------------------------------------------------------------------

void mainAcc1(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= N){ gl_FragColor = vec4(0.0); return; }

  float t = simClock();
  vec2 w = cellToWorld(tx, N);
  float T = terrain(w, t);

  float A = 0.0;
  for(int k = 0; k < NUTS; k++){
    vec3 nk = nutrient(k, t);
    vec2 dv = w - nk.xy;
    A += nk.z*exp(-dot(dv, dv)/(NUT_R*NUT_R));
  }

  vec4 me = fetchTexel(net, tx);
  float cellW = 2.0/N;

  vec2 cd = curlDir(w, t);
  vec2 h0 = heartPos(0, t), h1 = heartPos(1, t), h2 = heartPos(2, t);
  vec3 dSelf = vec3(length(w - h0), length(w - h1), length(w - h2));

  vec3 best = vec3(-1e9);
  float d0 = 8.0, d1 = 8.0, d2 = 8.0;            // 8 == no outflow (the heart)

  for(int k = 0; k < 8; k++){
    float fk = float(k);
    vec2 off = dirVec(fk);
    vec2 nx = tx + off;
    if(nx.x >= 0.0 && nx.y >= 0.0 && nx.x < N && nx.y < N){
      vec2 nw = cellToWorld(nx, N);
      float nT = terrain(nw, t);
      vec4 nn = fetchTexel(net, nx);
      float step = length(off)*cellW;
      // Downhill terrain slope and the pull of an existing channel decide which
      // of the admissible neighbours wins; the admissibility test decides that
      // the path terminates.
      float bias = meander*2.30*(T - nT)/step
                 - sealing*0.35*nn.a
                 + swirl*1.45*dot(off/length(off), cd)
                 + JITTER*(hash12(tx*3.1 + off*17.3) - 0.5);

      vec3 prog = (dSelf - vec3(length(nw - h0), length(nw - h1), length(nw - h2)))/step;
      if(prog.x > 0.0){ float s = PROGRESS_W*prog.x + bias + carve*2.4*nn.r;
                        if(s > best.x){ best.x = s; d0 = fk; } }
      if(prog.y > 0.0){ float s = PROGRESS_W*prog.y + bias + carve*2.4*nn.g;
                        if(s > best.y){ best.y = s; d1 = fk; } }
      if(prog.z > 0.0){ float s = PROGRESS_W*prog.z + bias + carve*2.4*nn.b;
                        if(s > best.z){ best.z = s; d2 = fk; } }
    }
  }

  gl_FragColor = vec4(packDirs(d0, d1, d2), T, A, 1.0);
}

// ---------------------------------------------------------------------------
// PASS 1 - acc: one flow-accumulation gather per frame
// ---------------------------------------------------------------------------

void mainAcc(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= N){ gl_FragColor = vec4(0.0); return; }

  float t = simClock();
  vec2 w = cellToWorld(tx, N);
  float mask = borderMask(w);

  // Rain. Scaled by 1/N^2 so that the fully accumulated trunk lands near 1.0:
  // half float has plenty of relative precision but a catchment summed as raw
  // cell counts would reach 25k, where the quantum is 16.
  vec3 vig = colonyVigour();
  vec4 hsw = heartShare(w, t, vig);
  vec3 share = hsw.xyz;
  // Baseline rain thins with distance from the heart, so a colony's catchment
  // fades out into dead stone instead of tiling the whole plain with capillaries.
  // Without this every cell drains somewhere and the plain is uniformly covered,
  // which is botanically correct and compositionally flat.
  //
  // The distance is the vigour-weighted one, so a colony that is winning the
  // nutrient share grows its whole catchment and a starving one pulls back into
  // its heart and lets its outer ground seal over.
  float reach = exp(-hsw.w*2.00);
  vec3 rain = share*(2.6/(N*N))*mask*(0.10 + 0.90*reach);

  // Nutrients inject flow at their own cell. Because routing is defined on the
  // terrain everywhere - not only where the network already exists - that flow
  // immediately has a drainage path to the heart, and the vein grows along it.
  // This is why there is no separate "exploration" rule: reaching for food is
  // what accumulation already does.
  for(int k = 0; k < NUTS; k++){
    vec3 nk = nutrient(k, t);
    vec2 dv = w - nk.xy;
    float blob = nk.z*exp(-dot(dv, dv)/(NUT_R*NUT_R));
    if(blob > 0.001){
      vec4 nd = fetchTexel(net, clamp(worldToGrid(nk.xy, N), vec2(0.0), vec2(N - 1.0)));
      vec3 hs = heartShare(nk.xy, t, vig).xyz;
      // Whoever already has conductance at the nutrient gets more of it, so a
      // colony that arrives first holds the supply until a rival grows thicker.
      vec3 claim = hs*(0.35 + 2.4*vec3(nd.r, nd.g, nd.b));
      claim /= max(claim.x + claim.y + claim.z, 1e-5);
      // Sized against the baseline rather than picked. Baseline rain totals
      // about 2.6 over a whole catchment; a nutrient blob covers only
      // pi*(NUT_R*N/2)^2 cells, so matching it takes a per-cell rate roughly
      // 0.6/NUT_R^2 times larger. At the first value tried, one nutrient was
      // worth a tenth of the baseline and nothing visibly reached for it.
      rain += claim*blob*nutrition*(55.0/(N*N))*mask;
    }
  }

  vec3 R = rain;
  for(int k = 0; k < 8; k++){
    vec2 off = dirVec(float(k));
    vec2 nx = tx + off;
    if(nx.x >= 0.0 && nx.y >= 0.0 && nx.x < N && nx.y < N){
      float dirBits = fetchTexel(acc1, nx).r;
      vec4 rn = fetchTexel(acc, nx);
      for(int c = 0; c < COLONIES; c++){
        float d = unpackDir(dirBits, c);
        if(d < 8.5){
          vec2 flow = dirVec(d) + off;          // neighbour drains into me when
          if(dot(flow, flow) < 0.1){            // its direction is exactly -off
            if(c == 0) R.x += rn.x;
            else if(c == 1) R.y += rn.y;
            else R.z += rn.z;
          }
        }
      }
    }
  }

  gl_FragColor = vec4(min(R, vec3(8.0)), 1.0);
}

// ---------------------------------------------------------------------------
// PASS 2 - net: conductance, crust, and the clock
// ---------------------------------------------------------------------------

void mainNet(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);

  // Control texel. Cold is detected from the buffer's own contents rather than
  // from FRAMEINDEX, which glslop derives from the requested poster time and
  // which therefore starts at ~480 on a shader's first ever frame.
  if(tx.x == 0.0 && tx.y == N){
    vec4 c = fetchTexel(net, vec2(0.0, N));
    if(c.b < 0.5){
      gl_FragColor = vec4(0.0, 0.0, 1.0, 1.0);
      return;
    }
    float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
    float dts = dtr*(1.0 + WARM_BOOST*c.a)*tempo;
    float frac = c.r + dts/CLOCK_UNIT;
    float epoch = c.g + floor(frac);
    gl_FragColor = vec4(fract(frac), mod(epoch, CLOCK_WRAP), 1.0, c.a*exp(-dts*WARM_RATE));
    return;
  }

  // Vigour texel. Everything physical in the piece is a function of the sim
  // clock, so scaling the clock with `tempo` scales the hearts, the nutrient
  // cycle, growth, decay and this - coherently, which no per-rate control can.
  if(tx.x == 1.0 && tx.y == N){
    vec4 v = fetchTexel(net, vec2(1.0, N));
    vec4 c0 = clockCell();
    if(c0.b < 0.5 || v.a < 0.5){
      gl_FragColor = vec4(1.0, 1.0, 1.0, 1.0);
      return;
    }
    float t = simClock();
    float dts = simDelta();
    vec3 vig = clamp(v.rgb, 0.35, 2.6);

    // Who is currently drinking from each nutrient. Proximity decides who can
    // reach it at all; conductance already present decides who holds it, so
    // arriving first is worth something and a rival has to out-grow you to take
    // it back.
    vec3 harvest = vec3(0.0);
    for(int k = 0; k < NUTS; k++){
      vec3 nk = nutrient(k, t);
      if(nk.z > 0.02){
        vec4 nd = fetchTexel(net, clamp(worldToGrid(nk.xy, N), vec2(0.0), vec2(N - 1.0)));
        // Deliberately UNWEIGHTED by vigour. Scoring the nutrient with the same
        // vigour-widened territory that vigour itself controls is a positive
        // feedback loop, and it ran away: whichever colony led erased the other
        // two inside a minute. Proximity is fixed by where the nutrient landed,
        // so the competition stays a lottery and vigour stays mean reverting.
        vec4 hs = heartShare(nk.xy, t, vec3(1.0));
        vec3 claim = hs.xyz*(0.45 + 1.2*vec3(nd.r, nd.g, nd.b));
        claim /= max(claim.x + claim.y + claim.z, 1e-5);
        harvest += claim*nk.z;
      }
    }

    // Normalised to shares, which makes the competition zero sum: one colony
    // can only boom at the expense of the others. Absolute nutrient abundance
    // deliberately does not raise everyone at once - that would be a brightness
    // control, not a conflict.
    float tot = harvest.x + harvest.y + harvest.z;
    vec3 frac2 = tot > 1e-4 ? harvest/tot : vec3(1.0/3.0);
    vec3 targetV = vec3(0.45) + (frac2 - vec3(1.0/3.0))*1.55*nutrition + vec3(0.55);
    vig += (targetV - vig)*(1.0 - exp(-dts*0.22));
    gl_FragColor = vec4(clamp(vig, 0.35, 2.6), 1.0);
    return;
  }
  if(tx.x >= N || tx.y >= N){ gl_FragColor = vec4(0.0); return; }

  vec4 c = clockCell();
  vec2 w = cellToWorld(tx, N);
  float t = simClock();
  float dts = simDelta();
  float mask = borderMask(w);

  vec4 me = fetchTexel(net, tx);
  vec3 D = vec3(me.r, me.g, me.b);
  float crust = me.a;

  if(c.b < 0.5){
    // Virgin ground: sealed stone, no network anywhere.
    gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
    return;
  }

  vec3 R = fetchTexel(acc, tx).rgb;

  // Thickness follows accumulated flow, compressed hard: a trunk carries
  // thousands of times the flow of a capillary and a linear map would show only
  // the trunk.
  // `density` sets how much accumulated flow a channel needs to survive, which
  // is a structural control: it changes what the network IS, not how fast it
  // gets there. Low leaves only trunks, high fills the territory with capillaries.
  float cut = mix(0.80, 0.22, clamp(density, 0.0, 1.0));
  vec3 target = clamp(vec3(2.40)*pow(max(R, vec3(0.0)) + 1e-6, vec3(0.385)) - vec3(cut), vec3(0.0), vec3(1.0));
  target *= mask;

  // The heart is always alive. This doubles as the seed: there is no init
  // branch anywhere in the simulation, so there is nothing for a mis-numbered
  // first frame to skip.
  for(int k = 0; k < COLONIES; k++){
    float dh = length(w - heartPos(k, t));
    float blob = exp(-dh*dh/(0.030*0.030));
    if(k == 0) target.x = max(target.x, blob*0.95);
    else if(k == 1) target.y = max(target.y, blob*0.95);
    else target.z = max(target.z, blob*0.95);
  }

  // Written as an exponential approach rather than target*rate*dt so that the
  // cold-start dt boost - which can hand this a 0.7s step - stays stable.
  vec3 rate = vec3(
    target.x > D.x ? 2.2 : wither*0.30,
    target.y > D.y ? 2.2 : wither*0.30,
    target.z > D.z ? 2.2 : wither*0.30);
  D += (target - D)*(vec3(1.0) - exp(-rate*dts));

  float activity = max(D.x, max(D.y, D.z));
  float open = smoothstep(0.04, 0.42, activity);
  // `sealing` sets how far dead ground heals, not how fast it heals. As a rate
  // it was invisible: the crust reached the same steady state either way and a
  // sweep across its whole range produced four identical frames. As a target it
  // is the difference between a plain that stays broken and dark and one that
  // closes over into pale sealed stone between the veins.
  float ct = (1.0 - open)*clamp(sealing*0.72, 0.0, 1.15);
  crust += (ct - crust)*(1.0 - exp(-(0.30 + open*1.6)*dts));

  gl_FragColor = vec4(clamp(D, 0.0, 1.0), clamp(crust, 0.0, 1.0));
}

// ---------------------------------------------------------------------------
// PASS 3 - mem: scar tint and the baked surface height
// ---------------------------------------------------------------------------

void mainMem(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= N){ gl_FragColor = vec4(0.0); return; }

  float dts = simDelta();
  vec4 me = fetchTexel(mem, tx);
  vec4 nt = fetchTexel(net, tx);
  float T = fetchTexel(acc1, tx).g;

  // Scar: where a colony has been, faintly and for a long time after it leaves.
  // Two timescales in one channel is impossible, which is why this is stored
  // separately from conductance rather than derived from a slow decay of it.
  vec3 scar = me.rgb;
  scar += (clamp(nt.rgb*1.25, 0.0, 1.0) - scar)*(1.0 - exp(-0.55*dts));
  scar *= exp(-0.0090*dts);

  // Height is baked from a tent-weighted 3x3 rather than from this cell alone.
  // A channel goes from full conductance to none across a single cell, so an
  // unsmoothed bake drops 1.4 relief units over one cell width - a vertical
  // cliff - and the interpolated march turns every such wall into a stair. The
  // blur turns the slot into a V-groove that the light can actually sit in.
  vec2 w = cellToWorld(tx, N);
  vec2 gq = worldToGrid(warpWorld(w), N);
  // 5x5 at one-cell spacing. The first attempt used 3x3 at 1.45 cells, which
  // undersamples a channel one cell wide: the blur then beads along every vein
  // at the aliasing period and the render grew a chain of glowing lumps down
  // each one. A smoothing kernel has to sample at or below the feature rate it
  // is smoothing.
  float wsum = 0.0;
  float aBlur = 0.0;
  float cBlur = 0.0;
  for(int j = 0; j < 25; j++){
    vec2 off = vec2(mod(float(j), 5.0) - 2.0, floor(float(j)/5.0) - 2.0);
    vec4 s = simBilerp(net, gq + off, N);
    float wgt = exp(-dot(off, off)*0.22);
    aBlur += max(s.r, max(s.g, s.b))*wgt;
    cBlur += s.a*wgt;
    wsum += wgt;
  }
  aBlur /= wsum;
  cBlur /= wsum;

  float open = smoothstep(0.04, 0.50, aBlur);
  // Crust stands proud, the fissure cuts down through it, and the terrain rolls
  // underneath both. The march reads this one channel per step, which is why it
  // is baked here instead of recomputed.
  float h = relief*(0.55*(T - 0.45)*2.0 + 0.62*cBlur) - relief*0.46*open;

  gl_FragColor = vec4(clamp(scar, 0.0, 1.0), h);
}

// ---------------------------------------------------------------------------
// PASS 4 - glow: blurred emission, at sim resolution
// ---------------------------------------------------------------------------

vec3 emissionAt(vec2 texel, float t){
  vec4 nt = fetchTexel(net, texel);
  vec3 D = nt.rgb;
  vec3 e = colonyColour(0)*D.x + colonyColour(1)*D.y + colonyColour(2)*D.z;
  return e*smoothstep(0.02, 0.30, max(D.x, max(D.y, D.z)));
}

// A conductance channel is one cell wide and runs diagonally, and a one-texel
// diagonal line put through bilinear upsampling reconstructs as a chain of
// diamonds - which is why every vein in the render was beaded. Nothing about
// the warp, the height bake, the interpolation weights or the fissure depth
// could fix that, because the line was already too thin in the field.
//
// So the render's emission is taken from here rather than from the raw network:
// this is a smooth field about three cells across, which upsamples cleanly.
//
// The blur is a 3x3 diffusion of the previous frame plus an injection, not a
// single wide disc. Repeated small kernels reach a wide smooth Gaussian for
// nine taps, and - more importantly - a 3x3 at one-cell spacing samples every
// cell, where the 17-tap disc of radius 6.5 that was here first sampled a
// 130-cell area at about a third of the channel rate and beaded on its own.
// The colour is injected fast and the density slowly, so one pass yields a
// tight core for the fissure and a broad halo for the air above it.
void mainGlow(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= N){ gl_FragColor = vec4(0.0); return; }

  float t = simClock();
  vec3 blurC = vec3(0.0);
  float blurD = 0.0;
  float wsum = 0.0;
  for(int j = 0; j < 9; j++){
    vec2 off = vec2(mod(float(j), 3.0) - 1.0, floor(float(j)/3.0) - 1.0);
    vec2 nx = clamp(tx + off, vec2(0.0), vec2(N - 1.0));
    vec4 g = fetchTexel(glow, nx);
    float wgt = exp(-dot(off, off)*0.62);
    blurC += g.rgb*wgt;
    blurD += g.a*wgt;
    wsum += wgt;
  }
  blurC /= wsum;
  blurD /= wsum;

  vec3 e = emissionAt(tx, t);
  float ed = max(e.x, max(e.y, e.z));
  gl_FragColor = vec4(mix(blurC, e, 0.21), mix(blurD, ed, 0.085));
}

// ---------------------------------------------------------------------------
// PASS 5 - scene: march the heightfield
// ---------------------------------------------------------------------------

float heightAt(vec2 w, float N){
  return simBilerp(mem, worldToGrid(w, N), N).a;
}

// Single tap. The baked height is already smoothed with a 5x5 kernel, so the
// difference from the bilinear read is far below the shadow test's own
// threshold - and the shadow march is a fixed cost on every lit pixel, where
// the interpolated read is four times the bandwidth for no visible return.
float heightFast(vec2 w, float N){
  vec2 g = clamp(floor(worldToGrid(w, N)), vec2(0.0), vec2(N - 1.0));
  return fetchTexel(mem, g).a;
}

vec3 sceneNormal(vec2 w, float N, float e){
  float hx = heightAt(w + vec2(e, 0.0), N) - heightAt(w - vec2(e, 0.0), N);
  float hz = heightAt(w + vec2(0.0, e), N) - heightAt(w - vec2(0.0, e), N);
  return normalize(vec3(-hx, 2.0*e, -hz));
}

// Slab intersection against the plain's bounding box. Marching only the span
// that can possibly contain the surface is what keeps the step count honest.
bool boxSpan(vec3 ro, vec3 rd, float lo, float hi, out float t0, out float t1){
  vec3 bmin = vec3(-1.55, lo, -1.55);
  vec3 bmax = vec3( 1.55, hi,  1.55);
  vec3 inv = 1.0/(rd + vec3(equal(rd, vec3(0.0)))*1e-6);
  vec3 a = (bmin - ro)*inv;
  vec3 b = (bmax - ro)*inv;
  vec3 lo3 = min(a, b), hi3 = max(a, b);
  t0 = max(max(lo3.x, lo3.y), lo3.z);
  t1 = min(min(hi3.x, hi3.y), hi3.z);
  t0 = max(t0, 0.0);
  return t1 > t0;
}

void mainScene(){
  float N = simEdge();
  vec2 uv = (gl_FragCoord.xy - 0.5*RENDERSIZE)/min(RENDERSIZE.x, RENDERSIZE.y);
  float t = simClock();

  // The camera is driven from TIME, not from the sim clock, deliberately: on a
  // poster the sim clock is capped by the step schedule but TIME is not, so
  // this is the one control that actually varies the composition with --time.
  // High and tilted, not low and grazing. A fissure is a few cells wide and
  // most of a relief unit deep, so from near the horizon the crust between two
  // channels presents edge-on and the plain reads as a forest of pillars. Above
  // about fifty degrees it reads as ground with cracks in it, which is the
  // whole idea, and there is still enough obliquity for parallax and for the
  // key light to graze.
  // A wide lens low over the plain buries the whole composition in perspective:
  // the nearest heart fills the bottom third and everything beyond it
  // compresses into the top edge. Higher, further back, and a longer lens
  // instead lays the three territories out across the frame with enough
  // obliquity left for the relief to read.
  // uv spans +/-0.5, not +/-1, so the half angle is atan(0.5/focal). Reading it
  // as atan(1/focal) put the camera at half the intended field of view and the
  // frame covered 0.7 of a two-unit plain - which is why single simulation
  // cells were seven pixels across and every field read as blocky. Verified
  // with a known-scale checkerboard rather than recomputed.
  float ca = TIME*0.045 + 2.20;
  float el = 1.90 + 0.09*sin(TIME*0.043);
  vec3 ro = vec3(cos(ca), 0.0, sin(ca))*1.25 + vec3(0.0, el, 0.0);
  vec3 ta = vec3(0.0, -0.02, 0.0);
  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(uv.x*rt + uv.y*up + 1.55*fw);

  float loH = -relief*1.10, hiH = relief*1.30;
  float t0, t1;
  bool span = boxSpan(ro, rd, loH, hiH, t0, t1);

  float tHit = -1.0;
  vec3 hitP = vec3(0.0);
  if(span){
    float tc = t0;
    float step = (t1 - t0)/36.0;
    float prevD = 1.0;
    float prevT = t0;
    for(int i = 0; i < 36; i++){
      vec3 p = ro + rd*tc;
      float d = p.y - heightAt(p.xz, N);
      if(d < 0.0){
        // Bisect between the last point above the surface and this one below.
        float lo = prevT, hi = tc;
        for(int j = 0; j < 8; j++){
          float mid = 0.5*(lo + hi);
          vec3 q = ro + rd*mid;
          if(q.y - heightAt(q.xz, N) < 0.0) hi = mid; else lo = mid;
        }
        tHit = hi;
        hitP = ro + rd*tHit;
        break;
      }
      prevD = d;
      prevT = tc;
      tc += step;
      if(tc > t1) break;
    }
  }

  vec3 col = vec3(0.0);
  float far = span ? t1 : 4.5;
  float tEnd = tHit > 0.0 ? tHit : far;

  if(tHit > 0.0){
    vec2 w = hitP.xz;
    // The network is read through the same warp the height was baked with, so
    // the emission stays inside the fissure it belongs to. The height itself is
    // read unwarped - it was already baked for this world position, and warping
    // it again would shift the surface out from under the shading.
    vec2 grid = worldToGrid(w, N);
    vec2 gridW = worldToGrid(warpWorld(w), N);
    vec4 nt = simBilerp(net, gridW, N);
    vec4 mm = simBilerp(mem, grid, N);
    vec3 scarW = simBilerp(mem, gridW, N).rgb;
    vec4 a1 = simBilerp(acc1, grid, N);
    vec3 gl = simBilerp(glow, gridW, N).rgb;

    vec3 n = sceneNormal(w, N, 1.4/N);
    // Fine mineral grain, evaluated once at the hit rather than per march step.
    // Kept weak and low frequency: at 0.55 and 103 cycles this perturbation was
    // steeper than the terrain it sat on and turned the plain into gravel.
    float fine = vnoise(w*19.0)*0.62 + vnoise(w*41.0)*0.38;
    float bare = 1.0 - smoothstep(0.05, 0.35, max(nt.r, max(nt.g, nt.b)));
    n = normalize(n + vec3(fine - 0.5, 0.0, vnoise(w*36.0 + 11.3) - 0.5)*0.085*bare);

    float activity = max(nt.r, max(nt.g, nt.b));
    float open = smoothstep(0.03, 0.45, activity);

    // Crust: cold, dark, faintly tinted by whoever used to own this ground.
    vec3 scar = scarW;
    vec3 rock = vec3(0.034, 0.032, 0.040)*(0.34 + 1.30*nt.a);
    rock *= 1.0 + fine*0.5;
    rock += (colonyColour(0)*scar.r + colonyColour(1)*scar.g + colonyColour(2)*scar.b)*0.055;

    // Key light: low and grazing, so the relief reads as stone.
    vec3 key = normalize(vec3(0.52, 0.46, -0.72));
    float ndl = max(dot(n, key), 0.0);
    float sh = 1.0;
    for(int i = 1; i < 7; i++){
      vec3 sp = hitP + key*(float(i)*0.042) + vec3(0.0, 0.004, 0.0);
      sh = min(sh, 1.0 - smoothstep(0.0, 0.016, heightFast(sp.xz, N) - sp.y));
    }
    col += rock*(0.115 + 0.62*ndl*sh)*vec3(0.70, 0.77, 1.0);

    // Cavity darkening from the baked height: the inside of a fissure is
    // shadowed before the emission fills it, which is what gives the crack a lip
    // instead of letting the glow sit flat on the surface.
    float cav = smoothstep(-1.15, 0.30, mm.a/max(relief, 1e-3));
    col *= mix(0.34, 1.0, cav);

    // Emission out of the fissure, with pulses running down the flow direction.
    float hd = min(length(w - heartPos(0, t)), min(length(w - heartPos(1, t)), length(w - heartPos(2, t))));
    float phase = hd*54.0 + a1.g*7.0 + t*3.60;
    float pulseK = 0.26 + 0.98*pow(0.5 + 0.5*sin(phase), 2.6);
    pulseK = mix(1.0, pulseK, pulse);
    // Emission from the smoothed field, sharpened with a power curve rather
    // than taken from the raw one-cell channel. Powering a smooth field tightens
    // the core without reintroducing the diagonal-line reconstruction artifact.
    float glLuma = clamp(max(gl.x, max(gl.y, gl.z)), 0.0, 1.0);
    vec3 emis = gl*(0.62 + 3.10*glLuma*glLuma)*pulseK*emissive*2.60;
    col += emis;

    // Bounce: the network is effectively the only light source in the piece, so
    // the crust beside a vein has to be lit by it or the whole thing reads flat.
    col += rock*gl*6.5*(0.35 + 0.65*ndl);

    // Distance fog toward black, so the far rim of the plain is never a visible
    // edge.
    col *= exp(-max(tHit - 1.4, 0.0)*0.42);
  } else {
    // Beyond the simulated square the plain carries on as bare stone. Without
    // this the domain reads as a floating slab with a hard rim against the sky,
    // which is the one thing that would give away that this is a field on a
    // grid. It is a plane rather than more marching because everything out here
    // is past the fog's half distance anyway.
    float tp = rd.y < -1e-4 ? (relief*0.45 - ro.y)/rd.y : -1.0;
    if(tp > 0.0){
      vec3 p = ro + rd*tp;
      float fine = vnoise(p.xz*19.0)*0.62 + vnoise(p.xz*41.0)*0.38;
      vec3 rock = vec3(0.034, 0.032, 0.040)*(1.05 + 0.5*fine);
      vec3 key = normalize(vec3(0.52, 0.46, -0.72));
      vec3 n = normalize(vec3((fine - 0.5)*0.16, 1.0, (vnoise(p.xz*36.0 + 11.3) - 0.5)*0.16));
      col = rock*(0.115 + 0.62*max(dot(n, key), 0.0))*vec3(0.70, 0.77, 1.0);
      col *= exp(-max(tp - 1.4, 0.0)*0.42);
      tEnd = min(tEnd, tp);
    } else {
      // Sky: almost nothing, a faint cool gradient so the horizon exists.
      float sky = max(rd.y, 0.0);
      col = mix(vec3(0.014, 0.016, 0.026), vec3(0.004, 0.005, 0.011), smoothstep(0.0, 0.5, sky));
    }
  }

  // Haze. Light leaving the plain and hanging in the air above it: this is the
  // layer that makes the piece read as a lit scene rather than as a field
  // visualisation, because it is the only thing in frame that is not a surface.
  float hs0 = span ? t0 : 0.0;
  float hs1 = min(tEnd, span ? t1 : 4.5);
  if(hs1 > hs0){
    float hstep = (hs1 - hs0)/8.0;
    float jit = hash12(gl_FragCoord.xy*0.37);
    vec3 hzc = vec3(0.0);
    for(int i = 0; i < 8; i++){
      float ht = hs0 + hstep*(float(i) + jit);
      vec3 p = ro + rd*ht;
      if(abs(p.x) < 1.0 && abs(p.z) < 1.0){
        vec4 g = fetchTexel(glow, clamp(floor(worldToGrid(p.xz, N)), vec2(0.0), vec2(N - 1.0)));
        float gm = max(g.x, max(g.y, g.z));
        float fall = exp(-max(p.y + relief*0.6, 0.0)*13.0);
        hzc += (g.rgb/max(gm, 1e-4))*g.a*fall*exp(-max(ht - 1.4, 0.0)*0.45);
      }
    }
    col += hzc*hstep*haze*1.55;
  }

  float lum = max(col.x, max(col.y, col.z));
  gl_FragColor = vec4(col, smoothstep(0.55, 1.6, lum));
}

// ---------------------------------------------------------------------------
// PASS 6 - composite: bloom, tonemap, grain
// ---------------------------------------------------------------------------

vec3 tonemap(vec3 x){
  // ACES-ish. Deep primaries matter here: the tonemap plus gamma desaturates,
  // and an emissive core that clips to white needs saturated colour in the
  // shoulder or the whole piece bleaches.
  return clamp((x*(2.51*x + 0.03))/(x*(2.43*x + 0.59) + 0.14), 0.0, 1.0);
}

void mainComposite(){
  vec2 uv = gl_FragCoord.xy/RENDERSIZE;
  vec4 s = texture2D(scene, uv);
  vec3 col = s.rgb;

  // Two-ring bloom gathered straight off the scene target, which keeps the
  // whole effect inside one pass and saves a texture unit - WebGL1 only
  // guarantees eight, and this shader already binds six.
  vec3 b = vec3(0.0);
  float wsum = 0.0;
  for(int i = 0; i < 22; i++){
    float fi = float(i);
    float ring = fi < 11.0 ? 0.014 : 0.052;
    float a = fi*2.39996 + (fi < 11.0 ? 0.0 : 0.7);
    vec2 off = vec2(cos(a), sin(a))*ring;
    off.x *= RENDERSIZE.y/RENDERSIZE.x;
    vec4 sp = texture2D(scene, clamp(uv + off, vec2(0.001), vec2(0.999)));
    float wgt = fi < 11.0 ? 1.0 : 0.55;
    b += sp.rgb*sp.a*wgt;
    wsum += wgt;
  }
  col += (b/wsum)*bloom*1.5;

  col = tonemap(col*1.15);

  float vig = 1.0 - 0.30*dot(uv - 0.5, uv - 0.5)*2.2;
  col *= clamp(vig, 0.0, 1.0);

  col = pow(col, vec3(1.0/2.2));

  // Dither before quantisation: large smooth gradients over near-black are
  // exactly where 8-bit banding shows, and this piece is mostly that.
  float grain = hash12(gl_FragCoord.xy*1.7 + 0.5) - 0.5;
  col += grain/255.0;

  gl_FragColor = vec4(col, 1.0);
}


void main(){
  if(PASSINDEX == 0) mainAcc1();
  else if(PASSINDEX == 1) mainAcc();
  else if(PASSINDEX == 2) mainNet();
  else if(PASSINDEX == 3) mainMem();
  else if(PASSINDEX == 4) mainGlow();
  else if(PASSINDEX == 5) mainScene();
  else mainComposite();
}