/*{
  "DESCRIPTION": "The Poo Rooms. The same infinite tiled natatorium, flooded to the ankles in something opaque. Structurally this is not a recolour of the Poolrooms, because the fluid stops being transparent and transparency was carrying most of that shader: with no transmission there is no refracted ray, so the refraction march, the Beer-Lambert absorption and the in-scattered turquoise all delete outright - and the caustics go with them, since a caustic is the Jacobian of the refracted map and there is no refracted map. That freed budget is what pays for the flies. Three things then have to be rebuilt rather than retuned. FIRST, VISCOSITY IS A LOW-PASS FILTER. Viscous damping goes as nu times k squared, so a thick fluid loses its short wavelengths first and keeps only long slow undulations; the four dispersive wave trains of the water version collapse to two, running at a fraction of the speed, and the surface is nearly static. Turn viscosity down and the top octaves come back and it starts behaving like water again, which is the honest way to expose that slider. SECOND, THE AMPLITUDE CROSSES THE PIXEL. The water was 13 mm of ripple - under a pixel at any distance past arm's length - which is the entire reason it could be solved as one plane hit plus a Newton step, with everything the viewer reacts to living in the normal. Lumps are 70 mm and they breach: the displacement is now several pixels, so it has to be traced, and the shader falls back to exactly what Seascape does for its 0.6 m waves - bracket the ray between a plane guaranteed above the surface and one guaranteed below, walk it with a step bounded by the field's slope, then refine the crossing by secant. Amplitude relative to a pixel, not taste, is what decides between the two methods. THIRD, IT IS NON-NEWTONIAN, so it holds shape: a crust field breaks the surface into plates that do NOT move with the undulation underneath, cracked along noise ridges, drier and matter where it has set. Bubbles surface and pop on a beat, and they are a shading feature rather than a geometric one on purpose - a 3 cm dome is smaller than the tracer's step and would simply be walked over, so it is applied as a normal perturbation and a ring in the albedo where it can never be missed. Reflection is a broad rough sheen rather than a march: an opaque surface at 0.035 Fresnel with this much roughness returns almost nothing legible, so a fifty-step reflection ray would be fifty map calls spent on mud, and a single wide specular lobe off the light gather is both cheaper and more correct. The flies are the other half. They are not geometry and never touch the distance field: each is a point tested against the primary ray by closest approach, which costs a dot product, and two things make them read as insects rather than as dots. They DART - the path is piecewise linear between hashed waypoints on a per-fly beat, so heading changes are discontinuous the way a real fly's are, where any smooth orbit reads as a planet. And they are MOTION BLURRED, because a fly covers several body lengths in one frame's exposure: the closest-approach distance is measured in an ellipse stretched along the velocity projected into the image, which turns each one into a streak for one extra normalise. Sub-pixel flies are the aliasing case, and fading them out with distance would be wrong because a swarm at ten metres is still visible - so a fly whose radius falls under a pixel is widened to a pixel and dimmed by the square of the widening, which conserves what it contributes and stops the whole swarm scintillating as the flies cross pixel centres. Clusters sit over the surface or up under a lit coffer, and the swarm is composited back to front so the near ones win. The tile, the arches, the plan, the coffers, the footprint filtering and the shafts are all as before; the accent course at the waterline is now buried under the tide line. Forty-one controls.",
  "CREDIT": "CC0. Plan generation, footprint filtering and the light gather follow the CC0 BACKROOMS shader.",
  "CATEGORIES": ["generator", "raymarch", "liminal", "horror", "gross", "flies", "procedural", "3d"],
  "INPUTS": [
    { "NAME": "uSpeed",      "TYPE": "float", "DEFAULT": 2.0,   "MIN": 0.0,   "MAX": 5.0 },
    { "NAME": "uLens",       "TYPE": "float", "DEFAULT": 1.50,  "MIN": 0.7,   "MAX": 3.0 },
    { "NAME": "uYaw",        "TYPE": "float", "DEFAULT": 0.0,   "MIN": -1.57, "MAX": 1.57 },
    { "NAME": "uPitch",      "TYPE": "float", "DEFAULT": -0.08, "MIN": -0.6,  "MAX": 0.6 },
    { "NAME": "uWander",     "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 1.6 },
    { "NAME": "uTurn",       "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uBob",        "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 2.0 },
    { "NAME": "uCell",       "TYPE": "float", "DEFAULT": 5.5,   "MIN": 3.5,   "MAX": 9.0 },
    { "NAME": "uDensity",    "TYPE": "float", "DEFAULT": 0.62,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uRooms",      "TYPE": "float", "DEFAULT": 0.25,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uArch",       "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uArchW",      "TYPE": "float", "DEFAULT": 0.66,  "MIN": 0.40,  "MAX": 1.10 },
    { "NAME": "uHeight",     "TYPE": "float", "DEFAULT": 3.60,  "MIN": 2.4,   "MAX": 6.0 },
    { "NAME": "uLedge",      "TYPE": "float", "DEFAULT": 0.26,  "MIN": 0.0,   "MAX": 0.60 },
    { "NAME": "uLevel",      "TYPE": "float", "DEFAULT": 0.30,  "MIN": 0.0,   "MAX": 1.00 },
    { "NAME": "uVisc",       "TYPE": "float", "DEFAULT": 0.80,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uLumps",      "TYPE": "float", "DEFAULT": 0.85,  "MIN": 0.0,   "MAX": 1.5 },
    { "NAME": "uCrust",      "TYPE": "float", "DEFAULT": 0.55,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uBubbles",    "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uSheen",      "TYPE": "float", "DEFAULT": 0.70,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uFilm",       "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uFlies",      "TYPE": "float", "DEFAULT": 0.65,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uSwarm",      "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uBlur",       "TYPE": "float", "DEFAULT": 0.55,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uLightSpace", "TYPE": "float", "DEFAULT": 3.05,  "MIN": 1.83,  "MAX": 7.3 },
    { "NAME": "uLightOn",    "TYPE": "float", "DEFAULT": 0.72,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uBright",     "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uZones",      "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uHaze",       "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.5 },
    { "NAME": "uFog",        "TYPE": "float", "DEFAULT": 0.045, "MIN": 0.0,   "MAX": 0.2 },
    { "NAME": "uPalette",    "TYPE": "float", "DEFAULT": 0.0,   "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uTileScale",  "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.5,   "MAX": 2.5 },
    { "NAME": "uAccent",     "TYPE": "float", "DEFAULT": 0.70,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uStain",      "TYPE": "float", "DEFAULT": 0.75,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uGrime",      "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uDamp",       "TYPE": "float", "DEFAULT": 0.50,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uAO",         "TYPE": "float", "DEFAULT": 0.75,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uShadow",     "TYPE": "float", "DEFAULT": 0.80,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uGloss",      "TYPE": "float", "DEFAULT": 0.75,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uVHS",        "TYPE": "float", "DEFAULT": 0.40,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uScan",       "TYPE": "float", "DEFAULT": 0.25,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uGrain",      "TYPE": "float", "DEFAULT": 0.24,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uVignette",   "TYPE": "float", "DEFAULT": 0.50,  "MIN": 0.0,   "MAX": 1.0 }
  ]
}*/

