/*{
  "DESCRIPTION": "THE ROAD AT NIGHT. A straight two-lane road through an empty grass field, sparse sodium streetlights, and around every lamp a swarm of several thousand moths. The swarm is the point, and it is not a particle system: a fragment shader has no particle list, it has a ray and a hash, so the moths ARE a three-dimensional LATTICE filling the volume around the lamp - one moth per cell - and the swarming is a FLOW FIELD applied to that lattice. To find which moth is at a sample point, the point is pushed BACKWARD through the flow and the lattice is read where it lands. Every moth in the volume is therefore carried coherently by the same vortex, for the cost of one inverse rotation per sample, and the population is set by the CELL SIZE rather than by any loop bound - which is the whole reason the count can be five thousand per lamp instead of the dozen a splat loop would afford. The flow has to be exactly invertible for this to be legal, so it is built only from operations that preserve the radius it is parameterised by: a rotation about the lamp's vertical axis whose rate is a function of that radius, and a vertical circulation that does not touch it. The rate falls with radius, so the inner moths lap the outer ones and the swarm SHEARS rather than turning like a plate, which is the difference between a swarm and a carousel. Once a moth is found in lattice space it is pushed FORWARD again into the world, so that it can be splatted against the true perpendicular distance from the ray - the along-ray separation is thrown away, because the extent of a blob along the view direction is invisible, and keeping it would force the march step down to the size of a moth and cost eighty samples where twenty-six will do. A moth smaller than a pixel widens to the pixel and DIMS by the ratio of the areas, so a distant swarm resolves to a haze around the lamp instead of aliasing into crawling sparkle. The grass is a marched height field, stepped by vertical clearance divided by the ray's descent rate rather than by a fixed increment, because a ray aimed at the horizon has to cross two hundred metres and a ray aimed at the verge has to resolve a leaf. Its detail is band-limited by distance from the camera - the blade octaves fade out with range while the mean height does not - so the marcher can never step over a blade it could not have resolved in the first place, which is the cheap way to buy an anti-aliased silhouette. The lamps are a z-repeated distance field with a bounding-box short circuit, and BOTH nearest periods are tested rather than one, because the alternating-sides mode puts consecutive lamps on opposite kerbs and the nearest lamp on the far side is then never the nearest period. Their light is a cobra-head distribution - a downward-weighted lobe rather than a point - gathered over the six nearest fittings, so the pools on the tarmac overlap the way real ones do and the gaps between them go properly black. The sky is a night gradient with a warm town glow sitting on the horizon, stars occluded by a domain-warped cloud layer whose coverage is a threshold rather than a fade, so the clouds have EDGES the way painted clouds do; wires between the utility poles are drawn analytically as closest-approach to a segment, because a three-millimetre catenary in a distance field is a hundred wasted steps. The grade is cel-banded and palette-driven - four anime night palettes - and the whole thing is aimed at one particular shot: a dark, empty road, lit at intervals, with somebody walking down it. Forty controls: walk, lens, yaw and pitch, position across the road, eye height, sway, road width, shoulder, markings, wear, grass length, patchiness, wind, lamp spacing, WHICH SIDES the lamps stand on, lamp height, brightness, colour temperature, dead lamps, light cone, MOTH COUNT, swarm radius, swarm spin, moth size, CLOUD COVER, cloud drift, stars, moon, town glow, utility poles, wires, treeline, four palettes, fog, cel banding, bloom, grain, vignette and contrast.",
  "CREDIT": "CC0. Original. Value noise constants after Dave Hoskins.",
  "CATEGORIES": ["generator", "raymarch", "night", "road", "anime", "moths", "swarm", "particles", "procedural", "3d"],
  "INPUTS": [
    { "NAME": "uSpeed",     "TYPE": "float", "DEFAULT": 1.10,  "MIN": 0.0,   "MAX": 6.0 },
    { "NAME": "uLens",      "TYPE": "float", "DEFAULT": 1.35,  "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.02,  "MIN": -0.5,  "MAX": 0.5 },
    { "NAME": "uCamX",      "TYPE": "float", "DEFAULT": 1.30,  "MIN": -6.0,  "MAX": 6.0 },
    { "NAME": "uCamY",      "TYPE": "float", "DEFAULT": 1.55,  "MIN": 0.30,  "MAX": 4.0 },
    { "NAME": "uBob",       "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 2.0 },
    { "NAME": "uRoadW",     "TYPE": "float", "DEFAULT": 3.10,  "MIN": 1.6,   "MAX": 7.0 },
    { "NAME": "uShoulder",  "TYPE": "float", "DEFAULT": 0.75,  "MIN": 0.10,  "MAX": 2.5 },
    { "NAME": "uMarks",     "TYPE": "float", "DEFAULT": 0.85,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uWear",      "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uGrass",     "TYPE": "float", "DEFAULT": 0.34,  "MIN": 0.0,   "MAX": 1.30 },
    { "NAME": "uPatch",     "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uWind",      "TYPE": "float", "DEFAULT": 0.70,  "MIN": 0.0,   "MAX": 2.0 },
    { "NAME": "uSpacing",   "TYPE": "float", "DEFAULT": 34.0,  "MIN": 10.0,  "MAX": 90.0 },
    { "NAME": "uSide",      "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uLampH",     "TYPE": "float", "DEFAULT": 8.20,  "MIN": 3.5,   "MAX": 13.0 },
    { "NAME": "uBright",    "TYPE": "float", "DEFAULT": 1.45,  "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uTemp",      "TYPE": "float", "DEFAULT": 2150.0,"MIN": 1700.0,"MAX": 5200.0 },
    { "NAME": "uLampOut",   "TYPE": "float", "DEFAULT": 0.12,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uCone",      "TYPE": "float", "DEFAULT": 0.85,  "MIN": 0.0,   "MAX": 2.5 },
    { "NAME": "uMothN",     "TYPE": "float", "DEFAULT": 0.55,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uSwarmR",    "TYPE": "float", "DEFAULT": 1.15,  "MIN": 0.30,  "MAX": 3.5 },
    { "NAME": "uSpin",      "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uMothSize",  "TYPE": "float", "DEFAULT": 0.034, "MIN": 0.010, "MAX": 0.090 },
    { "NAME": "uClouds",    "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uCloudSpd",  "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uStars",     "TYPE": "float", "DEFAULT": 0.70,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uMoon",      "TYPE": "float", "DEFAULT": 0.35,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uTown",      "TYPE": "float", "DEFAULT": 0.50,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uPoles",     "TYPE": "float", "DEFAULT": 0.85,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uWires",     "TYPE": "float", "DEFAULT": 0.90,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uTrees",     "TYPE": "float", "DEFAULT": 0.70,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uPalette",   "TYPE": "float", "DEFAULT": 0.0,   "MIN": 0.0,   "MAX": 3.0 },
    { "NAME": "uFog",       "TYPE": "float", "DEFAULT": 0.014, "MIN": 0.0,   "MAX": 0.070 },
    { "NAME": "uCel",       "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uBloom",     "TYPE": "float", "DEFAULT": 1.0,   "MIN": 0.0,   "MAX": 2.5 },
    { "NAME": "uGrain",     "TYPE": "float", "DEFAULT": 0.35,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uVignette",  "TYPE": "float", "DEFAULT": 0.40,  "MIN": 0.0,   "MAX": 1.0 },
    { "NAME": "uContrast",  "TYPE": "float", "DEFAULT": 0.80,  "MIN": 0.0,   "MAX": 1.5 }
  ]
}*/

