Organ

experimental backend WebGL2 · not opted in

Enable WebGPU controls in settings. WebGL2 remains unchanged.

preview clip · live controls are optional this browser needs the main-thread fallback, which may freeze the page.

params (15)

tempo 1.00

density 0.72

wither 1.20

nutrition 1.00

meander 0.62

swirl 1.25

fibrosis 0.45

turgor 0.55

breath 0.70

emissive 1.00

pulse 0.75

translucency 0.85

flesh 0.70

bloom 0.95

seed 7.00

.fs source · isf not running

/*{
  "DESCRIPTION": "Organ: a beating, three-lobed body of translucent flesh, lit from the inside by the vasculature growing through it. The shape is constructive solid geometry on spheres - three lobes fused with a smooth union and two spherical bites scooped out of the result - and turgor separates the lobes, so the silhouette runs from one round blob to a distinctly three-chambered mass. Three colonies drain the shell through a flow-accumulation network, and the light they make has to travel back out through skin: the vessels are read where the refracted ray reaches them rather than at the surface, so the network sits under the skin with real parallax, and the path through that layer goes as one over the cosine, so the limb is a warm scattered wash and a deep trunk arrives as dull red. A systolic front leaves each nucleus once a beat and runs out along the tree, dilating the vessels it passes and raising them out of the skin, while a chamber contracts on the same beat. Territory is a partition of the surface between three drifting nuclei, each carrying a vigour that tracks its share of the nutrients currently on the body. Forked from Vein Atlas.",
  "CATEGORIES": ["generative", "simulation", "organic"],
  "INPUTS": [
    { "NAME": "tempo", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.15, "MAX": 3.0 },
    { "NAME": "density", "TYPE": "float", "DEFAULT": 0.72, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "wither", "TYPE": "float", "DEFAULT": 1.20, "MIN": 0.05, "MAX": 2.0 },
    { "NAME": "nutrition", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "meander", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 2.4 },
    { "NAME": "swirl", "TYPE": "float", "DEFAULT": 1.25, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "fibrosis", "TYPE": "float", "DEFAULT": 0.45, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "turgor", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0, "MAX": 0.9 },
    { "NAME": "breath", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "emissive", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "pulse", "TYPE": "float", "DEFAULT": 0.75, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "translucency", "TYPE": "float", "DEFAULT": 0.85, "MIN": 0.10, "MAX": 2.2 },
    { "NAME": "flesh", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 1.4 },
    { "NAME": "bloom", "TYPE": "float", "DEFAULT": 0.95, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "seed", "TYPE": "float", "DEFAULT": 7.0, "MIN": 0.0, "MAX": 30.0 }
  ],
  "PASSES": [
    { "TARGET": "acc1", "FLOAT": true },
    { "TARGET": "acc", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "net", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "mem", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "glow", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "scene", "FLOAT": true },
    {}
  ]
}*/

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI  3.14159265359
#define TAU 6.28318530718

#define COLONIES 3
#define NUTS 9

#define CLOCK_UNIT 8.0
#define CLOCK_WRAP 1024.0

#define BODY_R 0.60

#define RBOUND 1.17

#define PROGRESS_W 0.28

#define JITTER 0.55

#define SKIN_H 0.052
#define SKIN_IOR 1.42

#define NUT_R 0.20
#define WARM_BOOST 16.0
#define WARM_RATE 0.25

float hash11(float p){ return fract(sin(p*127.1 + seed*19.73)*43758.5453); }
float hash12(vec2 p){ return fract(sin(dot(p, vec2(127.1, 311.7)) + seed*3.17)*43758.5453); }
float hash13(vec3 p){ return fract(sin(dot(p, vec3(127.1, 311.7, 74.7)) + seed*3.17)*43758.5453); }

float vnoise3(vec3 p){
  vec3 i = floor(p), f = fract(p);
  f = f*f*(3.0 - 2.0*f);
  float a = hash13(i + vec3(0.0, 0.0, 0.0));
  float b = hash13(i + vec3(1.0, 0.0, 0.0));
  float c = hash13(i + vec3(0.0, 1.0, 0.0));
  float d = hash13(i + vec3(1.0, 1.0, 0.0));
  float e = hash13(i + vec3(0.0, 0.0, 1.0));
  float g = hash13(i + vec3(1.0, 0.0, 1.0));
  float h = hash13(i + vec3(0.0, 1.0, 1.0));
  float k = hash13(i + vec3(1.0, 1.0, 1.0));
  return mix(mix(mix(a, b, f.x), mix(c, d, f.x), f.y),
             mix(mix(e, g, f.x), mix(h, k, f.x), f.y), f.z);
}

#define ASPECT 3.0
float simEdge(){
  float w = min(576.0, floor(RENDERSIZE.x) - 2.0);
  float h = min(192.0, floor(RENDERSIZE.y) - 2.0);
  return floor(min(w, h*ASPECT));
}
float simRows(float N){ return floor(N/ASPECT); }
vec2 simGrid(float N){ return vec2(N, simRows(N)); }

vec4 fetchTexel(sampler2D tex, vec2 texel){
  return texture2D(tex, (texel + 0.5)/RENDERSIZE);
}

vec3 uvToDir(vec2 uv){
  float phi = (uv.x - 0.5)*TAU;
  float sy = clamp(uv.y*2.0 - 1.0, -1.0, 1.0);
  float cy = sqrt(max(1.0 - sy*sy, 0.0));
  return vec3(cy*cos(phi), sy, cy*sin(phi));
}

vec2 dirToUV(vec3 d){
  return vec2(atan(d.z, d.x)/TAU + 0.5, d.y*0.5 + 0.5);
}

vec4 shellFetch(sampler2D tex, vec2 texel, float N){
  return fetchTexel(tex, vec2(mod(texel.x, N), clamp(texel.y, 0.0, simRows(N) - 1.0)));
}

vec4 shellBilerp(sampler2D tex, vec2 grid, float N){
  vec2 g = grid - 0.5;
  vec2 i = floor(g);
  vec2 f = g - i;
  vec4 a = shellFetch(tex, i, N);
  vec4 b = shellFetch(tex, i + vec2(1.0, 0.0), N);
  vec4 c = shellFetch(tex, i + vec2(0.0, 1.0), N);
  vec4 d = shellFetch(tex, i + vec2(1.0, 1.0), N);
  return mix(mix(a, b, f.x), mix(c, d, f.x), f.y);
}

vec2 dirVec(float d){
  if(d < 0.5) return vec2( 1.0,  0.0);
  if(d < 1.5) return vec2( 1.0,  1.0);
  if(d < 2.5) return vec2( 0.0,  1.0);
  if(d < 3.5) return vec2(-1.0,  1.0);
  if(d < 4.5) return vec2(-1.0,  0.0);
  if(d < 5.5) return vec2(-1.0, -1.0);
  if(d < 6.5) return vec2( 0.0, -1.0);
  return vec2( 1.0, -1.0);
}

float packDirs(float a, float b){ return a + 9.0*b; }
float unpackDir(float p, int c){
  if(c == 0) return mod(p, 9.0);
  return floor(p/9.0);
}

float polarMask(float sy){ return 1.0 - smoothstep(0.600, 0.780, abs(sy)); }

vec4 clockCell(){ return fetchTexel(net, vec2(0.0, simRows(simEdge()))); }

vec4 vigourCell(){ return fetchTexel(net, vec2(1.0, simRows(simEdge()))); }

vec3 colonyVigour(){
  vec4 v = vigourCell();
  if(v.a < 0.5) return vec3(1.0);
  return clamp(v.rgb, 0.35, 2.6);
}

float simClock(){
  vec4 c = clockCell();
  return (c.g + c.r)*CLOCK_UNIT;
}

float simDelta(){
  vec4 c = clockCell();
  float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
  return dtr*(1.0 + WARM_BOOST*c.a)*tempo;
}

float terrain3(vec3 d, float t){
  float s = 0.0, amp = 0.5, f = 2.00;
  vec3 dr = vec3(0.013, -0.009, 0.006);
  for(int o = 0; o < 3; o++){
    s += amp*vnoise3(d*f + dr*t*f);
    f *= 2.05;
    amp *= 0.44;
    dr = vec3(-dr.y, dr.z, dr.x)*1.31;
  }
  return s;
}

float swirlPot(vec3 d, float t){
  vec3 dr = vec3(0.0, t*0.008, t*0.005);

  return vnoise3(d*2.40 + dr) + 0.5*vnoise3(d*5.40 + dr*1.7)
       + 0.55*vnoise3(d*30.0 + dr*2.3);
}

vec3 curlTangent(vec3 d, float t){

  float e = 0.006;
  float p0 = swirlPot(d, t);
  vec3 g = vec3(swirlPot(d + vec3(e, 0.0, 0.0), t),
                swirlPot(d + vec3(0.0, e, 0.0), t),
                swirlPot(d + vec3(0.0, 0.0, e), t)) - vec3(p0);
  g -= d*dot(g, d);
  vec3 c = cross(d, g);
  return c/max(length(c), 1e-6);
}

vec2 warpUV(vec2 uv, float N){
  vec3 d = uvToDir(uv);
  vec2 a = vec2(vnoise3(d*3.80 + 3.1) - 0.5, vnoise3(d*3.80 + 19.7) - 0.5);
  vec2 b = vec2(vnoise3(d*8.50 + 7.7) - 0.5, vnoise3(d*8.50 + 31.3) - 0.5);

  vec2 c = vec2(vnoise3(d*17.0 + 5.3) - 0.5, vnoise3(d*17.0 + 41.9) - 0.5);

  return uv + (a*1.40 + b*0.50 + c*0.42)/simGrid(N);
}

vec3 heartDir(int c, float t){
  float fc = float(c);
  float a = TAU*fc/3.0 + t*0.019 + hash11(3.0)*TAU;
  float lat = (hash11(fc*5.1 + 1.0) - 0.5)*1.05;
  vec3 d = vec3(cos(lat)*cos(a), sin(lat), cos(lat)*sin(a));

  d += vec3(vnoise3(vec3(t*0.021, fc*7.3 + 2.0, 3.0)) - 0.5,
            vnoise3(vec3(fc*3.1 + 5.0, t*0.021 + 11.0, 8.0)) - 0.5,
            vnoise3(vec3(fc*1.7 + 9.0, 5.0, t*0.021 + 2.0)) - 0.5)*0.55;
  return normalize(d);
}

vec3 heartCore(int c, float t){ return heartDir(c, t)*BODY_R*0.30; }

vec4 nutrient(int k, float t){
  float fk = float(k);
  float period = 24.0 + 14.0*hash11(fk*3.7 + 2.0);
  float u = (t + hash11(fk*9.1 + 5.0)*period)/period;
  float e = mod(floor(u), 64.0);
  float ph = fract(u);
  float a = hash11(e*13.1 + fk*57.3)*TAU;
  float sy = hash11(e*29.7 + fk*17.9)*1.70 - 0.85;
  float cy = sqrt(max(1.0 - sy*sy, 0.0));
  float amt = smoothstep(0.0, 0.13, ph)*(1.0 - smoothstep(0.52, 0.98, ph));
  return vec4(cy*cos(a), sy, cy*sin(a), amt);
}

