/*{
  "DESCRIPTION": "Escher's regular division of the plane, carved in relief and set moving. One tessellation of interlocking birds and fish -- the interlock is guaranteed rather than drawn: the shapes come from warping a square lattice with a field that has the lattice's own period, so every piece stays congruent to its neighbours and the plane is covered with no gaps however far the shapes bend. Each piece is a solid with a rounded edge; birds carry a raised eye, beak and wing, fish an eye, gill, fin and rows of scales. Figure and ground are made literal in depth: where the pieces are birds the light ones stand proud and the dark ones sink, and along the spiral that swaps until the fish rise and the birds become the ground. The pieces are ray traced as a height field over a cone and lit with shadows and reflections. Three surfaces: PRINT GALLERY winds the Escher-Droste spiral down a funnel whose walls are mapped conformally, so going once round the centre is itself one whole step of the lattice; SMALLER AND SMALLER rings the creatures down to a point; SKY AND WATER lays them on a flat board, birds on the far side and fish on the near. The loop is exact: each LATTICE preset is one where a full turn and a full zoom both land on the lattice with the checkerboard intact, and every motion completes a whole number of cycles per period. FLAP beats the wings and nods the heads -- the silhouette itself moves, so the pieces still fit mid-stroke. RISE sends waves of creatures lifting out of the surface. MORPH sets how far the shapes bend, CREATURE which appear, DEPTH the funnel, TILT, ORBIT and LENS the camera, CENTRE where it looks, RELIEF the height of the pieces, CARVE the detail, LIGHT and SHADOW the lighting. Ten finishes: woodcut, sky and water, lithograph, three-colour print, gold leaf, delft, blueprint, stained glass, neon, and night and day.",
  "CREDIT": "CC0",
  "ISFVSN": "2",
  "CATEGORIES": ["tessellation", "escher", "conformal", "relief", "loop"],
  "INPUTS": [
    { "NAME": "uGeometry", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1, 2],
      "LABELS": ["print gallery", "smaller and smaller", "sky and water"] },
    { "NAME": "uLattice", "TYPE": "long", "DEFAULT": 1, "VALUES": [0, 1, 2, 3],
      "LABELS": ["8 arms", "18 arms", "16 arms", "24 arms"] },
    { "NAME": "uPeriod", "TYPE": "float", "DEFAULT": 18.0, "MIN": 6.0, "MAX": 60.0 },
    { "NAME": "uSpin", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1, -1, 2],
      "LABELS": ["still", "one turn", "one turn back", "two turns"] },
    { "NAME": "uDepth", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0, "MAX": 1.4 },

    { "NAME": "uTilt", "TYPE": "float", "DEFAULT": 58.0, "MIN": 25.0, "MAX": 85.0 },
    { "NAME": "uOrbit", "TYPE": "long", "DEFAULT": 1, "VALUES": [0, 1, 2, 3],
      "LABELS": ["fixed", "sway", "orbit", "orbit back"] },
    { "NAME": "uLens", "TYPE": "float", "DEFAULT": 2.0, "MIN": 1.0, "MAX": 5.0 },
    { "NAME": "uCentre", "TYPE": "point2D", "DEFAULT": [0.5, 0.5] },

    { "NAME": "uMorph", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 0.75 },
    { "NAME": "uCreature", "TYPE": "long", "DEFAULT": 0, "VALUES": [0, 1, 2, 3],
      "LABELS": ["birds and fish", "birds", "fish", "shapes only"] },
    { "NAME": "uFlap", "TYPE": "float", "DEFAULT": 0.7, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uRise", "TYPE": "float", "DEFAULT": 0.6, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uRelief", "TYPE": "float", "DEFAULT": 0.16, "MIN": 0.02, "MAX": 0.35 },
    { "NAME": "uCarve", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.0, "MAX": 2.0 },

    { "NAME": "uStyle", "TYPE": "long", "DEFAULT": 0,
      "VALUES": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
      "LABELS": ["woodcut", "sky and water", "lithograph", "three-colour print",
                 "gold leaf", "delft", "blueprint", "stained glass", "neon",
                 "night and day"] },
    { "NAME": "uLight", "TYPE": "float", "DEFAULT": 0.375, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uShadow", "TYPE": "float", "DEFAULT": 0.85, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "uGain", "TYPE": "float", "DEFAULT": 1.0, "MIN": 0.2, "MAX": 2.0 },
    { "NAME": "uVign", "TYPE": "float", "DEFAULT": 0.35, "MIN": 0.0, "MAX": 1.0 }
  ]
}*/