// THE POO ROOMS - ISF / GLSL ES 1.00 (gl_FragColor, TIME, RENDERSIZE).
//
// The Poolrooms plan, flooded with something opaque. The fluid is traced as a
// height field rather than solved analytically, because its lumps are seventy
// millimetres and no longer hide under a pixel; the flies are analytic points
// tested against the ray, because they are five millimetres and always will.

#define PI  3.14159265359
#define FAR 44.0

#define CTILE 0.6096
#define WTILE 0.1524
#define FTILE 0.2286

float gLG;
vec2  gLOF;
#define PANEL_C(lid) ((lid + 0.5)*gLG + gLOF)
#define PANEL_E      (vec2(CTILE, CTILE))

float gSL;      // still level of the sludge
float gWT;      // fluid clock
float gAMP;     // amplitude bound of the surface, for the trace bracket
vec2  gCam;
float gPX;

vec3 gTile, gGrout, gAcc, gPlas, gLite, gMuck, gDry, gStain;

// ---------------------------------------------------------------
// hashes and value noise
// ---------------------------------------------------------------
float hash11(float p){ p = fract(p*0.1031); p *= p + 33.33; p *= p + p; return fract(p); }

float hash12(vec2 p){
    vec3 q = fract(vec3(p.xyx)*0.1031);
    q += dot(q, q.yzx + 33.33);
    return fract((q.x + q.y)*q.z);
}

vec4 hash42(vec2 p){
    vec4 q = fract(vec4(p.xyxy)*vec4(0.1031, 0.1030, 0.0973, 0.1099));
    q += dot(q, q.wzxy + 33.33);
    return fract((q.xxyz + q.yzzw)*q.zywx);
}

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

// every octave rotated as well as scaled, or four sets of axis-aligned cell
// edges stack into a visible quilt on any large flat surface
float fbm(vec2 p){
    mat2 R = mat2(0.8040, 0.5946, -0.5946, 0.8040);
    float s = 0.0, a = 0.5;
    for (int i = 0; i < 4; i++){
        s += a*vnoise(p);
        p = R*p*2.07 + 19.3;
        a *= 0.5;
    }
    return s;
}

float fbm2(vec2 p){
    return 0.62*vnoise(p) + 0.38*vnoise(p*2.13 + 7.3);
}

// ---------------------------------------------------------------
// primitives and the footprint filter
// ---------------------------------------------------------------
float sdBox2(vec2 p, vec2 b){ vec2 q = abs(p) - b; return length(max(q, 0.0)) + min(max(q.x, q.y), 0.0); }

// Box-filter coverage of |x| < r by a pixel w wide. A grout line is 3 mm, so
// past a couple of metres - and immediately on any oblique wall - it is
// thinner than a pixel, and a hard-edged test keeps drawing it at full
// strength wherever it lands, which beats against the pixel grid as moire.
float band(float x, float r, float w){
    float h = 0.5*max(w, 1e-5);
    return clamp((min(x + h, r) - max(x - h, -r))/(2.0*h), 0.0, 1.0);
}

float fill(float d, float w){ return clamp(0.5 - d/max(w, 1e-5), 0.0, 1.0); }

// The transition has to be WIDE: everything keys off the footprint, and the
// footprint depends on dot(n, rd), so a narrow fade paints conics about the
// normal - concentric rings across every oblique wall.
float mip(float feature, float foot){ return smoothstep(1.10, 0.10, foot/feature); }

vec3 pick4(vec3 a, vec3 b, vec3 c, vec3 d, float f){
    float i = clamp(f, 0.0, 3.0);
    float n = floor(i), k = i - n;
    vec3 lo = a;
    lo = mix(lo, b, step(0.5, n)); lo = mix(lo, c, step(1.5, n)); lo = mix(lo, d, step(2.5, n));
    vec3 hi = b;
    hi = mix(hi, c, step(0.5, n)); hi = mix(hi, d, step(1.5, n));
    return mix(lo, hi, k);
}

// ---------------------------------------------------------------
// the plan
// ---------------------------------------------------------------
float pathX(float z){
    return uWander*(2.6*sin(z*0.0290) + 1.35*sin(z*0.0161 + 1.7))
         + uTurn*(11.0*sin(z*0.0473 + 2.1) + 7.0*sin(z*0.0281 + 5.3)
                + 4.5*sin(z*0.0917 + 0.7));
}

float pathSlope(float z){
    return uWander*(0.0754*cos(z*0.0290) + 0.0217*cos(z*0.0161 + 1.7))
         + uTurn*(0.5203*cos(z*0.0473 + 2.1) + 0.1967*cos(z*0.0281 + 5.3)
                + 0.4127*cos(z*0.0917 + 0.7));
}

// The arch as a CROSS-SECTION, not as a height-dependent width. The obvious
// way - half-width = W*sqrt(1 - u*u) above the springing, subtracted as a box
// - has infinite slope at the crown, so the field's Lipschitz constant blows
// up exactly there and the tracer steps straight through the arch head. Built
// in the plane of the wall as a jamb box unioned with a disc and then extruded
// through it, the same silhouette has a bounded gradient.
float archSection(float u, float v, float y, float thick){
    float W = uArchW;
    float S = 1.02 + 0.25*W;
    float jamb = sdBox2(vec2(u, y - 0.5*S), vec2(W, 0.5*S));
    float head = length(vec2(u, y - S)) - W;
    return max(min(jamb, head), abs(v) - thick);
}