float nutCellRate(float N){
  float cellChord = TAU/N;
  float rCells = NUT_R/cellChord;
  return 0.32/(PI*rCells*rCells);
}

vec4 heartShare(vec3 d, float t, vec3 vig){
  vec3 s = pow(max(vig, vec3(0.05)), vec3(0.45));
  vec3 dd;
  dd.x = length(d - heartDir(0, t))/s.x;
  dd.y = length(d - heartDir(1, t))/s.y;
  dd.z = length(d - heartDir(2, t))/s.z;
  float m = min(dd.x, min(dd.y, dd.z));

  vec3 e = exp(-(dd - vec3(m))*22.0);
  return vec4(e/max(e.x + e.y + e.z, 1e-5), m);
}

#define BEAT_RATE 0.62
#define BEAT_LAG  0.44

float beatPhase(float t){ return t*BEAT_RATE; }

float beatEnv(float ph){
  ph = fract(ph);
  float lub = exp(-ph*ph*140.0);
  float d = ph - 0.165;
  float dub = exp(-d*d*330.0)*0.58;

  float w = 1.0 - ph;
  lub += exp(-w*w*140.0);

  float tail = 0.34*exp(-ph*3.0);
  return lub + dub + tail;
}

vec3 colonyColour(int c){
  if(c == 0) return vec3(1.00, 0.30, 0.155);
  if(c == 1) return vec3(0.24, 1.00, 0.60);
  return vec3(0.50, 0.34, 1.00);
}

void mainAcc1(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

  float t = simClock();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float T = terrain3(d, t);

  float A = 0.0;
  for(int k = 0; k < NUTS; k++){
    vec4 nk = nutrient(k, t);
    float ch = length(d - nk.xyz);
    A += nk.w*exp(-ch*ch/(NUT_R*NUT_R));
  }

  vec3 cT = curlTangent(d, t);
  vec3 h0 = heartDir(0, t), h1 = heartDir(1, t), h2 = heartDir(2, t);
  vec3 dSelf = vec3(length(d - h0), length(d - h1), length(d - h2));

  vec3 best = vec3(-1e9);
  float d0 = 8.0, d1 = 8.0, d2 = 8.0;

  for(int k = 0; k < 8; k++){
    float fk = float(k);
    vec2 off = dirVec(fk);
    float ny = tx.y + off.y;
    if(ny >= 0.0 && ny < simRows(N)){
      float nxx = mod(tx.x + off.x, N);
      vec2 nuv = (vec2(nxx, ny) + 0.5)/simGrid(N);
      vec3 ndir = uvToDir(nuv);
      float nT = terrain3(ndir, t);
      vec4 nn = shellFetch(net, vec2(nxx, ny), N);

      vec3 stepV = ndir - d;
      float stepW = max(length(stepV), 1e-5);

      float bias = meander*2.76*(T - nT)/stepW
                 - fibrosis*0.35*nn.a
                 + swirl*3.60*dot(stepV/stepW, cT)
                 + JITTER*(hash12(tx*3.1 + off*17.3) - 0.5);

      vec3 prog = (dSelf - vec3(length(ndir - h0), length(ndir - h1), length(ndir - h2)))/stepW;

      if(prog.x > 0.0){ float s = PROGRESS_W*prog.x + bias + 1.32*nn.r;
                        if(s > best.x){ best.x = s; d0 = fk; } }
      if(prog.y > 0.0){ float s = PROGRESS_W*prog.y + bias + 1.32*nn.g;
                        if(s > best.y){ best.y = s; d1 = fk; } }
      if(prog.z > 0.0){ float s = PROGRESS_W*prog.z + bias + 1.32*nn.b;
                        if(s > best.z){ best.z = s; d2 = fk; } }
    }
  }

  gl_FragColor = vec4(packDirs(d0, d1), d2, T, A);
}

void mainAcc(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

  float t = simClock();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float mask = polarMask(d.y);

  vec3 vig = colonyVigour();
  vec4 hsw = heartShare(d, t, vig);
  vec3 share = hsw.xyz;

  float reach = exp(-hsw.w*1.15);
  vec3 rain = share*(2.6/(N*simRows(N)))*mask*(0.10 + 0.90*reach);

  float nrate = nutCellRate(N);
  for(int k = 0; k < NUTS; k++){
    vec4 nk = nutrient(k, t);
    float ch = length(d - nk.xyz);
    float blob = nk.w*exp(-ch*ch/(NUT_R*NUT_R));
    if(blob > 0.001){
      vec4 nd = shellBilerp(net, dirToUV(nk.xyz)*simGrid(N), N);
      vec3 hs = heartShare(nk.xyz, t, vig).xyz;

      vec3 claim = hs*(0.35 + 2.4*vec3(nd.r, nd.g, nd.b));
      claim /= max(claim.x + claim.y + claim.z, 1e-5);
      rain += claim*blob*nutrition*nrate*mask;
    }
  }

  vec3 R = rain;
  for(int k = 0; k < 8; k++){
    vec2 off = dirVec(float(k));
    float ny = tx.y + off.y;
    if(ny >= 0.0 && ny < simRows(N)){
      float nxx = mod(tx.x + off.x, N);
      vec4 a1 = shellFetch(acc1, vec2(nxx, ny), N);
      vec4 rn = shellFetch(acc, vec2(nxx, ny), N);
      for(int c = 0; c < COLONIES; c++){
        float dd = c == 2 ? a1.g : unpackDir(a1.r, c);
        if(dd < 7.5){
          vec2 flow = dirVec(dd) + off;
          if(dot(flow, flow) < 0.1){
            if(c == 0) R.x += rn.x;
            else if(c == 1) R.y += rn.y;
            else R.z += rn.z;
          }
        }
      }
    }
  }

  gl_FragColor = vec4(min(R, vec3(8.0)), 1.0);
}

void mainNet(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);

  if(tx.x == 0.0 && tx.y == simRows(N)){
    vec4 c = fetchTexel(net, vec2(0.0, simRows(N)));
    if(c.b < 0.5){
      gl_FragColor = vec4(0.0, 0.0, 1.0, 1.0);
      return;
    }
    float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
    float dts = dtr*(1.0 + WARM_BOOST*c.a)*tempo;
    float frac = c.r + dts/CLOCK_UNIT;
    float epoch = c.g + floor(frac);
    gl_FragColor = vec4(fract(frac), mod(epoch, CLOCK_WRAP), 1.0, c.a*exp(-dts*WARM_RATE));
    return;
  }

  if(tx.x == 1.0 && tx.y == simRows(N)){
    vec4 v = fetchTexel(net, vec2(1.0, simRows(N)));
    vec4 c0 = clockCell();
    if(c0.b < 0.5 || v.a < 0.5){
      gl_FragColor = vec4(1.0, 1.0, 1.0, 1.0);
      return;
    }
    float t = simClock();
    float dts = simDelta();
    vec3 vig = clamp(v.rgb, 0.35, 2.6);

    vec3 harvest = vec3(0.0);
    for(int k = 0; k < NUTS; k++){
      vec4 nk = nutrient(k, t);
      if(nk.w > 0.02){
        vec4 nd = shellBilerp(net, dirToUV(nk.xyz)*simGrid(N), N);

        vec4 hs = heartShare(nk.xyz, t, vec3(1.0));
        vec3 claim = hs.xyz*(0.45 + 1.2*vec3(nd.r, nd.g, nd.b));
        claim /= max(claim.x + claim.y + claim.z, 1e-5);
        harvest += claim*nk.w;
      }
    }

    float tot = harvest.x + harvest.y + harvest.z;
    vec3 frac2 = tot > 1e-4 ? harvest/tot : vec3(1.0/3.0);
    vec3 targetV = vec3(1.0) + (frac2 - vec3(1.0/3.0))*1.55*nutrition;
    vig += (targetV - vig)*(1.0 - exp(-dts*0.22));
    gl_FragColor = vec4(clamp(vig, 0.35, 2.6), 1.0);
    return;
  }
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

  vec4 c = clockCell();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float t = simClock();
  float dts = simDelta();
  float mask = polarMask(d.y);

  vec4 me = fetchTexel(net, tx);
  vec3 D = vec3(me.r, me.g, me.b);
  float scarTissue = me.a;

  if(c.b < 0.5){

    gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
    return;
  }

  vec3 R = fetchTexel(acc, tx).rgb;

  vec3 q = vec3(2.40)*pow(max(R, vec3(0.0)) + 1e-6, vec3(0.385));

  float cut = mix(0.55, 0.12, clamp(density, 0.0, 1.0));
  vec3 target = smoothstep(vec3(cut), vec3(cut*1.55 + 0.05), q);
  target *= mask;

  float outletR = 3.2*TAU/N;
  for(int k = 0; k < COLONIES; k++){
    float dh = length(d - heartDir(k, t));

    float blob = exp(-dh*dh/(outletR*outletR))*mask;
    if(k == 0) target.x = max(target.x, blob*0.95);
    else if(k == 1) target.y = max(target.y, blob*0.95);
    else target.z = max(target.z, blob*0.95);
  }

  float up = 2.40;
  vec3 rate3 = vec3(
    target.x > D.x ? up : up*clamp(wither, 0.02, 2.0)/1.20,
    target.y > D.y ? up : up*clamp(wither, 0.02, 2.0)/1.20,
    target.z > D.z ? up : up*clamp(wither, 0.02, 2.0)/1.20);
  D += (target - D)*(vec3(1.0) - exp(-rate3*dts));

  float activity = max(D.x, max(D.y, D.z));
  float open = smoothstep(0.04, 0.42, activity);

  float ct = (1.0 - open)*clamp(fibrosis*0.72, 0.0, 1.15);
  scarTissue += (ct - scarTissue)*(1.0 - exp(-(0.30 + open*1.6)*dts));

  gl_FragColor = vec4(clamp(D, 0.0, 1.0), clamp(scarTissue, 0.0, 1.0));
}

void mainMem(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

  float dts = simDelta();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec2 gq = warpUV(uv, N)*simGrid(N);
  vec4 me = fetchTexel(mem, tx);
  vec4 nt = shellBilerp(net, gq, N);

  vec3 scar = me.rgb;
  scar += (clamp(nt.rgb*1.25, 0.0, 1.0) - scar)*(1.0 - exp(-0.55*dts));
  scar *= exp(-0.0090*dts);

  float dcut = mix(0.55, 0.12, clamp(density, 0.0, 1.0));
  float depth = 0.0;
  for(int j = 0; j < 9; j++){
    vec2 off = vec2(mod(float(j), 3.0) - 1.0, floor(float(j)/3.0) - 1.0);

    vec3 R = shellFetch(acc, floor(gq) + off, N).rgb;
    float q = 2.40*pow(max(max(R.x, max(R.y, R.z)), 0.0) + 1e-6, 0.385);
    float wgt = exp(-dot(off, off)*0.30);
    depth = max(depth, smoothstep(dcut*1.03, dcut*2.45, q)*wgt);
  }

  gl_FragColor = vec4(clamp(scar, 0.0, 1.0), depth);
}

