/*{
  "DESCRIPTION": "A supersaturated vapour chamber. Charged particles rip through it and leave trails of condensation that bloom, fatten, drift and evaporate. The vocabulary is physical: fat dead-straight alpha stubs bursting from the needle source with their density piling up at the end of range; thin betas that wander and spiral tighter as they shed energy in the magnetic field; decay forks where a neutral gap opens and one track becomes two arms curving apart; and, rarely, a cosmic-ray shower fanning across the plate. Arrivals run on a jittered Poisson-ish clock, so it is quiet most of the time and the fogged plate is the picture until something comes. A persistent float target holds the droplet density: deposit along segments, diffuse, drift on the convection, evaporate.",
  "CREDIT": "",
  "CATEGORIES": ["physics", "simulation", "monochrome", "feedback"],
  "INPUTS": [
    { "NAME": "rate",       "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.05, "MAX": 3.0 },
    { "NAME": "field",      "TYPE": "float",   "DEFAULT": 1.0,  "MIN": -2.0, "MAX": 2.0 },
    { "NAME": "vapour",     "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.2,  "MAX": 2.0 },
    { "NAME": "diffusion",  "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.5 },
    { "NAME": "evaporate",  "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.25, "MAX": 3.0 },
    { "NAME": "drift",      "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 3.0 },
    { "NAME": "exposure",   "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.3,  "MAX": 2.5 },
    { "NAME": "plate",      "TYPE": "float",   "DEFAULT": 0.28, "MIN": 0.0,  "MAX": 1.0 },
    { "NAME": "grain",      "TYPE": "float",   "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.0 },
    { "NAME": "source",     "TYPE": "point2D", "DEFAULT": [0.34, 0.36] }
  ],
  "PASSES": [
    { "TARGET": "vap", "PERSISTENT": true, "FLOAT": true },
    { }
  ]
}*/

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

// ---------------------------------------------------------------------------
// Cloud Chamber
//
// pass 0 -> "vap"  RGBA16F persistent
//             r = droplet density (far tracks are deposited already defocused)
//             g = fresh/hot core  (fast decay: the head that is passing now)
//             b = scattered halo  (born wide, slow decay: the lit vapour)
//             a = CLOCK. Uniform across the buffer, holds mod(TIME,1) so the
//                 next frame can recover the true timestep.
// pass 1 -> screen: light sheet, plate tone, procedural fog, duotone, grain.
//
// -- host facts this shader is built around ---------------------------------
// * TIMEDELTA CANNOT BE TRUSTED. Measured on this renderer's poster path it is
//   a constant 1/60 while TIME steps by (--time)/120 per step, i.e. it is wrong
//   by 4x at --time 8. Using it made every track a string of beads (the deposit
//   covered a quarter of the arc the particle actually crossed) and left the
//   field aged 2 s no matter how long the schedule had run. So the real dt is
//   recovered from a clock the buffer carries in its own alpha channel.
// * FRAMEINDEX is derived from TIME, so it is never consulted; the field starts
//   at zero and zero is a correct empty chamber, so there is nothing to seed.
// * Every target is NEAREST, so the advection read is hand-rolled bilinear.
// * The update is dt-invariant by construction: deposits use max() (idempotent),
//   decay uses exp(-k*dt), and the diffusion stencil RADIUS grows as sqrt(dt) at
//   a fixed mix factor, which holds the diffusion coefficient constant for any
//   timestep and any resolution.
// ---------------------------------------------------------------------------

#define NEV  8       // arrival slots examined per frame
#define NBR  5       // max branches per event (primary + 4 secondaries)
#define NSUB 4       // sub-segments deposited per branch per frame
#define NARC 6       // constant-curvature pieces per track (spirals need these)
#define TAU  6.2831853

#define TAPI(ij) IMG_NORM_PIXEL(vap, ((ij) + 0.5) / RENDERSIZE)

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

vec4 hash41(float p){
  vec4 q = fract(vec4(p) * 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 hash21(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);
}

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

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-9), 0.0, 1.0);
  return length(pa - ba*h);
}