vec2 maze2(vec2 q, float y){
    float C = uCell;
    vec2 base = floor(q/C);
    float dW = C, dD = C;
    for (int j = 0; j < 3; j++){
        float rz = (base.y + float(j) - 0.5)*C;
        float px = pathX(rz);
        float sl = pathSlope(rz);
        // a fixed-width horizontal band is a corridor that narrows as the path
        // steepens; sqrt(1 + slope^2) keeps it constant measured ACROSS itself
        float pw = 0.80*C*sqrt(1.0 + sl*sl);
        for (int i = 0; i < 3; i++){
            vec2 id = base + vec2(float(i) - 1.0, float(j) - 1.0);
            vec2 c  = (id + 0.5)*C;
            vec4 h  = hash42(id);
            vec4 g  = hash42(id + 57.31);

            float on = step(h.x, uDensity)*step(pw, abs(c.x - px));

            float horiz = step(0.5, h.y);
            vec2 dirL = mix(vec2(0.0, 1.0), vec2(1.0, 0.0), horiz);
            vec2 dirT = mix(vec2(1.0, 0.0), vec2(0.0, 1.0), horiz);

            float pil = step(0.90, g.x);
            float len = mix((0.66 + 0.34*h.z)*0.5*C, 0.32, pil);
            float th  = mix(0.105, 0.32, pil);
            vec2  ctr = c + dirL*((h.w - 0.5)*(C - 2.0*len))
                          + dirT*((g.y - 0.5)*0.70*(C - 2.0*th));

            float isRoom = step(fract(h.y*31.7), uRooms);
            float rr = 0.5*C - 0.22;
            float ring = max(sdBox2(q - c, vec2(rr)), -sdBox2(q - c, vec2(rr - 0.10)));
            float slab = sdBox2(q - ctr, dirL*len + dirT*th);
            dW = min(dW, mix(slab, ring, isRoom) + (1.0 - on)*90.0);

            // a chamber always gets its arch, or it is a sealed box nobody can
            // ever be inside
            float dr = on*mix(step(g.z, uArch)*(1.0 - pil), 1.0, isRoom);
            vec2 dc = mix(ctr + dirL*((g.w - 0.5)*2.0*max(len - 1.15, 0.0)),
                          c + dirT*rr, isRoom);
            vec2 rel = q - dc;
            float u = dot(rel, dirL), v = dot(rel, dirT);
            dD = min(dD, archSection(u, v, y, mix(th, 0.10, isRoom) + 0.55) + (1.0 - dr)*90.0);

            if (dW < -0.5*C) break;
        }
        if (dW < -0.5*C) break;
    }
    return vec2(min(dW, C), min(dD, C));
}

// .x = distance, .y = material (1 floor, 2 ceiling, 3 wall, 4 ledge)
vec2 map(vec3 p){
    vec2 r = vec2(p.y, 1.0);
    float dC = uHeight - p.y;
    if (dC < r.x) r = vec2(dC, 2.0);

    vec2 mz = maze2(p.xz, p.y);
    float dW = max(mz.x, -mz.y);
    if (dW < r.x) r = vec2(dW, 3.0);

    // the ledge is defined relative to the LEVEL, not to the floor, because
    // its whole job is to be half in it
    if (uLedge > 0.02){
        float dl = max(max(mz.x - uLedge, -mz.y), p.y - (gSL + 0.22));
        if (dl < r.x) r = vec2(dl, 4.0);
    }
    return r;
}

vec3 calcNormal(vec3 p){
    vec2 e = vec2(1.0, -1.0)*0.0018;
    return normalize(e.xyy*map(p + e.xyy).x + e.yyx*map(p + e.yyx).x
                   + e.yxy*map(p + e.yxy).x + e.xxx*map(p + e.xxx).x);
}

float calcAO(vec3 p, vec3 n){
    float o = 0.0, s = 1.0;
    for (int i = 0; i < 5; i++){
        float h = 0.02 + 0.16*float(i);
        o += (h - map(p + n*h).x)*s;
        s *= 0.72;
    }
    return clamp(1.0 - 1.05*o, 0.0, 1.0);
}

float shadowRay(vec3 ro, vec3 rd, float maxt){
    float res = 1.0, t = 0.16, prev = 1e12;
    for (int i = 0; i < 24; i++){
        float h = map(ro + rd*t).x;
        float y = h*h/(2.0*prev);
        float d = sqrt(max(h*h - y*y, 0.0));
        res = min(res, d/max(0.0001, 7.0*(t - y)));
        prev = h;
        t += clamp(h, 0.05, 0.55);
        if (res < 0.02 || t > maxt) break;
    }
    return clamp(res, 0.0, 1.0);
}

// ---------------------------------------------------------------
// the coffers
// ---------------------------------------------------------------
float lightZone(vec2 xz){ return smoothstep(0.20, 0.80, vnoise(xz*0.026)); }

float panelOn(vec2 lid){
    vec4 h = hash42(lid + 11.7);
    float zn = lightZone(PANEL_C(lid));
    float on = step(h.x, uLightOn*mix(1.0, 0.05 + 1.90*zn, uZones));
    return on*(0.965 + 0.035*sin(TIME*47.0 + h.y*97.0 + h.z*13.0));
}

// Nine nearest coffers as RECTANGLES, plus one direction to fire a single
// shadow ray along - the WEIGHTED MEAN of the nine, not the strongest. The
// strongest changes as you slide along a wall, the loci where two tie are
// smooth curves, and the shadow term steps across them, painting the wall in
// concentric arcs. A mean has no seams.
float lighting(vec3 p, vec3 n, out vec3 keyDir, out float keyDist){
    float LG = gLG;
    vec2 lb = floor(p.xz/LG);
    float sum = 0.0, accD = 0.0;
    vec3 acc = vec3(0.0);
    for (int j = 0; j < 3; j++){
        for (int i = 0; i < 3; i++){
            vec2 lid = lb + vec2(float(i) - 1.0, float(j) - 1.0);
            vec2 lc  = PANEL_C(lid);
            float on = panelOn(lid);
            vec3 pl = vec3(clamp(p.x, lc.x - PANEL_E.x, lc.x + PANEL_E.x), uHeight - 0.04,
                           clamp(p.z, lc.y - PANEL_E.y, lc.y + PANEL_E.y));
            vec3 L = pl - p;
            float dl = max(length(L), 0.08);
            L /= dl;
            // a fixture must be worth NOTHING by the time the 3x3 window drops
            // it, or the sum steps at every cell edge
            vec2 ho = max(abs(p.xz - lc) - PANEL_E, 0.0);
            float w = on*max(dot(n, L), 0.0)*(0.30 + 0.70*max(L.y, 0.0))
                    * smoothstep(1.45*LG, 1.02*LG, length(ho))/(1.0 + 0.30*dl*dl);
            sum  += w;
            acc  += L*w;
            accD += dl*w;
        }
    }
    keyDir  = normalize(acc + n*0.001 + vec3(0.0, 0.0001, 0.0));
    keyDist = accD/max(sum, 0.0001);
    return sum;
}