vec3 emissionAt(vec2 uv, float N){
  vec3 d = uvToDir(uv);
  vec4 nt = shellBilerp(net, warpUV(uv, N)*simGrid(N), N);
  vec3 D = nt.rgb;
  vec3 e = colonyColour(0)*D.x + colonyColour(1)*D.y + colonyColour(2)*D.z;

  return e*smoothstep(0.02, 0.30, max(D.x, max(D.y, D.z)))*polarMask(d.y);
}

void mainGlow(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 blurC = vec3(0.0);
  float blurD = 0.0;
  float wsum = 0.0;
  for(int j = 0; j < 9; j++){
    vec2 off = vec2(mod(float(j), 3.0) - 1.0, floor(float(j)/3.0) - 1.0);
    vec4 g = shellFetch(glow, tx + off, N);
    float wgt = exp(-dot(off, off)*0.62);
    blurC += g.rgb*wgt;
    blurD += g.a*wgt;
    wsum += wgt;
  }
  blurC /= wsum;
  blurD /= wsum;

  vec3 e = emissionAt(uv, N);
  float ed = max(e.x, max(e.y, e.z));

  gl_FragColor = vec4(mix(blurC, e, 0.17), mix(blurD, ed, 0.085));
}

vec3 gL0, gL1, gL2;
float gK0, gK1, gK2;
vec3 gB0, gB1;
float gS0, gS1;
float gFuse, gScale;
vec3 gX0, gX1;
float gA0, gA1, gPer, gPerT;
float gRMax;
float gRFit;

float smin(float a, float b, float k){
  float h = clamp(0.5 + 0.5*(b - a)/k, 0.0, 1.0);
  return mix(b, a, h) - k*h*(1.0 - h);
}
float smax(float a, float b, float k){ return -smin(-a, -b, k); }

float lobeR(vec3 d, vec3 c, float k2){
  float b = dot(c, d);
  return b + sqrt(max(k2 + b*b, 0.0));
}

float biteT(vec3 d, vec3 q, float sq){
  float b = dot(q, d);
  float disc = b*b - sq;

  float away = step(disc, 0.0) + step(b, 0.0);
  return b - sqrt(max(disc, 0.0)) + away*8.0;
}

#define GUARD 1.10

void setBite(float a, float T, vec3 axis, out vec3 q, out float sq){
  float sum = T*T/a;
  float qm = 0.5*(sum + a);
  float Rq = 0.5*(sum - a);
  q = normalize(axis)*qm;
  sq = qm*qm - Rq*Rq;
}

void setupBody(float t){
  float w0 = t*0.037, w1 = t*0.029, w2 = t*0.043;

  float sys = beatEnv(beatPhase(t));
  float sq = clamp(pulse, 0.0, 1.5)*sys;

  float sep = 0.18 + 0.62*turgor;

  vec3 aL = normalize(vec3(cos(w0), 0.42, sin(w0)));
  vec3 aM = normalize(vec3(cos(w2*1.3 + 3.9), -0.62, sin(w2*1.3 + 3.9)));

  gL0 = vec3(0.03*sin(w1), 0.02, 0.03*cos(w1))*BODY_R;
  float R0 = (0.90 + 0.03*sin(t*0.21))*BODY_R;
  gK0 = R0*R0 - dot(gL0, gL0);

  gL1 = aL*(sep*0.88*BODY_R)*(1.0 - 0.10*sq);
  float R1 = (0.72 - 0.055*sq + 0.03*sin(t*0.17 + 1.7))*BODY_R;
  gK1 = R1*R1 - dot(gL1, gL1);

  gL2 = aM*(sep*1.02*BODY_R)*(1.0 - 0.05*sq);
  float R2 = (0.58 + 0.025*sin(t*0.26 + 3.1))*BODY_R;
  gK2 = R2*R2 - dot(gL2, gL2);

  gFuse = (0.075 + 0.035*sin(t*0.13))*BODY_R;

  setBite(0.445 + 0.022*sin(t*0.11 + 2.4), GUARD,
          vec3(cos(w1*1.7 + 2.0), 0.55, sin(w0*0.8 + 1.0)), gB0, gS0);
  setBite(0.560, GUARD,
          vec3(sin(w1*1.3 + 0.4), -0.28, cos(w2*0.9 + 2.2)), gB1, gS1);

  gX0 = normalize(vec3(cos(w0*1.4), 0.45, sin(w0*1.4)));
  gX1 = normalize(vec3(0.30, cos(w1*1.6), sin(w1*1.6)));
  gA0 = turgor*0.030;
  gA1 = turgor*0.045;
  gPer = breath*0.024;
  gPerT = -t*0.183;

  gScale = 1.0 + breath*(0.042*sin(t*0.95) + 0.015*sin(t*1.90 + 1.0)) - 0.026*sq;

  float lobeMax = max(length(gL0) + sqrt(max(gK0 + dot(gL0, gL0), 0.0)),
                  max(length(gL1) + sqrt(max(gK1 + dot(gL1, gL1), 0.0)),
                      length(gL2) + sqrt(max(gK2 + dot(gL2, gL2), 0.0))));
  gRMax = (lobeMax + 0.25*gFuse)*(1.0 + gA0 + 0.667*gA1 + gPer)*gScale*1.02;

  float fitMax = (0.966 + 0.30*sep)*BODY_R;
  gRFit = (fitMax + 0.25*gFuse)*(1.0 + gA0 + 0.667*gA1 + gPer)
        *(1.0 + 0.057*breath)*1.02;
}

float bodyRadius(vec3 d, float t){
  float r = lobeR(d, gL0, gK0);
  r = smax(r, lobeR(d, gL1, gK1), gFuse);
  r = smax(r, lobeR(d, gL2, gK2), gFuse);

  float ph = fract(d.y*1.03 + gPerT)*2.0 - 1.0;
  float ripple = ph*(1.0 - abs(ph))*4.0;
  float q0 = dot(d, gX0), q1 = dot(d, gX1);
  r *= 1.0 + gA0*q0 + gA1*(q1*q1 - 0.3333) + gPer*ripple;

  r = smin(r, min(biteT(d, gB0, gS0), biteT(d, gB1, gS1)), 0.160*BODY_R);

  return max(r*gScale, 0.10*BODY_R);
}

float bodyField(vec3 p, float t){
  float r = length(p);
  if(r < 1e-5) return -BODY_R;
  return r - bodyRadius(p/r, t);
}

vec3 bodyNormal(vec3 p, float t){
  vec3 d = normalize(p);
  vec3 ax = abs(d.y) < 0.9 ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0);
  vec3 tu = normalize(cross(ax, d));
  vec3 tv = cross(d, tu);
  float e = 0.004;
  float r0 = bodyRadius(d, t);
  float ru = bodyRadius(normalize(d + tu*e), t);
  float rv = bodyRadius(normalize(d + tv*e), t);
  return normalize(d - (tu*(ru - r0) + tv*(rv - r0))/(e*max(r0, 1e-4)));
}

bool sphereSpan(vec3 ro, vec3 rd, float R, out float t0, out float t1){
  float b = dot(ro, rd);
  float c = dot(ro, ro) - R*R;
  float h = b*b - c;
  if(h < 0.0){ t0 = 0.0; t1 = -1.0; return false; }
  h = sqrt(h);
  t0 = max(-b - h, 0.0);
  t1 = -b + h;
  return t1 > t0;
}

vec3 shellEmission(vec3 d, float t, vec3 sig, float N, float graze, vec3 rd){
  vec2 grid = dirToUV(d)*simGrid(N);

  vec3 Tv = rd - d*dot(rd, d);
  float tl = length(Tv);
  vec2 gstep = vec2(0.0);
  if(tl > 1e-4){
    Tv /= tl;

    float cl = max(sqrt(max(1.0 - d.y*d.y, 0.0)), 0.10);
    vec3 east = normalize(vec3(-d.z, 0.0, d.x) + vec3(1e-6));
    vec3 north = cross(d, east);
    vec2 duv = vec2(dot(Tv, east)/(TAU*cl), dot(Tv, north)*cl*0.5);
    vec2 gdir = duv*simGrid(N);
    float gl2 = length(gdir);

    if(gl2 > 1e-5) gstep = (gdir/gl2)*(graze*7.0*(N/288.0));
  }
  vec4 g = shellBilerp(glow, grid, N)*0.46
         + shellBilerp(glow, grid + gstep, N)*0.27
         + shellBilerp(glow, grid - gstep, N)*0.27;
  float depth = shellBilerp(mem, grid, N).a;

  float gm = max(g.x, max(g.y, g.z));
  vec3 broad = (g.rgb/max(gm, 1e-4))*g.a;

  float softD = clamp(depth, 0.0, 1.0);
  float softG = clamp(graze, 0.0, 1.0);

  float hd = min(length(d - heartDir(0, t)),
             min(length(d - heartDir(1, t)), length(d - heartDir(2, t))));
  float env = beatEnv(beatPhase(t) - hd*BEAT_LAG);
  float pk = clamp(pulse, 0.0, 1.5);

  float dilate = 0.75*env*pk*(1.0 - 0.60*softG);

  vec3 gl = g.rgb*(1.0 - 0.30*softD)*(1.0 - 0.92*softG)
          + broad*(0.26 + 0.85*softD + 0.70*softG + dilate);
  float glLuma = clamp(max(gl.x, max(gl.y, gl.z)), 0.0, 1.0);

  float pulseK = mix(1.0, 0.22 + 1.10*env, pk);

  vec3 e = gl*(0.70 + 2.10*glLuma*glLuma)*pulseK*emissive*2.05;

  e /= (1.0 + 0.55*max(e.x, max(e.y, e.z)));

  return e*exp(-sig*(0.020 + depth*0.290))*polarMask(d.y);
}

