/*{
  "DESCRIPTION": "Streamer III: ruled-strip silk in slow-motion divergence-free turbulence. Anisotropic drag lines each ribbon up with its own path so the tail retraces where the head has been; the tips steer as a bounded-turn-rate flock; and no hard constraint remains in the position solve, so nothing grips and releases between frames.",
  "INPUTS": [
    { "NAME": "speed", "TYPE": "float", "DEFAULT": 0.78, "MIN": 0.0, "MAX": 1.6 },
    { "NAME": "ribbonWidth", "TYPE": "float", "DEFAULT": 0.054, "MIN": 0.020, "MAX": 0.095 },
    { "NAME": "stiffness", "TYPE": "float", "DEFAULT": 0.64, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "wind", "TYPE": "float", "DEFAULT": 0.72, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "mess", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "twist", "TYPE": "float", "DEFAULT": 0.68, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "depth", "TYPE": "float", "DEFAULT": 0.76, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "sheen", "TYPE": "float", "DEFAULT": 0.74, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "seed", "TYPE": "float", "DEFAULT": 23.0, "MIN": 0.0, "MAX": 30.0 }
  ],
  "PASSES": [
    { "TARGET": "spines", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "spineDepth", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "ribbonLayer", "FLOAT": true },
    {}
  ]
}*/

// A conforming GLSL ES 1.00 host declares `precision mediump float;` for
// fragment shaders, since ES 1.00 gives fragment stages no default float
// precision and highp is optional there. mediump guarantees only about 10 bits
// of mantissa. This shader integrates node positions incrementally across
// thousands of frames in a persistent float target, so mediump would quantize
// the spine and reintroduce exactly the frame-to-frame stutter the motion work
// removed. Ask for highp where the hardware advertises it; the guard is
// required because referencing highp without it is an error on parts that lack
// it. Desktop drivers usually promote mediump to fp32 anyway, so the absence of
// this line would likely go unnoticed until the piece ran on a phone.
#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI 3.14159265359
#define RIBBONS 7
#define NODES 28

float hash11(float p){return fract(sin(p*127.1+seed*29.17)*43758.5453);}

vec3 spectrum(float x){
  x=clamp(x,0.0,1.0)*6.0;
  vec3 r=vec3(0.98,0.018,0.055),o=vec3(1.0,0.26,0.008),y=vec3(1.0,0.76,0.008);
  vec3 g=vec3(0.008,0.73,0.18),b=vec3(0.008,0.29,0.99),i=vec3(0.15,0.045,0.61),v=vec3(0.67,0.035,0.87);
  if(x<1.0)return mix(r,o,smoothstep(0.0,1.0,x));
  if(x<2.0)return mix(o,y,smoothstep(1.0,2.0,x));
  if(x<3.0)return mix(y,g,smoothstep(2.0,3.0,x));
  if(x<4.0)return mix(g,b,smoothstep(3.0,4.0,x));
  if(x<5.0)return mix(b,i,smoothstep(4.0,5.0,x));
  return mix(i,v,smoothstep(5.0,6.0,x));
}

vec2 nodeUv(int node,int ribbon){
  return (vec2(float(node),float(ribbon))+0.5)/RENDERSIZE;
}

// Limiting a 2D quantity has to bound its magnitude, not its components.
// clamp() on a vec2 bounds it to a square, so a saturating value gets pushed
// out to the corners of that box and direction snaps onto 45 degree diagonals
// - which is what was making the tails trace fast diamond shapes. Bounding the
// length leaves direction untouched. Used for velocities and for steering.
vec2 limitMagnitude(vec2 value,float limit){
  float magnitude=length(value);
  return magnitude>limit?value*(limit/max(magnitude,0.0001)):value;
}

// Starting layout only. The heads are free-flying boids once the simulation
// is running, so all this has to do is spread seven of them around a rosette
// at roughly the golden angle.
vec2 headSeedPosition(float fi){
  float angle=fi*2.39+0.6;
  return 0.5+vec2(cos(angle),sin(angle))*(0.135+0.085*hash11(fi+17.0));
}

// One octave of a divergence-free field: a unit vector perpendicular to the
// wavevector, modulated by cos(k.p). Because that vector is orthogonal to k,
// the divergence is identically zero however the amplitude is scaled, so
// octaves can be summed and weighted freely and the result still swirls
// without piling material into sinks. Rotating the wavevector over time keeps
// the pattern reorganising.
vec2 curlOctave(vec2 p,float angle,float freq,float phase,float amplitude){
  vec2 direction=vec2(cos(angle),sin(angle));
  return vec2(-direction.y,direction.x)*(cos(dot(direction,p)*freq+phase)*amplitude);
}