// THE ROAD AT NIGHT - ISF / GLSL ES 1.00 (gl_FragColor, TIME, RENDERSIZE).
//
// uSide: 0 = left kerb, 1 = right kerb, 2 = both, 3 = alternating.
// uPalette: 0 = prussian night, 1 = deep cyan, 2 = mercury green, 3 = violet dusk.

#define PI  3.14159265359
#define FAR 200.0

float gRoadH;      // road half width
float gShoW;       // shoulder width
float gSpc;        // lamp spacing
float gLampB;      // |x| of the lamp bases
float gPoleB;      // |x| of the utility poles
float gPoleS;      // utility pole spacing
float gLampH;
float gPX;         // world size of one pixel at unit distance
float gT;          // clock
vec2  gCam;        // camera plan position, for the grass level of detail

vec3 gLamp;        // lamp colour
vec3 pZen, pHor, pTown, pGrassA, pGrassB, pRoad, pLine, pCloud, pMoon;

// ---------------------------------------------------------------
// hashes and value noise (Dave Hoskins' constants)
// ---------------------------------------------------------------
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);
}

vec2 hash22(vec2 p){
    vec3 q = fract(vec3(p.xyx)*vec3(0.1031, 0.1030, 0.0973));
    q += dot(q, q.yzx + 33.33);
    return fract((q.xx + q.yz)*q.zy);
}

// The swarm lives on a 3D lattice, so the moth hash has to be three-dimensional
// and it has to be good: a weak 3D hash puts visible planes through the swarm,
// and a swarm with planes in it reads as a texture rather than as insects.
vec4 hash43(vec3 p){
    vec4 q = fract(vec4(p.xyzx)*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. Value noise has axis-aligned
// iso-contours, and unrotated octaves stack their cell edges into a quilt that
// is perfectly obvious on a flat field seen at a grazing angle.
float fbm(vec2 p, int oct){
    mat2 R = mat2(0.8040, 0.5946, -0.5946, 0.8040);
    float s = 0.0, a = 0.5;
    for (int i = 0; i < 6; i++){
        if (i >= oct) break;
        s += a*vnoise(p);
        p = R*p*2.03 + 17.7;
        a *= 0.5;
    }
    return s;
}

// ---------------------------------------------------------------
// blackbody, for the lamps
// ---------------------------------------------------------------
// Normalised by its own maximum channel and then pulled back to a fixed
// luminance, so the temperature slider changes the COLOUR of the sodium and
// not the exposure of the scene.
vec3 kelvin(float k){
    k = clamp(k, 1000.0, 15000.0)*0.01;
    float r, g, b;
    if (k <= 66.0){
        r = 1.0;
        g = clamp(0.39008158*log(max(k, 1.0)) - 0.63184144, 0.0, 1.0);
    } else {
        r = clamp(1.29293619*pow(max(k - 60.0, 0.001), -0.1332047), 0.0, 1.0);
        g = clamp(1.12989086*pow(max(k - 60.0, 0.001), -0.0755148), 0.0, 1.0);
    }
    if (k >= 66.0)      b = 1.0;
    else if (k <= 19.0) b = 0.0;
    else                b = clamp(0.54320679*log(max(k - 10.0, 1.0)) - 1.19625408, 0.0, 1.0);
    vec3 c = vec3(r, g, b)/max(max(r, g), max(b, 1e-3));
    return c*(0.92/max(dot(c, vec3(0.299, 0.587, 0.114)), 1e-3));
}

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);
}

// Coverage of the band |x| < r by a pixel w wide. A road line is 12 cm across
// and a lamp pole is 17 cm; both are under a pixel long before they are out of
// frame, and a hard test keeps drawing them at full strength wherever they
// happen to land, which crawls.
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);
}

// CEL SHADING, but band-limited. Quantising light into steps is what makes a
// painted frame read as painted; doing it with a bare floor() puts a hard
// contour across every smooth gradient, and those contours alias horribly once
// the camera moves. So the step edge is softened by the local derivative of the
// value itself, which is a one-line approximation to filtering the banding.
float cel(float x, float steps, float soft){
    float s = max(steps, 1.0);
    float g = x*s;
    float f = fract(g);
    float w = clamp(soft, 0.02, 0.9);
    return (floor(g) + smoothstep(0.5 - w, 0.5 + w, f))/s;
}

// ---------------------------------------------------------------
// primitives
// ---------------------------------------------------------------
float sdCapsule(vec3 p, vec3 a, vec3 b, float r){
    vec3 pa = p - a, ba = b - a;
    float h = clamp(dot(pa, ba)/dot(ba, ba), 0.0, 1.0);
    return length(pa - ba*h) - r;
}

float sdBox3(vec3 p, vec3 b){
    vec3 q = abs(p) - b;
    return length(max(q, 0.0)) + min(max(q.x, max(q.y, q.z)), 0.0);
}

float sdRoundBox(vec3 p, vec3 b, float r){ return sdBox3(p, b) - r; }

// Closest approach between a ray and a segment, in closed form. Returns the
// separation; `t` comes back as the distance along the ray, so the result can
// be depth-tested against whatever the ray actually hit.
float raySeg(vec3 ro, vec3 rd, vec3 a, vec3 b, out float t){
    vec3 ba = b - a, w0 = ro - a;
    float C = dot(ba, ba);
    float B = dot(rd, ba);
    float D = dot(rd, w0);
    float E = dot(ba, w0);
    float den = C - B*B;
    float s = clamp((E - B*D)/max(den, 1e-6), 0.0, 1.0);
    vec3 q = a + ba*s;
    t = max(dot(q - ro, rd), 0.0);
    return length(q - (ro + rd*t));
}