// ---------------------------------------------------------------
// the surface
// ---------------------------------------------------------------
// VISCOSITY IS A LOW-PASS FILTER. Viscous damping goes as nu*k^2, so a thick
// fluid loses its short wavelengths first: the four dispersive wave trains of
// the water version collapse to two long slow ones, and at full viscosity the
// surface barely moves at all. Take uVisc down and the octaves come back and
// it starts behaving like water again, which is the honest way to expose it.
float undulate(vec2 p, float t){
    mat2 R = mat2(0.8763, 0.4818, -0.4818, 0.8763);
    float sp = mix(2.4, 0.30, uVisc);          // and it creeps rather than runs
    float h = 0.0, a = 1.0, f = 1.0, s = 0.0;
    for (int i = 0; i < 3; i++){
        // the top octave is damped out entirely at high viscosity
        float vis = (i < 2) ? 1.0 : (1.0 - uVisc);
        vec2 q = p*f;
        float ph = 6.2831*vnoise(q*0.33 + float(i)*7.1 + t*0.02);
        h += a*vis*sin(q.x + 0.55*q.y + ph + t*sp*inversesqrt(f));
        s += a;
        p = R*p; a *= 0.55; f *= 2.15;
    }
    return h/s;
}

// The lumps: solids that have not dispersed, sitting proud of the level. A
// one-sided field - max against zero - because a lump is an addition to the
// surface and never a hole in it, which is also what stops the level reading
// as a noisy plane rather than as a thing with objects in it.
float lumps(vec2 p){
    if (uLumps < 0.005) return 0.0;
    float a = fbm2(p*1.30 + 4.0);
    float b = fbm2(p*2.90 + 21.0);
    float m = smoothstep(0.52, 0.86, a) + 0.55*smoothstep(0.62, 0.92, b);
    return uLumps*0.070*m;
}

// The crust. Non-Newtonian, so it HOLDS SHAPE: these plates do not move with
// the undulation underneath them, which is the whole difference between a
// liquid and this. Ridged noise gives the cracks between them.
float crust(vec2 p){
    if (uCrust < 0.005) return 0.0;
    float c = smoothstep(0.46, 0.72, fbm(p*0.85 + 13.0));
    float crack = smoothstep(0.030, 0.0, abs(fbm(p*1.9 + 40.0) - 0.5));
    return uCrust*c*(1.0 - 0.85*crack);
}

float surfaceH(vec2 xz){
    float t = gWT;
    float h = gSL + uLumps*0.0;                      // level
    float cr = crust(xz);
    // where it has set, the undulation is suppressed and the plate sits a
    // little proud of the liquid it is floating on
    h += mix(0.022, 0.004, uVisc)*undulate(xz*2.2, t)*(1.0 - 0.85*cr);
    h += lumps(xz);
    h += 0.012*cr;
    // a wake, because wading through it should disturb it - slower and
    // shorter-lived than water's, because it is thick
    float r = length(xz - gCam);
    h += mix(0.014, 0.004, uVisc)*sin(r*7.0 - t*2.2)*exp(-r*1.1)
       * smoothstep(0.05, 0.8, r)*clamp(uSpeed*0.4, 0.0, 1.2);
    return h;
}

// TRACED, not solved.
//
// The water version intersected a plane and took one Newton step, and it was
// right to: 13 mm of ripple is under a pixel past arm's length, so the
// displacement is invisible and only the normal matters. Lumps are 70 mm and
// they breach. Once the amplitude crosses the pixel the analytic shortcut is
// simply wrong, and the correct method is the one Seascape uses for its 0.6 m
// waves: bracket the ray between a plane guaranteed above the surface and one
// guaranteed below, walk it with a step bounded by the field's own slope, then
// refine the crossing by secant.
//
// The step is f/(|rd.y| + slopeBound) rather than a fixed increment because at
// grazing incidence a fixed increment either misses every lump or costs two
// hundred samples, and the slope bound is what makes the bounded version safe.
float traceSurface(vec3 ro, vec3 rd, float tmax, out float tHit){
    tHit = tmax;
    if (rd.y > -0.004 || ro.y < gSL) return 0.0;
    float t0 = (gSL + gAMP - ro.y)/rd.y;
    t0 = max(t0, 0.0);
    float t1 = min((gSL - 0.02 - ro.y)/rd.y, tmax);
    if (t1 <= t0) return 0.0;

    float t = t0, pf = 1.0, pt = t0;
    for (int i = 0; i < 26; i++){
        vec3 q = ro + rd*t;
        float f = q.y - surfaceH(q.xz);
        if (f < 0.0){
            // secant onto the crossing: two more evaluations buy an exact
            // silhouette on every lump, which matters because the lumps ARE
            // the silhouette
            float ta = pt, tb = t, fa = pf, fb = f;
            for (int k = 0; k < 3; k++){
                float tc = mix(ta, tb, clamp(fa/max(fa - fb, 1e-4), 0.0, 1.0));
                vec3 qc = ro + rd*tc;
                float fc = qc.y - surfaceH(qc.xz);
                if (fc < 0.0){ tb = tc; fb = fc; } else { ta = tc; fa = fc; }
            }
            tHit = mix(ta, tb, clamp(fa/max(fa - fb, 1e-4), 0.0, 1.0));
            return 1.0;
        }
        pf = f; pt = t;
        t += clamp(f*0.80/(abs(rd.y) + 0.55), 0.03, 1.40);
        if (t > t1) break;
    }
    return 0.0;
}

vec3 surfaceN(vec2 xz, float e){
    float h  = surfaceH(xz);
    float hx = surfaceH(xz + vec2(e, 0.0));
    float hz = surfaceH(xz + vec2(0.0, e));
    return normalize(vec3(h - hx, e, h - hz));
}

// ---------------------------------------------------------------
// tile, box-filtered
// ---------------------------------------------------------------
// A periodic grid seen edge-on averages to its DUTY CYCLE, not to nothing.
// Letting band() run to zero bleaches every distant wall to plain white;
// blending toward the average coverage as the footprint outgrows the tile is
// what makes a far wall read as tile you cannot resolve.
float grout(vec2 uv, float T, float foot, out vec2 tid){
    vec2 g = uv/T;
    tid = floor(g);
    vec2 d = abs(fract(g) - 0.5);
    float e = 0.5 - max(d.x, d.y);
    float w = foot/T;
    float gr = 0.034;
    return mix(band(e, gr, w), 4.0*gr, smoothstep(0.06, 0.40, w));
}