void mainScene(){
  float N = simEdge();
  vec2 uv = (gl_FragCoord.xy - 0.5*RENDERSIZE)/min(RENDERSIZE.x, RENDERSIZE.y);
  float t = simClock();
  setupBody(t);

  float ca = TIME*0.055 + 1.30;
  float el = 0.30 + 0.16*sin(TIME*0.037);
  float dist = gRFit/0.2895;
  vec3 ro = vec3(cos(ca), 0.0, sin(ca))*dist + vec3(0.0, el, 0.0);
  vec3 ta = vec3(0.0, 0.0, 0.0);
  vec3 fw = normalize(ta - ro);
  vec3 rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
  vec3 up = cross(rt, fw);
  vec3 rd = normalize(uv.x*rt + uv.y*up + 1.55*fw);

  vec3 sig = mix(vec3(1.0), vec3(0.50, 1.30, 1.75), flesh)*(3.2/translucency);

  vec3 col = vec3(0.0);
  float t0, t1;
  bool span = sphereSpan(ro, rd, gRMax, t0, t1);

  float tFront = -1.0, tBack = -1.0, inside = 0.0;
  float backLo = -1.0, backHi = -1.0;
  if(span){

    float stepT = (t1 - t0)/15.0;
    float prevT = t0;
    float prevS = bodyField(ro + rd*t0, t);
    for(int i = 1; i <= 15; i++){
      float tc = t0 + stepT*float(i);
      float s = bodyField(ro + rd*tc, t);
      if(s < 0.0) inside += stepT;
      if(prevS >= 0.0 && s < 0.0 && tFront < 0.0){
        float lo = prevT, hi = tc;
        for(int j = 0; j < 5; j++){
          float mid = 0.5*(lo + hi);
          if(bodyField(ro + rd*mid, t) < 0.0) hi = mid; else lo = mid;
        }
        tFront = hi;
      }
      if(prevS < 0.0 && s >= 0.0){ backLo = prevT; backHi = tc; }
      prevS = s;
      prevT = tc;
    }
    if(backHi > 0.0){
      float lo = backLo, hi = backHi;
      for(int j = 0; j < 5; j++){
        float mid = 0.5*(lo + hi);
        if(bodyField(ro + rd*mid, t) > 0.0) hi = mid; else lo = mid;
      }
      tBack = hi;
    } else if(tFront > 0.0){
      tBack = t1;
    }
  }

  float chord = max(inside, 0.0);

  vec3 through = exp(-sig*chord);

  if(tFront > 0.0){
    vec3 pF = ro + rd*tFront;
    vec3 dF = normalize(pF);
    vec3 nGeo = bodyNormal(pF, t);
    vec3 n = nGeo;
    vec2 guv = dirToUV(dF);

    float radial = clamp(dot(nGeo, dF), 0.0, 1.0);
    float steep = 1.0 - radial*radial;

    vec2 du = vec2(1.5)/simGrid(N);
    float gC = shellFetch(glow, guv*simGrid(N), N).a;
    float gU = shellFetch(glow, (guv + vec2(du.x, 0.0))*simGrid(N), N).a;
    float gV = shellFetch(glow, (guv + vec2(0.0, du.y))*simGrid(N), N).a;
    vec4 memF = shellBilerp(mem, guv*simGrid(N), N);
    vec3 tanU = normalize(cross(vec3(0.0, 1.0, 0.0), dF) + vec3(1e-5));
    vec3 tanV = cross(dF, tanU);

    float relief = 38.0*(1.0 + 0.85*clamp(pulse, 0.0, 1.5)*beatEnv(beatPhase(t)))
                 *(1.0 - 0.92*steep);
    n = normalize(n - (tanU*(gU - gC) + tanV*(gV - gC))*relief*(1.0 - 0.65*memF.a));

    float fres = pow(1.0 - clamp(dot(n, -rd), 0.0, 1.0), 3.4);

    vec3 key = normalize(vec3(0.78, 0.40, -0.48));

    float ndl = max(dot(n, key), 0.0);
    float wrap = clamp((dot(n, key) + 0.62)/1.62, 0.0, 1.0);
    vec4 ntF = shellBilerp(net, guv*simGrid(N), N);
    float micro = vnoise3(dF*24.0)*0.6 + vnoise3(dF*53.0)*0.4;

    float opaque = clamp(ntF.a, 0.0, 1.0);
    vec3 skin = mix(vec3(0.030, 0.013, 0.012), vec3(0.190, 0.140, 0.130), opaque);
    skin *= 0.90 + 0.20*micro;

    col += skin*(0.021 + 0.34*ndl + 0.26*wrap*wrap);

    col += vec3(0.62, 0.46, 0.40)*fres*0.055;

    col += (colonyColour(0)*memF.r + colonyColour(1)*memF.g + colonyColour(2)*memF.b)*0.030;

    vec3 windowT = exp(-sig*opaque*0.32);

    float facing = clamp(dot(nGeo, -rd), 0.0, 1.0);
    float graze = 1.0 - facing*facing;
    float foot = clamp(max(graze*graze, steep), 0.0, 1.0);

    float limbFade = smoothstep(0.05, 0.42, facing)*(1.0 - 0.80*steep);

    vec3 rIn = refract(rd, nGeo, 1.0/SKIN_IOR);
    if(dot(rIn, rIn) < 1e-6) rIn = rd;

    vec4 gW = shellBilerp(glow, guv*simGrid(N), N);
    vec3 dV = normalize(pF + rIn*(SKIN_H*(1.5 + 1.7*clamp(gW.a, 0.0, 1.0))));

    float skinPath = min(1.0/max(facing, 0.09), 7.0)*SKIN_H;
    vec3 skinAbs = exp(-sig*skinPath);

    col += shellEmission(dV, t, sig, N, foot, rd)*windowT*limbFade*skinAbs;

    float gWm = max(gW.r, max(gW.g, gW.b));
    vec3 washHue = gW.rgb/max(gWm, 1e-4);
    vec3 inScat = (vec3(1.0) - exp(-sig*0.60*skinPath))*exp(-sig*0.45*skinPath);
    col += washHue*gW.a*inScat*1.55*windowT*polarMask(dF.y);

    vec3 dB = normalize(ro + rd*max(tBack, tFront));
    if(tBack > tFront){
      col += shellEmission(dB, t, sig, N, max(foot, 0.55), rd)*through*windowT*limbFade;
    }

    vec4 gF = shellBilerp(glow, dirToUV(dF)*simGrid(N), N);
    vec4 gB = shellBilerp(glow, dirToUV(dB)*simGrid(N), N);
    float mF = max(gF.x, max(gF.y, gF.z));
    float mB = max(gB.x, max(gB.y, gB.z));
    vec3 medium = (gF.rgb/max(mF, 1e-4))*gF.a + (gB.rgb/max(mB, 1e-4))*gB.a;

    float medBeat = mix(1.0, 0.36 + 0.95*beatEnv(beatPhase(t) - 0.33),
                        clamp(pulse, 0.0, 1.5));
    col += medium*((vec3(1.0) - through)/max(sig, vec3(1e-4)))*1.9*windowT*medBeat;

    vec3 vg = colonyVigour();

    float coreBeat = mix(1.0, 0.30 + 1.55*beatEnv(beatPhase(t)),
                         clamp(pulse, 0.0, 1.5));
    for(int k = 0; k < COLONIES; k++){
      vec3 hp = heartCore(k, t);
      float tc = clamp(dot(hp - ro, rd), tFront, max(tBack, tFront));
      float dist = length(ro + rd*tc - hp);
      float core = 1.0/(1.0 + dist*dist/0.0030);
      float vk = k == 0 ? vg.x : (k == 1 ? vg.y : vg.z);
      vec3 cc = colonyColour(k)*core*0.45*vk*emissive*coreBeat*exp(-sig*max(tc - tFront, 0.0));

      col += cc/(1.0 + 0.60*max(cc.x, max(cc.y, cc.z)));
    }
  }

  float lum = max(col.x, max(col.y, col.z));
  gl_FragColor = vec4(col, smoothstep(0.55, 1.6, lum));
}

vec3 tonemap(vec3 x){

  return clamp((x*(2.51*x + 0.03))/(x*(2.43*x + 0.59) + 0.14), 0.0, 1.0);
}

void mainComposite(){
  vec2 uv = gl_FragCoord.xy/RENDERSIZE;
  vec4 s = texture2D(scene, uv);
  vec3 col = s.rgb;

  float ba = hash12(gl_FragCoord.xy*0.731)*0.18;
  vec2 rot = vec2(cos(ba), sin(ba));
  vec3 b = vec3(0.0);
  float wsum = 0.0;
  for(int i = 0; i < 40; i++){
    float fi = float(i) + 0.5;
    float rr = sqrt(fi/40.0);
    float a = fi*2.39996 + ba;
    vec2 c0 = vec2(cos(a), sin(a));

    vec2 off = vec2(c0.x*rot.x - c0.y*rot.y, c0.x*rot.y + c0.y*rot.x)*(rr*0.048);
    off.x *= RENDERSIZE.y/RENDERSIZE.x;
    vec4 sp = texture2D(scene, clamp(uv + off, vec2(0.001), vec2(0.999)));

    float wgt = exp(-rr*rr*1.9);
    b += sp.rgb*sp.a*wgt;
    wsum += wgt;
  }
  col += (b/wsum)*bloom*1.5;

  col = tonemap(col*1.04);

  float vig = 1.0 - 0.42*dot(uv - 0.5, uv - 0.5)*2.2;
  col *= clamp(vig, 0.0, 1.0);

  col = pow(col, vec3(1.0/2.2));

  float grain = hash12(gl_FragCoord.xy*1.7 + 0.5) - 0.5;
  col += (grain/255.0)*smoothstep(0.0, 0.006, max(col.r, max(col.g, col.b)));

  gl_FragColor = vec4(col, 1.0);
}

void main(){
  if(PASSINDEX == 0) mainAcc1();
  else if(PASSINDEX == 1) mainAcc();
  else if(PASSINDEX == 2) mainNet();
  else if(PASSINDEX == 3) mainMem();
  else if(PASSINDEX == 4) mainGlow();
  else if(PASSINDEX == 5) mainScene();
  else mainComposite();
}