// ---------------------------------------------------------------
// the ground
// ---------------------------------------------------------------
// A travelling shear rather than a scroll: grass does not slide, it leans, and
// the lean has to run across the field as a gust or the whole verge breathes
// in unison like a lung.
vec2 windOff(vec2 p){
    float a = sin(p.x*0.26 + p.y*0.17 - gT*0.85) + 0.7*sin(p.y*0.43 + gT*0.61)
            + 0.4*sin(p.x*0.9 - p.y*0.6 + gT*1.9);
    return vec2(a, a*0.55)*0.055*uWind;
}

// THE LEVEL OF DETAIL IS NOT COSMETIC. The blade octaves are faded out with
// distance from the camera while the mean height is kept, so a ray forty
// metres out is marching a smooth field and cannot step over a blade it could
// never have resolved. Leave the detail in and the march either misses grass
// and shows holes, or costs four times as many steps to stop missing it.
float grassHeight(vec2 p){
    float d = distance(p, gCam);
    float fb = exp(-d*0.115);
    float ft = exp(-d*0.011);

    vec2 q = p + windOff(p);
    float tuft  = fbm(q*0.55, 3);
    float blade = vnoise(q*21.0)*0.62 + vnoise(q*52.0 + 7.3)*0.38;
    float patch = mix(1.0, 0.30 + 1.25*fbm(p*0.105, 2), uPatch);
    return uGrass*patch*(0.34 + 0.66*mix(0.5, tuft, ft))*mix(0.55, blade, fb);
}

// 0 on the carriageway, 1 out in the field, with the gravel shoulder between
float grassMask(float x){
    float a = abs(x);
    return smoothstep(gRoadH + 0.05, gRoadH + gShoW, a);
}

float terrainH(vec2 p){
    float a = abs(p.x);
    // The camber. A road is crowned so it drains, and the crown is the reason
    // the lamp pools on the near lane are brighter than on the far one.
    float camber = 0.045*(1.0 - clamp(a/(gRoadH + gShoW), 0.0, 1.0));
    return camber + grassMask(p.x)*grassHeight(p);
}

vec3 terrainN(vec2 p, float e){
    e = max(e, 0.006);
    float h  = terrainH(p);
    float hx = terrainH(p + vec2(e, 0.0));
    float hz = terrainH(p + vec2(0.0, e));
    return normalize(vec3(h - hx, e, h - hz));
}

// ---------------------------------------------------------------
// the lamps
// ---------------------------------------------------------------
// Which kerbs carry a fitting at this index. Slot 0 and slot 1; a slot returns
// 0 when it is empty, otherwise the sign of the side it stands on.
float lampSlot(float zi, float k){
    float m = floor(uSide + 0.5);
    if (m < 0.5)      return (k < 0.5) ? -1.0 : 0.0;    // left kerb only
    else if (m < 1.5) return (k < 0.5) ?  1.0 : 0.0;    // right kerb only
    else if (m < 2.5) return (k < 0.5) ? -1.0 : 1.0;    // both
    float s = (mod(zi, 2.0) < 0.5) ? -1.0 : 1.0;        // alternating
    return (k < 0.5) ? s : 0.0;
}

float poleSide(){ return (floor(uSide + 0.5) < 0.5) ? 1.0 : -1.0; }

// where the light actually comes from: under the cobra head, out over the road
vec3 lampPos(float zi, float s){
    return vec3(s*(gLampB - 2.02), gLampH - 0.05, zi*gSpc);
}

// Dead lamps, and the handful that are on the way out. A flat "some are off"
// is worth having on its own; the ones caught mid-failure, strobing at a few
// hertz, are worth more than all the rest of the lighting model put together.
float lampLive(float zi, float s){
    float h = hash12(vec2(zi, s*7.3 + 2.1));
    float on = step(uLampOut, h);
    float dying = step(h, uLampOut + 0.09);
    float f = 1.0 - dying*0.75*smoothstep(0.30, 0.95, 0.5 + 0.5*sin(gT*(2.7 + 11.0*h) + h*40.0));
    return on*f;
}

// The fitting, in local coordinates: +x reaches in over the carriageway, y is
// up, z runs along the kerb.
float lampBody(vec3 q){
    float H = gLampH;
    // A bounding box that short-circuits the assembly. Without it every march
    // step pays for five capsules and a rounded box, and the lamp march
    // becomes the most expensive thing in the frame by a factor of four.
    float b = sdBox3(q - vec3(1.05, H*0.5, 0.0), vec3(1.30, H*0.5 + 0.42, 0.26));
    if (b > 0.26) return b;

    float d = sdCapsule(q, vec3(0.0, 0.0, 0.0), vec3(0.0, H - 0.75, 0.0), 0.086);
    d = min(d, sdCapsule(q, vec3(0.0, 0.0, 0.0), vec3(0.0, 0.42, 0.0), 0.125));
    d = min(d, sdCapsule(q, vec3(0.0, H - 0.75, 0.0), vec3(0.44, H - 0.03, 0.0), 0.074));
    d = min(d, sdCapsule(q, vec3(0.44, H - 0.03, 0.0), vec3(1.60, H + 0.10, 0.0), 0.064));
    d = min(d, sdRoundBox(q - vec3(2.00, H + 0.015, 0.0), vec3(0.26, 0.040, 0.090), 0.082));
    return d;
}

float poleBody(vec3 q){
    float H = gLampH*1.22;
    float b = sdBox3(q - vec3(0.0, H*0.5, 0.0), vec3(1.00, H*0.5 + 0.55, 0.24));
    if (b > 0.26) return b;

    float d = sdCapsule(q, vec3(0.0, 0.0, 0.0), vec3(0.05, H, 0.0), 0.100);
    d = min(d, sdCapsule(q, vec3(-0.88, H - 0.55, 0.0), vec3(0.88, H - 0.55, 0.0), 0.044));
    d = min(d, sdCapsule(q, vec3(-0.58, H - 1.42, 0.0), vec3(0.58, H - 1.42, 0.0), 0.038));
    d = min(d, sdRoundBox(q - vec3(0.24, H - 2.55, 0.0), vec3(0.075, 0.24, 0.075), 0.085));
    return d;
}

