/*{
  "DESCRIPTION": "Neon feedback fireworks, from the year 3000. A shell rises, decelerates to a dead stop, hangs -- then opens three times. First a perfect great circle, usually seen edge-on. Then the circle splits into two counter-rotating tilted rings and wires itself together with chords into a star polygon that holds for a beat. Then the whole figure unfolds along its own axis into a drifting helix and lets go. The physics is confidently not ours: sparks hang weightless long after they should have fallen, then fall along curves. Two persistent float targets hold the sky's memory -- a crisp trail and a slow luminous residue -- and both are advected every frame by the bursts' OWN velocity field, so each new shell shoves the wake of the last one around and re-lights whatever is still hanging there. Colour is driven by AGE and by which stage a spark is in, never by brightness, so a trail changes hue along its own lifetime instead of collapsing to one colour.",
  "CREDIT": "",
  "CATEGORIES": ["neon", "feedback", "fireworks", "particles", "sci-fi"],
  "INPUTS": [
    { "NAME": "cadence",   "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.35, "MAX": 2.5 },
    { "NAME": "order",     "TYPE": "long",  "DEFAULT": 5,    "MIN": 1,    "MAX": 7   },
    { "NAME": "wireframe", "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 1.5 },
    { "NAME": "gravity",   "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.5 },
    { "NAME": "hangtime",  "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.2,  "MAX": 2.5 },
    { "NAME": "memory",    "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.25, "MAX": 3.0 },
    { "NAME": "residue",   "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.5 },
    { "NAME": "wake",      "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.5 },
    { "NAME": "huedrift",  "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.0 },
    { "NAME": "aether",    "TYPE": "float", "DEFAULT": 0.45, "MIN": 0.0,  "MAX": 1.0 },
    { "NAME": "exposure",  "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.4,  "MAX": 2.5 }
  ],
  "PASSES": [
    { "TARGET": "sky",  "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "haze", "PERSISTENT": true, "FLOAT": true },
    { }
  ]
}*/

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

// ---------------------------------------------------------------------------
// pass 0 -> "sky"   RGBA16F persistent.  rgb = CRISP trail (fast decay ~0.5 s).
// pass 1 -> "haze"  RGBA16F persistent.  rgb = SLOW luminous residue (~5 s),
//                   fed by whatever "sky" is showing, advected harder, blurred.
//                   It is nearly black until a later burst re-lights it.
// pass 2 -> screen: composite, aether lattice, hue-preserving tone-map.
// Both buffers carry their own CLOCK in alpha: 1.0 + mod(TIME,1.0).
//
// Host facts this is built around (all of them cost somebody a build round):
//  * FRAMEINDEX is derived from the requested time, so it is never 0 on the
//    render path -> "uninitialised" is detected from alpha instead.
//  * TIMEDELTA lies on the render path (a constant 1/60 while TIME steps by
//    t/120) -> the real dt is recovered from the stored clock.
//  * the clock is stored mod 1 because RGBA16F cannot resolve one frame's
//    increment on top of a TIME near 100.
//  * every render target is NEAREST -> advection hand-rolls its bilinear read.
// ---------------------------------------------------------------------------

#define TAU 6.2831853
#define NSHELL 6
#define NP 16
#define NPf 16.0

float h11(float n){ return fract(sin(n*57.31 + 11.7)*43758.5453); }
vec3  h31(float n){ return fract(sin(n*vec3(127.1,311.7,74.7) + vec3(1.7,9.2,4.3))*43758.5453); }
float h21(vec2 q){ return fract(sin(dot(q, vec2(12.9898,78.233)))*43758.5453); }

float gau(float x){ return exp(-x*x); }
float lamp(float d, float w){ float x = d/w; return exp(-x*x); }

float sdSeg(vec2 p, vec2 a, vec2 b){
  vec2 pa = p-a, ba = b-a;
  float h = clamp(dot(pa,ba)/max(dot(ba,ba), 1e-7), 0.0, 1.0);
  return length(pa - ba*h);
}