// =============================================================================
//  regular division -- Escher's birds and fish, carved in relief and moving
// =============================================================================
//
//  world point q  --map-->  tile space T  --warp-->  W  --fract-->  one piece, local L
//  height(q) = funnel(q) + relief(piece, L)            ray traced, then lit
//
//  THE INTERLOCK IS GUARANTEED, NOT DRAWN. The pieces are the unit squares of
//  W = T + A*D(T), where D has the lattice's period. Every region floor(W) = n is
//  congruent to every other, and as long as the warp keeps a positive Jacobian
//  they cover the plane with nothing missing and nothing doubled. That is also
//  what makes the ANIMATION safe: the edge curves move every frame, and the
//  pieces cannot stop fitting, because the partition is re-derived from W each
//  time rather than being a shape that was drawn once.
//
//  WHY THE LATTICE PRESETS ARE WHAT THEY ARE. In the Print Gallery map, tile
//  space is T = (1 - i k)(u + i arg z) / h, with h = 2pi/N, k = M/N and u the
//  conformal coordinate down the funnel wall. Three things must land exactly on
//  the lattice or the picture tears:
//    going once around the centre   shifts T by (M, N)
//    zooming by one full repeat     shifts T by (M, -M*M/N)
//    and a two-colour checkerboard  needs both shifts to be EVEN in total.
//  So N | M*M, M + N even, M + M*M/N even: (4,8) (6,18) (8,16) (12,24). Every
//  phase that moves in time is built from T.x/M, which those two shifts change
//  by whole numbers, so wingbeats and waves also close the loop.
//
//  THE FUNNEL. A cone z = s*r is a surface of revolution, and on one the
//  conformal coordinate down the wall is u = sqrt(1 + s*s) * log r. Using that
//  instead of log r keeps the creatures unsheared on the wall itself, not merely
//  in a view from straight above. The relief is sized in proportion to the
//  local tile, so it is self-similar too, and the top of the relief shell is
//  just a steeper cone -- which is why the ray can be clipped to the shell with a
//  quadratic instead of being marched in from the camera.
// =============================================================================

#define TAU 6.283185307179586
#define FLAT_S 2.2        // tiles per world unit on the flat board
#define FLAT_SHIFT 8.0    // tiles the flat board scrolls per loop: even, and the colour period
#define BEV 0.10          // width of a piece's rounded edge, tile units
#define LIFT_MAX 2.3      // the largest relief() can be, in tile widths x RELIEF

// ---- set once per pixel in main ---------------------------------------------
float gM, gN, gH, gK, gHK, gCs, gShift, gSpinA, gPh, gFlat, gSlope, gWave, gNF, gNR;
float gBirdOn, gFishOn, gCreatureOn, gFade, gFine;

// ---- written by tile() and carve(): the piece under the last point asked -----
vec2  gT, gL, gCell;
float gTW, gDE, gA, gBand, gFig, gPar, gFlap, gNod, gRise, gCav, gPupil, gEye, gFy, gGy;

float h21(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.x + q.y) * q.z);
}

float latM(){
  if (uLattice == 0) return 4.0;
  if (uLattice == 1) return 6.0;
  if (uLattice == 2) return 8.0;
  return 12.0;
}
float latN(){
  if (uLattice == 0) return 8.0;
  if (uLattice == 1) return 18.0;
  if (uLattice == 2) return 16.0;
  return 24.0;
}

// A smooth bump and its slope. The slope is what turns a distance measured in
// warped W space back into a true distance: across a bent edge, W changes faster
// than the plane does, and an uncorrected bevel narrows wherever the edge leans.
vec2 bumpD(float x, float c, float w){
  float s = (x - c) / w;
  float t = clamp(1.0 - abs(s), 0.0, 1.0);
  return vec2(t * t * (3.0 - 2.0 * t), -6.0 * t * (1.0 - t) * sign(s) / w);
}