// Everything standing up out of the ground, as one field.
//
// BOTH nearest periods are tested, not one. Folding to the single nearest
// period is exact when every period is identical - but in the alternating mode
// consecutive lamps stand on opposite kerbs, so the nearest fitting on the far
// side is never in the nearest period, and a march that only looks at one
// walks straight through it.
float structures(vec3 p){
    float d = 1e6;

    float zi = floor(p.z/gSpc + 0.5);
    float zs = (p.z - zi*gSpc >= 0.0) ? 1.0 : -1.0;
    for (int n = 0; n < 2; n++){
        float zz = zi + float(n)*zs;
        float zl = p.z - zz*gSpc;
        for (int k = 0; k < 2; k++){
            float s = lampSlot(zz, float(k));
            if (abs(s) > 0.5) d = min(d, lampBody(vec3(-s*p.x + gLampB, p.y, zl)));
        }
    }

    if (uPoles > 0.02){
        float ps = poleSide();
        float pb = floor((p.z - gPoleS*0.5)/gPoleS + 0.5);
        float ss = (p.z - gPoleS*0.5 - pb*gPoleS >= 0.0) ? 1.0 : -1.0;
        for (int n = 0; n < 2; n++){
            float pz = (pb + float(n)*ss)*gPoleS + gPoleS*0.5;
            d = min(d, poleBody(vec3(-ps*p.x + gPoleB, p.y, p.z - pz)));
        }
    }
    return d;
}

vec3 structN(vec3 p){
    vec2 e = vec2(1.0, -1.0)*0.0025;
    return normalize(e.xyy*structures(p + e.xyy) + e.yyx*structures(p + e.yyx)
                   + e.yxy*structures(p + e.yxy) + e.xxx*structures(p + e.xxx));
}

// ---------------------------------------------------------------
// the sky
// ---------------------------------------------------------------
float starField(vec3 rd){
    if (uStars < 0.01 || rd.y < 0.0) return 0.0;
    vec2 s = vec2(atan(rd.z, rd.x)*3.0, acos(clamp(rd.y, -1.0, 1.0))*6.0)*11.0;
    vec2 id = floor(s), f = fract(s);
    vec2 h = hash22(id + 4.7);
    float d = length(f - (0.2 + 0.6*h));
    float mag = hash12(id*1.7 + 9.3);
    float pt = smoothstep(0.13, 0.0, d)*step(0.865, mag);
    // Scintillation, not a sine: a star twinkles because the column of air in
    // front of it is turbulent, so the rate has to differ per star or the whole
    // sky pulses together like a string of fairy lights.
    float tw = 0.35 + 0.65*(0.5 + 0.5*sin(gT*(1.3 + 5.0*h.x) + h.y*63.0));
    return pt*tw*(0.45 + 0.75*(mag - 0.865)*7.0);
}

// A THRESHOLD, not a fade. Painted clouds have edges; a raw fbm ramp reads as
// fog, and no amount of contrast afterwards puts the edge back in. The warp is
// what stops the threshold from looking like a contour map of noise.
float cloudField(vec3 rd){
    if (uClouds < 0.01 || rd.y < 0.004) return 0.0;
    vec2 p = rd.xz*(1.0/max(rd.y, 0.020))*0.0135 + vec2(gT*uCloudSpd*0.012, gT*uCloudSpd*0.004);
    vec2 w = vec2(fbm(p*0.75 + 1.3, 3), fbm(p*0.75 + 5.1, 3));
    float f = fbm(p + w*1.25, 4);
    float cov = mix(0.74, 0.235, clamp(uClouds, 0.0, 1.0));
    float c = smoothstep(cov, cov + 0.155, f);
    // fade the layer out at the horizon, where a plane projection goes singular
    return c*smoothstep(0.004, 0.075, rd.y);
}

vec3 skyCol(vec3 rd, float px){
    float up = clamp(rd.y, 0.0, 1.0);
    vec3 c = mix(pHor, pZen, pow(up, 0.52));

    // The town over the hill. This is most of what separates a night sky near
    // people from deep space, and it is the one warm thing in the frame that
    // is not a lamp.
    float hz = exp(-up*9.0);
    float az = 0.55 + 0.45*smoothstep(-0.5, 1.0, rd.z);
    c += pTown*hz*az*uTown*0.85;

    // the moon, low and small, with the haze it drags around itself
    if (uMoon > 0.01){
        vec3 md = normalize(vec3(-0.42, 0.30, 0.86));
        float ca = dot(rd, md);
        c += pMoon*uMoon*(0.55*pow(max(ca, 0.0), 900.0) + 0.10*pow(max(ca, 0.0), 14.0));
    }

    float cl = cloudField(rd);
    c += vec3(0.80, 0.86, 1.00)*starField(rd)*uStars*0.60*(1.0 - cl);

    // Cloud underlit by the town and topped by the moon. Anime clouds are two
    // flat values and a rim, not a gradient, so they get the cel treatment too.
    float shape = cel(clamp(cl, 0.0, 1.0), mix(2.0, 4.0, uCel), 0.22);
    vec3 cc = mix(pCloud*0.42, pCloud, smoothstep(0.02, 0.30, up));
    cc += pTown*uTown*1.15*hz*az;
    c = mix(c, cc, mix(cl, shape, uCel)*0.92);
    return c;
}

// A far treeline sitting on the horizon, which is what gives the road
// somewhere to vanish INTO rather than simply running out of resolution.
float treeLine(vec3 rd, float px){
    if (uTrees < 0.02) return 0.0;
    float D = 340.0;
    float a = atan(rd.x, max(abs(rd.z), 0.02)*sign(rd.z + 1e-6));
    float hh = 7.0 + 15.0*fbm(vec2(a*7.0, 0.0), 4)
             + 5.0*vnoise(vec2(a*33.0, 3.1)) + 3.0*vnoise(vec2(a*95.0, 8.7));
    float e = (hh - 1.6)/D;
    float w = max(px*1.2, 1.0e-4);
    return smoothstep(e + w, e - w, rd.y)*step(0.0, rd.y + 0.02)*uTrees;
}