// Shed vortices were the wrong medium for a constant flurry: a travelling bump
// with a sin(life*PI) envelope means every point in the frame spends most of
// its time between vortices, in dead air, which is exactly the lull that reads
// as stop-and-go. A spectrum of always-on octaves has shear everywhere at all
// times. The octaves beat against each other instead of switching on and off,
// so the field reorganises continuously without ever going quiet.
// `detail` low-passes the field in space. The fine octaves are what make a
// free end buzz: eddies a fraction of a ribbon long shake it at a rate the eye
// reads as noise rather than as flow. Fading them out along the chain leaves
// the tail riding broad currents only, so it drapes and sweeps instead of
// flailing, while the head keeps the full spectrum and stays lively.
vec2 windForce(vec2 p,float fi,float t,float detail){
  // Seven spines sharing one field move as a rigid sheet. Sampling each ribbon
  // at its own offset decorrelates them while keeping every ribbon's own flow
  // divergence-free and coherent along its length.
  p+=vec2(hash11(fi+31.0),hash11(fi+41.0))*mess*0.42;
  // The largest cell is deliberately about one ribbon long. Make it wider than
  // the frame and the whole cast sits inside a single eddy and orbits it
  // together, which drifts everything into one corner even though the field is
  // divergence-free - zero divergence forbids compression, not advection.
  // Slow phase rates. Turbulence photographed in slow motion still has all its
  // spatial structure - the eddies are the same size - it just evolves gently,
  // so the large cells are what move the silk and the small ones only texture
  // it. Speeding the phases up is what turns weather into vibration.
  vec2 flow=curlOctave(p,0.70+t*0.05,5.2,t*0.19,1.00);
  flow+=curlOctave(p,2.31-t*0.07,8.6,-t*0.29,0.55*detail);
  flow+=curlOctave(p,4.12+t*0.09,13.5,t*0.44,0.26*detail*detail);
  flow+=curlOctave(p,5.63-t*0.11,19.0,-t*0.60,0.10*detail*detail*detail);
  return flow*2.0;
}

vec2 quadraticBezier(vec2 p0,vec2 p1,vec2 p2,float u){
  vec2 a=mix(p0,p1,u);
  vec2 b=mix(p1,p2,u);
  return mix(a,b,u);
}

vec2 quadraticBezierD(vec2 p0,vec2 p1,vec2 p2,float u){
  return 2.0*mix(p1-p0,p2-p1,u);
}

vec2 quadraticBezierDD(vec2 p0,vec2 p1,vec2 p2){
  return 2.0*(p0-2.0*p1+p2);
}

float quadraticBezier1(float p0,float p1,float p2,float u){
  return mix(mix(p0,p1,u),mix(p1,p2,u),u);
}

float quadraticBezierD1(float p0,float p1,float p2,float u){
  return 2.0*mix(p1-p0,p2-p1,u);
}

vec2 projectDepth(vec2 point,float z){
  return 0.5+(point-0.5)*(1.0+z*0.20);
}

void considerQuadraticRoot(
  vec2 point,vec2 p0,vec2 p1,vec2 p2,float root,
  inout float bestDistanceSquared,inout float bestU
){
  float u=clamp(root,0.0,1.0);
  vec2 delta=quadraticBezier(p0,p1,p2,u)-point;
  float distanceSquared=dot(delta,delta);
  if(distanceSquared<bestDistanceSquared){
    bestDistanceSquared=distanceSquared;
    bestU=u;
  }
}

void refineQuadraticRoot(
  vec2 point,vec2 p0,vec2 p1,vec2 p2,float seed,
  inout float bestDistanceSquared,inout float bestU
){
  float u=clamp(seed,0.0,1.0);
  for(int solve=0;solve<4;solve++){
    vec2 curve=quadraticBezier(p0,p1,p2,u);
    vec2 firstDerivative=quadraticBezierD(p0,p1,p2,u);
    vec2 secondDerivative=quadraticBezierDD(p0,p1,p2);
    float numerator=dot(curve-point,firstDerivative);
    float denominator=dot(firstDerivative,firstDerivative)+dot(curve-point,secondDerivative);
    float safeDenominator=(denominator<0.0?-1.0:1.0)*max(abs(denominator),0.0000001);
    u=clamp(u-clamp(numerator/safeDenominator,-0.34,0.34),0.0,1.0);
  }
  considerQuadraticRoot(point,p0,p1,p2,u,bestDistanceSquared,bestU);
}

// Foot of the normal: the u whose ruling passes through this pixel, solving
// dot(p-B(u),B'(u))=0. Endpoints are deliberately NOT offered as candidates.
// An endpoint "hit" is a round cap, and unioning caps across sections is what
// turns a flat strip into a sweep of disks. Short convex pieces rarely have
// competing minima, so two safeguarded seeds cover both basins.
float footOfNormal(vec2 point,vec2 p0,vec2 p1,vec2 p2){
  vec2 chord=p2-p0;
  float chordSeed=dot(point-p0,chord)/max(dot(chord,chord),0.00000001);
  float bestDistanceSquared=1000000.0;
  float bestU=0.5;
  refineQuadraticRoot(point,p0,p1,p2,chordSeed,bestDistanceSquared,bestU);
  refineQuadraticRoot(point,p0,p1,p2,chordSeed<0.5?0.78:0.22,bestDistanceSquared,bestU);
  return bestU;
}