mat3 rotX(float a){ float c=cos(a), s=sin(a); return mat3(1.,0.,0., 0.,c,s, 0.,-s,c); }
mat3 rotY(float a){ float c=cos(a), s=sin(a); return mat3(c,0.,-s, 0.,1.,0., s,0.,c); }

// Cohesive neon ramp keyed to a hue PHASE. Saturation is pushed hard so that
// two bursts of different age never average into grey.
vec3 neon(float h){
  vec3 c = 0.5 + 0.5*cos(TAU*(h + vec3(0.0, 0.33, 0.67)));
  c = c*c*(3.0-2.0*c);
  return c*vec3(1.02, 0.96, 1.08);
}

// --------------------------- per-shell globals -----------------------------
vec2  gC;      // burst centre, p-units
mat3  gM;      // figure orientation
float gHueS;
float gTR, gT0, gT1, gT2, gLife;
float gV0, gDrag, gSpin, gTiltB, gGrav, gHang, gSwirl;
float gFire, gW, gDim, gRmax, gNP;

// Fills the globals for staggered shell slot s; returns that shell's age.
float shellSetup(float s){
  float P  = 1.30/max(cadence, 0.05);
  float sl = floor(TIME/P) - s;
  vec3  a  = h31(sl + 11.0);
  vec3  b  = h31(sl + 57.0);
  float c1 = h11(sl + 3.0);

  gFire = step(0.18, c1);                       // some slots stay dark
  float t0 = sl*P + a.x*P*0.5;

  float ar = min(RENDERSIZE.x/RENDERSIZE.y, 1.9);
  // sizes biased small so the sky has DEPTH: many far tight bursts, and only
  // occasionally one that owns the frame.
  gRmax = 0.048 + b.y*b.y*0.150;
  // Keep the whole figure in frame: the margin has to grow with the burst, or
  // the biggest and most structured shells are the ones that get cropped.
  float m  = 0.17 + 0.90*gRmax;
  float hx = max(0.5*ar - m, 0.04);
  float hy = max(0.5     - m, 0.04);
  gC    = vec2((a.y-0.5)*2.0*hx, mix(-hy, hy, 0.30 + 0.62*b.x));
  gDrag = 2.0 + b.z*1.4;
  gV0   = gRmax*gDrag;
  gW    = 0.0024 + 0.0085*gRmax;                // far bursts draw thinner
  // Brightness falls off with the SQUARE of the apparent size, like distance
  // does. A linear falloff leaves far shells bright enough that their own
  // trail fills their (small) disc solid, and they read as flat coloured
  // balls instead of faint sparks a long way off.
  gDim  = 0.08 + 21.7*gRmax*gRmax;

  gSpin  = (a.z-0.5)*3.2;
  gTiltB = 0.35 + b.x*0.75;
  gGrav  = gravity*(0.020 + 0.030*a.x);
  gHang  = hangtime*(0.65 + 0.70*b.z);
  gSwirl = (c1-0.5)*0.55;

  float hidx = floor(fract(a.z*13.7)*4.0);
  gHueS = fract(0.55 + hidx*0.25 + (a.y-0.5)*0.05);

  gTR   = 0.75 + a.x*0.35;
  gT0   = gTR + 0.22 + b.x*0.22;
  gT1   = gT0 + 0.85 + b.y*0.45;
  gT2   = gT1 + 0.90 + b.z*0.50;
  gLife = gT2 + 2.30;
  // Only the big near shells run the full three-act. A distant one is too far
  // away to resolve a lattice, so it stays a rotating ring and dies -- which
  // makes the wire-up a real event rather than something that happens
  // constantly all over the sky.
  float big = step(0.105, gRmax);
  gT1 = mix(300.0, gT1, big);
  gT2 = mix(300.0, gT2, big);
  // A far shell also carries FEWER sparks: at that radius 16 of them are closer
  // together than their own stroke width and merge into one saturated disc.
  gNP = mix(8.0, NPf, big);

  float age = TIME - t0;
  gM = rotY(a.z*TAU + 0.35*age)*rotX(0.25 + b.z*1.05);
  return age;
}