// ---- the piece under world point q --------------------------------------------
// The silhouette is two EDGE CURVES, Escher's own method: a curve for the left
// edge translated to the right, one for the bottom translated to the top. The
// head pushed out through the right edge is the notch cut into the left, and
// that notch is the tail. FLAP moves the curves: the wing hump beats and the
// head nods, with a phase that drifts slowly across the picture so the flock
// beats in a wave rather than in lockstep.
void tile(vec2 q){
  float r = max(length(q), 1e-4);
  float u = log(r) * gCs - gShift;
  float v = atan(q.y, q.x) + gSpinA;
  vec2 Ts = vec2(u + gK * v, v - gK * u) / gH;
  vec2 Tf = q * FLAT_S - vec2(gPh * FLAT_SHIFT, 0.0);
  gT  = mix(Ts, Tf, gFlat);
  gTW = mix(r * gHK, 1.0 / FLAT_S, gFlat);          // world units per tile
  gBand = mix(cos(TAU * gT.x / gM), clamp(q.y * 0.6, -1.0, 1.0), gFlat);
  // smoothstep has zero slope at zero, so there is no crease where the shapes
  // pass through a plain grid -- abs() alone would draw a line there
  gA = uMorph * smoothstep(0.0, 0.85, abs(gBand));

  float wave = TAU * (gT.x * gWave + gPh * gNF);
  gFlap = uFlap * sin(wave);
  gNod  = uFlap * sin(wave - 1.3);

  float y = fract(gT.y), x = fract(gT.x);
  vec2 F = -0.30 * bumpD(y, 0.56 + 0.07 * gNod, 0.36)          // head
           -0.17 * bumpD(y, 0.64 + 0.10 * gNod, 0.11)          // beak
           +0.13 * bumpD(y, 0.20, 0.17);                        // chin
  vec2 G = -(0.28 + 0.14 * gFlap) * bumpD(x, 0.40, 0.33)       // wing hump
           -(0.10 + 0.14 * gFlap) * bumpD(x, 0.16, 0.09)       // wingtip
           +0.10 * bumpD(x, 0.80, 0.14);                        // dip behind the head
  vec2 W = gT + gA * vec2(F.x, G.x);
  gFy = F.y; gGy = G.y;
  gCell = floor(W);
  gL = W - gCell;
  gDE = min(min(gL.x, 1.0 - gL.x) / sqrt(1.0 + gA * gA * F.y * F.y),
            min(gL.y, 1.0 - gL.y) / sqrt(1.0 + gA * gA * G.y * G.y));
  gPar = mod(gCell.x + gCell.y, 2.0);                            // 0 light, 1 dark

  // which colour is the figure here; the other is the ground
  float bird = smoothstep(-0.15, 0.75, gBand) * gBirdOn;
  float fish = (1.0 - smoothstep(-0.75, 0.15, gBand)) * gFishOn;
  gFig = mix(fish, bird, 1.0 - gPar);
  float c = 0.5 + 0.5 * cos(TAU * (gT.x * gWave * 2.0 - gPh * gNR));
  c *= c;
  gRise = uRise * c * c * gFig;
}

// Height of the piece above the funnel, world units. A quarter-round edge and a
// swell toward the middle; figures stand proud, ground sinks, RISE lifts figures
// further, and FLAP rocks each figure about its body so it banks as it beats.
// Everything is multiplied down to zero at the piece's edge, so neighbouring
// pieces meet in a continuous groove whatever each of them is doing.
float relief(){
  float e = clamp(gDE / BEV, 0.0, 1.0);
  float pillow = sqrt(e * (2.0 - e));
  float dome = smoothstep(0.0, 0.45, gDE);
  float lift = 0.30 + 0.70 * gFig + 0.90 * gRise;
  float rock = 0.34 * gFig * (gFlap * (gL.y - 0.5) + gNod * (gL.x - 0.5));
  return gTW * uRelief * ((0.85 * pillow + 0.15 * dome) * lift + rock * pillow);
}

float surfZ(vec2 q){
  tile(q);
  return gSlope * length(q) + relief();
}

// ---- carved detail --------------------------------------------------------------
// Tile units of height, for the normal only -- the ray is traced against the
// pieces without it, and the shadow ignores it too. Marks are placed against the
// OUTLINE: in L space the piece is a unit square, so a mark at L = (0.80, 0.585)
// sits just inside the head however the edge curves are bent this frame.
//
// Bird and fish share one silhouette, so they share most of their anatomy: the
// beak is the mouth, the wing is the dorsal fin, the forked tail is the tail fin.
// Those are drawn once, with positions blended by which colour the piece is.
//
// The normal needs this height at three points. Written as three calls -- or as
// one call in a loop of three, which the Windows compiler unrolls anyway -- it is
// three copies of the most expensive function in the file to link (measured 287
// ms for one copy, 359 for the loop). So every operation here is done on vec3s
// instead: x, y, z are the height at L + (e,0), at L + (0,e), and at L itself.
// One copy of the code, three answers.
vec3 ridge3(vec3 d, float w){
  vec3 t = clamp(1.0 - d / w, 0.0, 1.0);
  return t * t * (3.0 - 2.0 * t);
}
vec3 len3(vec3 x, vec3 y){ return sqrt(x * x + y * y); }
vec3 seg3(vec3 X, vec3 Y, vec2 a, vec2 b){
  vec3 px = X - a.x, py = Y - a.y;
  vec2 ba = b - a;
  vec3 h = clamp((px * ba.x + py * ba.y) / dot(ba, ba), 0.0, 1.0);
  return len3(px - ba.x * h, py - ba.y * h);
}