void updateSpines(vec2 uv,float t){
  ivec2 cell=ivec2(floor(uv*RENDERSIZE));
  int node=cell.x;
  int ribbon=cell.y;
  if(node>=NODES||ribbon>=RIBBONS){gl_FragColor=vec4(0.0);return;}
  float fi=float(ribbon);
  float dt=min(TIMEDELTA,1.0/30.0);
  vec2 head=headSeedPosition(fi);
  float segmentLength=0.036+0.004*sin(fi*1.7);
  vec2 initialDirection=normalize(vec2(cos(fi*2.17+0.4),sin(fi*1.63+0.8)));
  float cruise=0.10+0.17*speed;
  // Seed on an empty buffer, not on a frame number. FRAMEINDEX is not portable
  // as an "is this the first frame" test: a host that renders a poster at t=8
  // may hand the shader FRAMEINDEX 480 on the very first frame it ever draws,
  // so a FRAMEINDEX<2 branch never runs, the persistent target keeps its
  // cleared contents, and every node reads position (0,0) - which puts the
  // whole piece in the bottom-left corner. That is exactly what glslop's poster
  // renderer did. A cleared float target is exactly zero, and a live node is
  // never at the origin with precisely zero velocity, so the buffer's own
  // contents are the reliable signal. FRAMEINDEX is kept as a fast path.
  vec4 stored=IMG_NORM_PIXEL(spines,nodeUv(node,ribbon));
  if(FRAMEINDEX<2||dot(stored,stored)<=0.0){
    vec2 initial=head-initialDirection*segmentLength*float(node);
    // Seed the tips already at cruise and tangential, so the flock is milling
    // from the first frame instead of accelerating out of a dead start.
    vec2 seedVelocity=vec2(-initialDirection.y,initialDirection.x)*cruise;
    gl_FragColor=vec4(clamp(initial,vec2(-0.35),vec2(1.35)),seedVelocity);
    return;
  }
  if(node==0){
    // The tips are boids. Seven of them is small enough that all-pairs costs
    // six texture reads at one node per ribbon, so no spatial structure is
    // needed. The reason it is worth doing at all is that steering separates
    // heading from pace: speed is renormalised to a fixed cruise every step,
    // so a tip only ever changes direction, never rate. That gets a decided,
    // constant lead without the uniformity of putting every head on a rail.
    vec4 headState=IMG_NORM_PIXEL(spines,nodeUv(0,ribbon));
    vec2 headPos=headState.xy;
    vec2 headVel=headState.zw;

    vec2 separation=vec2(0.0);
    vec2 alignment=vec2(0.0);
    vec2 cohesion=vec2(0.0);
    float neighbours=0.0;
    for(int k=0;k<RIBBONS;k++){
      // Masked rather than broken out of, to stay inside the ES 1.00 loop form.
      float notSelf=step(0.5,abs(float(k)-fi));
      vec4 otherState=IMG_NORM_PIXEL(spines,nodeUv(0,k));
      vec2 offset=headPos-otherState.xy;
      float gap=max(length(offset),0.0001);
      separation+=offset/gap*(notSelf*max(0.42-gap,0.0));
      alignment+=otherState.zw*notSelf;
      cohesion+=otherState.xy*notSelf;
      neighbours+=notSelf;
    }
    alignment/=max(neighbours,1.0);
    cohesion/=max(neighbours,1.0);

    // Steering authority scales with cruise. A path's turn radius is
    // cruise^2/acceleration, so holding these gains fixed while the tips speed
    // up flattens every trajectory - and since the body is towed along the
    // path the head traces, flat paths mean straight ribbons that sprawl
    // across the frame. Scaling keeps the shape of the motion invariant to
    // pace, including when the speed input is moved.
    float authority=cruise/0.35;
    vec2 steer=separation*(13.5*authority);
    // Alignment kept low. Strong alignment is what makes a flock dart as one
    // body, and a synchronised dart is the opposite of weather - it reads as
    // seven ribbons on one string. Separation and the wind do the work; this
    // only stops heads from crossing head-on.
    steer+=(alignment-headVel)*0.45;
    steer+=(cohesion-headPos)*(0.08*authority);
    // Soft containment. Boids in an unbounded field commit to a heading and
    // leave; a turn-back that engages only near the rim keeps them milling in
    // frame without flattening the flock against a wall.
    vec2 fromCentre=headPos-0.5;
    float outward=length(fromCentre);
    steer-=fromCentre/max(outward,0.0001)*smoothstep(0.36,0.54,outward)*(6.0*authority);
    // The wind leads the tips now rather than merely nudging them, so a head
    // is steered by the same weather its body is hanging in.
    steer+=windForce(headPos,fi,t,1.0)*(0.40+0.85*wind)*authority;
    // A slow per-tip bias perpendicular to travel. Once speed is fixed, cross
    // acceleration is pure turning, and this term is written as perp(velocity)
    // so its magnitude is already proportional to cruise - which makes the
    // coefficient an angular rate in rad/s directly, independent of pace. At
    // this rate a tip curls through a radius of roughly a sixth of the frame,
    // so the ribbon behind it lays down loops rather than a straight wake. The
    // rate is well under the wind's, so it sets the size of the loops without
    // competing with the weather for who is steering.
    steer+=vec2(-headVel.y,headVel.x)*sin(t*(0.29+0.11*hash11(fi+3.0))+fi*2.2)*1.25;

    // Bound the steering. This is what classic boids clamp and I had not.
    // Because speed is renormalised below, acceleration is converted entirely
    // into turning, at a rate of |steer|/cruise - so the gains above are really
    // angular rates, and unclamped they are enormous: separation alone reaches
    // about 38 rad/s and containment about 17, either of which reverses a
    // heading in a handful of frames. Most of the time no term is near its
    // peak and the motion is calm, but when a head passes close to another or
    // crosses the containment rim, it snaps to a new direction - and the body,
    // which now faithfully follows the head's path, is dragged around that
    // corner. That is the occasional spaz between long chill stretches. The
    // curl term above never caused it because writing it as perp(velocity)
    // already made its coefficient an angular rate.
    steer=limitMagnitude(steer,cruise*2.2);

    headVel+=steer*dt;
    float pace=length(headVel);
    headVel=pace>0.0001?headVel/pace*cruise:initialDirection*cruise;
    headPos=clamp(headPos+headVel*dt,vec2(-0.30),vec2(1.30));
    gl_FragColor=vec4(headPos,headVel);
    return;
  }

  vec4 self=IMG_NORM_PIXEL(spines,nodeUv(node,ribbon));
  vec2 pos=self.xy;
  vec2 vel=self.zw;
  vec2 prev=IMG_NORM_PIXEL(spines,nodeUv(node-1,ribbon)).xy;
  vec2 toPrev=prev-pos;
  float prevDistance=max(length(toPrev),0.0001);
  float springStrength=34.0+stiffness*54.0;
  vec2 force=toPrev/prevDistance*(prevDistance-segmentLength)*springStrength;
  if(node<NODES-1){
    vec2 next=IMG_NORM_PIXEL(spines,nodeUv(node+1,ribbon)).xy;
    vec2 toNext=next-pos;
    float nextDistance=max(length(toNext),0.0001);
    force+=toNext/nextDistance*(nextDistance-segmentLength)*springStrength*0.62;
    // Light straightening only. Curvature is limited in the position solve
    // below, where it can be enforced hard without being stiff, and the
    // transverse wave is dissipated by drag rather than by bending stiffness.
    force+=(prev+next-2.0*pos)*(5.0+stiffness*7.0);
  }
  float tail=float(node)/float(NODES-1);
  float detail=mix(1.0,0.12,tail);
  // Wind authority is now near-uniform along the chain. It used to more than
  // double from head to tail, which is backwards: a ribbon is smaller than the
  // eddies carrying it, so every part of it sits in much the same flow. Giving
  // the free end the most forcing, the least damping and the largest ripple
  // stacked three amplifiers at exactly the point that already reflects waves.
  force+=windForce(pos,fi,t,detail)*(0.22+0.50*wind)*(0.90+0.10*tail);

  vec2 tangent=toPrev/prevDistance;
  vec2 lateral=vec2(-tangent.y,tangent.x);
  // A slight, slow undulation - the "variance" on top of following, not a
  // driver of it. Kept small deliberately: with anisotropic drag below, any
  // across-path force is strongly resisted, so this now reads as the ribbon
  // breathing rather than as something shaking it.
  float ripple=sin(t*(0.78+0.16*hash11(fi+13.0))-tail*(3.4+2.2*twist)+fi*1.77);
  force+=lateral*ripple*(0.03+0.07*tail*(1.30-tail))*(0.55+0.85*mess);

  vel+=force*dt;

  // Anisotropic drag - the piece that was actually missing. A ribbon sliding
  // along its own length pushes almost no air; moving broadside it pushes a
  // whole face of it. Damping isotropically models a bead on a string, free to
  // swing sideways exactly as easily as it trails, so a small turn at the head
  // launches a transverse wave that travels down the chain, and because the
  // free end has nothing beyond it to carry the wave onward, it reflects there
  // and doubles - the thrash. Resisting across-path motion several times
  // harder than along-path motion dissipates that wave as it travels, which is
  // why a real streamer lines itself up with its own path and the tail simply
  // follows where the head has already been.
  vec2 alongVelocity=tangent*dot(vel,tangent);
  vec2 acrossVelocity=vel-alongVelocity;
  // Both grades now rise toward the free end rather than falling. The far end
  // is where the wave arrives with nowhere to go, so it wants the most
  // absorption, not the least.
  float alongDrag=mix(2.2,2.9,tail)*(1.0+stiffness*0.30);
  float acrossDrag=mix(7.0,17.0,tail)*(1.0+stiffness*0.30);
  vel=alongVelocity*exp(-dt*alongDrag)+acrossVelocity*exp(-dt*acrossDrag);
  vel=limitMagnitude(vel,0.85);
  vec2 predicted=pos+vel*dt;

  // Length constraint, projected the way position-based dynamics requires:
  // move the predicted position onto the segment-length circle around the
  // predecessor, then derive velocity from the displacement that actually
  // happened. Correcting position while leaving the old velocity in place
  // leaves the node still driving into its neighbour, so the chain keeps
  // collapsing and the ribbon renders as a balloon.
  vec2 toPrevSettled=predicted-prev;
  float settledDistance=max(length(toPrevSettled),0.0001);
  vec2 corrected=prev+toPrevSettled/settledDistance*mix(settledDistance,segmentLength,0.72);

  // There was a curvature-limiting projection here, capping the turn per joint
  // so a hairpin could not close tighter than the strip is wide and render as
  // a solid paddle. It has been removed, and its removal is the fix for the
  // stick-slip stutter: a hard `if (turn < limit) snap to limit` grips on one
  // frame and releases the next for any joint sitting near the threshold,
  // which is the same mechanism as a rubber band dragged across a table.
  // Smoothing the engagement did not help - it only spread a partial
  // correction across many more joints and measured worse.
  //
  // It is safe to delete rather than repair because the anisotropic drag added
  // above already prevents the folding it was catching: a fold IS across-path
  // motion, which is now the heavily resisted direction, so the coils never
  // form in the first place. Verified over 100 s at defaults and at maximum
  // wind and mess with zero stiffness and maximum speed - no coils, no bulbs.
  // If an extreme setting ever reintroduces them, the lever is drag, not a
  // constraint.

  vel=limitMagnitude((corrected-pos)/max(dt,0.001),1.05);
  pos=corrected;
  pos=clamp(pos,vec2(-0.30),vec2(1.30));
  gl_FragColor=vec4(pos,vel);
}