// ---------------------------------------------------------------
// the wires
// ---------------------------------------------------------------
// Drawn analytically. A three-millimetre catenary in a distance field is a
// hundred wasted steps and it still crawls; the closest approach between a ray
// and a segment is a closed form. A wire thinner than a pixel is widened to the
// pixel and DIMMED by the ratio, which is what a real thin line does and what
// keeps it from strobing as the camera walks.
void wires(vec3 ro, vec3 rd, float tmax, out vec3 wcol, out float wa){
    wcol = vec3(0.0); wa = 0.0;
    if (uWires < 0.02 || uPoles < 0.02) return;

    float ps = poleSide();
    float H  = gLampH*1.22;
    float bx = ps*gPoleB;
    float base = floor((ro.z - gPoleS*0.5)/gPoleS) ;

    for (int i = 0; i < 5; i++){
        float z0 = (base + float(i))*gPoleS + gPoleS*0.5;
        float z1 = z0 + gPoleS;
        for (int k = 0; k < 3; k++){
            float fk = float(k);
            // three conductors: two on the top arm, one on the lower
            float ox = (fk < 2.0) ? (fk*1.52 - 0.76) : 0.0;
            float hy = (fk < 2.0) ? (H - 0.60) : (H - 1.47);
            float sag = ((fk < 2.0) ? 0.55 : 0.42)*(gPoleS/34.0);
            vec3 a = vec3(bx - ps*ox, hy, z0);
            vec3 b = vec3(bx - ps*ox, hy, z1);
            vec3 m = 0.5*(a + b) - vec3(0.0, sag, 0.0);
            // the catenary as two chords; a third would not be visible at any
            // span this shader can put on screen
            for (int q = 0; q < 2; q++){
                vec3 s0 = (q == 0) ? a : m;
                vec3 s1 = (q == 0) ? m : b;
                float t;
                float dd = raySeg(ro, rd, s0, s1, t);
                if (t > 0.5 && t < tmax){
                    float rw = 0.019;
                    float pr = max(rw, t*gPX*0.75);
                    float cov = clamp(1.0 - dd/pr, 0.0, 1.0);
                    cov = cov*cov*(3.0 - 2.0*cov)*clamp(rw/pr, 0.0, 1.0)*uWires;
                    // a wire is a silhouette with one lit edge, from whatever
                    // lamp happens to be under it
                    vec3 wc = pZen*0.55 + gLamp*0.10;
                    wcol = wcol + wc*cov*(1.0 - wa);
                    wa   = wa + cov*(1.0 - wa);
                }
            }
        }
    }
}

// ---------------------------------------------------------------
// the light
// ---------------------------------------------------------------
// A cobra head is not a point source. It throws a lobe DOWNWARD and along the
// kerb, and the difference matters: a point source lights the far verge as
// brightly as the tarmac and the road stops reading as a road. `L` points from
// the lamp toward the surface.
float lampGain(vec3 L){
    float dn = clamp(-L.y, 0.0, 1.0);
    float along = 0.55 + 0.45*abs(L.z);
    return (0.10 + 0.90*pow(dn, 0.70))*along;
}

// The six nearest fittings, gathered as points, plus one direction to key the
// specular off. Three periods each side covers a lamp behind, the one you are
// under, and the one ahead - which is exactly the set that contributes
// anything at a spacing where the pools do not quite meet.
vec3 lampLight(vec3 p, vec3 n, out vec3 keyDir, out float keyDist){
    float zi0 = floor(p.z/gSpc + 0.5);
    vec3 sum = vec3(0.0), acc = vec3(0.0);
    float sw = 0.0, accD = 0.0;
    for (int j = 0; j < 3; j++){
        float zz = zi0 + float(j) - 1.0;
        for (int k = 0; k < 2; k++){
            float s = lampSlot(zz, float(k));
            if (abs(s) > 0.5){
                float live = lampLive(zz, s);
                if (live > 0.002){
                    vec3 lp = lampPos(zz, s);
                    vec3 dv = lp - p;
                    float dl = max(length(dv), 0.30);
                    vec3 L = dv/dl;
                    float w = live*lampGain(-L)/(1.0 + 0.020*dl*dl);
                    float ndl = max(dot(n, L), 0.0);
                    sum  += gLamp*w*ndl;
                    acc  += L*w*ndl;
                    accD += dl*w*ndl;
                    sw   += w*ndl;
                }
            }
        }
    }
    keyDir  = normalize(acc + n*0.001 + vec3(0.0, 0.0001, 0.0));
    keyDist = accD/max(sw, 1e-4);
    return sum;
}

// ---------------------------------------------------------------
// the road surface
// ---------------------------------------------------------------
float roadMarks(vec2 p, float foot){
    if (uMarks < 0.01) return 0.0;
    // Japanese rural centre line: 5 m of paint, 5 m of nothing.
    float dash = step(0.5, fract(p.y/10.0));
    float m = band(p.x, 0.075, foot)*dash;
    m = max(m, band(abs(p.x) - (gRoadH - 0.30), 0.070, foot));
    // worn, and worn WHERE THE TYRES GO: the paint survives at the edges of
    // the lane and is polished away in the wheel tracks
    float wear = 0.55 + 0.45*fbm(p*vec2(2.2, 0.28) + 3.0, 3);
    return m*uMarks*clamp(wear, 0.0, 1.0);
}

vec3 shadeGround(vec3 p, vec3 rd, vec3 n, float foot, float dist){
    float ax = abs(p.x);
    float g  = grassMask(p.x);
    vec3 alb; float gloss; float rough;

    // ---------------- carriageway ----------------
    float agg = vnoise(p.xz*95.0)*0.55 + vnoise(p.xz*310.0)*0.45;    // aggregate
    float aggF = smoothstep(1.10, 0.10, foot/0.012);
    vec3 tar = pRoad*(0.86 + 0.30*mix(0.5, agg, aggF));
    // the polished wheel tracks, and the coarser strip between them
    float track = exp(-pow((ax - gRoadH*0.52)/0.55, 2.0));
    tar *= 1.0 - 0.16*track;
    // patches, cracks and old repairs
    float rep = smoothstep(0.56, 0.72, fbm(p.xz*vec2(0.28, 0.10) + 11.0, 3));
    tar = mix(tar, pRoad*0.72, rep*uWear*smoothstep(1.10, 0.10, foot/0.30));
    float mk = roadMarks(p.xz, foot);
    tar = mix(tar, pLine, mk);
    float tarGloss = mix(0.30, 0.12, uWear) + 0.35*track + 0.45*mk;

    // ---------------- gravel shoulder ----------------
    float sh = smoothstep(gRoadH - 0.02, gRoadH + 0.22, ax)*(1.0 - g);
    vec3 grav = pRoad*vec3(1.35, 1.30, 1.20)*(0.72 + 0.55*vnoise(p.xz*46.0));
    tar = mix(tar, grav, sh);
    tarGloss = mix(tarGloss, 0.10, sh);

    // ---------------- the field ----------------
    float blade = vnoise(p.xz*23.0 + windOff(p.xz)*3.0);
    float clump = fbm(p.xz*0.62 + 5.0, 3);
    vec3 gr = mix(pGrassA, pGrassB, clamp(0.25 + 0.9*clump + 0.5*(blade - 0.5), 0.0, 1.0));
    // Grass is dark at its own root. The height field already knows where the
    // roots are, so the occlusion is free: low ground is shadowed by the blades
    // standing over it.
    float rel = clamp(terrainH(p.xz)/max(uGrass, 1e-3), 0.0, 1.0);
    float gao = mix(0.30, 1.0, pow(rel, 0.65));
    gr *= mix(0.55, 1.0, gao);
    // dead straw showing through, and the trodden verge next to the gravel
    gr = mix(gr, pGrassB*vec3(1.25, 1.10, 0.72),
             smoothstep(0.62, 0.86, fbm(p.xz*0.9 + 21.0, 3))*0.45);

    alb   = mix(tar, gr, g);
    gloss = mix(tarGloss, 0.08, g);
    float ao = mix(1.0, gao, g);

    // ---- light ----
    vec3 kd; float kdist;
    vec3 lit = lampLight(p, n, kd, kdist);

    // Cel banding on the LIGHT, never on the albedo. Quantising colour posters
    // the texture; quantising illumination is what a painted cel actually is -
    // a flat body tone and a flat shadow tone with a hard edge between them.
    float lum = dot(lit, vec3(0.299, 0.587, 0.114));
    float ce = cel(clamp(lum*1.5, 0.0, 1.0), mix(3.0, 5.0, 1.0 - uCel),
                   mix(0.45, 0.10, uCel));
    lit = mix(lit, gLamp*ce/1.5, uCel*0.80);

    vec3 sky = mix(pHor, pZen, 0.40)*(0.30 + 0.55*max(n.y, 0.0));
    vec3 col = alb*(lit*(2.6 + 3.6*uBright) + sky*ao + pTown*uTown*0.075*ao);

    // the sheen off wet-looking tarmac and off paint, which is most of what
    // makes the road read as a surface at all at this light level
    if (gloss > 0.02){
        vec3 hv = normalize(kd - rd);
        float sp = pow(max(dot(n, hv), 0.0), mix(12.0, 90.0, gloss));
        col += gLamp*sp*gloss*length(lit)*(0.5 + 0.8*uBright)*0.9;
    }
    return col;
}