vec3 carve3(vec2 L, float e){
  vec3 X = L.x + vec3(e, 0.0, 0.0);
  vec3 Y = L.y + vec3(0.0, e, 0.0);
  float light = 1.0 - gPar;                 // 1 bird, 0 fish
  float fig = gFig * gCreatureOn;
  // eye, and the pupil just forward of centre
  vec2 ec = mix(vec2(0.78, 0.570), vec2(0.80, 0.585), light);
  vec3 ex = X - ec.x, ey = Y - ec.y;
  vec3 eyeM = sqrt(max(1.0 - (ex * ex + ey * ey) / mix(0.0030, 0.0036, light), 0.0));
  vec3 pup = ridge3(len3(ex - 0.011, ey - 0.002), 0.025);
  // beak or mouth, wing or dorsal fin, the forked tail
  vec3 beak = ridge3(seg3(X, Y, mix(vec2(0.90, 0.62), vec2(0.86, 0.612), light),
                                mix(vec2(0.99, 0.628), vec2(0.99, 0.64), light)), 0.035);
  vec2 wc = mix(vec2(0.40, 0.50), vec2(0.40, 0.555 - 0.05 * gFlap), light);
  vec3 wing = ridge3(abs(len3(X - wc.x, Y - wc.y) - mix(0.36, 0.33, light)), 0.048)
            * smoothstep(0.12, 0.20, X) * (1.0 - smoothstep(0.58, 0.66, X))
            * smoothstep(mix(0.76, 0.64, light), mix(0.80, 0.72, light), Y);
  vec3 tail = ridge3(seg3(X, Y, vec2(0.02, 0.71), vec2(0.16, 0.60)), 0.030)
            + ridge3(seg3(X, Y, vec2(0.02, 0.47), vec2(0.16, 0.555)), 0.030);
  // bird only: feathers, a row of cuts across the wing
  vec3 fea = ridge3(abs(fract(X * 8.0 - Y * 2.5) - 0.5), 0.09)
           * smoothstep(0.22, 0.28, X) * (1.0 - smoothstep(0.52, 0.58, X))
           * smoothstep(0.69, 0.73, Y) * (1.0 - smoothstep(0.80, 0.84, Y));
  // fish only: a gill arc, and scales as a pattern -- sparse little curves read as
  // letters ("3e1o" in the flat version), dense staggered rows read as scales
  vec3 gill = ridge3(abs(len3(X - 0.83, Y - 0.55) - 0.20), 0.04)
            * (1.0 - smoothstep(0.64, 0.69, X)) * smoothstep(0.35, 0.40, Y)
            * (1.0 - smoothstep(0.71, 0.76, Y));
  vec3 sx = X * 10.0;
  vec3 scy = fract(Y * 8.0 + 0.5 * mod(floor(sx), 2.0)) - 0.5;
  vec3 scale = ridge3(abs(len3(fract(sx), scy) - 0.46), 0.16) * gFine
             * smoothstep(0.18, 0.26, X) * (1.0 - smoothstep(0.60, 0.66, X))
             * smoothstep(0.30, 0.37, Y) * (1.0 - smoothstep(0.70, 0.76, Y));

  vec3 raised = mix(0.26, 0.30, light) * eyeM + mix(0.20, 0.22, light) * wing
              + light * 0.18 * beak + (1.0 - light) * (0.16 * gill + 0.08 * scale);
  vec3 cut = min(tail + mix(beak, fea * gFine, light), 1.0);
  gEye = step(0.001, eyeM.z) * fig;
  gPupil = pup.z * fig;
  gCav = cut.z * fig * min(uCarve, 1.0);
  // the ground keeps the gouge marks of a block cleared by hand
  vec3 gouge = 0.5 + 0.5 * cos(TAU * 5.0 * (X * 0.8 + Y));
  return uCarve * fig * (raised - 0.14 * cut - 0.18 * pup)
       - (1.0 - gFig) * 0.022 * gouge * gFine;
}