void updateSpineDepth(vec2 uv,float t){
  ivec2 cell=ivec2(floor(uv*RENDERSIZE));
  int node=cell.x;
  int ribbon=cell.y;
  if(node>=NODES||ribbon>=RIBBONS){gl_FragColor=vec4(0.0);return;}
  float fi=float(ribbon);
  float along=float(node)/float(NODES-1);
  float dt=min(TIMEDELTA,1.0/30.0);
  float targetHead=0.42*sin(t*0.31+fi*1.37)+0.10*cos(t*0.17-fi*0.77);

  // Cumulative arc length from the head, carried in the spare .z channel. Twist
  // parameterized by node index instead kinks at every joint, because the
  // springs stretch and sections are not equal length; the kink shows directly
  // in a ruled silhouette. This pass is 28x7 pixels, so the prefix sum is free.
  // Masked rather than broken out of, to stay inside the ES 1.00 loop form.
  float arc=0.0;
  vec2 walkPoint=IMG_NORM_PIXEL(spines,nodeUv(0,ribbon)).xy;
  float walkDepth=IMG_NORM_PIXEL(spineDepth,nodeUv(0,ribbon)).x;
  for(int k=1;k<NODES;k++){
    vec2 nextPoint=IMG_NORM_PIXEL(spines,nodeUv(k,ribbon)).xy;
    float nextDepth=IMG_NORM_PIXEL(spineDepth,nodeUv(k,ribbon)).x;
    arc+=step(float(k),float(node))*length(vec3(nextPoint-walkPoint,(nextDepth-walkDepth)*0.32));
    walkPoint=nextPoint;
    walkDepth=nextDepth;
  }

  // Same portability problem as the spine pass, but this buffer cannot use an
  // all-zero test: node 0 legitimately stores arc 0, and its depth and velocity
  // both pass through zero in normal motion, so a zero test would re-seed it
  // mid-flight. The w channel was unused, so it carries an explicit
  // initialised flag instead. Every write below sets it.
  vec4 state=IMG_NORM_PIXEL(spineDepth,nodeUv(node,ribbon));
  if(FRAMEINDEX<2||state.w<0.5){
    float initial=targetHead*0.45+0.20*sin(fi*1.31+along*5.4);
    gl_FragColor=vec4(initial,0.0,arc,1.0);
    return;
  }
  float z=state.x;
  float velocity=state.y;
  if(node==0){
    velocity=clamp((targetHead-z)/max(dt,0.001),-1.4,1.4);
    gl_FragColor=vec4(targetHead,velocity,arc,1.0);
    return;
  }
  float previous=IMG_NORM_PIXEL(spineDepth,nodeUv(node-1,ribbon)).x;
  float force=(previous-z)*(24.0+stiffness*34.0);
  if(node<NODES-1){
    float next=IMG_NORM_PIXEL(spineDepth,nodeUv(node+1,ribbon)).x;
    force+=(next-z)*(13.0+stiffness*18.0);
    force+=(previous+next-2.0*z)*(5.0+stiffness*8.0);
  }
  // Long and slow, matching the in-plane undulation. Thirteen wavelengths of
  // depth wave along a ribbon is a corrugation, not a swell, and it was adding
  // fast apparent motion through the width projection even where the ribbon
  // was barely moving in plane.
  force+=sin(t*(0.34+fi*0.012)-along*(2.8+twist*2.2)+fi*1.11)*(0.20+0.44*mess);
  velocity+=force*dt;
  velocity*=exp(-dt*(3.4+stiffness*1.5));
  velocity=clamp(velocity,-1.15,1.15);
  z=clamp(z+velocity*dt,-0.62,0.62);
  gl_FragColor=vec4(z,velocity,arc,1.0);
}