// ---------------------------------------------------------------
// the ironmongery
// ---------------------------------------------------------------
vec3 shadeStruct(vec3 p, vec3 rd, vec3 n, float foot){
    // Is this the underside of a cobra head? If so it is not a surface, it is
    // the source, and it has to be blown out or the brightest thing in the
    // frame ends up being the tarmac it is lighting.
    float zi = floor(p.z/gSpc + 0.5);
    float best = 1e6, live = 0.0;
    for (int j = 0; j < 3; j++){
        float zz = zi + float(j) - 1.0;
        for (int k = 0; k < 2; k++){
            float s = lampSlot(zz, float(k));
            if (abs(s) > 0.5){
                float d = length(p - lampPos(zz, s));
                if (d < best){ best = d; live = lampLive(zz, s); }
            }
        }
    }
    float head = smoothstep(0.34, 0.14, best)*clamp(-n.y*1.6 + 0.35, 0.0, 1.0);

    vec3 alb = mix(vec3(0.19, 0.20, 0.22), vec3(0.26, 0.23, 0.19), 0.5);
    alb *= 0.85 + 0.30*vnoise(vec2(p.y*30.0, p.x*30.0 + p.z*30.0));
    vec3 kd; float kdist;
    vec3 lit = lampLight(p, n, kd, kdist);
    vec3 sky = mix(pHor, pZen, 0.5)*(0.16 + 0.26*max(n.y, 0.0));
    vec3 col = alb*(lit*(2.0 + 2.6*uBright) + sky + pTown*uTown*0.08);
    col += gLamp*head*live*(3.0 + 9.0*uBright);
    return col;
}

// ---------------------------------------------------------------
// THE SWARM
// ---------------------------------------------------------------
// The flow field the whole swarm is carried by. It must be EXACTLY invertible,
// because the lattice is read by pushing a sample backward through it, so it is
// built only from operations that preserve the radius it is parameterised by:
// a rotation about the lamp's axis whose rate is a function of that radius, and
// a vertical circulation which does not touch it.
//
// The rate FALLS with radius, so the inner moths lap the outer ones and the
// swarm shears. That is the whole difference between a swarm and a carousel.
vec3 flow(vec3 q, float sgn){
    float r = length(q.xz) + 0.04;
    float a = sgn*uSpin*(1.45/(0.40 + r*1.30))*gT;
    float yb = sgn*0.15*uSpin*sin(gT*0.63 + r*1.85);
    float c = cos(a), s = sin(a);
    return vec3(c*q.x - s*q.z, q.y + yb, s*q.x + c*q.z);
}