// First crossing of the ray with the cone z = a*r + b (a plane when a = 0), from
// a camera above it. Squaring admits the lower nappe, so roots are kept only
// where the ray is actually above the apex.
float hitCone(vec3 ro, vec3 rd, float a, float b){
  float w0 = ro.z - b;
  if (a < 1e-5) return rd.z < 0.0 ? w0 / -rd.z : 1e4;
  float a2 = a * a;
  float A = rd.z * rd.z - a2 * dot(rd.xy, rd.xy);
  float B = 2.0 * (rd.z * w0 - a2 * dot(ro.xy, rd.xy));
  float C = w0 * w0 - a2 * dot(ro.xy, ro.xy);
  if (abs(A) < 1e-7) {
    float tl = -C / B;
    return (tl > 0.0 && w0 + rd.z * tl > -1e-4) ? tl : 1e4;
  }
  float disc = B * B - 4.0 * A * C;
  if (disc < 0.0) return 1e4;
  float s = sqrt(disc);
  float t1 = (-B - s) / (2.0 * A), t2 = (-B + s) / (2.0 * A);
  float lo = min(t1, t2), hi = max(t1, t2);
  float res = 1e4;
  if (hi > 0.0 && w0 + rd.z * hi > -1e-4) res = hi;
  if (lo > 0.0 && w0 + rd.z * lo > -1e-4) res = lo;
  return res;
}

// A plain studio: dim floor, bright sky, one soft window where the light is.
vec3 env(vec3 d, vec3 ld){
  vec3 c = mix(vec3(0.16, 0.15, 0.14), vec3(0.95, 0.92, 0.86), smoothstep(-0.3, 0.9, d.z));
  return c + vec3(1.3, 1.2, 1.0) * pow(max(dot(d, ld), 0.0), 40.0);
}

// Cell colour for the stained glass. A cell's NUMBER changes at the atan cut and
// again every loop -- both renumber cells by a lattice step -- so hashing the
// number would put a seam down the picture and a jump at the loop. Hash its
// class modulo those steps instead.
float cellHash(vec2 c){
  if (gFlat > 0.5) return h21(vec2(mod(c.x, FLAT_SHIFT), c.y) * 0.37 + 5.0);
  float p = floor(gK * gN + 0.5), q = floor(gK * gM + 0.5);
  float D = abs(p * q + gM * gN);
  float A = mod(-q * c.x - gM * c.y, D);
  float B = mod(-gN * c.x + p * c.y, D);
  return h21(vec2(A, B) * 0.37 + 5.0);
}