vec3 shadeCloth(
  float fi,float facing,float edge,vec3 normal,float along,float sceneDepth,float t
){
  vec3 ink=spectrum(fi/6.0);
  vec3 back=mix(ink,vec3(0.020,0.012,0.038),0.46);
  bool frontFacing=facing>=0.0;
  vec3 cloth=frontFacing?ink:back;
  vec3 lightDirection=normalize(vec3(-0.35,0.55,0.76));
  cloth*=0.68+0.27*abs(dot(normal,lightDirection));
  cloth*=mix(0.72,1.10,smoothstep(-0.75,0.75,sceneDepth));
  float broadGlint=pow(max(0.0,abs(dot(normal,normalize(vec3(0.25,0.28,0.93))))),18.0);
  cloth+=mix(ink,vec3(1.0),0.78)*broadGlint*sheen*0.17;
  vec3 edgeInk=mix(ink,vec3(0.006,0.008,0.018),frontFacing?0.42:0.55);
  cloth=mix(cloth,edgeInk,edge*0.55);
  // Real silk flares as it turns edge-on. Paying the vanishing moment out as a
  // grazing highlight keeps the fold reading as material rather than a hole.
  cloth+=mix(ink,vec3(1.0),0.60)*pow(1.0-abs(facing),3.0)*sheen*0.50;
  cloth*=0.97+0.03*sin(along*18.0-t*1.1+fi);
  return cloth;
}