// Screen position of spark idx of the current shell at absolute age.
vec2 sparkScreen(float idx, float age, out float bri){
  float te = max(age - gT0, 0.0);
  float w1 = smoothstep(gT1, gT1+0.38, age);
  float w2 = smoothstep(gT2, gT2+0.50, age);
  float rr = gV0*(1.0 - exp(-gDrag*te))/gDrag;

  float a  = TAU*idx/gNP;
  float ca = cos(a), sa = sin(a);

  // A -- a great circle. Too perfect to be an explosion, which is the point.
  vec3 pA = vec3(ca, sa, 0.0)*rr;

  // B -- it splits into two counter-rotating tilted rings.
  float side = 1.0 - 2.0*mod(idx, 2.0);
  float bb = a + gSpin*te*side;
  float ct = cos(gTiltB), st = sin(gTiltB);
  vec3 pB = vec3(cos(bb), sin(bb)*ct, sin(bb)*st*side)*rr*1.06;

  // C -- and unfolds along its own axis instead of radially.
  float u  = idx/gNP - 0.5 + 0.5/gNP;
  float gg = a*0.5 + 2.4*te;
  // radii are staggered per spark: a helix whose sparks share one radius packs
  // them all into a few texels and clips to a flat disc
  float rc = rr*(0.34 + 0.30*fract(idx*0.618));
  vec3 pC = vec3(cos(gg), sin(gg), 0.0)*rc + vec3(0.0, 0.0, u*rr*2.8);

  vec3 wp = gM*mix(mix(pA, pB, w1), pC, w2);

  // Gravity that is not ours: a dead weightless hang, then a curved fall.
  float tf = max(te - gHang, 0.0);
  wp.y -= gGrav*tf*tf;
  wp.x += gSwirl*sin(1.9*te + a)*tf*0.6;

  float pz = 1.0/(1.0 - 0.32*clamp(wp.z, -2.0, 1.4));
  bri = pz*pz;
  return gC + wp.xy*pz;
}

// --------- the sky's velocity field: the bursts' own expanding fronts -------
vec2 velField(vec2 p){
  // a slow thermal so the sky breathes even when nothing is firing
  vec2 v = vec2(0.012*sin(p.y*3.1 + TIME*0.13), 0.010 + 0.008*cos(p.x*2.7 - TIME*0.11));
  for (int s=0; s<NSHELL; s++){
    float age = shellSetup(float(s));
    if (gFire < 0.5) continue;
    if (age < gT0 || age > gLife) continue;
    float te = age - gT0;
    float rr = gV0*(1.0 - exp(-gDrag*te))/gDrag;
    vec2  dv = p - gC;
    float rd = length(dv) + 1e-4;
    // the pressure front, plus the slower spark front behind it
    // Keep both fronts NARROW and short-lived. A wide or long-lived tangential
    // push winds the residue into a permanent vortex within a few seconds --
    // the classic way feedback turns to soup.
    float q1 = (rd - (0.20 + 1.2*gRmax)*(1.0 - exp(-2.6*te)))/0.10;
    float q2 = (rd - rr)/0.09;
    float en = exp(-te*1.6);
    float sp = gV0*exp(-gDrag*te) + 0.35*en;
    vec2  rn = dv/rd;
    float ff = 0.28*exp(-q1*q1)*exp(-te*3.2) + exp(-q2*q2)*en;
    // rn is a SINGULARITY at the burst centre: an arbitrary direction at full
    // magnitude. Left alone it blows the bright hang point outward along every
    // direction at once and, because the trail blends with max(), the buffer
    // keeps the whole history -- which is a solid filled disc, not a trail.
    // Outflow must vanish at r = 0 by symmetry anyway.
    ff *= smoothstep(0.0, 0.05, rd);
    v += (rn*sp*0.30 + vec2(-rn.y, rn.x)*sp*0.05*sign(gSpin))*ff;
  }
  return v;
}