// ---------------------------------------------------------------
// shading. qual 1 = full, 0 = the cheap version (unused here but kept so the
// reflection path can be reinstated if anyone wants it back)
// ---------------------------------------------------------------
vec3 shade(vec3 p, vec3 rd, vec3 n, float mat, float foot, float qual, out vec3 emis){
    emis = vec3(0.0);
    float det  = mip(0.008, foot);
    float detT = mip(0.030, foot);
    float detM = mip(0.45,  foot);
    float above = p.y - gSL;
    float sub  = smoothstep(0.012, -0.012, above);

    vec3 alb = gTile;
    float gloss = 0.0;
    vec2 tid;

    if (mat < 1.5){
        float T = FTILE*uTileScale;
        float gm = grout(p.xz, T, foot, tid);
        vec4 th = hash42(tid + 4.3);
        alb = gTile*(0.90 + 0.14*th.x)*(1.0 + 0.10*(vnoise(p.xz*40.0) - 0.5)*det);
        alb = mix(alb, gGrout*(0.85 + 0.3*th.y), gm);
        alb *= mix(1.0, 0.66 + 0.40*fbm(p.xz*0.55 + 9.0), uGrime*detM);
        gloss = 0.35;
    } else if (mat < 2.5){
        vec2 lid = floor(p.xz/gLG);
        vec2 dp  = abs(p.xz - PANEL_C(lid)) - PANEL_E;
        if (max(dp.x, dp.y) < 0.0){
            float on = panelOn(lid);
            vec2 pg = fract((p.xz - PANEL_C(lid))/(CTILE*0.25));
            float pr = smoothstep(0.0, 0.16, min(min(pg.x, 1.0 - pg.x), min(pg.y, 1.0 - pg.y)));
            pr = mix(1.0, 0.74 + 0.38*pr, det);
            // and the diffuser is full of forty years of dead insects
            float dead = smoothstep(0.62, 0.90, fbm(p.xz*7.0 + 3.0));
            emis = gLite*(1.5 + 6.5*uBright)*on*pr*(1.0 - 0.45*dead*detT);
            alb  = gLite*0.55;
        } else {
            float gm = grout(p.xz, CTILE, foot, tid);
            vec4 th = hash42(tid + 3.7);
            alb = gPlas*(0.88 + 0.18*th.x)*(1.0 + 0.16*(vnoise(p.xz*90.0) - 0.5)*det);
            alb = mix(alb, gPlas*0.62, gm*0.8);
            float ring = smoothstep(0.54, 0.82, fbm(p.xz*1.3 + th.yz*11.0))*uDamp*step(th.w, 0.5)*detM;
            alb = mix(alb, gStain*0.7, ring*0.75);
            alb *= 1.0 - 0.50*smoothstep(-0.05, 0.0, max(dp.x, dp.y));
        }
    } else {
        float T = WTILE*uTileScale;
        float top = step(0.55, abs(n.y));
        vec2 uv;
        if (top > 0.5)           uv = p.xz;
        else if (abs(n.x) > 0.5) uv = vec2(p.z, p.y);
        else                     uv = vec2(p.x, p.y);
        if (mat > 3.5 && top > 0.5) T = FTILE*uTileScale;

        float gm = grout(uv, T, foot, tid);
        vec4 th = hash42(tid + 17.1);
        alb = gTile*(0.92 + 0.13*th.x);
        alb *= 1.0 + 0.13*(vnoise(uv*220.0) - 0.5)*mip(0.010, foot)
                   + 0.07*(vnoise(uv*55.0)  - 0.5)*mip(0.050, foot);
        alb *= mix(1.0, 0.55, step(0.988, th.z)*detT);

        // the accent course is still there, it is just under the tide line now
        if (top < 0.5 && uAccent > 0.01){
            float cov = band(uv.y - (gSL + 0.14), 0.152, foot);
            alb = mix(alb, gAcc*(0.85 + 0.3*th.y), cov*uAccent);
        }
        alb = mix(alb, gGrout*(0.88 + 0.24*th.y), gm);

        if (top < 0.5){
            // THE TIDE LINE. Two separate things and they read differently: a
            // hard band at the high mark where it stood and dried, and streaks
            // running down from it, which are what say the level has DROPPED.
            float gr = fbm(uv*1.25 + 7.0);
            float mark = gSL + 0.30 + 0.10*gr;
            float wall = smoothstep(mark + 0.04, mark - 0.10, uv.y);
            float run  = smoothstep(0.55, 0.95, vnoise(vec2(uv.x*22.0, uv.y*0.6)))
                       * smoothstep(mark + 0.85, mark, uv.y);
            float st = uStain*detM*max(wall*(0.55 + 0.45*gr), run*0.75);
            alb = mix(alb, gStain*(0.7 + 0.5*gr), clamp(st, 0.0, 0.92));
            alb *= mix(1.0, 0.66 + 0.42*gr, smoothstep(1.4, 0.0, above)*uGrime*detM);
        }
        gloss = 0.55;
    }

    // wet is not a colour, it is a darkening plus a gloss, and it climbs a
    // little way up everything the level touches
    float wet = sub + (1.0 - sub)*exp(-max(above, 0.0)*9.0);
    alb *= 1.0 - 0.42*wet;
    alb = mix(alb, gMuck*0.8, wet*0.55);
    gloss = mix(gloss, 0.92, wet);

    float ao = 1.0;
    if (qual > 0.5) ao = mix(1.0, calcAO(p, n), uAO);
    vec3 kd; float kdist;
    float lit = lighting(p, n, kd, kdist);
    float sh = 1.0;
    if (qual > 0.5 && uShadow > 0.01) sh = mix(1.0, shadowRay(p + n*0.06, kd, kdist - 0.20), uShadow);

    // No caustics anywhere in this shader. A caustic is the Jacobian of the
    // REFRACTED map, and an opaque fluid has no refracted map - so the light
    // that danced on the walls of the Poolrooms is simply gone, and the walls
    // above the level are lit by the coffers alone. That absence is most of
    // why the room reads as dead rather than as a pool.
    vec3 diff = gLite*lit*(0.52 + 0.48*sh)*(2.4 + 3.8*uBright);

    float dn = 0.5 - 0.5*n.y;
    float zdim = mix(1.0, 0.25 + 1.05*lightZone(p.xz), uZones);
    // the bounce off the level is dark and brown, where the pool's was bright
    // and cyan - one line, and it recolours the entire room
    vec3 amb = ao*zdim*(mix(gTile, gLite, 0.5)*(0.075 + 0.09*uBright)
                      + gMuck*gLite*(0.10 + 0.14*uBright)*dn);

    vec3 col = alb*(diff*ao*0.85 + amb);

    if (uGloss > 0.01 && lit > 0.001){
        vec3 hv = normalize(kd - rd);
        float sp = pow(max(dot(n, hv), 0.0), mix(28.0, 150.0, gloss));
        col += gLite*sp*lit*sh*gloss*uGloss*(0.35 + 0.5*uBright);
    }
    return col;
}