void renderRibbons(vec2 uv,float t){
  float pixelSize=1.0/min(RENDERSIZE.x,RENDERSIZE.y);
  float frontZ=-1000.0;
  float frontMask=0.0;
  vec3 frontColor=vec3(0.0);
  float backZ=-1000.0;
  float backMask=0.0;
  vec3 backColor=vec3(0.0);
  float shadow=0.0;

  for(int r=0;r<RIBBONS;r++){
    float fi=float(r);
    float ribbonArcLength=max(IMG_NORM_PIXEL(spineDepth,nodeUv(NODES-1,r)).z,0.0001);
    float ribbonScore=-1000.0;
    float ribbonSurfaceDistance=-20.0;
    float ribbonSignedSide=0.0;
    float ribbonBaseWidth=0.0;
    float ribbonEffectiveProjection=1.0;
    float ribbonFacing=1.0;
    float ribbonAlong=0.0;
    float ribbonPhase=0.0;
    vec2 ribbonTangent=vec2(1.0,0.0);
    float ribbonSceneDepth=0.0;
    float ribbonDepthTangent=0.0;
    for(int j=0;j<NODES;j++){
      // Ternaries, not min/max: GLSL ES 1.00 has no integer overload of
      // min/max, so max(j-1,0) silently resolves to the float version and then
      // fails to assign back to an int. Our local renderer compiles as ES 3.00,
      // where the integer overloads do exist, so it accepts this and the error
      // only appears on a conforming ES 1.00 host.
      int previousIndex=j>0?j-1:0;
      int nextIndex=j<NODES-1?j+1:NODES-1;
      vec4 previousState=IMG_NORM_PIXEL(spineDepth,nodeUv(previousIndex,r));
      vec4 controlState=IMG_NORM_PIXEL(spineDepth,nodeUv(j,r));
      vec4 nextState=IMG_NORM_PIXEL(spineDepth,nodeUv(nextIndex,r));
      float previousZ=previousState.x;
      float controlZ=controlState.x;
      float nextZ=nextState.x;
      vec2 previous=projectDepth(IMG_NORM_PIXEL(spines,nodeUv(previousIndex,r)).xy,previousZ);
      vec2 control=projectDepth(IMG_NORM_PIXEL(spines,nodeUv(j,r)).xy,controlZ);
      vec2 next=projectDepth(IMG_NORM_PIXEL(spines,nodeUv(nextIndex,r)).xy,nextZ);

      // Midpoint-to-midpoint quadratic B-spline section. Neighboring pieces
      // share the same endpoint and first derivative exactly. Repeated endpoint
      // controls create clean first/last caps without a special raster path.
      vec2 p0=j==0?control:0.5*(previous+control);
      vec2 p1=control;
      vec2 p2=j==NODES-1?control:0.5*(control+next);
      float z0=j==0?controlZ:0.5*(previousZ+controlZ);
      float z1=controlZ;
      float z2=j==NODES-1?controlZ:0.5*(controlZ+nextZ);
      // Arc length rides the same midpoint construction as position, so the
      // twist phase advances uniformly in space instead of per node index.
      float a0=j==0?controlState.z:0.5*(previousState.z+controlState.z);
      float a1=controlState.z;
      float a2=j==NODES-1?controlState.z:0.5*(controlState.z+nextState.z);
      float boundsPadding=ribbonWidth+0.040;
      vec2 boundsMin=min(min(p0,p1),p2)-vec2(boundsPadding);
      vec2 boundsMax=max(max(p0,p1),p2)+vec2(boundsPadding);
      if(any(lessThan(uv,boundsMin))||any(greaterThan(uv,boundsMax)))continue;

      float u=footOfNormal(uv,p0,p1,p2);
      vec2 nearest=quadraticBezier(p0,p1,p2,u);
      vec2 tangentVector=quadraticBezierD(p0,p1,p2,u);
      if(dot(tangentVector,tangentVector)<0.00000001)tangentVector=p2-p0;
      vec2 tangent=normalize(tangentVector+vec2(0.00001,0.0));
      vec2 normal=vec2(-tangent.y,tangent.x);
      vec2 fromCurve=uv-nearest;

      // Stationarity means the longitudinal residual vanishes; a non-zero one
      // says the solve clamped and this pixel's ruling lies past the piece.
      float overshoot=max(abs(dot(fromCurve,tangent))-1.5*pixelSize,0.0);

      float along=quadraticBezier1(a0,a1,a2,u)/ribbonArcLength;
      float phase=along*(3.2+twist*7.2)-t*(0.76+fi*0.028)+fi*1.71;
      phase+=0.30*sin(along*8.2-t*0.48+fi*0.73);
      float facing=cos(phase);
      float taper=1.0-smoothstep(0.74,1.0,along)*0.55;
      // Perspective has to reach the width too. projectDepth spreads the
      // centreline by the same factor, so leaving width alone drew near ribbons
      // too narrow and far ones too fat.
      float curveDepth=quadraticBezier1(z0,z1,z2,u);
      float baseWidth=ribbonWidth*taper*(1.0+curveDepth*0.20);
      float projection=abs(facing);
      // Projected half extent of a rectangular cross-section: the width edge
      // foreshortens as |cos|, the thickness edge fills in as |sin|. A ribbon
      // turning edge-on becomes a thin square-ended sliver, which is what silk
      // does, instead of collapsing to a resolution-dependent minimum line.
      float width=max(baseWidth*projection+ribbonWidth*0.060*abs(sin(phase)),pixelSize);
      float effectiveProjection=width/max(baseWidth,0.001);

      // Ruled-strip coverage: the lateral material coordinate against the local
      // half width, unioned as square-ended rectangles in each piece's own
      // frame. Testing distance-to-centreline instead unions circles, which
      // inflates the silhouette near width maxima and can never form the sharp
      // crossing where a twisting ribbon passes edge-on. Neighbouring pieces
      // share their end ruling exactly, so rectangles butt together with no gap
      // and no cap bulge, and an interior w-|v| always outranks a neighbour's
      // clamped exterior distance.
      float signedSide=dot(fromCurve,normal);
      float lateralOut=abs(signedSide)-width;
      float surfaceDistance=overshoot>0.0
        ? -length(vec2(overshoot,max(lateralOut,0.0)))
        : -lateralOut;

      float physicalSide=clamp(signedSide/effectiveProjection,-baseWidth,baseWidth);
      float sectionZ=curveDepth*1.65
        +physicalSide/max(baseWidth,0.001)*sin(phase)*0.16*depth;

      // Coverage is the union of every piece; which piece gets *shaded* is a
      // separate question, answered by depth. Taking the deepest-inside sample
      // instead makes a self-crossing show whichever layer happens to be more
      // interior, and lets the winner flip wherever two pieces nearly tie.
      // Interior samples all score above 100, so they always outrank exterior
      // ones and the boundary still finds a frame to shade.
      ribbonSurfaceDistance=max(ribbonSurfaceDistance,surfaceDistance);
      float sectionScore=surfaceDistance>0.0?100.0+sectionZ:surfaceDistance;
      if(sectionScore>ribbonScore){
        ribbonScore=sectionScore;
        ribbonSignedSide=signedSide;
        ribbonBaseWidth=baseWidth;
        ribbonEffectiveProjection=effectiveProjection;
        ribbonFacing=facing;
        ribbonAlong=along;
        ribbonPhase=phase;
        ribbonTangent=tangent;
        ribbonSceneDepth=sectionZ;
        ribbonDepthTangent=quadraticBezierD1(z0,z1,z2,u);
      }
    }

    if(ribbonSurfaceDistance>-10.0){
      // No global cap planes. They were infinite half-space clips, so a ribbon
      // that curled back past its own head or tail plane got sliced straight
      // through mid-body, amputating whole spans. Square ends come free from
      // the rectangle union: the end pieces already terminate flat, normal to
      // the tangent at the end node.
      // A fixed pixel footprint stays stable across the bounding-box branches;
      // screen derivatives of the rejected sentinel would draw the boxes.
      float aa=max(1.15/min(RENDERSIZE.x,RENDERSIZE.y),0.0008);
      float mask=smoothstep(-aa,aa,ribbonSurfaceDistance);
      float width=ribbonBaseWidth*ribbonEffectiveProjection;
      // The trim is a fraction of the *current* half width, never a fixed
      // radius. A fixed band swallows the whole cross-section once a ribbon
      // turns edge-on, which reads as the material fading out instead of
      // narrowing. Cloth has a crisp cut edge, not a shaded falloff.
      float halfWidth=max(width,aa);
      float edge=1.0-smoothstep(0.0,min(0.20*halfWidth,aa*1.6),abs(ribbonSurfaceDistance));

      // Width rotates out of the image plane, so an edge-on face becomes a line
      // while its two physical edges separate in depth rather than swelling.
      vec3 centerTangent=normalize(vec3(ribbonTangent,ribbonDepthTangent*0.32));
      vec3 widthAxis=normalize(
        vec3(-ribbonTangent.y,ribbonTangent.x,sin(ribbonPhase)*0.75)
      );
      vec3 planeNormal=normalize(cross(centerTangent,widthAxis));
      if(ribbonFacing<0.0)planeNormal=-planeNormal;
      shadow=max(shadow,smoothstep(-0.032,0.004,ribbonSurfaceDistance)*0.20);

      if(mask>0.004){
        vec3 cloth=shadeCloth(
          fi,ribbonFacing,edge,planeNormal,ribbonAlong,ribbonSceneDepth,t
        );
        // Keep the two frontmost layers rather than one argmax winner, so a
        // crossing edge antialiases against the ribbon behind it. A single
        // winner picked by depth-plus-coverage lets a fully covered far ribbon
        // beat a half-covered near one, which flips hard along the seam.
        if(ribbonSceneDepth>frontZ){
          backZ=frontZ;backMask=frontMask;backColor=frontColor;
          frontZ=ribbonSceneDepth;frontMask=mask;frontColor=cloth;
        }else if(ribbonSceneDepth>backZ){
          backZ=ribbonSceneDepth;backMask=mask;backColor=cloth;
        }
      }
    }
  }

  float alpha=frontMask+(1.0-frontMask)*backMask;
  if(alpha<=0.0){gl_FragColor=vec4(0.0,0.0,0.0,shadow);return;}
  vec3 premultiplied=frontColor*frontMask+(1.0-frontMask)*backColor*backMask;
  gl_FragColor=vec4(premultiplied/alpha,alpha);
}

