// Pass 1: eight sheets, one atlas tile each; sheet 5 holds a shard frozen in the ice.
// Alpha stores linear depth, so pass 2 absorbs it like any other sheet.
precision highp float;
uniform vec2 u_res; // size of one tile
uniform float u_time;

float hash(vec2 p) { return fract(sin(dot(p, vec2(127.1, 311.7))) * 43758.5453); }
float noise(vec2 p) {
  vec2 i = floor(p), f = fract(p);
  f = f*f*(3.0-2.0*f);
  return mix(mix(hash(i), hash(i+vec2(1.0,0.0)), f.x),
             mix(hash(i+vec2(0.0,1.0)), hash(i+vec2(1.0,1.0)), f.x), f.y);
}
float fbm(vec2 p) {
  float sum = 0.0, amp = 0.5;
  for (int i = 0; i < 4; i++) {
    sum += amp * noise(p);
    p = mat2(1.62, 1.17, -1.17, 1.62) * p;
    amp *= 0.5;
  }
  return sum;
}
float fbm2(vec2 p) {
  return 0.5 * noise(p) + 0.25 * noise(mat2(1.62, 1.17, -1.17, 1.62) * p);
}
const float EXTENT = 3.8; // half-width of a sheet in world units; must match the ice pass

// Signed distance to an elongated octahedron turning once per cycle. Sharp edges,
// cheap to march: this is the one thing in the slab that is a real 3D object.
vec3 shardLocal(vec3 p, float phase) {
  float c1 = cos(phase), s1 = sin(phase);
  p.xz = mat2(c1, -s1, s1, c1) * p.xz;
  float c2 = cos(phase * 2.0 + 1.0), s2 = sin(phase * 2.0 + 1.0);
  p.yz = mat2(c2, -s2, s2, c2) * p.yz;
  p.yz /= vec2(1.0, 1.45);
  return abs(p);
}
const float SHARD = 0.95; // how big the shard is in sheet units
// Cross two octahedra: the intersection is an eight-sided gem, so the shard still
// reads as a cut stone from any angle instead of a flat kite.
vec3 gemAxes(vec3 l, out vec2 r) {
  r = mat2(0.70711, -0.70711, 0.70711, 0.70711) * l.xz;
  return vec3(abs(r.x), l.y, abs(r.y));
}
float sdShard(vec3 p, float phase) {
  vec3 l = shardLocal(p / SHARD, phase);
  vec2 r;
  vec3 m = gemAxes(l, r);
  float a = (l.x + l.y * 0.86 + l.z - 1.0) * 0.45;
  float b = (m.x * 1.34 + m.y * 0.86 + m.z * 1.34 - 1.0) * 0.40;
  return max(a, b) * SHARD;
}