/*{
  "DESCRIPTION": "Organ: a beating, three-lobed body of translucent flesh, lit from the inside by the vasculature growing through it. The shape is constructive solid geometry on spheres - three lobes fused with a smooth union and two spherical bites scooped out of the result - and turgor separates the lobes, so the silhouette runs from one round blob to a distinctly three-chambered mass. Three colonies drain the shell through a flow-accumulation network, and the light they make has to travel back out through skin: the vessels are read where the refracted ray reaches them rather than at the surface, so the network sits under the skin with real parallax, and the path through that layer goes as one over the cosine, so the limb is a warm scattered wash and a deep trunk arrives as dull red. A systolic front leaves each nucleus once a beat and runs out along the tree, dilating the vessels it passes and raising them out of the skin, while a chamber contracts on the same beat. Territory is a partition of the surface between three drifting nuclei, each carrying a vigour that tracks its share of the nutrients currently on the body. Forked from Vein Atlas.",
  "CATEGORIES": ["generative", "simulation", "organic"],
  "INPUTS": [
    { "NAME": "tempo", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.15, "MAX": 3.0 },
    { "NAME": "density", "TYPE": "float", "DEFAULT": 0.72, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "wither", "TYPE": "float", "DEFAULT": 1.20, "MIN": 0.05, "MAX": 2.0 },
    { "NAME": "nutrition", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "meander", "TYPE": "float", "DEFAULT": 0.62, "MIN": 0.0, "MAX": 2.4 },
    { "NAME": "swirl", "TYPE": "float", "DEFAULT": 1.25, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "fibrosis", "TYPE": "float", "DEFAULT": 0.45, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "turgor", "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0, "MAX": 0.9 },
    { "NAME": "breath", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "emissive", "TYPE": "float", "DEFAULT": 1.00, "MIN": 0.0, "MAX": 2.5 },
    { "NAME": "pulse", "TYPE": "float", "DEFAULT": 0.75, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "translucency", "TYPE": "float", "DEFAULT": 0.85, "MIN": 0.10, "MAX": 2.2 },
    { "NAME": "flesh", "TYPE": "float", "DEFAULT": 0.70, "MIN": 0.0, "MAX": 1.4 },
    { "NAME": "bloom", "TYPE": "float", "DEFAULT": 0.95, "MIN": 0.0, "MAX": 2.0 },
    { "NAME": "seed", "TYPE": "float", "DEFAULT": 7.0, "MIN": 0.0, "MAX": 30.0 }
  ],
  "PASSES": [
    { "TARGET": "acc1", "FLOAT": true },
    { "TARGET": "acc", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "net", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "mem", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "glow", "PERSISTENT": true, "FLOAT": true },
    { "TARGET": "scene", "FLOAT": true },
    {}
  ]
}*/

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI  3.14159265359
#define TAU 6.28318530718

#define COLONIES 3
#define NUTS 9

#define CLOCK_UNIT 8.0
#define CLOCK_WRAP 1024.0

#define BODY_R 0.60

#define RBOUND 1.17

#define PROGRESS_W 0.28

#define JITTER 0.55

#define SKIN_H 0.052
#define SKIN_IOR 1.42

#define NUT_R 0.20
#define WARM_BOOST 16.0
#define WARM_RATE 0.25

float hash11(float p){ return fract(sin(p*127.1 + seed*19.73)*43758.5453); }
float hash12(vec2 p){ return fract(sin(dot(p, vec2(127.1, 311.7)) + seed*3.17)*43758.5453); }
float hash13(vec3 p){ return fract(sin(dot(p, vec3(127.1, 311.7, 74.7)) + seed*3.17)*43758.5453); }

float vnoise3(vec3 p){
  vec3 i = floor(p), f = fract(p);
  f = f*f*(3.0 - 2.0*f);
  float a = hash13(i + vec3(0.0, 0.0, 0.0));
  float b = hash13(i + vec3(1.0, 0.0, 0.0));
  float c = hash13(i + vec3(0.0, 1.0, 0.0));
  float d = hash13(i + vec3(1.0, 1.0, 0.0));
  float e = hash13(i + vec3(0.0, 0.0, 1.0));
  float g = hash13(i + vec3(1.0, 0.0, 1.0));
  float h = hash13(i + vec3(0.0, 1.0, 1.0));
  float k = hash13(i + vec3(1.0, 1.0, 1.0));
  return mix(mix(mix(a, b, f.x), mix(c, d, f.x), f.y),
             mix(mix(e, g, f.x), mix(h, k, f.x), f.y), f.z);
}

#define ASPECT 3.0
float simEdge(){
  float w = min(576.0, floor(RENDERSIZE.x) - 2.0);
  float h = min(192.0, floor(RENDERSIZE.y) - 2.0);
  return floor(min(w, h*ASPECT));
}
float simRows(float N){ return floor(N/ASPECT); }
vec2 simGrid(float N){ return vec2(N, simRows(N)); }

vec4 fetchTexel(sampler2D tex, vec2 texel){
  return texture2D(tex, (texel + 0.5)/RENDERSIZE);
}

vec3 uvToDir(vec2 uv){
  float phi = (uv.x - 0.5)*TAU;
  float sy = clamp(uv.y*2.0 - 1.0, -1.0, 1.0);
  float cy = sqrt(max(1.0 - sy*sy, 0.0));
  return vec3(cy*cos(phi), sy, cy*sin(phi));
}

vec2 dirToUV(vec3 d){
  return vec2(atan(d.z, d.x)/TAU + 0.5, d.y*0.5 + 0.5);
}

vec4 shellFetch(sampler2D tex, vec2 texel, float N){
  return fetchTexel(tex, vec2(mod(texel.x, N), clamp(texel.y, 0.0, simRows(N) - 1.0)));
}

vec4 shellBilerp(sampler2D tex, vec2 grid, float N){
  vec2 g = grid - 0.5;
  vec2 i = floor(g);
  vec2 f = g - i;
  vec4 a = shellFetch(tex, i, N);
  vec4 b = shellFetch(tex, i + vec2(1.0, 0.0), N);
  vec4 c = shellFetch(tex, i + vec2(0.0, 1.0), N);
  vec4 d = shellFetch(tex, i + vec2(1.0, 1.0), N);
  return mix(mix(a, b, f.x), mix(c, d, f.x), f.y);
}

vec2 dirVec(float d){
  if(d < 0.5) return vec2( 1.0,  0.0);
  if(d < 1.5) return vec2( 1.0,  1.0);
  if(d < 2.5) return vec2( 0.0,  1.0);
  if(d < 3.5) return vec2(-1.0,  1.0);
  if(d < 4.5) return vec2(-1.0,  0.0);
  if(d < 5.5) return vec2(-1.0, -1.0);
  if(d < 6.5) return vec2( 0.0, -1.0);
  return vec2( 1.0, -1.0);
}

float packDirs(float a, float b){ return a + 9.0*b; }
float unpackDir(float p, int c){
  if(c == 0) return mod(p, 9.0);
  return floor(p/9.0);
}

float polarMask(float sy){ return 1.0 - smoothstep(0.600, 0.780, abs(sy)); }

vec4 clockCell(){ return fetchTexel(net, vec2(0.0, simRows(simEdge()))); }

vec4 vigourCell(){ return fetchTexel(net, vec2(1.0, simRows(simEdge()))); }

vec3 colonyVigour(){
  vec4 v = vigourCell();
  if(v.a < 0.5) return vec3(1.0);
  return clamp(v.rgb, 0.35, 2.6);
}

float simClock(){
  vec4 c = clockCell();
  return (c.g + c.r)*CLOCK_UNIT;
}

float simDelta(){
  vec4 c = clockCell();
  float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
  return dtr*(1.0 + WARM_BOOST*c.a)*tempo;
}

float terrain3(vec3 d, float t){
  float s = 0.0, amp = 0.5, f = 2.00;
  vec3 dr = vec3(0.013, -0.009, 0.006);
  for(int o = 0; o < 3; o++){
    s += amp*vnoise3(d*f + dr*t*f);
    f *= 2.05;
    amp *= 0.44;
    dr = vec3(-dr.y, dr.z, dr.x)*1.31;
  }
  return s;
}

float swirlPot(vec3 d, float t){
  vec3 dr = vec3(0.0, t*0.008, t*0.005);

return vnoise3(d*2.40 + dr) + 0.5*vnoise3(d*5.40 + dr*1.7)
       + 0.55*vnoise3(d*30.0 + dr*2.3);
}

vec3 curlTangent(vec3 d, float t){

float e = 0.006;
  float p0 = swirlPot(d, t);
  vec3 g = vec3(swirlPot(d + vec3(e, 0.0, 0.0), t),
                swirlPot(d + vec3(0.0, e, 0.0), t),
                swirlPot(d + vec3(0.0, 0.0, e), t)) - vec3(p0);
  g -= d*dot(g, d);
  vec3 c = cross(d, g);
  return c/max(length(c), 1e-6);
}

vec2 warpUV(vec2 uv, float N){
  vec3 d = uvToDir(uv);
  vec2 a = vec2(vnoise3(d*3.80 + 3.1) - 0.5, vnoise3(d*3.80 + 19.7) - 0.5);
  vec2 b = vec2(vnoise3(d*8.50 + 7.7) - 0.5, vnoise3(d*8.50 + 31.3) - 0.5);

vec2 c = vec2(vnoise3(d*17.0 + 5.3) - 0.5, vnoise3(d*17.0 + 41.9) - 0.5);

return uv + (a*1.40 + b*0.50 + c*0.42)/simGrid(N);
}

vec3 heartDir(int c, float t){
  float fc = float(c);
  float a = TAU*fc/3.0 + t*0.019 + hash11(3.0)*TAU;
  float lat = (hash11(fc*5.1 + 1.0) - 0.5)*1.05;
  vec3 d = vec3(cos(lat)*cos(a), sin(lat), cos(lat)*sin(a));

d += vec3(vnoise3(vec3(t*0.021, fc*7.3 + 2.0, 3.0)) - 0.5,
            vnoise3(vec3(fc*3.1 + 5.0, t*0.021 + 11.0, 8.0)) - 0.5,
            vnoise3(vec3(fc*1.7 + 9.0, 5.0, t*0.021 + 2.0)) - 0.5)*0.55;
  return normalize(d);
}

vec3 heartCore(int c, float t){ return heartDir(c, t)*BODY_R*0.30; }

vec4 nutrient(int k, float t){
  float fk = float(k);
  float period = 24.0 + 14.0*hash11(fk*3.7 + 2.0);
  float u = (t + hash11(fk*9.1 + 5.0)*period)/period;
  float e = mod(floor(u), 64.0);
  float ph = fract(u);
  float a = hash11(e*13.1 + fk*57.3)*TAU;
  float sy = hash11(e*29.7 + fk*17.9)*1.70 - 0.85;
  float cy = sqrt(max(1.0 - sy*sy, 0.0));
  float amt = smoothstep(0.0, 0.13, ph)*(1.0 - smoothstep(0.52, 0.98, ph));
  return vec4(cy*cos(a), sy, cy*sin(a), amt);
}

float nutCellRate(float N){
  float cellChord = TAU/N;
  float rCells = NUT_R/cellChord;
  return 0.32/(PI*rCells*rCells);
}

vec4 heartShare(vec3 d, float t, vec3 vig){
  vec3 s = pow(max(vig, vec3(0.05)), vec3(0.45));
  vec3 dd;
  dd.x = length(d - heartDir(0, t))/s.x;
  dd.y = length(d - heartDir(1, t))/s.y;
  dd.z = length(d - heartDir(2, t))/s.z;
  float m = min(dd.x, min(dd.y, dd.z));

vec3 e = exp(-(dd - vec3(m))*22.0);
  return vec4(e/max(e.x + e.y + e.z, 1e-5), m);
}

#define BEAT_RATE 0.62
#define BEAT_LAG  0.44

float beatPhase(float t){ return t*BEAT_RATE; }

float beatEnv(float ph){
  ph = fract(ph);
  float lub = exp(-ph*ph*140.0);
  float d = ph - 0.165;
  float dub = exp(-d*d*330.0)*0.58;

float w = 1.0 - ph;
  lub += exp(-w*w*140.0);

float tail = 0.34*exp(-ph*3.0);
  return lub + dub + tail;
}