// ------------------------- this frame's emission ---------------------------
void emit(vec2 p, float dt, out vec3 core, out vec3 halo){
  core = vec3(0.0); halo = vec3(0.0);
  float KK = float(order);

  for (int s=0; s<NSHELL; s++){
    float age = shellSetup(float(s));
    if (gFire < 0.5) continue;
    if (age < 0.0 || age > gLife) continue;

    float te = max(age - gT0, 0.0);
    float rr = gV0*(1.0 - exp(-gDrag*te))/gDrag;
    vec2  dv = p - gC;
    float rd = length(dv) + 1e-4;
    vec3  hb = neon(gHueS);

    // staging flashes -- also what re-lights the residue in pass 2
    float fl = gau((age-gT0)/0.10) + 0.80*gau((age-gT1)/0.12) + 0.70*gau((age-gT2)/0.15);
    // hot and SMALL plus a wide faint lobe. One broad bright lobe clips to a
    // flat featureless disc at the instant of the break.
    halo += hb*fl*gDim*(0.60/(1.0 + 1600.0*rd*rd) + 0.10/(1.0 + 190.0*rd*rd));

    // NOTE: there is deliberately NO emissive pressure ring here. An expanding
    // ring drawn into an accumulating buffer is max-blended in at every radius
    // it crosses, so it never reads as a moving ring: it leaves a filled disc
    // while it is small and a hanging soap bubble once it is big, and both own
    // the frame. The pressure front survives where it actually earns its keep,
    // in velField() -- shoving the residue outward without drawing anything.


    // ------------------------------ rise + hang ------------------------------
    if (age < gT0){
      float u  = clamp(age/gTR, 0.0, 1.0);
      float u0 = clamp((age - dt*1.15)/gTR, 0.0, 1.0);
      vec2  lp = vec2(gC.x*0.92, -0.70);
      vec2  A  = mix(lp, gC, 1.0-(1.0-u0)*(1.0-u0));
      vec2  B  = mix(lp, gC, 1.0-(1.0-u )*(1.0-u ));
      // Deposit the ascent DIM. The residue inherits whatever brightness the
      // head had when it passed, so a bright ascent segment leaves a hard laser
      // stripe across the whole frame that the wake then smears sideways.
      float d  = sdSeg(p, A, B);
      float ld = 0.45 + 0.75*gDim;
      core += hb*ld*0.75*lamp(d, 0.0034);
      halo += hb*ld*0.05*lamp(d, 0.014);
      core += hb*ld*1.0*(1.0-u*0.5)*lamp(length(p-B), 0.0042);
      float hh = smoothstep(gTR-0.10, gT0, age);
      float dh = length(p - gC);
      core += hb*ld*(0.2 + 1.25*hh)*lamp(dh, 0.0055);
      halo += hb*hh*0.30*gDim*lamp(dh, 0.016);
      continue;
    }

    // cheap bound: past this radius the spark loop has nothing to add
    if (rd > rr*1.55 + 0.10) continue;

    float w1   = smoothstep(gT1, gT1+0.38, age);
    float w2   = smoothstep(gT2, gT2+0.50, age);
    // keyed off gLife, NOT gT2 -- a shell that never reaches stage B has gT2
    // parked past the end of time and a gT2-keyed fade would blank it entirely
    float fade = 1.0 - smoothstep(gLife-1.60, gLife, age);
    // Hue is driven by AGE and by which stage we are in -- never by density.
    float hA   = gHueS + huedrift*(0.030*min(te,3.0) + 0.055*w1 + 0.095*w2);
    vec3  cA   = neon(hA);
    vec3  cB   = neon(hA + 0.085);
    // The lattice is a MOMENT, not the whole life: it wires up at the second
    // break, holds, and is gone by the time the figure unfolds.
    float chord = w1*(1.0 - w2)*wireframe*(1.0 + 1.4*gau((age-gT1-0.30)/0.45));

    // Sparks resolve out of the fireball rather than existing from frame one:
    // at the break they are all coincident, and drawn at full strength they sum
    // to ~40 in one texel, which clips to a flat featureless disc.
    float spread = smoothstep(0.18, 0.55, 1.0 - exp(-gDrag*te));
    vec3 sc = vec3(0.0), sh = vec3(0.0);

    for (int i=0; i<NP; i++){
      float fi = float(i);
      if (fi >= gNP) break;      // also keeps the driver from fully unrolling
      float b1, b0, b2;
      vec2  P1 = sparkScreen(fi, age, b1);
      vec2  P0 = sparkScreen(fi, age - dt*1.15, b0);
      float d  = sdSeg(p, P0, P1);
      vec3  cc = mix(cA, cB, fract(fi*0.618));
      float tw = 0.72 + 0.45*sin(TIME*9.0 + fi*7.31 + gHueS*61.0);
      float amp = min(b1, 2.2)*fade*gDim;
      sc += cc*amp*(1.10*lamp(d, gW) + 0.26*lamp(d, gW*3.1));
      sc += cc*amp*tw*1.35*lamp(length(p-P1), gW*0.85);   // the glittering head
      // The wide halo is what FILLS a burst's swept area in the trail buffer.
      // On a far shell it is fat compared with the burst's own radius, so the
      // whole thing accumulates into one solid ball. Near shells only.
      sh += cc*amp*0.15*smoothstep(0.05, 0.15, gRmax)*lamp(d, gW*7.5);
      if (chord > 0.012){
        vec2  P2 = sparkScreen(mod(fi + KK, gNP), age, b2);
        float d2 = sdSeg(p, P1, P2);
        sc += cc*amp*chord*(0.26*lamp(d2, gW*0.90));
        sh += cc*amp*chord*0.06*lamp(d2, gW*3.8);
      }
    }
    // Per-shell soft knee, applied BEFORE the sum ever reaches the display
    // tone-map. The display knee is 1-exp(-l), which is already flat white by
    // l=3, so a sum of coincident samples is past the shoulder before the
    // tone-map sees it. This one has a hard ceiling (T+S) and leaves anything
    // below T untouched, so single sparks keep their brightness and only piles
    // get compressed. A ring seen edge-on stacks every spark onto one line and
    // needs this just as much as the fireball at the break does.
    float m  = max(max(sc.r, sc.g), sc.b);
    float e  = max(m - 0.95, 0.0);
    core += sc*spread*((min(m, 0.95) + e/(1.0 + e/1.00))/max(m, 1e-4));
    float mh = max(max(sh.r, sh.g), sh.b);
    float eh = max(mh - 0.55, 0.0);
    halo += sh*spread*((min(mh, 0.55) + eh/(1.0 + eh/0.70))/max(mh, 1e-4));
  }
}