void main() {
  vec2 tile = floor(gl_FragCoord.xy / u_res);
  float layer = tile.x + tile.y * 4.0;
  vec2 uv = fract(gl_FragCoord.xy / u_res);
  vec2 p = (uv - 0.5) * 2.0 * vec2(EXTENT * u_res.x / u_res.y, EXTENT);
  p *= 0.85 + 0.10 * layer; // far sheets carry finer detail, which reads as distance
  float angle = hash(vec2(layer, 3.1)) * 6.2831853;
  p = mat2(cos(angle), sin(angle), -sin(angle), cos(angle)) * p;

  // Periodic time: every field returns exactly after 30 seconds.
  float phase = u_time * 6.2831853 / 30.0;
  vec2 o = vec2(cos(phase), sin(phase));
  vec2 seed = vec2(layer * 17.3, layer * 9.1);

  // Two nested warp stages, the marbled-plasma motion.
  vec2 q = vec2(fbm(p * 0.50 + seed + o * 0.55),
                fbm(p * 0.50 + seed + 23.1 - o * 0.55));
  vec2 r = vec2(fbm(p * 0.85 + 2.0 * q + seed + vec2(o.y, -o.x) * 0.40),
                fbm(p * 0.85 + 2.0 * q + seed + 31.7 - vec2(o.y, -o.x) * 0.40));
  vec2 w = p * 1.10 + 2.6 * r + seed;   // warped coords of the density field

  // The silhouette is shared by all sheets (slid a little per depth), so the gaps
  // stay open right through the stack instead of filling in with overlap.
  float mass = smoothstep(0.520, 0.548, fbm((p + o * 0.30 + vec2(0.10, -0.07) * layer) * 0.62 + r * 1.1 + 13.7));

  // Light the density as a height surface: $cheap 2-octave gradient, strong terminator.
  float h = 0.05;
  vec2 slope = vec2(fbm2(w + vec2(h, 0.0)) - fbm2(w - vec2(h, 0.0)),
                    fbm2(w + vec2(0.0, h)) - fbm2(w - vec2(0.0, h))) * (10.0 / (2.0 * h));
  vec3 nrm = normalize(vec3(-slope, 1.0));
  float lit = clamp(dot(nrm, normalize(vec3(-0.38, 0.52, 0.76))), 0.0, 1.0);

  // Folded planes across the volume, feathered so they never read as a grid.
  vec2 g = (p + r * 2.0) * 1.4;
  vec2 cell = floor(g), f = fract(g);
  float feather = smoothstep(0.0, 0.10, min(min(f.x, 1.0 - f.x), min(f.y, 1.0 - f.y)));
  float facet = mix(1.0, mix(0.70, 1.25, hash(cell + seed + step(1.0, f.x + f.y) * 37.0)), feather);

  float shade = 0.05 + 1.10 * lit * lit;                // lit tops, near-black undersides
  float e = mass * facet * shade * 0.085
          + mass * (1.0 - mass) * 4.0 * 0.040           // thin torn silhouette
          + mass * pow(lit, 14.0) * 0.110;              // glint where a fold faces the light
  vec3 col = e * mix(vec3(0.02, 0.42, 0.60), vec3(0.08, 0.64, 0.88), facet * 0.5);
  // Faint dye in the open gaps so darkness reads as depth, not as a hole.
  col += smoothstep(0.40, 0.85, fbm(p * 0.30 + seed * 0.4 + o * 0.3)) * vec3(0.0003, 0.0011, 0.0015);

  float solid = 0.45; // sheet 1 keeps its plasma, but its alpha becomes coverage
  if (layer > 0.5 && layer < 1.5) {
    vec3 ro = vec3(0.0, 0.0, 6.0);
    vec3 lit = vec3(0.0);
    float hits = 0.0;
    // 3x3 supersample: pass 2 magnifies this sheet where the surface steepens at the
    // cracks, and a one-pixel silhouette would show up there as big hard steps.
    for (int sy = 0; sy < 3; sy++) {
      for (int sx = 0; sx < 3; sx++) {
        vec2 jitter = (vec2(float(sx), float(sy)) / 3.0 + 1.0 / 6.0 - 0.5) / u_res;
        vec3 rd = normalize(vec3((uv + jitter - 0.5) * 2.0 * vec2(u_res.x / u_res.y, 1.0), -1.0));
        float b = dot(ro, rd);
        float disc = b * b - dot(ro, ro) + 2.4 * 2.4;
        if (disc <= 0.0) continue;
        float near = max(-b - sqrt(disc), 0.0);
        float t = near, d = 1.0;
        for (int k = 0; k < 40; k++) {
          d = sdShard(ro + rd * t, phase);
          if (d < 0.003 || t > near + 5.0) break;
          t += d;
        }
        if (d >= 0.003) continue;
        float e2 = 0.004;
        vec3 pos = ro + rd * t;
        vec3 nrm = normalize(vec3(
          sdShard(pos + vec3(e2, 0.0, 0.0), phase) - sdShard(pos - vec3(e2, 0.0, 0.0), phase),
          sdShard(pos + vec3(0.0, e2, 0.0), phase) - sdShard(pos - vec3(0.0, e2, 0.0), phase),
          sdShard(pos + vec3(0.0, 0.0, e2), phase) - sdShard(pos - vec3(0.0, 0.0, e2), phase)));
        // A key light and a cool fill from the plasma side: flat faces take one tone
        // each, which is what makes the shard read as a cut stone rather than a blob.
        float key = max(dot(nrm, normalize(vec3(-0.45, 0.72, 0.52))), 0.0);
        float fill = max(dot(nrm, normalize(vec3(0.60, -0.20, -0.75))), 0.0);
        float rim = pow(1.0 - max(dot(nrm, -rd), 0.0), 4.0);
        vec3 c = mix(vec3(0.004, 0.016, 0.026), vec3(0.060, 0.290, 0.340), key);
        c += vec3(0.010, 0.060, 0.080) * fill;
        c += vec3(0.20, 0.62, 0.72) * rim * 0.55;
        lit += c;
        hits += 1.0;
      }
    }
    float cov = hits / 9.0;
    col = mix(col, lit / max(hits, 1.0), cov);
    solid = mix(0.45, 1.0, cov);
  }
  float alpha = (layer > 0.5 && layer < 1.5) ? solid : (0.42 + layer * 0.26) / 2.6;
  gl_FragColor = vec4(col, alpha);
}