// Position + heading at arc length s along a chain of constant-curvature arcs.
// Curvature is multiplied by kg each piece (kg > 1 is the electron spiral, its
// radius shrinking as it loses energy) and perturbed per piece by wob, which is
// multiple scattering: the reason a beta wanders and an alpha does not.
void trackAt(float s, vec2 p0, float th0, float k0, float kg, float wob,
             float sd, float segL, out vec2 pos, out float th){
  pos = p0; th = th0;
  float kc = k0;
  float rem = s;
  for (int m = 0; m < NARC; m++){
    float take = min(segL, rem);
    if (take <= 0.0) break;
    float kk = kc + wob*(hash11(sd + float(m)*7.31)*2.0 - 1.0);
    float ak = max(abs(kk), 3e-3);            // huge-radius arc stands in for straight
    kk = kk < 0.0 ? -ak : ak;
    float th1 = th + kk*take;
    pos += vec2(sin(th1) - sin(th), cos(th) - cos(th1)) / kk;
    th = th1;
    rem -= take;
    kc *= kg;
  }
  pos += vec2(cos(th), sin(th)) * max(rem, 0.0);
}

void main(){
  vec2 res = RENDERSIZE;
  float asp = res.x / max(res.y, 1.0);
  vec2 P = (gl_FragCoord.xy - 0.5*res) / max(res.y, 1.0);   // y in [-.5,.5], square units

  // =========================================================================
  // PASS 0 - advance the droplet field
  // =========================================================================
  if (PASSINDEX == 0){
    vec4 self = IMG_THIS_PIXEL(vap);

    // --- recover the real timestep from our own clock ----------------------
    float clk = mod(TIME, 1.0);
    float dt = clk - self.a;
    if (dt < 0.0) dt += 1.0;                  // the clock wrapped
    dt = clamp(dt, 0.0002, 0.22);

    // --- slow convection of the vapour (semi-Lagrangian back-trace) --------
    vec2 velF = vec2(0.0, -0.030);
    velF += 0.026*vec2(sin(P.y*4.3 + TIME*0.19 + 1.1), cos(P.x*3.6 - TIME*0.15));
    velF *= drift;
    vec2 spx = gl_FragCoord.xy - velF*res.y*dt;

    vec2 g  = spx - 0.5;
    vec2 gi = floor(g);
    vec2 gf = g - gi;
    vec4 t00 = TAPI(gi);
    vec4 t10 = TAPI(gi + vec2(1.0, 0.0));
    vec4 t01 = TAPI(gi + vec2(0.0, 1.0));
    vec4 t11 = TAPI(gi + vec2(1.0, 1.0));
    vec3 st  = mix(mix(t00, t10, gf.x), mix(t01, t11, gf.x), gf.y).rgb;

    // --- diffusion: droplets growing and the trail spreading ---------------
    // The stencil is an 8-point RING, not a 4-point cross. A sparse cross does
    // not blur a thin line, it REPLICATES it into a bead lattice.
    float Dm  = 5.5e-5 * diffusion;
    float oxf = clamp(floor(sqrt(4.0*Dm*dt/0.42)*res.y + 0.5), 1.0, 5.0);
    float hA  = oxf / res.y;
    vec3  Dch = vec3(1.4e-5, 1.0e-5, 5.5e-5) * diffusion;
    vec3  aCh = min(4.0*Dch*dt/(hA*hA), vec3(0.55));

    vec2 ci = floor(spx);
    float od = oxf*0.7071;
    vec3 nb = TAPI(ci + vec2( oxf, 0.0)).rgb + TAPI(ci + vec2(-oxf, 0.0)).rgb
            + TAPI(ci + vec2(0.0,  oxf)).rgb + TAPI(ci + vec2(0.0, -oxf)).rgb
            + TAPI(ci + vec2( od,  od)).rgb  + TAPI(ci + vec2(-od,  od)).rgb
            + TAPI(ci + vec2( od, -od)).rgb  + TAPI(ci + vec2(-od, -od)).rgb;
    st = mix(st, nb*0.125, aCh);

    // --- evaporation -------------------------------------------------------
    st *= exp(-vec3(0.30, 2.20, 0.58)*evaporate*dt);
    st = max(st - vec3(0.0030, 0.020, 0.0016)*evaporate*dt, 0.0);

    // =======================================================================
    // arrivals
    // =======================================================================
    float rr    = clamp(rate, 0.05, 3.0);
    float slot  = 1.75 / clamp(rr, 0.25, 1.6);          // seconds per lattice slot
    float skipP = mix(0.62, 0.04, clamp(rr/3.0, 0.0, 1.0));
    float nNow  = floor(TIME / slot);

    vec2 srcP = source - 0.5;
    srcP.x *= asp;

    vec3 dep = vec3(0.0);

    for (int i = 0; i < NEV; i++){
      float n = nNow - float(i);
      if (n < -0.5) break;
      vec4 H = hash41(n*1.1943 + 0.66);
      if (H.x < skipP && n > 1.5) continue;             // the quiet stretches

      float t0  = (n + 0.12 + 0.72*H.y) * slot;
      float age = TIME - t0;
      if (age <= 0.0 || age > 3.0) continue;            // outside the depositing window

      vec4 G = hash41(n*1.9330 + 71.3);                 // species parameters
      vec4 Q = hash41(n*3.7710 + 19.1);                 // geometry

      // ---- species ------------------------------------------------------
      float nbr   = 1.0;      // branch count
      float len   = 2.1;      // arc length of the primary
      float vel0  = 0.95;     // units / second
      float tau   = 0.0;      // >0 : decelerating, s = len*(1-exp(-age/tau))
      float wid   = 0.0032;
      float amp   = 0.60;
      float kB    = 0.30;     // curvature, rad per unit
      float kG    = 1.0;      // curvature growth per piece (the spiral)
      float wobA  = 0.05;     // multiple scattering
      float dpth  = 0.26 + 0.68*Q.z;
      float bLen  = 0.0;      // secondary length
      float sv    = 0.0;      // vertex arc
      float fan   = 0.0;
      float gap   = 0.0;      // neutral (invisible) flight before a fork
      float kArm  = 0.0;
      float bragg = 0.0;      // density piling up at end of range
      float fadeE = 0.0;      // the track thinning out as energy is spent
      float dly   = 0.0;      // stagger between burst members
      float fSrc  = 0.0;      // launched from the needle source

      float ty = H.z;
      if (ty < 0.44){                       // ---- muon: long, thin, near-straight
        len = 2.2; vel0 = 0.72 + 0.55*G.x; wid = 0.0038 + 0.0022*G.y;
        amp = 0.80 + 0.45*G.z; kB = 0.20 + 0.55*G.w; wobA = 0.05;
      } else if (ty < 0.605){               // ---- alpha burst: short, fat, dead straight
        nbr = 2.0 + floor(G.x*2.0);
        fSrc = 1.0;
        len = 0.150 + 0.070*G.y; vel0 = 0.52; wid = 0.0072 + 0.0028*G.z;
        amp = 1.05; kB = 0.02; wobA = 0.015; bragg = 1.0;
        dly = 0.24 + 0.30*G.w; dpth = 0.62 + 0.34*Q.z;
      } else if (ty < 0.785){               // ---- beta: thin, wandering, curls tight
        // it comes in on a wide arc and winds up as energy goes; the last coil
        // is the densest, because a slow electron ionizes hardest
        len = 0.95 + 0.65*G.x; tau = 0.70; wid = 0.0028 + 0.0013*G.y;
        amp = 0.68; kB = 1.9 + 1.3*G.z; kG = 1.56; wobA = 0.50; bragg = 0.55;
        dpth = 0.24 + 0.58*Q.z;
      } else if (ty < 0.90){                // ---- decay fork: neutral gap, then a V
        nbr = 3.0; len = 0.50 + 0.35*G.x; vel0 = 0.85;
        wid = 0.0044; amp = 1.00; kB = 0.14; wobA = 0.04;
        sv = len; gap = 0.055 + 0.090*G.y; bLen = 0.62 + 0.50*G.z;
        fan = 0.28 + 0.20*G.w; kArm = 0.60;
      } else {                              // ---- cosmic-ray shower
        nbr = 5.0; len = 0.30 + 0.22*G.x; vel0 = 1.25;
        wid = 0.0038; amp = 0.90; kB = 0.10; wobA = 0.03;
        sv = len; bLen = 1.55; fan = 0.40 + 0.28*G.y; kArm = 0.26;
        dpth = 0.30 + 0.62*Q.z;
      }

      // depth: far events are smaller, slower, dimmer and defocused
      float persp = mix(0.68, 1.0, dpth);
      float defoc = 1.0 + 2.4*(1.0 - dpth);
      len  *= persp;  vel0 *= persp;  bLen *= persp;  sv *= persp;  gap *= persp;
      wid  *= persp;

      float chg = Q.w < 0.5 ? -1.0 : 1.0;   // sign of the charge: which way it bends
      kB   *= chg * field / persp;
      kArm *= field / persp;

      // ---- entry ---------------------------------------------------------
      float ea = TAU * H.w;
      vec2 p0 = vec2((0.5*asp + 0.12)*cos(ea), 0.62*sin(ea)) * persp;
      float th0 = atan(-p0.y, -p0.x) + (Q.x - 0.5)*1.30;
      if (fSrc > 0.5){ p0 = srcP; th0 = TAU*Q.x; }

      // ---- cheap rejection -------------------------------------------------
      float reach = max(len, sv + gap + bLen);
      float mgn = wid*defoc*11.0 + 0.05;
      if (distance(P, p0) > reach + mgn) continue;
      if (nbr < 1.5){
        // near-straight single track: reject against the chord plus its sagitta
        float kt = abs(kB) + wobA;
        if (kt*len < 1.2){
          vec2 pe = p0 + vec2(cos(th0), sin(th0))*len;
          if (sdSeg(P, p0, pe) > kt*len*len*0.175 + mgn) continue;
        }
      }

      // ---- fork / shower vertex -------------------------------------------
      vec2 vp = p0; float vth = th0;
      if (nbr > 1.5 && fSrc < 0.5){
        trackAt(sv, p0, th0, kB, kG, wobA, n, sv/float(NARC), vp, vth);
        vp += vec2(cos(vth), sin(vth))*gap;
      }
      float tv = (sv + gap) / max(vel0, 1e-3);

      // ---- branches ---------------------------------------------------------
      for (int b = 0; b < NBR; b++){
        if (float(b) >= nbr) break;

        vec4 Bh = hash41(n*11.317 + float(b)*3.71 + 0.9);
        vec2 bp0; float bth, bk, bl, bw, bam, bage, bwob, bkg, btau, bbr, bfd;
        bkg = 1.0; btau = 0.0; bbr = 0.0; bfd = 0.0;

        if (b == 0){
          bp0 = p0;  bth = th0;  bk = kB;  bl = len;  bw = wid;  bam = amp;
          bage = age; bwob = wobA; bkg = kG; btau = tau; bbr = bragg; bfd = fadeE;
        } else if (fSrc > 0.5){
          // another alpha out of the same speck of source, a moment later
          bp0 = p0;  bth = TAU*Bh.x;
          bk  = kB * (Bh.y < 0.5 ? -1.0 : 1.0);
          bl  = len*(0.72 + 0.55*Bh.z);
          bw  = wid*(0.82 + 0.34*Bh.w);
          bam = amp*(0.80 + 0.35*Bh.x);
          bage = age - dly*float(b);
          bwob = wobA; bbr = bragg;
        } else {
          // secondaries out of the vertex, fanned. A fork's two arms take
          // opposite charge and therefore bend apart.
          float j  = float(b) - 1.0;
          float hb = max(0.5*(nbr - 2.0), 0.5);
          bp0 = vp;
          bth = vth + fan*((j - 0.5*(nbr - 2.0))/hb) + (Bh.x - 0.5)*0.38;
          bk  = kArm * (mod(float(b), 2.0) < 0.5 ? 1.0 : -1.0) * chg;
          bl  = bLen*(0.55 + 0.95*Bh.y);
          bw  = wid*(0.86 + 0.32*Bh.z);
          bam = amp*(0.82 + 0.30*Bh.w);
          bage = age - tv;
          bwob = wobA + 0.25;
          bfd = 0.28;
        }

        if (bage <= 0.0) continue;

        float s1, s0;
        if (btau > 0.0){
          s1 = bl*(1.0 - exp(-bage/btau));
          s0 = bl*(1.0 - exp(-max(bage - dt, 0.0)/btau));
        } else {
          s1 = min(bage*vel0, bl);
          s0 = min(max(bage - dt, 0.0)*vel0, bl);
        }
        if (s1 <= s0 + 1e-7) continue;      // not started yet, or already complete

        float segL = bl / float(NARC);
        float sd = n*3.1 + float(b)*17.7;
        vec2 pa; float tha;
        trackAt(s0, bp0, bth, bk, bkg, bwob, sd, segL, pa, tha);

        for (int k = 0; k < NSUB; k++){
          float sk = mix(s0, s1, float(k + 1)/float(NSUB));
          vec2 pb; float thb;
          trackAt(sk, bp0, bth, bk, bkg, bwob, sd, segL, pb, thb);
          float d = sdSeg(P, pa, pb);
          pa = pb;

          float sMid = mix(s0, s1, (float(k) + 0.5)/float(NSUB));
          float u = sMid / max(bl, 1e-4);
          float endr = smoothstep(0.55, 1.0, u);
          float A = bam * (1.0 + 0.9*bbr*endr) * (1.0 - bfd*smoothstep(0.15, 1.0, u));
          A *= 0.82 + 0.36*hash11(sd*3.7 + floor(sMid*70.0));  // droplet granularity
          float w = max(bw*(1.0 + 0.30*bbr*endr), 2.6/res.y);

          // defocus: a far track is laid down already soft, and dimmer for it
          float wd = w*defoc;
          A *= (0.45 + 0.55*dpth) / sqrt(defoc);

          float q = d/wd;
          float core = A * exp(-min(2.3*q*q, 40.0));
          float hq = d/(wd*6.0);
          // the halo is capped: an alpha is 3x brighter than a muon but must not
          // flood the plate with a bloom 3x as strong
          float hal = min(A, 1.5) * 0.34 * exp(-min(1.5*hq*hq, 40.0));

          dep.r = max(dep.r, core);
          dep.g = max(dep.g, core);
          dep.b = max(dep.b, hal);
        }
      }
    }

    st = max(st, dep*vapour);
    gl_FragColor = vec4(clamp(st, 0.0, 6.0), clk);
    return;
  }

  // =========================================================================
  // PASS 1 - photograph the chamber
  // =========================================================================
  vec2 uv = gl_FragCoord.xy / res;
  vec4 V = IMG_THIS_PIXEL(vap);

  // one wide ring: the halo channel is already born broad, this only veils it
  vec3 veil = vec3(0.0);
  for (int i = 0; i < 8; i++){
    float a = float(i)*0.7854 + 0.31;
    vec2 off = vec2(cos(a), sin(a))*0.040;
    off.x /= max(asp, 0.35);
    veil += IMG_NORM_PIXEL(vap, uv + off).rgb;
  }
  veil *= 0.125;

  // --- illumination: a slab of light across the chamber --------------------
  float sy = P.y + 0.06;
  float sheet = exp(-1.9*sy*sy);
  float lit = mix(1.0, 0.78, smoothstep(-0.62, 0.72, P.x/max(asp, 0.6)));
  float lowGrad = smoothstep(0.50, -0.48, P.y);

  // --- vapour texture: without this a quiet frame is an empty grey square ---
  vec2 fp = P*vec2(1.7, 3.4);          // vapour stratifies in horizontal layers
  float fog = 0.50*vnoise(fp + vec2(TIME*0.013, -TIME*0.021))
            + 0.29*vnoise(fp*2.3 + vec2(-TIME*0.030, TIME*0.018))
            + 0.21*vnoise(fp*5.3 + vec2(TIME*0.045, TIME*0.037));

  float base = 0.090 + 0.084*lit*sheet + 0.050*lowGrad*sheet
             + (fog - 0.48)*0.225*sheet*lit;

  // --- droplets -------------------------------------------------------------
  float hal = V.b + veil.b*0.85 + veil.r*0.30;
  float lum = base
            + V.r*1.85*sheet*lit
            + hal*0.78*lit
            + V.g*0.60*lit;

  lum = 1.0 - exp(-lum*1.55*exposure);
  lum = mix(lum, lum*lum*(3.0 - 2.0*lum), 0.32);

  // --- grain: live emulsion plus a fixed speck of plate ---------------------
  float g1 = hash21(gl_FragCoord.xy + vec2(fract(TIME*13.0)*331.0, fract(TIME*7.3)*197.0));
  float g2 = hash21(floor(gl_FragCoord.xy*0.5) + 11.3);
  float ga = grain*(0.022 + 0.040*(1.0 - lum));
  lum += (g1 - 0.5)*ga + (g2 - 0.5)*ga*0.40;
  lum = clamp(lum, 0.0, 1.0);

  // --- plate tone: cool blue-grey <-> warm bromide, never a hue ramp --------
  vec3 shC = mix(vec3(0.050, 0.078, 0.132), vec3(0.118, 0.090, 0.064), plate);
  vec3 mdC = mix(vec3(0.352, 0.432, 0.535), vec3(0.532, 0.462, 0.352), plate);
  vec3 hiC = mix(vec3(0.955, 0.962, 0.955), vec3(1.000, 0.955, 0.858), plate);
  vec3 col = lum < 0.5 ? mix(shC, mdC, lum*2.0) : mix(mdC, hiC, (lum - 0.5)*2.0);

  col *= 1.0 - 0.22*clamp(dot(P, P)*1.9, 0.0, 1.0);

  gl_FragColor = vec4(col, 1.0);
}