// -------- what is burning right now, as a light source for the residue -----
vec3 relight(vec2 p){
  vec3 L = vec3(0.0);
  for (int s=0; s<NSHELL; s++){
    float age = shellSetup(float(s));
    if (gFire < 0.5) continue;
    if (age < 0.0 || age > gLife) continue;
    float fl = 0.70*gau((age-gT0)/0.16) + 0.55*gau((age-gT1)/0.20) + 0.50*gau((age-gT2)/0.24);
    fl += 0.26*smoothstep(gT0, gT0+0.5, age)*(1.0 - smoothstep(gLife-1.6, gLife, age));
    float d = length(p - gC);
    L += neon(gHueS)*fl*gDim/(1.0 + 20.0*d*d);
  }
  return L;
}

// hand-rolled bilinear -- every render target is allocated NEAREST
vec3 bilSky(vec2 uv){
  vec2 ts = 1.0/RENDERSIZE, lo = ts*0.5, hi = 1.0 - ts*0.5;
  vec2 px = uv*RENDERSIZE - 0.5;
  vec2 f  = fract(px);
  vec2 b  = (floor(px) + 0.5)*ts;
  vec3 a00 = IMG_NORM_PIXEL(sky, clamp(b,                 lo, hi)).rgb;
  vec3 a10 = IMG_NORM_PIXEL(sky, clamp(b+vec2(ts.x,0.0),  lo, hi)).rgb;
  vec3 a01 = IMG_NORM_PIXEL(sky, clamp(b+vec2(0.0,ts.y),  lo, hi)).rgb;
  vec3 a11 = IMG_NORM_PIXEL(sky, clamp(b+ts,              lo, hi)).rgb;
  return mix(mix(a00,a10,f.x), mix(a01,a11,f.x), f.y);
}
vec3 bilHaze(vec2 uv){
  vec2 ts = 1.0/RENDERSIZE, lo = ts*0.5, hi = 1.0 - ts*0.5;
  vec2 px = uv*RENDERSIZE - 0.5;
  vec2 f  = fract(px);
  vec2 b  = (floor(px) + 0.5)*ts;
  vec3 a00 = IMG_NORM_PIXEL(haze, clamp(b,                lo, hi)).rgb;
  vec3 a10 = IMG_NORM_PIXEL(haze, clamp(b+vec2(ts.x,0.0), lo, hi)).rgb;
  vec3 a01 = IMG_NORM_PIXEL(haze, clamp(b+vec2(0.0,ts.y), lo, hi)).rgb;
  vec3 a11 = IMG_NORM_PIXEL(haze, clamp(b+ts,             lo, hi)).rgb;
  return mix(mix(a00,a10,f.x), mix(a01,a11,f.x), f.y);
}