vec3 colonyColour(int c){
  if(c == 0) return vec3(1.00, 0.30, 0.155);
  if(c == 1) return vec3(0.24, 1.00, 0.60);
  return vec3(0.50, 0.34, 1.00);
}

void mainAcc1(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

float t = simClock();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float T = terrain3(d, t);

float A = 0.0;
  for(int k = 0; k < NUTS; k++){
    vec4 nk = nutrient(k, t);
    float ch = length(d - nk.xyz);
    A += nk.w*exp(-ch*ch/(NUT_R*NUT_R));
  }

vec3 cT = curlTangent(d, t);
  vec3 h0 = heartDir(0, t), h1 = heartDir(1, t), h2 = heartDir(2, t);
  vec3 dSelf = vec3(length(d - h0), length(d - h1), length(d - h2));

vec3 best = vec3(-1e9);
  float d0 = 8.0, d1 = 8.0, d2 = 8.0;

for(int k = 0; k < 8; k++){
    float fk = float(k);
    vec2 off = dirVec(fk);
    float ny = tx.y + off.y;
    if(ny >= 0.0 && ny < simRows(N)){
      float nxx = mod(tx.x + off.x, N);
      vec2 nuv = (vec2(nxx, ny) + 0.5)/simGrid(N);
      vec3 ndir = uvToDir(nuv);
      float nT = terrain3(ndir, t);
      vec4 nn = shellFetch(net, vec2(nxx, ny), N);

vec3 stepV = ndir - d;
      float stepW = max(length(stepV), 1e-5);

float bias = meander*2.76*(T - nT)/stepW
                 - fibrosis*0.35*nn.a
                 + swirl*3.60*dot(stepV/stepW, cT)
                 + JITTER*(hash12(tx*3.1 + off*17.3) - 0.5);

vec3 prog = (dSelf - vec3(length(ndir - h0), length(ndir - h1), length(ndir - h2)))/stepW;

if(prog.x > 0.0){ float s = PROGRESS_W*prog.x + bias + 1.32*nn.r;
                        if(s > best.x){ best.x = s; d0 = fk; } }
      if(prog.y > 0.0){ float s = PROGRESS_W*prog.y + bias + 1.32*nn.g;
                        if(s > best.y){ best.y = s; d1 = fk; } }
      if(prog.z > 0.0){ float s = PROGRESS_W*prog.z + bias + 1.32*nn.b;
                        if(s > best.z){ best.z = s; d2 = fk; } }
    }
  }

gl_FragColor = vec4(packDirs(d0, d1), d2, T, A);
}

void mainAcc(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

float t = simClock();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float mask = polarMask(d.y);

vec3 vig = colonyVigour();
  vec4 hsw = heartShare(d, t, vig);
  vec3 share = hsw.xyz;

float reach = exp(-hsw.w*1.15);
  vec3 rain = share*(2.6/(N*simRows(N)))*mask*(0.10 + 0.90*reach);

float nrate = nutCellRate(N);
  for(int k = 0; k < NUTS; k++){
    vec4 nk = nutrient(k, t);
    float ch = length(d - nk.xyz);
    float blob = nk.w*exp(-ch*ch/(NUT_R*NUT_R));
    if(blob > 0.001){
      vec4 nd = shellBilerp(net, dirToUV(nk.xyz)*simGrid(N), N);
      vec3 hs = heartShare(nk.xyz, t, vig).xyz;

vec3 claim = hs*(0.35 + 2.4*vec3(nd.r, nd.g, nd.b));
      claim /= max(claim.x + claim.y + claim.z, 1e-5);
      rain += claim*blob*nutrition*nrate*mask;
    }
  }

vec3 R = rain;
  for(int k = 0; k < 8; k++){
    vec2 off = dirVec(float(k));
    float ny = tx.y + off.y;
    if(ny >= 0.0 && ny < simRows(N)){
      float nxx = mod(tx.x + off.x, N);
      vec4 a1 = shellFetch(acc1, vec2(nxx, ny), N);
      vec4 rn = shellFetch(acc, vec2(nxx, ny), N);
      for(int c = 0; c < COLONIES; c++){
        float dd = c == 2 ? a1.g : unpackDir(a1.r, c);
        if(dd < 7.5){
          vec2 flow = dirVec(dd) + off;
          if(dot(flow, flow) < 0.1){
            if(c == 0) R.x += rn.x;
            else if(c == 1) R.y += rn.y;
            else R.z += rn.z;
          }
        }
      }
    }
  }

gl_FragColor = vec4(min(R, vec3(8.0)), 1.0);
}

void mainNet(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);

if(tx.x == 0.0 && tx.y == simRows(N)){
    vec4 c = fetchTexel(net, vec2(0.0, simRows(N)));
    if(c.b < 0.5){
      gl_FragColor = vec4(0.0, 0.0, 1.0, 1.0);
      return;
    }
    float dtr = clamp(TIMEDELTA, 0.0, 1.0/24.0);
    float dts = dtr*(1.0 + WARM_BOOST*c.a)*tempo;
    float frac = c.r + dts/CLOCK_UNIT;
    float epoch = c.g + floor(frac);
    gl_FragColor = vec4(fract(frac), mod(epoch, CLOCK_WRAP), 1.0, c.a*exp(-dts*WARM_RATE));
    return;
  }

if(tx.x == 1.0 && tx.y == simRows(N)){
    vec4 v = fetchTexel(net, vec2(1.0, simRows(N)));
    vec4 c0 = clockCell();
    if(c0.b < 0.5 || v.a < 0.5){
      gl_FragColor = vec4(1.0, 1.0, 1.0, 1.0);
      return;
    }
    float t = simClock();
    float dts = simDelta();
    vec3 vig = clamp(v.rgb, 0.35, 2.6);

vec3 harvest = vec3(0.0);
    for(int k = 0; k < NUTS; k++){
      vec4 nk = nutrient(k, t);
      if(nk.w > 0.02){
        vec4 nd = shellBilerp(net, dirToUV(nk.xyz)*simGrid(N), N);

vec4 hs = heartShare(nk.xyz, t, vec3(1.0));
        vec3 claim = hs.xyz*(0.45 + 1.2*vec3(nd.r, nd.g, nd.b));
        claim /= max(claim.x + claim.y + claim.z, 1e-5);
        harvest += claim*nk.w;
      }
    }

float tot = harvest.x + harvest.y + harvest.z;
    vec3 frac2 = tot > 1e-4 ? harvest/tot : vec3(1.0/3.0);
    vec3 targetV = vec3(1.0) + (frac2 - vec3(1.0/3.0))*1.55*nutrition;
    vig += (targetV - vig)*(1.0 - exp(-dts*0.22));
    gl_FragColor = vec4(clamp(vig, 0.35, 2.6), 1.0);
    return;
  }
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

vec4 c = clockCell();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 d = uvToDir(uv);
  float t = simClock();
  float dts = simDelta();
  float mask = polarMask(d.y);

vec4 me = fetchTexel(net, tx);
  vec3 D = vec3(me.r, me.g, me.b);
  float scarTissue = me.a;

if(c.b < 0.5){

gl_FragColor = vec4(0.0, 0.0, 0.0, 1.0);
    return;
  }

vec3 R = fetchTexel(acc, tx).rgb;

vec3 q = vec3(2.40)*pow(max(R, vec3(0.0)) + 1e-6, vec3(0.385));

float cut = mix(0.55, 0.12, clamp(density, 0.0, 1.0));
  vec3 target = smoothstep(vec3(cut), vec3(cut*1.55 + 0.05), q);
  target *= mask;

float outletR = 3.2*TAU/N;
  for(int k = 0; k < COLONIES; k++){
    float dh = length(d - heartDir(k, t));

float blob = exp(-dh*dh/(outletR*outletR))*mask;
    if(k == 0) target.x = max(target.x, blob*0.95);
    else if(k == 1) target.y = max(target.y, blob*0.95);
    else target.z = max(target.z, blob*0.95);
  }

float up = 2.40;
  vec3 rate3 = vec3(
    target.x > D.x ? up : up*clamp(wither, 0.02, 2.0)/1.20,
    target.y > D.y ? up : up*clamp(wither, 0.02, 2.0)/1.20,
    target.z > D.z ? up : up*clamp(wither, 0.02, 2.0)/1.20);
  D += (target - D)*(vec3(1.0) - exp(-rate3*dts));

float activity = max(D.x, max(D.y, D.z));
  float open = smoothstep(0.04, 0.42, activity);

float ct = (1.0 - open)*clamp(fibrosis*0.72, 0.0, 1.15);
  scarTissue += (ct - scarTissue)*(1.0 - exp(-(0.30 + open*1.6)*dts));

gl_FragColor = vec4(clamp(D, 0.0, 1.0), clamp(scarTissue, 0.0, 1.0));
}

void mainMem(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

float dts = simDelta();
  vec2 uv = (tx + 0.5)/simGrid(N);
  vec2 gq = warpUV(uv, N)*simGrid(N);
  vec4 me = fetchTexel(mem, tx);
  vec4 nt = shellBilerp(net, gq, N);

vec3 scar = me.rgb;
  scar += (clamp(nt.rgb*1.25, 0.0, 1.0) - scar)*(1.0 - exp(-0.55*dts));
  scar *= exp(-0.0090*dts);

float dcut = mix(0.55, 0.12, clamp(density, 0.0, 1.0));
  float depth = 0.0;
  for(int j = 0; j < 9; j++){
    vec2 off = vec2(mod(float(j), 3.0) - 1.0, floor(float(j)/3.0) - 1.0);

vec3 R = shellFetch(acc, floor(gq) + off, N).rgb;
    float q = 2.40*pow(max(max(R.x, max(R.y, R.z)), 0.0) + 1e-6, 0.385);
    float wgt = exp(-dot(off, off)*0.30);
    depth = max(depth, smoothstep(dcut*1.03, dcut*2.45, q)*wgt);
  }

gl_FragColor = vec4(clamp(scar, 0.0, 1.0), depth);
}

vec3 emissionAt(vec2 uv, float N){
  vec3 d = uvToDir(uv);
  vec4 nt = shellBilerp(net, warpUV(uv, N)*simGrid(N), N);
  vec3 D = nt.rgb;
  vec3 e = colonyColour(0)*D.x + colonyColour(1)*D.y + colonyColour(2)*D.z;

return e*smoothstep(0.02, 0.30, max(D.x, max(D.y, D.z)))*polarMask(d.y);
}