void swarm(vec3 ro, vec3 rd, float tmax, float jt, out vec3 mcol, out float ma){
    mcol = vec3(0.0); ma = 0.0;
    if (uMothN < 0.01 || uSwarmR < 0.05) return;

    float R  = uSwarmR;
    // The cell size sets the POPULATION - about fifteen cells across the
    // sphere is three and a half thousand cells, and the occupancy hash keeps
    // roughly a third of them, so this is a few thousand moths per lamp with
    // no loop anywhere that knows how many there are.
    float cs = R/7.5;
    float zi0 = floor(ro.z/gSpc + 0.5);

    for (int j = 0; j < 3; j++){
        float zz = zi0 + float(j);
        for (int k = 0; k < 2; k++){
            float s = lampSlot(zz, float(k));
            if (abs(s) > 0.5){
                float live = lampLive(zz, s);
                if (live > 0.02){
                    vec3 LP = lampPos(zz, s);
                    vec3 C  = LP + vec3(0.0, -0.22*R, 0.0);
                    vec3 oc = ro - C;
                    float b = dot(oc, rd);
                    float cq = dot(oc, oc) - R*R;
                    float hh = b*b - cq;
                    if (hh > 0.0){
                        hh = sqrt(hh);
                        float t0 = max(-b - hh, 0.10);
                        float t1 = min(-b + hh, tmax);
                        if (t1 > t0){
                            float dt = max(cs, (t1 - t0)/26.0);
                            for (int i = 0; i < 26; i++){
                                float t = t0 + (float(i) + jt)*dt;
                                if (t > t1 || ma > 0.985) break;
                                vec3 q2 = flow(ro + rd*t - C, -1.0);   // into lattice space
                                vec3 ci = floor(q2/cs);
                                vec4 h4 = hash43(ci + 0.5);
                                vec3 cc = (ci + 0.5)*cs;
                                float rr = length(cc*vec3(1.0, 1.25, 1.0))/R;
                                // moths pile up ON the lamp, not around it
                                float dens = 0.92*uMothN*smoothstep(1.04, 0.04, rr);
                                if (h4.w < dens){
                                    vec3 mp = cc + 0.80*cs*(h4.xyz - 0.5);
                                    mp += 0.32*cs*vec3(sin(gT*9.1 + h4.x*41.0),
                                                       sin(gT*12.7 + h4.y*33.0),
                                                       sin(gT*10.3 + h4.z*27.0));
                                    vec3 mw = flow(mp, 1.0) + C;       // back to the world
                                    vec3 w = mw - ro;
                                    float md = dot(w, rd);
                                    if (md > 0.12 && md < tmax){
                                        // The along-ray separation is DISCARDED.
                                        // A blob's extent along the view
                                        // direction is invisible, and keeping it
                                        // would force the step down to the size
                                        // of a moth - eighty samples instead of
                                        // twenty-six for a picture nobody can
                                        // tell apart.
                                        float perp = length(w - rd*md);
                                        float wb = 0.30 + 0.70*abs(sin(gT*12.0 + h4.x*50.0));
                                        float ms = uMothSize*(0.70 + 0.65*h4.y)*(0.72 + 0.55*wb);
                                        float pr = max(ms, md*gPX*0.80);
                                        float cov = clamp(1.0 - perp/pr, 0.0, 1.0);
                                        cov = cov*cov*clamp((ms*ms)/(pr*pr), 0.0, 1.0);
                                        if (cov > 0.001){
                                            float dl = length(mw - LP);
                                            float li = live/(1.0 + 1.9*dl*dl);
                                            vec3 mc = gLamp*(0.20 + 3.4*li)*wb*(0.45 + 0.70*uBright)
                                                    + pZen*0.55;
                                            mc *= exp(-uFog*md);
                                            mcol += mc*cov*(1.0 - ma);
                                            ma   += cov*(1.0 - ma);
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

// ---------------------------------------------------------------
// the light in the air
// ---------------------------------------------------------------
// The cone under each lamp, which on a damp night is half the picture. The
// window on the 3x3 set has to reach zero before a lamp leaves it, or the sum
// steps at every lattice plane - and those planes project to straight vertical
// lines through the vanishing point.
vec3 volGlow(vec3 ro, vec3 rd, float tmax, float jt){
    if (uCone < 0.01) return vec3(0.0);
    float lim = min(tmax, 140.0);
    vec3 g = vec3(0.0);
    for (int i = 0; i < 20; i++){
        float u = (float(i) + jt)/20.0;
        float t = lim*u*u;
        float ds = lim*2.0*u/20.0;
        vec3 q = ro + rd*t;
        float zi = floor(q.z/gSpc + 0.5);
        float zs = (q.z - zi*gSpc >= 0.0) ? 1.0 : -1.0;
        for (int n = 0; n < 2; n++){
            float zz = zi + float(n)*zs;
            for (int k = 0; k < 2; k++){
                float s = lampSlot(zz, float(k));
                if (abs(s) > 0.5){
                    vec3 lp = lampPos(zz, s);
                    vec3 dv = lp - q;
                    float dl = max(length(dv), 0.35);
                    float win = smoothstep(1.55*gSpc, 1.05*gSpc, abs(q.z - lp.z));
                    // damp air lies low, and the cone is only visible in it
                    float mist = 0.55 + 0.85*exp(-max(q.y, 0.0)*0.30);
                    g += gLamp*lampLive(zz, s)*lampGain(-dv/dl)*win*mist
                       * ds*exp(-uFog*1.6*t)/(1.0 + 0.60*dl*dl);
                }
            }
        }
    }
    return g*uCone*0.030;
}

// The halo. A lamp seen directly is a hard little rectangle of geometry; the
// bloom around it is not a post effect here but a real integral along the ray,
// so it is occluded correctly by anything standing in front of it.
vec3 lampHalo(vec3 ro, vec3 rd, float tmax){
    if (uBloom < 0.01) return vec3(0.0);
    vec3 g = vec3(0.0);
    float zi0 = floor(ro.z/gSpc + 0.5);
    for (int j = 0; j < 5; j++){
        float zz = zi0 + float(j) - 1.0;
        for (int k = 0; k < 2; k++){
            float s = lampSlot(zz, float(k));
            if (abs(s) > 0.5){
                float live = lampLive(zz, s);
                if (live > 0.02){
                    vec3 lp = lampPos(zz, s);
                    vec3 w = lp - ro;
                    float md = dot(w, rd);
                    if (md > 0.3 && md < tmax){
                        float perp = length(w - rd*md);
                        float r0 = 0.30 + 0.55*uBloom;
                        g += gLamp*live*(r0*r0)/(perp*perp + r0*r0)
                           * exp(-uFog*md)/(1.0 + 0.010*md*md);
                    }
                }
            }
        }
    }
    return g*uBloom*(0.16 + 0.24*uBright);
}

// ---------------------------------------------------------------
void main(){
    vec2 fc = gl_FragCoord.xy;
    vec2 R  = RENDERSIZE.xy;
    float T = TIME;
    gT = T;

    // Four anime night palettes. Every one of them is built the same way: a
    // near-black zenith, a horizon two stops brighter, and exactly one warm
    // thing in the frame that is not a lamp - the town over the hill.
    pZen    = pick4(vec3(0.018,0.030,0.070), vec3(0.020,0.045,0.095), vec3(0.022,0.038,0.055), vec3(0.040,0.028,0.080), uPalette);
    pHor    = pick4(vec3(0.055,0.085,0.150), vec3(0.045,0.130,0.190), vec3(0.055,0.095,0.105), vec3(0.130,0.075,0.145), uPalette);
    pTown   = pick4(vec3(0.115,0.055,0.030), vec3(0.130,0.080,0.045), vec3(0.070,0.105,0.060), vec3(0.165,0.075,0.085), uPalette);
    pGrassA = pick4(vec3(0.130,0.190,0.115), vec3(0.105,0.205,0.160), vec3(0.120,0.180,0.125), vec3(0.150,0.165,0.140), uPalette);
    pGrassB = pick4(vec3(0.265,0.310,0.165), vec3(0.215,0.330,0.225), vec3(0.220,0.290,0.180), vec3(0.290,0.265,0.190), uPalette);
    pRoad   = pick4(vec3(0.200,0.212,0.245), vec3(0.170,0.220,0.262), vec3(0.196,0.216,0.208), vec3(0.222,0.196,0.250), uPalette);
    pLine   = pick4(vec3(0.78,0.77,0.72),    vec3(0.82,0.84,0.80),    vec3(0.76,0.79,0.72),    vec3(0.80,0.76,0.78),    uPalette);
    pCloud  = pick4(vec3(0.085,0.105,0.155), vec3(0.075,0.140,0.185), vec3(0.080,0.110,0.110), vec3(0.130,0.090,0.150), uPalette);
    pMoon   = pick4(vec3(0.85,0.88,1.00),    vec3(0.80,0.90,1.00),    vec3(0.85,0.95,0.95),    vec3(0.95,0.88,1.00),    uPalette);
    gLamp   = kelvin(uTemp);

    gRoadH = uRoadW;
    gShoW  = uShoulder;
    gSpc   = uSpacing;
    gLampH = uLampH;
    gLampB = gRoadH + gShoW*0.85 + 0.45;
    gPoleB = gRoadH + gShoW + 1.60;
    gPoleS = max(gSpc*0.80, 22.0);

    // ---- the walk ----
    float z = T*uSpeed;
    vec3 ro = vec3(uCamX + uBob*0.038*sin(z*1.95 + 0.6),
                   uCamY + uBob*0.028*sin(z*3.90),
                   z);
    gCam = ro.xz;

    float yaw  = uYaw   + uBob*0.018*sin(z*0.87 + 1.0);
    float pit  = uPitch + uBob*0.011*sin(z*1.71);
    vec3 fw = normalize(vec3(sin(yaw), pit, cos(yaw)));
    vec3 rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
    vec3 up = cross(rt, fw);
    float roll = uBob*0.009*sin(z*1.30 + 0.4);
    vec3 r2 = rt*cos(roll) + up*sin(roll);
    up = up*cos(roll) - rt*sin(roll);
    rt = r2;

    vec2 uv = (fc - 0.5*R)/R.y;
    vec3 rd = normalize(uv.x*rt + uv.y*up + uLens*fw);
    gPX = 1.0/(uLens*R.y);

    float jt  = hash12(fc + floor(T*29.0)*vec2(7.1, 3.7));
    float jt2 = hash12(fc*1.37 + floor(T*29.0)*vec2(11.3, 19.1));

    // ---- the ground, as a marched height field ----
    //
    // Stepped by vertical clearance divided by the ray's rate of descent
    // rather than by a fixed increment: a ray aimed at the horizon has to
    // cross two hundred metres and a ray aimed at the verge has to resolve a
    // blade, and one loop has to do both. The clamp at the bottom is what
    // stops a nearly-level ray from taking an infinite first step.
    float maxH = 0.06 + uGrass*1.40;
    float tG = -1.0;
    if (!(rd.y > 0.0 && ro.y > maxH)){
        float t = 0.02;
        if (ro.y > maxH && rd.y < -1e-4) t = max(t, (maxH - ro.y)/rd.y);
        float pt = t, ph = terrainH((ro + rd*t).xz);
        for (int i = 0; i < 150; i++){
            vec3 p = ro + rd*t;
            float h = terrainH(p.xz);
            if (p.y < h){
                float d0 = (ro + rd*pt).y - ph;
                float d1 = p.y - h;
                tG = mix(pt, t, clamp(d0/max(d0 - d1, 1e-5), 0.0, 1.0));
                break;
            }
            pt = t; ph = h;
            t += clamp((p.y - h)*0.55/max(-rd.y, 0.045), 0.030 + 0.009*t, 2.5 + 0.10*t);
            if (t > FAR) break;
        }
    }

    // ---- the ironmongery ----
    float tS = -1.0;
    {
        float t = 0.05;
        for (int i = 0; i < 96; i++){
            float d = structures(ro + rd*t);
            if (d < 0.0018*t + 0.0016){ tS = t; break; }
            t += d*0.94;
            if (t > 175.0) break;
        }
    }

    float tHit = FAR;
    int mat = 0;
    if (tG > 0.0){ tHit = tG; mat = 1; }
    if (tS > 0.0 && tS < tHit){ tHit = tS; mat = 2; }

    vec3 col;
    if (mat == 1){
        vec3 p = ro + rd*tG;
        // The normal is taken at the PIXEL's epsilon, not at a fixed one, so
        // the surface it reports is the surface the pixel can actually see -
        // which is the cheapest anti-aliasing there is on a field of grass.
        float e0 = max(0.005, tG*gPX*1.3);
        vec3 n = terrainN(p.xz, e0);
        float ndv = min(dot(n, rd), -0.02);
        float foot = tG*gPX/sqrt(ndv*ndv + 0.004);
        col = shadeGround(p, rd, n, foot, tG);
    } else if (mat == 2){
        vec3 p = ro + rd*tS;
        vec3 n = structN(p);
        col = shadeStruct(p, rd, n, tS*gPX);
    } else {
        col = skyCol(rd, gPX);
        float tl = treeLine(rd, gPX);
        vec3 tc = mix(pZen, pHor, 0.30)*0.42 + pTown*uTown*0.10;
        col = mix(col, tc, tl);
    }

    // ---- distance ----
    vec3 fogc = mix(pHor, pHor + pTown*uTown*0.8, 0.45)*0.92;
    if (tHit < FAR){
        float f = 1.0 - exp(-uFog*tHit);
        // the fade has to finish before rays start missing, or the last few
        // metres of road keep a fifth of their contrast while the ray beside
        // them keeps none, and that edge is a hard line at the vanishing point
        col = mix(col, fogc, max(f, smoothstep(FAR*0.55, FAR*0.93, tHit)));
    }

    // ---- the light in the air, then what is flying through it ----
    col += volGlow(ro, rd, tHit, jt);
    col += lampHalo(ro, rd, tHit);

    vec3 wc; float wa;
    wires(ro, rd, tHit, wc, wa);
    col = col*(1.0 - wa) + wc;

    vec3 mc; float ma;
    swarm(ro, rd, tHit, jt2, mc, ma);
    col = col*(1.0 - ma) + mc;

    // ---- the print ----
    col = max(col, 0.0);
    col = col/(1.0 + col*0.55);                    // roll the lamps off
    col = mix(col, col*col*(3.0 - 2.0*col), 0.55*uContrast);
    float lum = dot(col, vec3(0.299, 0.587, 0.114));
    col = mix(vec3(lum), col, 1.08 + 0.30*uContrast);
    // night stock leans cold in the shadows, which is what leaves the sodium
    // room to be the only warm thing on screen
    float sh = (1.0 - lum)*(1.0 - lum);
    // Weighted by the shadow term and by nothing else. An unconditional cool
    // multiply also cools the sodium, and a sodium lamp that has come out cyan
    // is the one mistake this palette cannot survive.
    col = mix(col, col*vec3(0.88, 0.96, 1.14) + vec3(0.0, 0.003, 0.011)*sh, 0.80*sh);

    vec2 sc = fc/R;
    col *= 1.0 - uVignette*1.15*dot(sc - 0.5, sc - 0.5);

    // Grain, weighted into the shadows. Film grain lives in the toe; spreading
    // it evenly is the tell that it was added by a slider.
    float nz = hash12(fc + floor(T*24.0)*vec2(13.7, 91.3));
    col += (nz - 0.5)*uGrain*0.075*(0.30 + 1.10*sh);

    // half an LSB, because a night sky is one enormous smooth gradient and
    // that is exactly where eight bits shows its seams
    col += (hash12(fc + 0.5) - 0.5)/255.0;

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