// a standing lattice in the air. Invisible until residue drifts through it.
float aetherField(vec2 q){
  float k = 95.0;
  return ( cos(dot(q, vec2( 1.0, 0.0   ))*k)
         + cos(dot(q, vec2( 0.5, 0.8660))*k)
         + cos(dot(q, vec2(-0.5, 0.8660))*k) )/3.0;
}

// recover the true timestep from a clock stored in a buffer's alpha
float clockDt(float pa, float cur, out float fresh){
  fresh = step(pa, 0.5);                       // 1.0 when the buffer is cleared
  float dt = mod(cur - (pa - 1.0) + 1.0, 1.0);
  return mix(clamp(dt, 0.0004, 0.30), 1.0/60.0, fresh);
}

void main(){
  vec2 uv  = isf_FragNormCoord;
  vec2 p   = (gl_FragCoord.xy - 0.5*RENDERSIZE)/RENDERSIZE.y;
  vec2 asp = RENDERSIZE.y/RENDERSIZE;
  float cur = mod(TIME, 1.0);

  // ------------------------------- pass 0: trail --------------------------
  if (PASSINDEX == 0){
    float fresh;
    float dt = clockDt(IMG_THIS_PIXEL(sky).a, cur, fresh);

    vec3 core, halo;
    emit(p, dt, core, halo);
    vec2 vel = velField(p);

    // CFL-ish guard: never move the field more than a few texels in one step,
    // so a long frame (or the poster path's fat dt) cannot teleport it.
    vec2 duv = -vel*dt*wake*asp;
    duv *= min(1.0, 4.0/max(length(duv*RENDERSIZE), 1e-5));
    vec3 prev = bilSky(clamp(uv + duv, vec2(0.0), vec2(1.0)));

    // Multiplicative decay only shrinks a smear; the subtractive floor is what
    // ERODES it, so faint wash reaches true black instead of hanging forever.
    float lum = max(max(prev.r, prev.g), prev.b);
    prev *= exp(-dt*(1.90/max(memory,0.05) + 0.30*lum));
    prev  = max(prev - dt*0.34, vec3(0.0));

    vec3 outc = max(prev, core) + halo*dt*0.55;
    outc = mix(outc, core, fresh);
    vec3 hot = max(outc - 1.5, vec3(0.0));
    outc = min(outc, vec3(1.5)) + hot/(1.0 + 0.5*hot);
    gl_FragColor = vec4(clamp(outc, 0.0, 8.0), 1.0 + cur);
    return;
  }

  // ------------------------ pass 1: the slow residue ----------------------
  if (PASSINDEX == 1){
    float fresh;
    float dt = clockDt(IMG_THIS_PIXEL(haze).a, cur, fresh);

    // the residue is free-floating, so it rides the wake harder than the trail
    vec2 vel = velField(p);
    vec2 dvv = -vel*dt*wake*1.7*asp;
    dvv *= min(1.0, 6.0/max(length(dvv*RENDERSIZE), 1e-5));
    vec2 src = clamp(uv + dvv, vec2(0.0), vec2(1.0));
    vec3 hz  = bilHaze(src);

    // soften it: the residue is a cloud, not a wire
    vec2 ts = 2.2/RENDERSIZE;
    vec3 bl = ( IMG_NORM_PIXEL(haze, clamp(src+vec2( ts.x,0.0), vec2(0.0), vec2(1.0))).rgb
              + IMG_NORM_PIXEL(haze, clamp(src+vec2(-ts.x,0.0), vec2(0.0), vec2(1.0))).rgb
              + IMG_NORM_PIXEL(haze, clamp(src+vec2(0.0, ts.y), vec2(0.0), vec2(1.0))).rgb
              + IMG_NORM_PIXEL(haze, clamp(src+vec2(0.0,-ts.y), vec2(0.0), vec2(1.0))).rgb )*0.25;
    hz = mix(hz, bl, clamp(dt*10.0, 0.0, 0.30));

    hz *= exp(-dt*0.55);
    hz  = max(hz - dt*0.012, vec3(0.0));
    hz += IMG_THIS_PIXEL(sky).rgb*dt*0.45;      // fed by whatever just burned

    hz = mix(hz, vec3(0.0), fresh);
    vec3 hot = max(hz - 1.2, vec3(0.0));
    hz = min(hz, vec3(1.2)) + hot/(1.0 + 0.8*hot);
    gl_FragColor = vec4(clamp(hz, 0.0, 6.0), 1.0 + cur);
    return;
  }

  // ------------------------------ pass 2: screen --------------------------
  vec3 tr = IMG_THIS_PIXEL(sky).rgb;
  vec3 hz = IMG_THIS_PIXEL(haze).rgb;
  vec3 L  = relight(p);

  float lat = mix(1.0, 0.36 + 1.32*pow(0.5 + 0.5*aetherField(p), 1.3), aether);

  // Residue accumulated over many bursts averages toward grey. Pull the
  // saturation back up and tint what is NOT currently re-lit cold, so the
  // hanging field reads as cooling plasma rather than smoke.
  float hzl = dot(hz, vec3(0.30, 0.5, 0.2));
  hz = max(mix(vec3(hzl), hz, 1.45), vec3(0.0));

  vec3 col = tr*(1.0 + 0.9*L)
           + hz*residue*(0.26*vec3(0.30,0.55,1.10) + 1.15*L)*lat
           + L*0.018;

  col += mix(vec3(0.012,0.016,0.034), vec3(0.001,0.002,0.009), smoothstep(0.0,0.8,uv.y));

  col *= 1.50*exposure;
  col *= 1.0 - 0.50*clamp(pow(length(p)*0.92, 2.4), 0.0, 1.0);

  // hue-preserving tone-map: only the very hottest cores desaturate to white
  float l = max(max(col.r, col.g), col.b);
  col = mix(col, vec3(l), smoothstep(3.0, 12.0, l)*0.50);
  col *= (1.0 - exp(-l))/max(l, 1e-4);
  col = pow(max(col, 0.0), vec3(0.88));
  col += (h21(gl_FragCoord.xy) - 0.5)*(1.6/255.0);

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