// ---------------------------------------------------------------
// the level itself
// ---------------------------------------------------------------
vec3 shadeSurface(vec3 p, vec3 rd, float dist, vec3 fogCol){
    float foot = dist*gPX/max(-rd.y, 0.02);
    float e = max(0.010, 0.55*foot);
    vec3 n = surfaceN(p.xz, e);
    float cr = crust(p.xz);

    // BUBBLES ARE SHADING, NOT GEOMETRY, on purpose. A 3 cm dome is smaller
    // than the tracer's step and would be walked straight over half the time;
    // applied as a normal perturbation and a ring in the albedo it cannot be
    // missed, and nobody can tell from outside which one it was.
    float bub = 0.0;
    if (uBubbles > 0.005){
        vec2 bg = p.xz*3.4;
        vec2 bid = floor(bg);
        vec2 bf = fract(bg) - 0.5;
        vec4 bh = hash42(bid + 5.9);
        float life = fract(gWT*0.22*(0.6 + bh.x) + bh.y);      // swell, then pop
        float rad = 0.16*smoothstep(0.0, 0.55, life)*smoothstep(1.0, 0.72, life)*(0.5 + bh.z);
        float d = length(bf - (bh.xy - 0.5)*0.5);
        bub = step(bh.w, uBubbles*0.55)*smoothstep(rad, rad*0.55, d);
        n = normalize(n + vec3(bf.x, 0.0, bf.y)*bub*2.4);
    }

    // Albedo. The wet fraction is dark and specular, the set fraction is
    // lighter, drier and matt, and the cracks between the plates are darker
    // than either because you are seeing down into it.
    float v = fbm(p.xz*1.7 + 2.0);
    vec3 alb = mix(gMuck*(0.70 + 0.60*v), gDry*(0.80 + 0.40*v), cr);
    alb *= 1.0 - 0.35*smoothstep(0.60, 0.95, lumps(p.xz)/max(uLumps*0.070, 1e-4));
    alb = mix(alb, alb*0.55, bub*0.4);

    float gloss = mix(0.95, 0.25, cr);

    vec3 kd; float kdist;
    float lit = lighting(p + vec3(0.0, 0.02, 0.0), n, kd, kdist);
    float sh = 1.0;
    if (uShadow > 0.01) sh = mix(1.0, shadowRay(p + vec3(0.0, 0.05, 0.0), kd, kdist - 0.20), uShadow);
    float zdim = mix(1.0, 0.25 + 1.05*lightZone(p.xz), uZones);

    vec3 col = alb*(gLite*lit*(0.55 + 0.45*sh)*(1.9 + 3.0*uBright)
                  + mix(gTile, gLite, 0.5)*(0.085 + 0.10*uBright)*zdim);

    // NO REFLECTION MARCH. An opaque surface at 0.035 Fresnel with this much
    // roughness returns nothing legible, so fifty-odd map calls would be spent
    // resolving a reflection nobody can read; one wide specular lobe off the
    // light gather is cheaper AND more correct. It is only worth marching a
    // reflection when the surface is smooth enough for the result to survive.
    float ndv = max(dot(-rd, n), 0.0);
    float F = 0.035 + 0.55*pow(1.0 - ndv, 5.0);
    vec3 hv = normalize(kd - rd);
    float sp = pow(max(dot(n, hv), 0.0), mix(12.0, 220.0, gloss*uSheen));
    col += gLite*sp*lit*sh*uSheen*(0.6 + 0.7*uBright)*(0.4 + 0.6*F);
    col += fogCol*F*1.4*(1.0 - cr);
    // the bright meniscus around a bubble, which is what makes it read as a
    // bubble rather than as a stain
    col += gLite*bub*smoothstep(0.35, 0.9, 1.0 - abs(dot(n, -rd)))*0.35*(0.4 + 0.6*uBright);

    // A thin biofilm at grazing angles. Iridescence needs the shift to depend
    // on the ANGLE, not on position, or it is just a coloured stain: the phase
    // runs on 1/ndv, which is the optical path through a film of fixed
    // thickness, so the bands sit in the surface and swim as the camera moves.
    if (uFilm > 0.01){
        float ph = 7.0/max(ndv, 0.12) + 3.0*v;
        vec3 ir = 0.5 + 0.5*cos(vec3(ph, ph + 2.1, ph + 4.2));
        col += ir*pow(1.0 - ndv, 3.0)*uFilm*0.14*(0.3 + 0.7*lit)*(1.0 - cr);
    }
    return col;
}