void compose(vec2 uv){
  vec4 layer=IMG_NORM_PIXEL(ribbonLayer,uv);
  vec2 px=1.0/RENDERSIZE;
  float softShadow=0.0;
  softShadow=max(softShadow,IMG_NORM_PIXEL(ribbonLayer,uv+px*vec2(-6.0,7.0)).a*0.24);
  softShadow=max(softShadow,IMG_NORM_PIXEL(ribbonLayer,uv+px*vec2(-10.0,12.0)).a*0.13);
  softShadow=max(softShadow,IMG_NORM_PIXEL(ribbonLayer,uv+px*vec2(-15.0,18.0)).a*0.07);
  softShadow=max(softShadow,layer.a*0.0);
  vec3 bg=mix(vec3(0.006,0.008,0.021),vec3(0.035,0.016,0.050),uv.y);
  bg+=vec3(0.025,0.03,0.067)*pow(max(0.0,1.0-length(uv-0.5)*1.42),4.0);
  bg*=1.0-softShadow*(1.0-layer.a);
  vec3 color=mix(bg,layer.rgb,layer.a);
  color=1.0-exp(-max(color,0.0)*1.22);
  gl_FragColor=vec4(pow(color,vec3(0.93)),1.0);
}

void main(){
  vec2 uv=isf_FragNormCoord;
  float t=TIME*speed;
  if(PASSINDEX==0){updateSpines(uv,t);return;}
  if(PASSINDEX==1){updateSpineDepth(uv,t);return;}
  if(PASSINDEX==2){renderRibbons(uv,t);return;}
  compose(uv);
}