void main(){
  vec2 R = RENDERSIZE.xy;
  vec2 p = (gl_FragCoord.xy - 0.5 * R) / (0.5 * R.y);

  // ---- constants for this frame ---------------------------------------------
  gPh = fract(TIME / uPeriod);
  gM = latM(); gN = latN(); gH = TAU / gN; gK = gM / gN;
  gFlat = uGeometry == 2 ? 1.0 : 0.0;
  if (uGeometry == 1) gK = 0.0;
  gSlope = uDepth * (1.0 - gFlat);
  gCs = sqrt(1.0 + gSlope * gSlope);
  gHK = gH / sqrt(1.0 + gK * gK);
  gShift = gPh * gM * gH;                        // one full zoom repeat per loop
  gSpinA = gPh * TAU * float(uSpin);             // whole turns only
  gWave = mix(1.0 / gM, 1.0 / FLAT_SHIFT, gFlat);
  gNF = floor(uPeriod * 0.7 + 0.5);              // wingbeats per loop, a whole number
  gNR = floor(uPeriod / 4.5 + 0.5);              // rising waves per loop
  gBirdOn = uCreature == 2 ? 0.0 : 1.0;
  gFishOn = uCreature == 1 ? 0.0 : 1.0;
  gCreatureOn = uCreature == 3 ? 0.0 : 1.0;

  // ---- camera ---------------------------------------------------------------
  float aTop = (gSlope + gHK * uRelief * LIFT_MAX) * (1.0 - gFlat);
  float bTop = gFlat * uRelief * LIFT_MAX / FLAT_S;
  // never let the camera sink below the funnel's rim line: it has to see the centre
  float el = radians(max(uTilt, degrees(atan(aTop * 1.4))));
  float orbit = 0.0;
  if (uOrbit == 1) orbit = 0.30 * sin(TAU * gPh);
  if (uOrbit == 2) orbit = TAU * gPh;
  if (uOrbit == 3) orbit = -TAU * gPh;
  float az = -0.25 * TAU + orbit;
  float dist = 1.9 * uLens;                       // longer lens, further back: same framing
  vec3 ro = dist * vec3(cos(el) * cos(az), cos(el) * sin(az), sin(el));
  vec3 look = vec3((uCentre - 0.5) * 2.4, 0.0);
  vec3 ww = normalize(look - ro);
  vec3 uu = normalize(cross(ww, vec3(0.0, 0.0, 1.0)));
  vec3 vv = cross(uu, ww);
  vec3 rd = normalize(p.x * uu + p.y * vv + uLens * ww);

  float lAz = TAU * uLight;
  vec3 ld = normalize(vec3(cos(lAz) * 0.77, sin(lAz) * 0.77, 0.64));

  // ---- trace, shadow and coarse normal: ONE loop ----------------------------------
  // Link time on Windows is paid per inlined COPY of the surface function, not per
  // iteration -- a loop body compiles once however many times it runs. So the march
  // (clip to the relief shell, 16 steps, 6 bisections), the soft shadow (6 samples
  // toward the light) and the normal (2 offsets, then the hit itself LAST, which
  // leaves the piece globals describing the hit point) all share one call.
  float tT = hitCone(ro, rd, aTop, bTop);
  float tB = hitCone(ro, rd, gSlope, -1e-3);
  float tE = min(tB, tT + 12.0);
  float dt = (tE - tT) / 16.0;
  float lo = tT, hi = tE, t = tT, tm = tT, tHit = tT;
  bool hit = false;
  vec3 P = ro;
  float rP = 1.0, twP = 1.0, foot = 1.0, eps = 1.0;
  float sh = 1.0, hx = 0.0, hy = 0.0, h0 = 0.0;
  for (int i = 0; i < 31; i++) {
    vec3 at;
    float st = 0.0;
    int k = i - 22;
    if (i < 22) {
      tm = hit ? 0.5 * (lo + hi) : t + dt;
      at = ro + rd * tm;
    } else {
      if (i == 22) {
        tHit = 0.5 * (lo + hi);
        P = ro + rd * tHit;
        rP = length(P.xy);
        twP = mix(rP * gHK, 1.0 / FLAT_S, gFlat);
        foot = tHit * 2.0 / (R.y * uLens);              // one pixel, world units, here
        eps = max(foot, twP * 0.003);
      }
      st = twP * (0.04 + 0.2 * float(k));
      at = k < 6 ? P + ld * st + vec3(0.0, 0.0, twP * 0.01)
                 : vec3(P.xy + (k == 6 ? vec2(eps, 0.0) : (k == 7 ? vec2(0.0, eps) : vec2(0.0))), 0.0);
    }
    float hh = surfZ(at.xy);
    if (i < 22) {
      bool under = at.z < hh;
      if (hit) { if (under) hi = tm; else lo = tm; }
      else if (under) { hit = true; lo = t; hi = tm; }
      else t = tm;
    }
    else if (k < 6) sh = min(sh, clamp(8.0 * (at.z - hh) / st, 0.0, 1.0));
    else if (k == 6) hx = hh;
    else if (k == 7) hy = hh;
    else h0 = hh;
  }
  if (tT > 1e3) hit = false;

  // ---- style ------------------------------------------------------------------
  vec3 cL = vec3(0.86, 0.76, 0.58), cD = vec3(0.07, 0.062, 0.055), cG = vec3(0.50, 0.36, 0.22);
  vec3 bgA = vec3(0.90, 0.86, 0.78), bgB = vec3(0.50, 0.45, 0.38);
  float spL = 0.10, spD = 0.16, shL = 16.0, shD = 30.0, mtL = 0.0;
  if (uStyle == 1) { cL = vec3(0.82, 0.86, 0.88); cD = vec3(0.16, 0.24, 0.33); cG = vec3(0.10, 0.14, 0.18);
                     bgA = vec3(0.78, 0.84, 0.88); bgB = vec3(0.30, 0.38, 0.46); spL = 0.2; spD = 0.3; }
  if (uStyle == 2) { cL = vec3(0.80, 0.78, 0.74); cD = vec3(0.34, 0.33, 0.32); cG = vec3(0.20, 0.19, 0.18);
                     bgA = vec3(0.86, 0.85, 0.82); bgB = vec3(0.42, 0.41, 0.40); spL = 0.05; spD = 0.08; }
  if (uStyle == 3) { cL = vec3(0.90, 0.72, 0.42); cD = vec3(0.56, 0.19, 0.11); cG = vec3(0.14, 0.09, 0.07);
                     bgA = vec3(0.92, 0.84, 0.66); bgB = vec3(0.45, 0.30, 0.20); spL = 0.45; spD = 0.45; shL = 70.0; shD = 70.0; }
  if (uStyle == 4) { cL = vec3(1.00, 0.76, 0.33); cD = vec3(0.025, 0.025, 0.03); cG = vec3(0.04, 0.035, 0.03);
                     bgA = vec3(0.16, 0.14, 0.12); bgB = vec3(0.02); spD = 0.7; shD = 120.0; mtL = 1.0; }
  if (uStyle == 5) { cL = vec3(0.94, 0.94, 0.91); cD = vec3(0.10, 0.20, 0.56); cG = vec3(0.08, 0.14, 0.40);
                     bgA = vec3(0.95, 0.95, 0.93); bgB = vec3(0.60, 0.64, 0.72); spL = 0.8; spD = 0.8; shL = 140.0; shD = 140.0; }
  if (uStyle == 6) { cL = vec3(0.22, 0.40, 0.68); cD = vec3(0.13, 0.27, 0.54); cG = vec3(0.08, 0.17, 0.36);
                     bgA = vec3(0.16, 0.30, 0.56); bgB = vec3(0.06, 0.13, 0.30); spL = 0.1; spD = 0.1; }
  if (uStyle == 8) { cL = vec3(0.05, 0.05, 0.07); cD = vec3(0.02, 0.02, 0.03); cG = vec3(0.0);
                     bgA = vec3(0.03, 0.02, 0.05); bgB = vec3(0.0); spL = 0.6; spD = 0.6; shL = 80.0; shD = 80.0; }
  if (uStyle == 9) {
    float s = smoothstep(-0.6, 0.6, p.x);
    vec3 a = vec3(0.92, 0.88, 0.78), b = vec3(0.08, 0.10, 0.16);
    cL = mix(a, b, s); cD = mix(b, a, s); cG = mix(vec3(0.35, 0.30, 0.25), vec3(0.03), s);
    bgA = mix(vec3(0.90, 0.86, 0.78), vec3(0.10, 0.12, 0.20), s); bgB = bgA * 0.5;
    spL = 0.3; spD = 0.3; shL = 40.0; shD = 40.0;
  }
  vec3 bg = mix(bgB, bgA, smoothstep(-1.0, 1.0, p.y));

  vec3 col = bg;
  if (hit) {
    float pixT = foot / max(twP, 1e-6);               // one pixel, tile widths
    gFade = 1.0 - smoothstep(0.02, 0.07, pixT);
    gFine = 1.0 - smoothstep(0.006, 0.02, pixT);
    float tiny = smoothstep(0.08, 0.30, pixT);         // pieces too small to draw

    // The carving's slope is taken in the piece's own coordinates -- no need to re-run
    // the map for it -- and carried to the world through the Jacobian of the map,
    // which is known in closed form because the map is conformal and the warp is two
    // one-dimensional curves.
    float eL = max(pixT, 0.004);
    vec3 cv = carve3(gL, eL);                          // +x, +y, and the hit: gEye etc. end up here
    vec2 dC = vec2(cv.x - cv.z, cv.y - cv.z) / eL;
    float r2 = max(rP * rP, 1e-8);
    vec2 qn = vec2(-P.y, P.x);
    vec2 dTx = mix((gCs * P.xy + gK * qn) / (gH * r2), vec2(FLAT_S, 0.0), gFlat);
    vec2 dTy = mix((qn - gK * gCs * P.xy) / (gH * r2), vec2(0.0, FLAT_S), gFlat);
    vec2 dLx = dTx + gA * gFy * dTy;
    vec2 dLy = gA * gGy * dTx + dTy;
    vec2 grad = vec2(hx - h0, hy - h0) / eps + twP * uRelief * gFade * (dC.x * dLx + dC.y * dLy);
    vec3 nrm = normalize(vec3(-grad, 1.0));
    vec3 nBase = normalize(vec3(-gSlope * P.xy / max(rP, 1e-4), 1.0));
    nrm = normalize(mix(nrm, nBase, tiny));

    // ---- material ------------------------------------------------------------
    float light = 1.0 - gPar;
    float jewel = cellHash(gCell);
    vec3 emi = vec3(0.0);
    // fully the groove colour AT the seam: the two pieces either side have different
    // colours, and any of that difference left at the seam draws a jagged pixel line
    float groove = 1.0 - smoothstep(0.0, max(BEV * 0.4, pixT * 1.5), gDE);
    if (uStyle == 7) {
      vec3 j = 0.5 + 0.45 * cos(TAU * (jewel + vec3(0.0, 0.33, 0.67)));
      cL = mix(j, vec3(0.95, 0.85, 0.55), 0.2); cD = j * 0.55; cG = vec3(0.06);
      spL = 0.9; spD = 0.9; shL = 160.0; shD = 160.0;
      emi = mix(cD, cL, light) * 0.55 * (1.0 - groove);           // lit from behind
    }
    vec3 alb = mix(cD, cL, light);
    float spK = mix(spD, spL, light), shin = mix(shD, shL, light), met = mtL * light;
    if (uStyle == 0) alb *= 1.0 - 0.07 * light * gFine * (0.5 + 0.5 * sin(gL.y * 55.0 + 3.0 * sin(gL.x * 7.0)));
    if (uStyle == 2) alb *= 0.9 + 0.2 * h21(floor(gL * 60.0) + gCell) * gFine;
    // grooves and carved cuts fill with the groove colour, eyes get a white and a pupil
    alb = mix(alb, cG, max(groove, gCav * 0.7));
    alb = mix(alb, mix(cL, vec3(0.95), 0.5), gEye * gFade * (1.0 - light) * 0.8);
    alb = mix(alb, vec3(0.02), gPupil * gFade);
    if (uStyle == 6)
      emi = vec3(0.85, 0.93, 1.0) * max(1.0 - smoothstep(pixT, pixT * 2.5 + 0.004, gDE), gCav * 0.6) * (1.0 - tiny);
    if (uStyle == 8)
      emi = mix(vec3(1.0, 0.25, 0.75), vec3(0.15, 0.95, 0.90), light)
          * (1.6 * (1.0 - smoothstep(0.0, BEV * 0.5, gDE)) + 0.25 * exp(-gDE * 14.0) + 0.8 * gCav);
    vec3 avgAlb = mix(mix(cD, cL, 0.5), cG, 0.25);
    alb = mix(alb, avgAlb, tiny);
    emi *= 1.0 - 0.5 * tiny;

    // ---- light ----------------------------------------------------------------
    float e = clamp(gDE / BEV, 0.0, 1.0);
    float ao = (0.68 + 0.32 * e) * (0.82 + 0.18 * (0.3 + 0.7 * gFig)) * (1.0 - 0.35 * gCav);
    ao = mix(1.0, mix(ao, 0.8, tiny), uShadow);
    float shade = mix(1.0, sh, uShadow);
    vec3 V = -rd;
    float ndl = max(dot(nrm, ld), 0.0);
    float spec = pow(max(dot(nrm, normalize(ld + V)), 0.0), shin) * (shin + 8.0) / 25.0;
    float fres = 0.04 + 0.96 * pow(1.0 - max(dot(nrm, V), 0.0), 5.0);
    vec3 envc = env(reflect(rd, nrm), ld);
    vec3 diel = alb * (ndl * shade * vec3(1.05, 1.0, 0.92) + (0.55 + 0.45 * nrm.z) * ao * vec3(0.30, 0.32, 0.36))
              + spK * (spec * shade * vec3(1.0, 0.97, 0.9) + envc * fres * ao * 0.5);
    vec3 metal = alb * (envc * 0.75 * ao + ndl * shade * 0.35 + spec * shade * 0.8);
    col = mix(diel, metal, met) + emi;

    // the funnel darkens as it falls away toward the centre
    col *= mix(1.0, 0.2 + 0.8 * smoothstep(0.0, 0.9, rP), min(uDepth, 1.0) * (1.0 - gFlat));
    col = mix(col, bg, 1.0 - exp(-max(tHit - dist * 0.8, 0.0) * 0.05));
  }

  // ---- gain, vignette ----------------------------------------------------------
  col *= uGain;
  vec2 vg = p * vec2(R.y / R.x, 1.0);
  col *= 1.0 - uVign * smoothstep(0.35, 1.45, dot(vg, vg));
  col += (h21(gl_FragCoord.xy) - 0.5) * 0.012;      // still, so the loop still closes
  gl_FragColor = vec4(clamp(col, 0.0, 1.0), 1.0);
}