// ---------------------------------------------------------------
void main(){
    vec2 fragCoord = gl_FragCoord.xy;
    vec2 R = RENDERSIZE.xy;
    float T = TIME;
    float t = T*uSpeed;

    gTile  = pick4(vec3(0.88, 0.90, 0.86), vec3(0.86, 0.84, 0.74), vec3(0.72, 0.76, 0.78), vec3(0.84, 0.87, 0.72), uPalette);
    gGrout = pick4(vec3(0.58, 0.58, 0.53), vec3(0.57, 0.52, 0.43), vec3(0.46, 0.50, 0.52), vec3(0.54, 0.58, 0.46), uPalette);
    gAcc   = pick4(vec3(0.10, 0.30, 0.36), vec3(0.18, 0.28, 0.26), vec3(0.06, 0.16, 0.28), vec3(0.16, 0.33, 0.20), uPalette);
    gPlas  = pick4(vec3(0.84, 0.86, 0.82), vec3(0.84, 0.80, 0.70), vec3(0.68, 0.72, 0.76), vec3(0.80, 0.85, 0.71), uPalette);
    gLite  = pick4(vec3(0.95, 0.97, 0.94), vec3(1.00, 0.94, 0.78), vec3(0.76, 0.86, 0.98), vec3(0.88, 1.00, 0.82), uPalette);
    gMuck  = pick4(vec3(0.20, 0.12, 0.055), vec3(0.26, 0.16, 0.06), vec3(0.16, 0.13, 0.09), vec3(0.18, 0.17, 0.06), uPalette);
    gDry   = pick4(vec3(0.38, 0.27, 0.15), vec3(0.46, 0.34, 0.18), vec3(0.32, 0.29, 0.23), vec3(0.36, 0.34, 0.16), uPalette);
    gStain = pick4(vec3(0.34, 0.24, 0.13), vec3(0.42, 0.30, 0.15), vec3(0.28, 0.26, 0.21), vec3(0.32, 0.31, 0.14), uPalette);

    float ltn = max(3.0, floor(uLightSpace/CTILE + 0.5));
    gLG  = CTILE*ltn;
    // the coffer is two tiles across, so its centre must land on a tile
    // BOUNDARY: free when the lattice is even, a half-tile nudge when it is odd
    gLOF = vec2(0.5*CTILE*mod(ltn, 2.0));
    gSL  = uLevel;
    gWT  = T;
    gAMP = uLumps*0.075 + 0.05;                 // bound for the trace bracket

    // ---- the wade ----
    float z = t*1.05;
    vec3 ro = vec3(pathX(z) + uBob*0.026*sin(t*2.4), 1.50 + uBob*0.034*sin(t*4.8), z);
    gCam = ro.xz;

    float dz = 0.35;
    vec3 fw = normalize(vec3(pathX(z + dz) - pathX(z - dz), 0.0, 2.0*dz));
    float yaw = uYaw + uBob*0.05*sin(t*0.79 + 1.0);
    float cy = cos(yaw), sy = sin(yaw);
    fw = normalize(vec3(fw.x*cy + fw.z*sy, 0.0, fw.z*cy - fw.x*sy));
    vec3 rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
    vec3 up = cross(rt, fw);
    fw = normalize(fw + up*(uPitch + uBob*0.020*sin(t*1.6)));
    rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
    up = cross(rt, fw);
    float roll = uBob*0.013*sin(t*2.0 + 0.5);
    vec3 rr2 = rt*cos(roll) + up*sin(roll);
    up = up*cos(roll) - rt*sin(roll);
    rt = rr2;

    // ---- the lens, bent before anything is traced ----
    vec2 uv = (fragCoord - 0.5*R)/R.y;
    float vhs = uVHS;
    float ln = hash11(floor(fragCoord.y) + floor(T*24.0)*17.0);
    uv.x += vhs*(0.0024*(ln - 0.5) + 0.0028*sin(T*1.7 + fragCoord.y*0.031));
    uv *= 1.0 + vhs*0.045*dot(uv, uv);
    vec3 rd = normalize(uv.x*rt + uv.y*up + uLens*fw);
    gPX = 1.0/(uLens*R.y);

    // ---- the room ----
    float tt = 0.0, tHit = FAR, mat = 0.0;
    for (int i = 0; i < 144; i++){
        vec2 m = map(ro + rd*tt);
        // the hit test opens up with distance or a grazing ray down a 40 m
        // hall runs out of steps, but not too far: stopping short by e along
        // the ray is a slide of e/cos ALONG a grazing surface, which bands
        // every lighting term computed from the hit point
        if (m.x < 0.0024*tt + 0.0009){ tHit = tt; mat = m.y; break; }
        tt += m.x*0.88;
        if (tt > FAR) break;
    }

    float zf = lightZone((ro + rd*18.0).xz);
    vec3 fogCol = mix(gTile, gLite, 0.5)*(0.05 + 0.10*uBright)
                * mix(1.0, 0.10 + 1.50*zf, uZones);
    fogCol = mix(fogCol, gDry*0.5, 0.25);        // the air in here is not clean
    vec3 col = fogCol;
    float tShow = FAR;

    // ---- the level, traced between the room and the camera ----
    float tS;
    float hitS = traceSurface(ro, rd, tHit, tS);

    if (hitS > 0.5){
        vec3 p = ro + rd*tS;
        col = shadeSurface(p, rd, tS, fogCol);
        tShow = tS;
    } else if (mat > 0.5){
        vec3 p = ro + rd*tHit;
        vec3 n = calcNormal(p);
        // one Newton step onto the true intersection: moving dt along the ray
        // changes the signed distance by dt*dot(n, rd), so the residual maps
        // straight onto a correction in t
        float ndv0 = min(dot(n, rd), -0.06);
        tHit += clamp(-map(p).x/ndv0, 0.0, 0.75);
        p = ro + rd*tHit;
        n = calcNormal(p);
        float ndv = dot(n, rd);
        // sqrt(c*c + k), not max(|c|, k): a hard clamp has an EDGE where it
        // engages, and on a flat wall that edge is a cone about the normal - a
        // ring painted across the surface
        float foot = tHit*gPX/sqrt(ndv*ndv + 0.004);
        vec3 em;
        col = shade(p, rd, n, mat, foot, 1.0, em) + em;
        tShow = tHit;
    }

    if (tShow < FAR){
        float fg = 1.0 - exp(-uFog*tShow);
        // the fade must COMPLETE before rays start missing, or a hit at 42 m
        // keeps 15 % of its surface while the ray beside it, which reached 44
        // and missed, keeps none - a hard stripe at the vanishing point
        col = mix(col, fogCol, max(fg, smoothstep(FAR*0.55, FAR*0.92, tShow)));
    }

    // ---- in-scattering: the shafts ----
    if (uHaze > 0.005){
        float tm = min(tShow, FAR);
        float jt = hash12(fragCoord + floor(T*30.0));
        float glow = 0.0;
        for (int i = 0; i < 18; i++){
            // sample positions depend ONLY on i, never on where the ray landed:
            // spacing them tm/18 apart ties the integral to tHit, which jumps
            // several metres between the ray that reaches the far plane and its
            // neighbour that clips a wall, and that is a hard vertical stripe
            // through the vanishing point
            float u = (float(i) + jt)/18.0;
            float s = FAR*u*u;
            if (s > tm) break;
            float ds = FAR*2.0*u/18.0;
            ds *= min(1.0, (tm - s)/max(ds, 1e-3));
            vec3 q = ro + rd*s;
            // nearest coffer by FOLDING the lattice, and the local density of
            // working tubes rather than one fixture's on/off: floor() plus
            // panelOn() are both piecewise constant, so the glow steps at every
            // lattice plane - and the camera stands ON the plane x = 0
            vec2 f  = mod(q.xz - gLOF, gLG) - 0.5*gLG;
            vec2 d2 = max(abs(f) - PANEL_E, 0.0);
            float dl = length(vec3(d2.x, uHeight - 0.04 - q.y, d2.y));
            float on = uLightOn*mix(1.0, 0.05 + 1.90*lightZone(q.xz), uZones);
            float humid = 1.0 + 1.2*exp(-max(q.y - gSL, 0.0)*0.7);   // it rises
            glow += on*humid*ds*exp(-uFog*2.0*s)/(1.0 + 1.5*dl*dl);
        }
        col += mix(gLite, gDry, 0.20)*glow*uHaze*0.30;
    }

    // ---- flies ----
    //
    // Not geometry, and they never touch the distance field: each is a point
    // tested against the primary ray by closest approach, which is a dot
    // product. Two things make them read as insects rather than as dots.
    //
    // THEY DART. The path is piecewise linear between hashed waypoints on a
    // per-fly beat, so the heading changes discontinuously the way a real
    // fly's does. Any smooth orbit - a Lissajous, a helix - reads as a planet,
    // and no amount of speed fixes that.
    //
    // AND THEY ARE MOTION BLURRED. A fly crosses several body lengths in one
    // frame's exposure, so a sharp dot is the one thing it can never be. The
    // closest-approach distance is measured in an ellipse stretched along the
    // velocity projected into the image plane, which is a streak for the cost
    // of one normalise.
    if (uFlies > 0.005){
        float segL = 7.0;
        float nF = floor(1.0 + 11.0*uFlies);
        // back to front, so the near flies win the alpha
        for (int gi = 0; gi < 6; gi++){
            float zi = floor(z/segL) + (5.0 - float(gi));
            float ez = (zi + 0.5)*segL + 3.0*(hash11(zi*1.7 + 0.3) - 0.5);
            float dA = ez - z;
            if (dA < 0.5 || dA > 30.0) continue;

            vec4 ch = hash42(vec2(zi, 3.7));
            // a swarm sits over the level, or up under a coffer where the
            // warmth is - both are where you would actually find them
            float hi = step(0.74, ch.x);
            vec3 cc = vec3(pathX(ez) + 5.5*(ch.y - 0.5),
                           mix(gSL + 0.16 + 0.55*ch.z, uHeight - 0.60, hi), ez);
            float rad = mix(0.50 + 0.70*ch.w, 0.95, hi);

            for (int k = 0; k < 12; k++){
                if (float(k) >= nF) break;
                float id = zi*13.3 + float(k)*7.77;
                float rate = mix(0.9, 3.4, hash11(id + 2.0))*(0.4 + 1.2*uSwarm);
                float ph = T*rate + hash11(id)*13.0;
                float kb = floor(ph), fr = fract(ph);
                vec3 w0 = vec3(hash11(id + kb*1.31), hash11(id + kb*2.71 + 4.0),
                               hash11(id + kb*3.91 + 8.0)) - 0.5;
                vec3 w1 = vec3(hash11(id + (kb + 1.0)*1.31), hash11(id + (kb + 1.0)*2.71 + 4.0),
                               hash11(id + (kb + 1.0)*3.91 + 8.0)) - 0.5;
                vec3 stretch = vec3(2.0, 1.1, 2.0)*rad;
                vec3 fpos = cc + mix(w0, w1, smoothstep(0.0, 1.0, fr))*stretch;
                vec3 vel  = (w1 - w0)*stretch*rate;

                vec3 oc = fpos - ro;
                float tca = dot(oc, rd);
                if (tca < 0.25 || tca > tShow) continue;

                vec3 perp = oc - rd*tca;                     // ray to fly, perpendicular
                vec3 vp = vel - rd*dot(vel, rd);             // velocity in the image plane
                float vl = length(vp);
                float along = 0.0, across = length(perp);
                if (vl > 1e-4){
                    vec3 vd = vp/vl;
                    along = dot(perp, vd);
                    across = length(perp - vd*along);
                }
                float smear = 1.0 + uBlur*clamp(vl*0.05, 0.0, 6.0);
                float d = sqrt((along/smear)*(along/smear) + across*across);

                // A SUB-PIXEL FLY IS THE ALIASING CASE. Fading it out with
                // distance would be wrong - a swarm at ten metres is plainly
                // visible - so a fly whose radius falls under a pixel is
                // widened to a pixel and dimmed by the SQUARE of the widening.
                // That conserves what it contributes and stops the whole swarm
                // scintillating as the flies cross pixel centres.
                float r0 = 0.0048;
                float px = tca*gPX*1.15;
                float rr3 = max(r0, px);
                float energy = (r0*r0)/(rr3*rr3);
                float a = smoothstep(rr3, rr3*0.35, d)*clamp(energy, 0.0, 1.0);
                a *= smoothstep(30.0, 22.0, dA)*exp(-uFog*0.8*tca);
                if (a < 0.002) continue;

                // near-black, with the one green-blue glint off a wing that is
                // the only colour a fly ever shows
                float glint = pow(max(dot(normalize(vp + rd*0.001), rt), 0.0), 8.0);
                vec3 fc = mix(vec3(0.045, 0.042, 0.040),
                              vec3(0.16, 0.30, 0.22), glint*0.7)*(0.4 + 0.8*uBright);
                fc = mix(fc, fogCol, 1.0 - exp(-uFog*tca));
                col = mix(col, fc, clamp(a, 0.0, 1.0));
            }
        }
    }

    // ---- and then it is put on tape ----
    vec2 sc = fragCoord/R;
    col = col/(1.0 + col*0.80);
    col *= 1.0 - uScan*vhs*0.18*(0.5 + 0.5*sin(fragCoord.y*PI));

    float nz = hash12(fragCoord + floor(T*24.0)*vec2(13.7, 91.3));
    float hb = smoothstep(0.026, 0.0, sc.y);
    col = mix(col, vec3(nz)*0.5 + col*0.45, hb*vhs*0.65);

    float drow = hash11(floor(sc.y*180.0) + floor(T*12.0)*7.3);
    float dseg = step(0.996, hash11(floor(sc.x*40.0)*3.1 + floor(sc.y*180.0)*11.7 + floor(T*12.0)*29.0));
    col += vhs*0.18*dseg*step(0.989, drow)*vec3(1.0, 0.98, 0.93);

    col += (nz - 0.5)*uGrain*0.13;
    float r2 = dot(uv, uv);
    col.r *= 1.0 - vhs*0.07*r2;
    col.b *= 1.0 - vhs*0.14*r2;
    col *= 1.0 - uVignette*1.25*dot(sc - 0.5, sc - 0.5);
    float lum = dot(col, vec3(0.299, 0.587, 0.114));
    col = mix(vec3(lum), col, 1.40 - vhs*0.20);
    col = col*(1.0 - vhs*0.05) + vhs*0.020;

    // half an LSB before the 8-bit write: a dim brown room is almost entirely
    // smooth low-amplitude gradient, which is exactly where quantisation shows
    // up as contour bands - and the grain slider cannot be relied on, because
    // a viewer is entitled to set it to zero
    col += (hash12(fragCoord + 0.5) - 0.5)/255.0;

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