void mainGlow(){
  float N = simEdge();
  vec2 tx = floor(gl_FragCoord.xy);
  if(tx.x >= N || tx.y >= simRows(N)){ gl_FragColor = vec4(0.0); return; }

vec2 uv = (tx + 0.5)/simGrid(N);
  vec3 blurC = vec3(0.0);
  float blurD = 0.0;
  float wsum = 0.0;
  for(int j = 0; j < 9; j++){
    vec2 off = vec2(mod(float(j), 3.0) - 1.0, floor(float(j)/3.0) - 1.0);
    vec4 g = shellFetch(glow, tx + off, N);
    float wgt = exp(-dot(off, off)*0.62);
    blurC += g.rgb*wgt;
    blurD += g.a*wgt;
    wsum += wgt;
  }
  blurC /= wsum;
  blurD /= wsum;

vec3 e = emissionAt(uv, N);
  float ed = max(e.x, max(e.y, e.z));

gl_FragColor = vec4(mix(blurC, e, 0.17), mix(blurD, ed, 0.085));
}

vec3 gL0, gL1, gL2;
float gK0, gK1, gK2;
vec3 gB0, gB1;
float gS0, gS1;
float gFuse, gScale;
vec3 gX0, gX1;
float gA0, gA1, gPer, gPerT;
float gRMax;
float gRFit;

float smin(float a, float b, float k){
  float h = clamp(0.5 + 0.5*(b - a)/k, 0.0, 1.0);
  return mix(b, a, h) - k*h*(1.0 - h);
}
float smax(float a, float b, float k){ return -smin(-a, -b, k); }

float lobeR(vec3 d, vec3 c, float k2){
  float b = dot(c, d);
  return b + sqrt(max(k2 + b*b, 0.0));
}

float biteT(vec3 d, vec3 q, float sq){
  float b = dot(q, d);
  float disc = b*b - sq;

float away = step(disc, 0.0) + step(b, 0.0);
  return b - sqrt(max(disc, 0.0)) + away*8.0;
}

#define GUARD 1.10

void setBite(float a, float T, vec3 axis, out vec3 q, out float sq){
  float sum = T*T/a;
  float qm = 0.5*(sum + a);
  float Rq = 0.5*(sum - a);
  q = normalize(axis)*qm;
  sq = qm*qm - Rq*Rq;
}

void setupBody(float t){
  float w0 = t*0.037, w1 = t*0.029, w2 = t*0.043;

float sys = beatEnv(beatPhase(t));
  float sq = clamp(pulse, 0.0, 1.5)*sys;

float sep = 0.18 + 0.62*turgor;

vec3 aL = normalize(vec3(cos(w0), 0.42, sin(w0)));
  vec3 aM = normalize(vec3(cos(w2*1.3 + 3.9), -0.62, sin(w2*1.3 + 3.9)));

gL0 = vec3(0.03*sin(w1), 0.02, 0.03*cos(w1))*BODY_R;
  float R0 = (0.90 + 0.03*sin(t*0.21))*BODY_R;
  gK0 = R0*R0 - dot(gL0, gL0);

gL1 = aL*(sep*0.88*BODY_R)*(1.0 - 0.10*sq);
  float R1 = (0.72 - 0.055*sq + 0.03*sin(t*0.17 + 1.7))*BODY_R;
  gK1 = R1*R1 - dot(gL1, gL1);

gL2 = aM*(sep*1.02*BODY_R)*(1.0 - 0.05*sq);
  float R2 = (0.58 + 0.025*sin(t*0.26 + 3.1))*BODY_R;
  gK2 = R2*R2 - dot(gL2, gL2);

gFuse = (0.075 + 0.035*sin(t*0.13))*BODY_R;

setBite(0.445 + 0.022*sin(t*0.11 + 2.4), GUARD,
          vec3(cos(w1*1.7 + 2.0), 0.55, sin(w0*0.8 + 1.0)), gB0, gS0);
  setBite(0.560, GUARD,
          vec3(sin(w1*1.3 + 0.4), -0.28, cos(w2*0.9 + 2.2)), gB1, gS1);

gX0 = normalize(vec3(cos(w0*1.4), 0.45, sin(w0*1.4)));
  gX1 = normalize(vec3(0.30, cos(w1*1.6), sin(w1*1.6)));
  gA0 = turgor*0.030;
  gA1 = turgor*0.045;
  gPer = breath*0.024;
  gPerT = -t*0.183;

gScale = 1.0 + breath*(0.042*sin(t*0.95) + 0.015*sin(t*1.90 + 1.0)) - 0.026*sq;

float lobeMax = max(length(gL0) + sqrt(max(gK0 + dot(gL0, gL0), 0.0)),
                  max(length(gL1) + sqrt(max(gK1 + dot(gL1, gL1), 0.0)),
                      length(gL2) + sqrt(max(gK2 + dot(gL2, gL2), 0.0))));
  gRMax = (lobeMax + 0.25*gFuse)*(1.0 + gA0 + 0.667*gA1 + gPer)*gScale*1.02;

float fitMax = (0.966 + 0.30*sep)*BODY_R;
  gRFit = (fitMax + 0.25*gFuse)*(1.0 + gA0 + 0.667*gA1 + gPer)
        *(1.0 + 0.057*breath)*1.02;
}

float bodyRadius(vec3 d, float t){
  float r = lobeR(d, gL0, gK0);
  r = smax(r, lobeR(d, gL1, gK1), gFuse);
  r = smax(r, lobeR(d, gL2, gK2), gFuse);

float ph = fract(d.y*1.03 + gPerT)*2.0 - 1.0;
  float ripple = ph*(1.0 - abs(ph))*4.0;
  float q0 = dot(d, gX0), q1 = dot(d, gX1);
  r *= 1.0 + gA0*q0 + gA1*(q1*q1 - 0.3333) + gPer*ripple;

r = smin(r, min(biteT(d, gB0, gS0), biteT(d, gB1, gS1)), 0.160*BODY_R);

return max(r*gScale, 0.10*BODY_R);
}

float bodyField(vec3 p, float t){
  float r = length(p);
  if(r < 1e-5) return -BODY_R;
  return r - bodyRadius(p/r, t);
}

vec3 bodyNormal(vec3 p, float t){
  vec3 d = normalize(p);
  vec3 ax = abs(d.y) < 0.9 ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0);
  vec3 tu = normalize(cross(ax, d));
  vec3 tv = cross(d, tu);
  float e = 0.004;
  float r0 = bodyRadius(d, t);
  float ru = bodyRadius(normalize(d + tu*e), t);
  float rv = bodyRadius(normalize(d + tv*e), t);
  return normalize(d - (tu*(ru - r0) + tv*(rv - r0))/(e*max(r0, 1e-4)));
}

bool sphereSpan(vec3 ro, vec3 rd, float R, out float t0, out float t1){
  float b = dot(ro, rd);
  float c = dot(ro, ro) - R*R;
  float h = b*b - c;
  if(h < 0.0){ t0 = 0.0; t1 = -1.0; return false; }
  h = sqrt(h);
  t0 = max(-b - h, 0.0);
  t1 = -b + h;
  return t1 > t0;
}

vec3 shellEmission(vec3 d, float t, vec3 sig, float N, float graze, vec3 rd){
  vec2 grid = dirToUV(d)*simGrid(N);

vec3 Tv = rd - d*dot(rd, d);
  float tl = length(Tv);
  vec2 gstep = vec2(0.0);
  if(tl > 1e-4){
    Tv /= tl;

float cl = max(sqrt(max(1.0 - d.y*d.y, 0.0)), 0.10);
    vec3 east = normalize(vec3(-d.z, 0.0, d.x) + vec3(1e-6));
    vec3 north = cross(d, east);
    vec2 duv = vec2(dot(Tv, east)/(TAU*cl), dot(Tv, north)*cl*0.5);
    vec2 gdir = duv*simGrid(N);
    float gl2 = length(gdir);

if(gl2 > 1e-5) gstep = (gdir/gl2)*(graze*7.0*(N/288.0));
  }
  vec4 g = shellBilerp(glow, grid, N)*0.46
         + shellBilerp(glow, grid + gstep, N)*0.27
         + shellBilerp(glow, grid - gstep, N)*0.27;
  float depth = shellBilerp(mem, grid, N).a;

float gm = max(g.x, max(g.y, g.z));
  vec3 broad = (g.rgb/max(gm, 1e-4))*g.a;

float softD = clamp(depth, 0.0, 1.0);
  float softG = clamp(graze, 0.0, 1.0);

float hd = min(length(d - heartDir(0, t)),
             min(length(d - heartDir(1, t)), length(d - heartDir(2, t))));
  float env = beatEnv(beatPhase(t) - hd*BEAT_LAG);
  float pk = clamp(pulse, 0.0, 1.5);

float dilate = 0.75*env*pk*(1.0 - 0.60*softG);

vec3 gl = g.rgb*(1.0 - 0.30*softD)*(1.0 - 0.92*softG)
          + broad*(0.26 + 0.85*softD + 0.70*softG + dilate);
  float glLuma = clamp(max(gl.x, max(gl.y, gl.z)), 0.0, 1.0);

float pulseK = mix(1.0, 0.22 + 1.10*env, pk);

vec3 e = gl*(0.70 + 2.10*glLuma*glLuma)*pulseK*emissive*2.05;

e /= (1.0 + 0.55*max(e.x, max(e.y, e.z)));

return e*exp(-sig*(0.020 + depth*0.290))*polarMask(d.y);
}

void mainScene(){
  float N = simEdge();
  vec2 uv = (gl_FragCoord.xy - 0.5*RENDERSIZE)/min(RENDERSIZE.x, RENDERSIZE.y);
  float t = simClock();
  setupBody(t);

float ca = TIME*0.055 + 1.30;
  float el = 0.30 + 0.16*sin(TIME*0.037);
  float dist = gRFit/0.2895;
  vec3 ro = vec3(cos(ca), 0.0, sin(ca))*dist + vec3(0.0, el, 0.0);
  vec3 ta = vec3(0.0, 0.0, 0.0);
  vec3 fw = normalize(ta - ro);
  vec3 rt = normalize(cross(fw, vec3(0.0, 1.0, 0.0)));
  vec3 up = cross(rt, fw);
  vec3 rd = normalize(uv.x*rt + uv.y*up + 1.55*fw);

vec3 sig = mix(vec3(1.0), vec3(0.50, 1.30, 1.75), flesh)*(3.2/translucency);

vec3 col = vec3(0.0);
  float t0, t1;
  bool span = sphereSpan(ro, rd, gRMax, t0, t1);

float tFront = -1.0, tBack = -1.0, inside = 0.0;
  float backLo = -1.0, backHi = -1.0;
  if(span){

float stepT = (t1 - t0)/15.0;
    float prevT = t0;
    float prevS = bodyField(ro + rd*t0, t);
    for(int i = 1; i <= 15; i++){
      float tc = t0 + stepT*float(i);
      float s = bodyField(ro + rd*tc, t);
      if(s < 0.0) inside += stepT;
      if(prevS >= 0.0 && s < 0.0 && tFront < 0.0){
        float lo = prevT, hi = tc;
        for(int j = 0; j < 5; j++){
          float mid = 0.5*(lo + hi);
          if(bodyField(ro + rd*mid, t) < 0.0) hi = mid; else lo = mid;
        }
        tFront = hi;
      }
      if(prevS < 0.0 && s >= 0.0){ backLo = prevT; backHi = tc; }
      prevS = s;
      prevT = tc;
    }
    if(backHi > 0.0){
      float lo = backLo, hi = backHi;
      for(int j = 0; j < 5; j++){
        float mid = 0.5*(lo + hi);
        if(bodyField(ro + rd*mid, t) > 0.0) hi = mid; else lo = mid;
      }
      tBack = hi;
    } else if(tFront > 0.0){
      tBack = t1;
    }
  }

float chord = max(inside, 0.0);

vec3 through = exp(-sig*chord);

if(tFront > 0.0){
    vec3 pF = ro + rd*tFront;
    vec3 dF = normalize(pF);
    vec3 nGeo = bodyNormal(pF, t);
    vec3 n = nGeo;
    vec2 guv = dirToUV(dF);

float radial = clamp(dot(nGeo, dF), 0.0, 1.0);
    float steep = 1.0 - radial*radial;

vec2 du = vec2(1.5)/simGrid(N);
    float gC = shellFetch(glow, guv*simGrid(N), N).a;
    float gU = shellFetch(glow, (guv + vec2(du.x, 0.0))*simGrid(N), N).a;
    float gV = shellFetch(glow, (guv + vec2(0.0, du.y))*simGrid(N), N).a;
    vec4 memF = shellBilerp(mem, guv*simGrid(N), N);
    vec3 tanU = normalize(cross(vec3(0.0, 1.0, 0.0), dF) + vec3(1e-5));
    vec3 tanV = cross(dF, tanU);

float relief = 38.0*(1.0 + 0.85*clamp(pulse, 0.0, 1.5)*beatEnv(beatPhase(t)))
                 *(1.0 - 0.92*steep);
    n = normalize(n - (tanU*(gU - gC) + tanV*(gV - gC))*relief*(1.0 - 0.65*memF.a));

float fres = pow(1.0 - clamp(dot(n, -rd), 0.0, 1.0), 3.4);

vec3 key = normalize(vec3(0.78, 0.40, -0.48));

float ndl = max(dot(n, key), 0.0);
    float wrap = clamp((dot(n, key) + 0.62)/1.62, 0.0, 1.0);
    vec4 ntF = shellBilerp(net, guv*simGrid(N), N);
    float micro = vnoise3(dF*24.0)*0.6 + vnoise3(dF*53.0)*0.4;

float opaque = clamp(ntF.a, 0.0, 1.0);
    vec3 skin = mix(vec3(0.030, 0.013, 0.012), vec3(0.190, 0.140, 0.130), opaque);
    skin *= 0.90 + 0.20*micro;

col += skin*(0.021 + 0.34*ndl + 0.26*wrap*wrap);

col += vec3(0.62, 0.46, 0.40)*fres*0.055;

col += (colonyColour(0)*memF.r + colonyColour(1)*memF.g + colonyColour(2)*memF.b)*0.030;

vec3 windowT = exp(-sig*opaque*0.32);

float facing = clamp(dot(nGeo, -rd), 0.0, 1.0);
    float graze = 1.0 - facing*facing;
    float foot = clamp(max(graze*graze, steep), 0.0, 1.0);

float limbFade = smoothstep(0.05, 0.42, facing)*(1.0 - 0.80*steep);

vec3 rIn = refract(rd, nGeo, 1.0/SKIN_IOR);
    if(dot(rIn, rIn) < 1e-6) rIn = rd;

vec4 gW = shellBilerp(glow, guv*simGrid(N), N);
    vec3 dV = normalize(pF + rIn*(SKIN_H*(1.5 + 1.7*clamp(gW.a, 0.0, 1.0))));

float skinPath = min(1.0/max(facing, 0.09), 7.0)*SKIN_H;
    vec3 skinAbs = exp(-sig*skinPath);

col += shellEmission(dV, t, sig, N, foot, rd)*windowT*limbFade*skinAbs;

float gWm = max(gW.r, max(gW.g, gW.b));
    vec3 washHue = gW.rgb/max(gWm, 1e-4);
    vec3 inScat = (vec3(1.0) - exp(-sig*0.60*skinPath))*exp(-sig*0.45*skinPath);
    col += washHue*gW.a*inScat*1.55*windowT*polarMask(dF.y);

vec3 dB = normalize(ro + rd*max(tBack, tFront));
    if(tBack > tFront){
      col += shellEmission(dB, t, sig, N, max(foot, 0.55), rd)*through*windowT*limbFade;
    }

vec4 gF = shellBilerp(glow, dirToUV(dF)*simGrid(N), N);
    vec4 gB = shellBilerp(glow, dirToUV(dB)*simGrid(N), N);
    float mF = max(gF.x, max(gF.y, gF.z));
    float mB = max(gB.x, max(gB.y, gB.z));
    vec3 medium = (gF.rgb/max(mF, 1e-4))*gF.a + (gB.rgb/max(mB, 1e-4))*gB.a;

float medBeat = mix(1.0, 0.36 + 0.95*beatEnv(beatPhase(t) - 0.33),
                        clamp(pulse, 0.0, 1.5));
    col += medium*((vec3(1.0) - through)/max(sig, vec3(1e-4)))*1.9*windowT*medBeat;

vec3 vg = colonyVigour();

float coreBeat = mix(1.0, 0.30 + 1.55*beatEnv(beatPhase(t)),
                         clamp(pulse, 0.0, 1.5));
    for(int k = 0; k < COLONIES; k++){
      vec3 hp = heartCore(k, t);
      float tc = clamp(dot(hp - ro, rd), tFront, max(tBack, tFront));
      float dist = length(ro + rd*tc - hp);
      float core = 1.0/(1.0 + dist*dist/0.0030);
      float vk = k == 0 ? vg.x : (k == 1 ? vg.y : vg.z);
      vec3 cc = colonyColour(k)*core*0.45*vk*emissive*coreBeat*exp(-sig*max(tc - tFront, 0.0));

col += cc/(1.0 + 0.60*max(cc.x, max(cc.y, cc.z)));
    }
  }

float lum = max(col.x, max(col.y, col.z));
  gl_FragColor = vec4(col, smoothstep(0.55, 1.6, lum));
}

vec3 tonemap(vec3 x){

return clamp((x*(2.51*x + 0.03))/(x*(2.43*x + 0.59) + 0.14), 0.0, 1.0);
}

void mainComposite(){
  vec2 uv = gl_FragCoord.xy/RENDERSIZE;
  vec4 s = texture2D(scene, uv);
  vec3 col = s.rgb;

float ba = hash12(gl_FragCoord.xy*0.731)*0.18;
  vec2 rot = vec2(cos(ba), sin(ba));
  vec3 b = vec3(0.0);
  float wsum = 0.0;
  for(int i = 0; i < 40; i++){
    float fi = float(i) + 0.5;
    float rr = sqrt(fi/40.0);
    float a = fi*2.39996 + ba;
    vec2 c0 = vec2(cos(a), sin(a));

vec2 off = vec2(c0.x*rot.x - c0.y*rot.y, c0.x*rot.y + c0.y*rot.x)*(rr*0.048);
    off.x *= RENDERSIZE.y/RENDERSIZE.x;
    vec4 sp = texture2D(scene, clamp(uv + off, vec2(0.001), vec2(0.999)));

float wgt = exp(-rr*rr*1.9);
    b += sp.rgb*sp.a*wgt;
    wsum += wgt;
  }
  col += (b/wsum)*bloom*1.5;

col = tonemap(col*1.04);

float vig = 1.0 - 0.42*dot(uv - 0.5, uv - 0.5)*2.2;
  col *= clamp(vig, 0.0, 1.0);

col = pow(col, vec3(1.0/2.2));

float grain = hash12(gl_FragCoord.xy*1.7 + 0.5) - 0.5;
  col += (grain/255.0)*smoothstep(0.0, 0.006, max(col.r, max(col.g, col.b)));

gl_FragColor = vec4(col, 1.0);
}

void main(){
  if(PASSINDEX == 0) mainAcc1();
  else if(PASSINDEX == 1) mainAcc();
  else if(PASSINDEX == 2) mainNet();
  else if(PASSINDEX == 3) mainMem();
  else if(PASSINDEX == 4) mainGlow();
  else if(PASSINDEX == 5) mainScene();
  else mainComposite();
}

#ay8z9jlq

⑂ fork it { } json

Organ: a beating, three-lobed body of translucent flesh, lit from the inside by the vasculature growing through it. The shape is constructive solid geometry on spheres - three lobes fused with a smooth union and two spherical bites scooped out of the result - and turgor separates the lobes, so the silhouette runs from one round blob to a distinctly three-chambered mass. Three colonies drain the shell through a flow-accumulation network, and the light they make has to travel back out through skin: the vessels are read where the refracted ray reaches them rather than at the surface, so the network sits under the skin with real parallax, and the path through that layer goes as one over the cosine, so the limb is a warm scattered wash and a deep trunk arrives as dull red. A systolic front leaves each nucleus once a beat and runs out along the tree, dilating the vessels it passes and raising them out of the skin, while a chamber contracts on the same beat. Territory is a partition of the surface between three drifting nuclei, each carrying a vigour that tracks its share of the nutrients currently on the body. Forked from Vein Atlas.

#generative #simulation #organic #csg #subsurface #vasculature #3dgenerativesimulationorganic

by ufffdwitnessed_generationmade with glslopCC0-1.0activecreated 2026-08-23 18:58Z116 views0 forks15 paramsmultipass ×7claude-opus-5 · anthropic

generation prompt

A translucent beating organ: three colonies drain the shell of a booleaned body through a flow-accumulation network, and the light they make travels back out through flesh.

lineage

forked from

Vein Atlas #y2o38gx3

export clip · render a video in your browser

Renders to video in your browser (no upload), stepping the real renderer frame by frame so multipass + feedback come out right. motion auto-animates the shader's own knobs — each export rolls a fresh take (🎲 re-roll, 🔒 to pin one you like). GIF is slower with a tighter length cap; MP4/WebM are faster. Some devices record in real time (noted below if so) — keep the tab visible.

format length max 300 sec at this fps start frame motion fps size quality

caps: ≤1920 per side, ≤60fps, ≤18000 frames (GIF ≤300)

embed · drop this shader into any page

glslop backlink play/pause + fullscreen size

<iframe src="https://glslop.com/embed/ay8z9jlq" width="512" height="512" style="border:0;max-width:100%;aspect-ratio:1/1" loading="lazy" allowfullscreen title="Organ — a shader on glslop"></iframe>

preview ↗

Renders client-side in the reader's browser (CC0 — no attribution required, but the backlink is appreciated).