Crystal Bloom

// Stage 1 · One piece
// A mesh shader makes its own geometry. The host asks how many corners and
// triangles each piece wants, then for each corner and each triangle by
// number. Here there is one piece: a flat hexagon, a corner in the middle
// and six round it, joined by six triangles, seen as the camera goes round.
//
// @mesh 1

// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot

constant float TAU = 6.2831853;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, and how hot it is.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float heat;
};

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many corners and triangles the piece asks for: 7 and 6.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    return uint2(7, 6);
}

// Corner 0 is the middle; corners 1 to 6 go round it, standing upright.
Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    float3 local = float3(0.0);
    if (i > 0u) {
        float a = float(i - 1u) * TAU / 6.0;
        local = float3(cos(a), sin(a), 0.0) * 0.6;
    }
    Vertex o;
    o.world = local + float3(0.0, 0.25, 0.0);
    o.position = bloomScreen(o.world, u.time, u.resolution.x / u.resolution.y);
    o.heat = i == 0u ? 1.0 : 0.3;
    return o;
}

// Triangle i joins the middle to two corners next to each other.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    return uint3(0u, 1u + i, 1u + (i + 1u) % 6u);
}

// Once for every pixel of the triangles: the heat, blended across each
// triangle from its corners, through the palette.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    return float4(kiln(v.heat + (p.heat - 0.5) * 0.4), 1.0);
}

// Behind them: the dark of the kiln, with a glow round the middle.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// Stage 2 · A crystal
// One piece becomes a crystal: rings of corners stacked up a prism, a
// point on top, and two triangles for each face of each band between rings.
// A piece may have 64 corners and 64 triangles at most, so the more faces
// go round it, the fewer rings it can have. Each face is lit on its own.
//
// @mesh 1

// @param sides 6 3 8 Faces around each crystal
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot

constant float TAU = 6.2831853;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, and how far up its crystal it is.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float along;
};

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many faces go round a crystal, and how many rings it may have:
// as many as keep it to 64 corners and 64 triangles.
uint bloomSides(constant Params& p) {
    return uint(clamp(floor(p.sides + 0.5), 3.0, 8.0));
}
uint bloomRings(uint s) {
    return min((64u - s) / (2u * s), 63u / s - 1u);
}

// s corners on each ring, one more ring than bands, and the point.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    uint s = bloomSides(p), rings = bloomRings(s);
    return uint2(s * (rings + 1u) + 1u, 2u * s * rings + s);
}

// Corner i is on ring i / s, side i % s; the last corner is the point.
Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = bloomSides(p), rings = bloomRings(s);
    float len = 1.1, radius = 0.15, tip = radius * 1.6;
    float3 local;
    float up;
    if (i < s * (rings + 1u)) {
        uint j = i / s, side = i % s;
        up = float(j) * len / float(rings);
        float a = float(side) * TAU / float(s);
        local = float3(cos(a) * radius, up, sin(a) * radius);
    } else {
        up = len + tip;
        local = float3(0.0, up, 0.0);
    }
    Vertex o;
    o.world = local + float3(0.0, -0.45, 0.0);
    o.position = bloomScreen(o.world, u.time, u.resolution.x / u.resolution.y);
    o.along = up / (len + tip);
    return o;
}

// The triangles: two for each face of each band between rings, then a fan
// of s up to the tip.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = bloomSides(p), bands = bloomRings(s);
    if (i < 2u * s * bands) {
        uint j = i / (2u * s), q = i % (2u * s), a = q / 2u;
        uint b = (a + 1u) % s;
        uint c00 = j * s + a, c01 = j * s + b;
        uint c10 = c00 + s, c11 = c01 + s;
        return q % 2u == 0u ? uint3(c00, c01, c11) : uint3(c00, c11, c10);
    }
    uint a = (i - 2u * s * bands) % s, b = (a + 1u) % s;
    uint top = bands * s;
    return uint3(top + a, top + b, top + s);
}

// Each face is flat: its normal is the cross of how the world point changes
// across the screen, so the faces catch the light one by one, as cut stones
// do. The colour runs from ember at the root to gold up the tip.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;
    float3 body = kiln(0.38 + 0.22 * v.along + shift);
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float3 c = body * (0.3 + 0.7 * diffuse)
             + float3(1.0, 0.95, 0.87) * glint * 0.9;
    return float4(c, 1.0);
}

// Behind them: the dark of the kiln, with a glow round the middle.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// Stage 3 · A cluster
// 96 pieces: piece 0 is the rock, a lumpy ball of four rings and two poles,
// 64 triangles; the other 95 are crystals. Their roots are spread over the
// rock by the golden angle, and each is turned to lean out from it by
// spread, a little off true, with a length of its own.
//
// @mesh 96

// @param size 1 0.4 1.4 How long the largest crystals grow
// @param sides 6 3 8 Faces around each crystal
// @param spread 0.75 0 1 How far the crystals lean out from the rock: at 0 they all stand up
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot

constant float TAU = 6.2831853;
// The rock's size and middle.
constant float3 ROCK = float3(0.55, 0.36, 0.55);
constant float3 ROCK_AT = float3(0.0, -0.1, 0.0);
// Every piece but the rock is a crystal.
constant uint CRYSTALS = 95;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, how far up its crystal it is, and 1 on the rock.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float along;
    float rock;
};

// Cracked Clay's hash: a whole number in, a random 0 to 1 out.
float hashU(uint n) {
    n = (n << 13u) ^ n;
    n = n * (n * n * 15731u + 789221u) + 1376312589u;
    return float(n & 0x7fffffffu) / float(0x7fffffff);
}

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many faces go round a crystal, and how many rings it may have:
// as many as keep it to 64 corners and 64 triangles.
uint bloomSides(constant Params& p) {
    return uint(clamp(floor(p.sides + 0.5), 3.0, 8.0));
}
uint bloomRings(uint s) {
    return min((64u - s) / (2u * s), 63u / s - 1u);
}

// Where crystal k points: over the rock by the golden angle.
float3 bloomDir(uint k) {
    float y = 0.95 - float(k) / float(CRYSTALS - 1) * 1.25;
    float r = sqrt(max(1.0 - y * y, 0.0));
    float a = float(k) * 2.39996;
    return float3(cos(a) * r, y, sin(a) * r);
}

// Everything about crystal k: where its root is, which way it grows and
// how long it is.
struct Crystal {
    float3 root;
    float3 axis;
    float full;
};

Crystal bloomCrystal(uint k, constant Params& p) {
    Crystal c;
    float3 d = bloomDir(k);
    c.root = ROCK_AT + d * ROCK * 0.8;
    float3 out = normalize(d / ROCK);
    float3 wob = float3(hashU(k * 5u + 1u), hashU(k * 5u + 2u),
                        hashU(k * 5u + 3u)) - 0.5;
    c.axis = normalize(mix(float3(0, 1, 0), out, p.spread) + wob * 0.3);

    // A length of its own: a few long ones, many short.
    float h = hashU(k * 7u + 9u);
    c.full = p.size * (0.2 + 0.8 * h * h) * (0.65 + 0.35 * max(d.y, 0.0));
    return c;
}

// How many corners and triangles a piece asks for. The rock is a fixed 34
// and 64; a crystal all the rings it may have.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    if (piece == 0u) return uint2(34, 64);
    uint s = bloomSides(p), rings = bloomRings(s);
    return uint2(s * (rings + 1u) + 1u, 2u * s * rings + s);
}

// A rock corner: four rings of eight round a lumpy ball, then its two poles.
float3 rockCorner(uint i) {
    if (i >= 32u) {
        float y = i == 32u ? -1.0 : 1.0;
        return ROCK_AT + float3(0.0, y, 0.0) * ROCK * (0.9 + 0.2 * hashU(i));
    }
    uint ring = i / 8u, j = i % 8u;
    float lat = mix(-1.0, 0.95, float(ring) / 3.0);
    float a = (float(j) + 0.5 * float(ring)) * TAU / 8.0;
    float3 d = float3(cos(lat) * cos(a), sin(lat), cos(lat) * sin(a));
    return ROCK_AT + d * ROCK * (0.82 + 0.3 * hashU(i + 50u));
}

Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    Vertex o;
    float aspect = u.resolution.x / u.resolution.y;
    if (piece == 0u) {
        o.world = rockCorner(i);
        o.position = bloomScreen(o.world, u.time, aspect);
        o.along = 0.0; o.rock = 1.0;
        return o;
    }

    uint k = piece - 1u;
    Crystal c = bloomCrystal(k, p);
    uint s = bloomSides(p), most = bloomRings(s), rings = most;
    float len = c.full;
    float radius = 0.02 + 0.06 * c.full;
    float tip = radius * 1.6;

    // Ring j sits at j steps up the crystal.
    float3 local;
    float up;
    if (i < s * (rings + 1u)) {
        uint j = i / s, side = i % s;
        up = float(j) * c.full / float(most);
        float a = float(side) * TAU / float(s) + hashU(k + 900u) * TAU;
        local = float3(cos(a) * radius, up, sin(a) * radius);
    } else {
        up = len + tip;
        local = float3(0.0, up, 0.0);
    }

    // Turned so that its y runs along the crystal's axis.
    float3 d = c.axis;
    float3 t = normalize(cross(abs(d.y) < 0.99 ? float3(0, 1, 0)
                                               : float3(1, 0, 0), d));
    float3 b = cross(d, t);
    o.world = c.root + t * local.x + d * local.y + b * local.z;
    o.position = bloomScreen(o.world, u.time, aspect);

    o.along = clamp(up / max(len + tip, 0.001), 0.0, 1.0);
    o.rock = 0.0;
    return o;
}

// The triangles: two for each face of each band between rings, then a fan
// of s up to the tip. The rock's are bands between its rings and fans to
// its poles.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = piece == 0u ? 8u : bloomSides(p);
    uint bands = piece == 0u ? 3u : (size(piece, u, p).y - s) / (2u * s);
    if (i < 2u * s * bands) {
        uint j = i / (2u * s), q = i % (2u * s), a = q / 2u;
        uint b = (a + 1u) % s;
        uint c00 = j * s + a, c01 = j * s + b;
        uint c10 = c00 + s, c11 = c01 + s;
        return q % 2u == 0u ? uint3(c00, c01, c11) : uint3(c00, c11, c10);
    }
    uint a = (i - 2u * s * bands) % s, b = (a + 1u) % s;
    if (piece == 0u) {
        if (i - 2u * s * bands < s) return uint3(a, b, 32u);
        return uint3(24u + a, 24u + b, 33u);
    }
    uint top = bands * s;
    return uint3(top + a, top + b, top + s);
}

// Once for every pixel of the triangles. Each face is flat: its normal is
// the cross of how the world point changes across the screen, so the faces
// catch the light one by one, as cut stones do.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;

    // The rock: dark, lit by the light above and by the crystals' glow,
    // warmest near its top where they grow.
    if (v.rock > 0.5) {
        float warm = smoothstep(-0.2, 0.3, v.world.y);
        float3 c = kiln(0.16 + shift) * (0.35 + 0.9 * diffuse)
                 + kiln(0.4 + shift) * warm * 0.12;
        return float4(c, 1.0);
    }

    // A crystal: ember at the root to gold up the tip, and a glint on faces
    // turned to the light and the eye.
    float3 body = kiln(0.38 + 0.22 * v.along + shift);
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float3 c = body * (0.3 + 0.7 * diffuse)
             + float3(1.0, 0.95, 0.87) * glint * 0.9;
    return float4(c, 1.0);
}

// Behind them: the dark of the kiln, with a glow round the cluster.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// Stage 4 · Growing
// Now each crystal lives. It starts after a delay, the top ones first, grows
// over the first 30% of a life, holds, then melts back into the rock and
// grows again, a new length each life. As it grows it asks for a ring more,
// so its corners and triangles change from frame to frame, and before it
// starts it asks for none. A new ring starts on top of the last and slides up.
//
// @mesh 96

// @param speed 1 0 3 How fast the crystals grow and melt back
// @param size 1 0.4 1.4 How long the largest crystals grow
// @param sides 6 3 8 Faces around each crystal
// @param spread 0.75 0 1 How far the crystals lean out from the rock: at 0 they all stand up
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot

constant float TAU = 6.2831853;
// The rock's size and middle.
constant float3 ROCK = float3(0.55, 0.36, 0.55);
constant float3 ROCK_AT = float3(0.0, -0.1, 0.0);
// Every piece but the rock is a crystal.
constant uint CRYSTALS = 95;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, how far up its crystal it is, and 1 on the rock.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float along;
    float rock;
};

// Cracked Clay's hash: a whole number in, a random 0 to 1 out.
float hashU(uint n) {
    n = (n << 13u) ^ n;
    n = n * (n * n * 15731u + 789221u) + 1376312589u;
    return float(n & 0x7fffffffu) / float(0x7fffffff);
}

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many faces go round a crystal, and how many rings it may have:
// as many as keep it to 64 corners and 64 triangles.
uint bloomSides(constant Params& p) {
    return uint(clamp(floor(p.sides + 0.5), 3.0, 8.0));
}
uint bloomRings(uint s) {
    return min((64u - s) / (2u * s), 63u / s - 1u);
}

// Where crystal k points: over the rock by the golden angle.
float3 bloomDir(uint k) {
    float y = 0.95 - float(k) / float(CRYSTALS - 1) * 1.25;
    float r = sqrt(max(1.0 - y * y, 0.0));
    float a = float(k) * 2.39996;
    return float3(cos(a) * r, y, sin(a) * r);
}

// When crystal k first starts, and how long each of its lives lasts.
float2 bloomLife(uint k) {
    float delay = 0.4 + 7.0 * float(k) / float(CRYSTALS);
    return float2(delay + 1.2 * hashU(k + 300u), 18.0 + 8.0 * hashU(k + 600u));
}

// Everything about crystal k at shader time T: where its root is, which way
// it grows, how long it is when grown, how far grown it is (0 to 1) and
// which life it is in.
struct Crystal {
    float3 root;
    float3 axis;
    float full;
    float grown;
    float life;
};

Crystal bloomCrystal(uint k, float T, constant Params& p) {
    Crystal c;
    float3 d = bloomDir(k);
    c.root = ROCK_AT + d * ROCK * 0.8;
    float3 out = normalize(d / ROCK);
    float3 wob = float3(hashU(k * 5u + 1u), hashU(k * 5u + 2u),
                        hashU(k * 5u + 3u)) - 0.5;
    c.axis = normalize(mix(float3(0, 1, 0), out, p.spread) + wob * 0.3);

    // Its lives: it starts after a delay, top crystals first, and each life
    // lasts 18 to 26 s of shader time.
    float2 life = bloomLife(k);
    float span = life.y;
    float local = T - life.x;
    c.life = floor(max(local, 0.0) / span);
    float f = fract(max(local, 0.0) / span);

    // Grows over the first 30% of a life, melts back over the last 12%.
    c.grown = local < 0.0 ? 0.0
            : smoothstep(0.0, 0.3, f) * (1.0 - smoothstep(0.88, 1.0, f));

    // A new length every life: a few long ones, many short.
    float h = hashU(k * 7u + uint(c.life) * 131u + 9u);
    c.full = p.size * (0.2 + 0.8 * h * h) * (0.65 + 0.35 * max(d.y, 0.0));
    return c;
}

// How many corners and triangles a piece asks for this frame. The rock is
// a fixed 34 and 64; a crystal grows a ring at a time, and asks for none
// before it has started.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    if (piece == 0u) return uint2(34, 64);
    Crystal c = bloomCrystal(piece - 1u, u.time * p.speed, p);
    if (c.grown <= 0.001) return uint2(0, 0);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    return uint2(s * (rings + 1u) + 1u, 2u * s * rings + s);
}

// A rock corner: four rings of eight round a lumpy ball, then its two poles.
float3 rockCorner(uint i) {
    if (i >= 32u) {
        float y = i == 32u ? -1.0 : 1.0;
        return ROCK_AT + float3(0.0, y, 0.0) * ROCK * (0.9 + 0.2 * hashU(i));
    }
    uint ring = i / 8u, j = i % 8u;
    float lat = mix(-1.0, 0.95, float(ring) / 3.0);
    float a = (float(j) + 0.5 * float(ring)) * TAU / 8.0;
    float3 d = float3(cos(lat) * cos(a), sin(lat), cos(lat) * sin(a));
    return ROCK_AT + d * ROCK * (0.82 + 0.3 * hashU(i + 50u));
}

Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    Vertex o;
    float aspect = u.resolution.x / u.resolution.y;
    if (piece == 0u) {
        o.world = rockCorner(i);
        o.position = bloomScreen(o.world, u.time, aspect);
        o.along = 0.0; o.rock = 1.0;
        return o;
    }

    uint k = piece - 1u;
    Crystal c = bloomCrystal(k, u.time * p.speed, p);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    float len = c.full * c.grown;
    float radius = (0.02 + 0.06 * c.full) * (0.35 + 0.65 * c.grown);
    float tip = radius * 1.6;

    // Ring j sits at j steps up, but never past the crystal's length, so a
    // new ring starts on top of the last one and slides up out of it.
    float3 local;
    float up;
    if (i < s * (rings + 1u)) {
        uint j = i / s, side = i % s;
        up = min(float(j) * c.full / float(most), len);
        float a = float(side) * TAU / float(s) + hashU(k + 900u) * TAU;
        local = float3(cos(a) * radius, up, sin(a) * radius);
    } else {
        up = len + tip;
        local = float3(0.0, up, 0.0);
    }

    // Turned so that its y runs along the crystal's axis.
    float3 d = c.axis;
    float3 t = normalize(cross(abs(d.y) < 0.99 ? float3(0, 1, 0)
                                               : float3(1, 0, 0), d));
    float3 b = cross(d, t);
    o.world = c.root + t * local.x + d * local.y + b * local.z;
    o.position = bloomScreen(o.world, u.time, aspect);

    o.along = clamp(up / max(len + tip, 0.001), 0.0, 1.0);
    o.rock = 0.0;
    return o;
}

// The triangles: two for each face of each band between rings, then a fan
// of s up to the tip. The rock's are bands between its rings and fans to
// its poles.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = piece == 0u ? 8u : bloomSides(p);
    uint bands = piece == 0u ? 3u : (size(piece, u, p).y - s) / (2u * s);
    if (i < 2u * s * bands) {
        uint j = i / (2u * s), q = i % (2u * s), a = q / 2u;
        uint b = (a + 1u) % s;
        uint c00 = j * s + a, c01 = j * s + b;
        uint c10 = c00 + s, c11 = c01 + s;
        return q % 2u == 0u ? uint3(c00, c01, c11) : uint3(c00, c11, c10);
    }
    uint a = (i - 2u * s * bands) % s, b = (a + 1u) % s;
    if (piece == 0u) {
        if (i - 2u * s * bands < s) return uint3(a, b, 32u);
        return uint3(24u + a, 24u + b, 33u);
    }
    uint top = bands * s;
    return uint3(top + a, top + b, top + s);
}

// Once for every pixel of the triangles. Each face is flat: its normal is
// the cross of how the world point changes across the screen, so the faces
// catch the light one by one, as cut stones do.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;

    // The rock: dark, lit by the light above and by the crystals' glow,
    // warmest near its top where they grow.
    if (v.rock > 0.5) {
        float warm = smoothstep(-0.2, 0.3, v.world.y);
        float3 c = kiln(0.16 + shift) * (0.35 + 0.9 * diffuse)
                 + kiln(0.4 + shift) * warm * 0.12;
        return float4(c, 1.0);
    }

    // A crystal: ember at the root to gold up the tip, and a glint on faces
    // turned to the light and the eye.
    float3 body = kiln(0.38 + 0.22 * v.along + shift);
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float3 c = body * (0.3 + 0.7 * diffuse)
             + float3(1.0, 0.95, 0.87) * glint * 0.9;
    return float4(c, 1.0);
}

// Behind them: the dark of the kiln, with a glow round the cluster.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// Once for every pixel of the triangles. Each face is flat: its normal is
// the cross of how the world point changes across the screen, so the faces
// catch the light one by one, as cut stones do.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;

    // The rock: dark, lit by the light above and by the crystals' glow,
    // warmest near its top where they grow.
    if (v.rock > 0.5) {
        float warm = smoothstep(-0.2, 0.3, v.world.y);
        float3 c = kiln(0.16 + shift) * (0.35 + 0.9 * diffuse)
                 + kiln(0.4 + shift) * warm * 0.12;
        return float4(c, 1.0);
    }

    // A crystal: its body colour from ember near the root to gold up the
    // tip, a glow where it is hot, a glint on faces turned to the light and
    // the eye, and a rim where it is seen edge on.
    float h = v.heat * p.glow;
    float3 body = kiln(0.34 + 0.18 * v.tint + 0.22 * v.along + shift);
    float3 burn = kiln(0.55 + 0.45 * min(h, 1.0) + shift) * h;
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float rim = pow(1.0 - max(dot(n, toEye), 0.0), 3.0);
    float3 c = body * (0.3 + 0.7 * diffuse) + burn * 0.9
             + float3(1.0, 0.95, 0.87) * glint * 0.9
             + kiln(0.7 + shift) * rim * 0.45;
    return float4(c, 1.0);
}

// Ray Kiln · Crystal Bloom
//
// @episode Crystal Bloom
// @length 60
// @variation 30
// @variations 1,3,4,5,8,9,11,12
// @still 12.0
// @teaches Geometry made on the GPU with mesh shaders
// @category Mesh
// @tags mesh shaders, crystals, flat shading, golden angle, emission, fresnel, procedural sound
// @final Crystal bloom | Crystals push out of a dark rock ring by ring, their growing tips white-hot, cool to gold once grown, and melt back to grow again, each one chiming as it finishes.
//
// A cluster of crystals growing out of a rock, made by mesh shaders. The
// geometry is not a fixed grid: it is made fresh every frame, in 96 pieces.
// Piece 0 is the rock, a lumpy ball of 64 triangles. Every other piece is one
// crystal, a prism with a pointed tip, and it is built in rings: as it grows
// it adds rings, so the number of corners and triangles it asks for changes
// from frame to frame, which only a mesh shader can do. A crystal that has not
// started yet asks for none and costs nothing. The growing tip burns
// white-hot, a grown crystal cools to gold, and after a while each one melts
// back into the rock, glowing again, and grows anew.
// Each @param line below becomes a slider in the host app and a row in the page's parameters table:
//   // @param name default min max description
//
// @mesh 96

// @param speed 1 0 3 How fast the crystals grow and melt back
// @param size 1 0.4 1.4 How long the largest crystals grow
// @param sides 6 3 8 Faces around each crystal
// @param spread 0.75 0 1 How far the crystals lean out from the rock: at 0 they all stand up
// @param glow 1 0 2 How hot the growing tips burn
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// The rock's size and middle.
constant float3 ROCK = float3(0.55, 0.36, 0.55);
constant float3 ROCK_AT = float3(0.0, -0.1, 0.0);
// Every piece but the rock is a crystal.
constant uint CRYSTALS = 95;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, how hot it burns, how far up its crystal it is,
// the crystal's own tint, and 1 on the rock.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float heat;
    float along;
    float tint;
    float rock;
};

// Cracked Clay's hash: a whole number in, a random 0 to 1 out.
float hashU(uint n) {
    n = (n << 13u) ^ n;
    n = n * (n * n * 15731u + 789221u) + 1376312589u;
    return float(n & 0x7fffffffu) / float(0x7fffffff);
}

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many faces go round a crystal, and how many rings it may have:
// as many as keep it to 64 corners and 64 triangles.
uint bloomSides(constant Params& p) {
    return uint(clamp(floor(p.sides + 0.5), 3.0, 8.0));
}
uint bloomRings(uint s) {
    return min((64u - s) / (2u * s), 63u / s - 1u);
}

// Where crystal k points: over the rock by the golden angle.
float3 bloomDir(uint k) {
    float y = 0.95 - float(k) / float(CRYSTALS - 1) * 1.25;
    float r = sqrt(max(1.0 - y * y, 0.0));
    float a = float(k) * 2.39996;
    return float3(cos(a) * r, y, sin(a) * r);
}

// When crystal k first starts, and how long each of its lives lasts.
float2 bloomLife(uint k) {
    float delay = 0.4 + 7.0 * float(k) / float(CRYSTALS);
    return float2(delay + 1.2 * hashU(k + 300u), 18.0 + 8.0 * hashU(k + 600u));
}

// Everything about crystal k at shader time T: where its root is, which way
// it grows, how long it is when grown, how far grown it is (0 to 1), how
// hot it burns and which life it is in.
struct Crystal {
    float3 root;
    float3 axis;
    float full;
    float grown;
    float glow;
    float life;
};

Crystal bloomCrystal(uint k, float T, constant Params& p) {
    Crystal c;
    float3 d = bloomDir(k);
    c.root = ROCK_AT + d * ROCK * 0.8;
    float3 out = normalize(d / ROCK);
    float3 wob = float3(hashU(k * 5u + 1u), hashU(k * 5u + 2u),
                        hashU(k * 5u + 3u)) - 0.5;
    c.axis = normalize(mix(float3(0, 1, 0), out, p.spread) + wob * 0.3);

    // Its lives: it starts after a delay, top crystals first, and each life
    // lasts 18 to 26 s of shader time.
    float2 life = bloomLife(k);
    float span = life.y;
    float local = T - life.x;
    c.life = floor(max(local, 0.0) / span);
    float f = fract(max(local, 0.0) / span);

    // Grows over the first 30% of a life, melts back over the last 12%.
    c.grown = local < 0.0 ? 0.0
            : smoothstep(0.0, 0.3, f) * (1.0 - smoothstep(0.88, 1.0, f));
    float cool = exp(-max(f - 0.3, 0.0) * span / 2.5);
    c.glow = f < 0.3 ? 1.0 : max(cool, smoothstep(0.86, 0.95, f));

    // A new length every life: a few long ones, many short.
    float h = hashU(k * 7u + uint(c.life) * 131u + 9u);
    c.full = p.size * (0.2 + 0.8 * h * h) * (0.65 + 0.35 * max(d.y, 0.0));
    return c;
}

// How many corners and triangles a piece asks for this frame. The rock is
// a fixed 34 and 64; a crystal grows a ring at a time, and asks for none
// before it has started.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    if (piece == 0u) return uint2(34, 64);
    Crystal c = bloomCrystal(piece - 1u, u.time * p.speed, p);
    if (c.grown <= 0.001) return uint2(0, 0);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    return uint2(s * (rings + 1u) + 1u, 2u * s * rings + s);
}

// A rock corner: four rings of eight round a lumpy ball, then its two poles.
float3 rockCorner(uint i) {
    if (i >= 32u) {
        float y = i == 32u ? -1.0 : 1.0;
        return ROCK_AT + float3(0.0, y, 0.0) * ROCK * (0.9 + 0.2 * hashU(i));
    }
    uint ring = i / 8u, j = i % 8u;
    float lat = mix(-1.0, 0.95, float(ring) / 3.0);
    float a = (float(j) + 0.5 * float(ring)) * TAU / 8.0;
    float3 d = float3(cos(lat) * cos(a), sin(lat), cos(lat) * sin(a));
    return ROCK_AT + d * ROCK * (0.82 + 0.3 * hashU(i + 50u));
}

Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    Vertex o;
    float aspect = u.resolution.x / u.resolution.y;
    if (piece == 0u) {
        o.world = rockCorner(i);
        o.position = bloomScreen(o.world, u.time, aspect);
        o.heat = 0.0; o.along = 0.0; o.tint = 0.0; o.rock = 1.0;
        return o;
    }

    uint k = piece - 1u;
    Crystal c = bloomCrystal(k, u.time * p.speed, p);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    float len = c.full * c.grown;
    float radius = (0.02 + 0.06 * c.full) * (0.35 + 0.65 * c.grown);
    float tip = radius * 1.6;

    // Ring j sits at j steps up, but never past the crystal's length, so a
    // new ring starts on top of the last one and slides up out of it.
    float3 local;
    float up;
    if (i < s * (rings + 1u)) {
        uint j = i / s, side = i % s;
        up = min(float(j) * c.full / float(most), len);
        float a = float(side) * TAU / float(s) + hashU(k + 900u) * TAU;
        local = float3(cos(a) * radius, up, sin(a) * radius);
    } else {
        up = len + tip;
        local = float3(0.0, up, 0.0);
    }

    // Turned so that its y runs along the crystal's axis.
    float3 d = c.axis;
    float3 t = normalize(cross(abs(d.y) < 0.99 ? float3(0, 1, 0)
                                               : float3(1, 0, 0), d));
    float3 b = cross(d, t);
    o.world = c.root + t * local.x + d * local.y + b * local.z;
    o.position = bloomScreen(o.world, u.time, aspect);

    // Hot where it is growing: the newest part of the crystal, near its tip.
    o.heat = c.glow * exp(-max(len - up, 0.0) / 0.12);
    o.along = clamp(up / max(len + tip, 0.001), 0.0, 1.0);
    o.tint = hashU(k + 1200u);
    o.rock = 0.0;
    return o;
}

// The triangles: two for each face of each band between rings, then a fan
// of s up to the tip. The rock's are bands between its rings and fans to
// its poles.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = piece == 0u ? 8u : bloomSides(p);
    uint bands = piece == 0u ? 3u : (size(piece, u, p).y - s) / (2u * s);
    if (i < 2u * s * bands) {
        uint j = i / (2u * s), q = i % (2u * s), a = q / 2u;
        uint b = (a + 1u) % s;
        uint c00 = j * s + a, c01 = j * s + b;
        uint c10 = c00 + s, c11 = c01 + s;
        return q % 2u == 0u ? uint3(c00, c01, c11) : uint3(c00, c11, c10);
    }
    uint a = (i - 2u * s * bands) % s, b = (a + 1u) % s;
    if (piece == 0u) {
        if (i - 2u * s * bands < s) return uint3(a, b, 32u);
        return uint3(24u + a, 24u + b, 33u);
    }
    uint top = bands * s;
    return uint3(top + a, top + b, top + s);
}

// Once for every pixel of the triangles. Each face is flat: its normal is
// the cross of how the world point changes across the screen, so the faces
// catch the light one by one, as cut stones do.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;

    // The rock: dark, lit by the light above and by the crystals' glow,
    // warmest near its top where they grow.
    if (v.rock > 0.5) {
        float warm = smoothstep(-0.2, 0.3, v.world.y);
        float3 c = kiln(0.16 + shift) * (0.35 + 0.9 * diffuse)
                 + kiln(0.4 + shift) * warm * 0.12;
        return float4(c, 1.0);
    }

    // A crystal: its body colour from ember near the root to gold up the
    // tip, a glow where it is hot, a glint on faces turned to the light and
    // the eye, and a rim where it is seen edge on.
    float h = v.heat * p.glow;
    float3 body = kiln(0.34 + 0.18 * v.tint + 0.22 * v.along + shift);
    float3 burn = kiln(0.55 + 0.45 * min(h, 1.0) + shift) * h;
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float rim = pow(1.0 - max(dot(n, toEye), 0.0), 3.0);
    float3 c = body * (0.3 + 0.7 * diffuse) + burn * 0.9
             + float3(1.0, 0.95, 0.87) * glint * 0.9
             + kiln(0.7 + shift) * rim * 0.45;
    return float4(c, 1.0);
}

// Behind them: the dark of the kiln, with a glow round the cluster.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// The sound: each crystal chimes as it finishes growing, a struck glass
// note whose pitch falls the longer the crystal, from where it is on screen,
// over a low hum that swells while many are growing.

// Kaleidoscope's scale: D minor pentatonic over three octaves.
float scaleNote(int i) {
    const float notes[16] = {
        293.66, 349.23, 392.00, 440.00, 523.25, 587.33, 698.46, 783.99,
        880.00, 1046.50, 1174.66, 1396.91, 1567.98, 1760.00, 2093.00, 2349.32};
    return notes[clamp(i, 0, 15)];
}

// A struck glass note `since` seconds after it was struck: four partials,
// the higher ones dying sooner, as a glass's do, and softer on high notes.
float chime(float f, float since) {
    float on = 1.0 - exp(-since * 500.0);
    float soft = clamp(900.0 / f, 0.0, 1.0);
    float s = sin(TAU * f * since) * exp(-since * 2.4)
            + 0.5 * sin(TAU * f * 2.32 * since) * exp(-since * 4.0)
            + 0.25 * soft * sin(TAU * f * 4.25 * since) * exp(-since * 7.0)
            + 0.12 * soft * sin(TAU * f * 6.63 * since) * exp(-since * 11.0);
    return s * on;
}

float2 sound(float t, constant Params& p) {
    float T = t * p.speed;
    float rate = max(p.speed, 0.01);
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));

    float2 out = float2(0.0);
    float growing = 0.0;
    for (uint k = 0u; k < CRYSTALS; k++) {
        Crystal c = bloomCrystal(k, T, p);
        float2 life = bloomLife(k);
        float span = life.y;
        float local = T - life.x;
        if (local < 0.0) continue;
        float f0 = fract(local / span);
        growing += f0 < 0.3 ? 1.0 : 0.0;

        // When this crystal last finished growing, in seconds.
        float done = (local - 0.3 * span) / span;
        if (done < 0.0) continue;
        float since = fract(done) * span / rate;
        if (since > 3.0) continue;
        float h = hashU(k * 7u + uint(floor(done)) * 131u + 9u);
        int note = 13 - int(h * h * 9.0) - int(hashU(k + 1500u) * 3.0);
        float pan = clamp(dot(c.root - eye, r) * 0.9, -1.0, 1.0) * 0.5 + 0.5;
        float s = chime(scaleNote(note), since) * 0.1;
        out += s * float2(sqrt(1.0 - pan), sqrt(pan));
    }

    // The hum: a low D and the A above it, a little apart in each ear so
    // they beat slowly, louder while more crystals grow.
    float swell = 0.05 + 0.11 * min(growing / 25.0, 1.0);
    float hum = sin(TAU * 73.42 * t) + 0.6 * sin(TAU * 110.0 * t);
    float hum2 = sin(TAU * 73.62 * t) + 0.6 * sin(TAU * 110.3 * t);
    out += swell * float2(hum, hum2) * 0.8;

    out = tanh(out * 1.2) / 1.2;
    return out * p.volume * smoothstep(0.0, 0.5, t);
}

// Ray Kiln · Crystal Bloom
//
// @episode Crystal Bloom
// @length 60
// @variation 30
// @variations 1,3,4,5,8,9,11,12
// @still 12.0
// @teaches Geometry made on the GPU with mesh shaders
// @category Mesh
// @tags mesh shaders, crystals, flat shading, golden angle, emission, fresnel, procedural sound
// @final Crystal bloom | Crystals push out of a dark rock ring by ring, their growing tips white-hot, cool to gold once grown, and melt back to grow again, each one chiming as it finishes.
//
// A cluster of crystals growing out of a rock, made by mesh shaders. The
// geometry is not a fixed grid: it is made fresh every frame, in 96 pieces.
// Piece 0 is the rock, a lumpy ball of 64 triangles. Every other piece is one
// crystal, a prism with a pointed tip, and it is built in rings: as it grows
// it adds rings, so the number of corners and triangles it asks for changes
// from frame to frame, which only a mesh shader can do. A crystal that has not
// started yet asks for none and costs nothing. The growing tip burns
// white-hot, a grown crystal cools to gold, and after a while each one melts
// back into the rock, glowing again, and grows anew.
// Each @param line below becomes a slider in the host app and a row in the page's parameters table:
//   // @param name default min max description
//
// @mesh 96

// @param speed 1 0 3 How fast the crystals grow and melt back
// @param size 1 0.4 1.4 How long the largest crystals grow
// @param sides 6 3 8 Faces around each crystal
// @param spread 0.75 0 1 How far the crystals lean out from the rock: at 0 they all stand up
// @param glow 1 0 2 How hot the growing tips burn
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// The rock's size and middle.
constant float3 ROCK = float3(0.55, 0.36, 0.55);
constant float3 ROCK_AT = float3(0.0, -0.1, 0.0);
// Every piece but the rock is a crystal.
constant uint CRYSTALS = 95;

// What each corner carries to the pixels between: where it lands on screen,
// where it is in the world, how hot it burns, how far up its crystal it is,
// the crystal's own tint, and 1 on the rock.
struct Vertex {
    float4 position [[position]];
    float3 world;
    float heat;
    float along;
    float tint;
    float rock;
};

// Cracked Clay's hash: a whole number in, a random 0 to 1 out.
float hashU(uint n) {
    n = (n << 13u) ^ n;
    n = n * (n * n * 15731u + 789221u) + 1376312589u;
    return float(n & 0x7fffffffu) / float(0x7fffffff);
}

// Flowing Color Field's palette: near-black, ember, orange, gold, white-hot.
float3 kiln(float t) {
    float3 c0 = float3(0.020, 0.012, 0.020);
    float3 c1 = float3(0.280, 0.040, 0.030);
    float3 c2 = float3(0.880, 0.260, 0.050);
    float3 c3 = float3(1.000, 0.680, 0.200);
    float3 c4 = float3(1.000, 0.970, 0.840);
    t = clamp(t, 0.0, 1.0);
    float3 c = mix(c0, c1, smoothstep(0.00, 0.22, t));
    c = mix(c, c2, smoothstep(0.22, 0.50, t));
    c = mix(c, c3, smoothstep(0.50, 0.76, t));
    return mix(c, c4, smoothstep(0.76, 1.00, t));
}

// The eye, going slowly round the rock and rising and falling a little.
float3 bloomEye(float t) {
    float a = t * 0.12 + 0.6;
    return float3(sin(a) * 2.5, 1.05 + 0.25 * sin(t * 0.09), cos(a) * 2.5);
}

// A point in the world as the host wants it: x and y on screen (times the
// depth), the depth from 0 to 1 for the depth test, and the depth itself.
float4 bloomScreen(float3 world, float t, float aspect) {
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));
    float3 up = cross(f, r);
    float3 v = world - eye;
    float3 view = float3(dot(v, r), dot(v, up), dot(v, f));
    float zoom = 1.95, near = 0.1, far = 20.0;
    float depth = (view.z - near) / (far - near) * view.z;
    return float4(view.x * zoom / aspect, view.y * zoom, depth, view.z);
}

// How many faces go round a crystal, and how many rings it may have:
// as many as keep it to 64 corners and 64 triangles.
uint bloomSides(constant Params& p) {
    return uint(clamp(floor(p.sides + 0.5), 3.0, 8.0));
}
uint bloomRings(uint s) {
    return min((64u - s) / (2u * s), 63u / s - 1u);
}

// Where crystal k points: over the rock by the golden angle.
float3 bloomDir(uint k) {
    float y = 0.95 - float(k) / float(CRYSTALS - 1) * 1.25;
    float r = sqrt(max(1.0 - y * y, 0.0));
    float a = float(k) * 2.39996;
    return float3(cos(a) * r, y, sin(a) * r);
}

// When crystal k first starts, and how long each of its lives lasts.
float2 bloomLife(uint k) {
    float delay = 0.4 + 7.0 * float(k) / float(CRYSTALS);
    return float2(delay + 1.2 * hashU(k + 300u), 18.0 + 8.0 * hashU(k + 600u));
}

// Everything about crystal k at shader time T: where its root is, which way
// it grows, how long it is when grown, how far grown it is (0 to 1), how
// hot it burns and which life it is in.
struct Crystal {
    float3 root;
    float3 axis;
    float full;
    float grown;
    float glow;
    float life;
};

Crystal bloomCrystal(uint k, float T, constant Params& p) {
    Crystal c;
    float3 d = bloomDir(k);
    c.root = ROCK_AT + d * ROCK * 0.8;
    float3 out = normalize(d / ROCK);
    float3 wob = float3(hashU(k * 5u + 1u), hashU(k * 5u + 2u),
                        hashU(k * 5u + 3u)) - 0.5;
    c.axis = normalize(mix(float3(0, 1, 0), out, p.spread) + wob * 0.3);

    // Its lives: it starts after a delay, top crystals first, and each life
    // lasts 18 to 26 s of shader time.
    float2 life = bloomLife(k);
    float span = life.y;
    float local = T - life.x;
    c.life = floor(max(local, 0.0) / span);
    float f = fract(max(local, 0.0) / span);

    // Grows over the first 30% of a life, melts back over the last 12%.
    c.grown = local < 0.0 ? 0.0
            : smoothstep(0.0, 0.3, f) * (1.0 - smoothstep(0.88, 1.0, f));
    float cool = exp(-max(f - 0.3, 0.0) * span / 2.5);
    c.glow = f < 0.3 ? 1.0 : max(cool, smoothstep(0.86, 0.95, f));

    // A new length every life: a few long ones, many short.
    float h = hashU(k * 7u + uint(c.life) * 131u + 9u);
    c.full = p.size * (0.2 + 0.8 * h * h) * (0.65 + 0.35 * max(d.y, 0.0));
    return c;
}

// How many corners and triangles a piece asks for this frame. The rock is
// a fixed 34 and 64; a crystal grows a ring at a time, and asks for none
// before it has started.
uint2 size(uint piece, constant Uniforms& u, constant Params& p) {
    if (piece == 0u) return uint2(34, 64);
    Crystal c = bloomCrystal(piece - 1u, u.time * p.speed, p);
    if (c.grown <= 0.001) return uint2(0, 0);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    return uint2(s * (rings + 1u) + 1u, 2u * s * rings + s);
}

// A rock corner: four rings of eight round a lumpy ball, then its two poles.
float3 rockCorner(uint i) {
    if (i >= 32u) {
        float y = i == 32u ? -1.0 : 1.0;
        return ROCK_AT + float3(0.0, y, 0.0) * ROCK * (0.9 + 0.2 * hashU(i));
    }
    uint ring = i / 8u, j = i % 8u;
    float lat = mix(-1.0, 0.95, float(ring) / 3.0);
    float a = (float(j) + 0.5 * float(ring)) * TAU / 8.0;
    float3 d = float3(cos(lat) * cos(a), sin(lat), cos(lat) * sin(a));
    return ROCK_AT + d * ROCK * (0.82 + 0.3 * hashU(i + 50u));
}

Vertex corner(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    Vertex o;
    float aspect = u.resolution.x / u.resolution.y;
    if (piece == 0u) {
        o.world = rockCorner(i);
        o.position = bloomScreen(o.world, u.time, aspect);
        o.heat = 0.0; o.along = 0.0; o.tint = 0.0; o.rock = 1.0;
        return o;
    }

    uint k = piece - 1u;
    Crystal c = bloomCrystal(k, u.time * p.speed, p);
    uint s = bloomSides(p), most = bloomRings(s);
    uint rings = clamp(uint(ceil(c.grown * float(most))), 1u, most);
    float len = c.full * c.grown;
    float radius = (0.02 + 0.06 * c.full) * (0.35 + 0.65 * c.grown);
    float tip = radius * 1.6;

    // Ring j sits at j steps up, but never past the crystal's length, so a
    // new ring starts on top of the last one and slides up out of it.
    float3 local;
    float up;
    if (i < s * (rings + 1u)) {
        uint j = i / s, side = i % s;
        up = min(float(j) * c.full / float(most), len);
        float a = float(side) * TAU / float(s) + hashU(k + 900u) * TAU;
        local = float3(cos(a) * radius, up, sin(a) * radius);
    } else {
        up = len + tip;
        local = float3(0.0, up, 0.0);
    }

    // Turned so that its y runs along the crystal's axis.
    float3 d = c.axis;
    float3 t = normalize(cross(abs(d.y) < 0.99 ? float3(0, 1, 0)
                                               : float3(1, 0, 0), d));
    float3 b = cross(d, t);
    o.world = c.root + t * local.x + d * local.y + b * local.z;
    o.position = bloomScreen(o.world, u.time, aspect);

    // Hot where it is growing: the newest part of the crystal, near its tip.
    o.heat = c.glow * exp(-max(len - up, 0.0) / 0.12);
    o.along = clamp(up / max(len + tip, 0.001), 0.0, 1.0);
    o.tint = hashU(k + 1200u);
    o.rock = 0.0;
    return o;
}

// The triangles: two for each face of each band between rings, then a fan
// of s up to the tip. The rock's are bands between its rings and fans to
// its poles.
uint3 triangle(uint piece, uint i, constant Uniforms& u, constant Params& p) {
    uint s = piece == 0u ? 8u : bloomSides(p);
    uint bands = piece == 0u ? 3u : (size(piece, u, p).y - s) / (2u * s);
    if (i < 2u * s * bands) {
        uint j = i / (2u * s), q = i % (2u * s), a = q / 2u;
        uint b = (a + 1u) % s;
        uint c00 = j * s + a, c01 = j * s + b;
        uint c10 = c00 + s, c11 = c01 + s;
        return q % 2u == 0u ? uint3(c00, c01, c11) : uint3(c00, c11, c10);
    }
    uint a = (i - 2u * s * bands) % s, b = (a + 1u) % s;
    if (piece == 0u) {
        if (i - 2u * s * bands < s) return uint3(a, b, 32u);
        return uint3(24u + a, 24u + b, 33u);
    }
    uint top = bands * s;
    return uint3(top + a, top + b, top + s);
}

// Once for every pixel of the triangles. Each face is flat: its normal is
// the cross of how the world point changes across the screen, so the faces
// catch the light one by one, as cut stones do.
float4 shade(Vertex v, constant Uniforms& u, constant Params& p) {
    float3 n = normalize(cross(dfdx(v.world), dfdy(v.world)));
    float3 toEye = normalize(bloomEye(u.time) - v.world);
    if (dot(n, toEye) < 0.0) n = -n;
    float3 light = normalize(float3(0.4, 1.0, 0.5));
    float diffuse = max(dot(n, light), 0.0);
    float shift = (p.heat - 0.5) * 0.4;

    // The rock: dark, lit by the light above and by the crystals' glow,
    // warmest near its top where they grow.
    if (v.rock > 0.5) {
        float warm = smoothstep(-0.2, 0.3, v.world.y);
        float3 c = kiln(0.16 + shift) * (0.35 + 0.9 * diffuse)
                 + kiln(0.4 + shift) * warm * 0.12;
        return float4(c, 1.0);
    }

    // A crystal: its body colour from ember near the root to gold up the
    // tip, a glow where it is hot, a glint on faces turned to the light and
    // the eye, and a rim where it is seen edge on.
    float h = v.heat * p.glow;
    float3 body = kiln(0.34 + 0.18 * v.tint + 0.22 * v.along + shift);
    float3 burn = kiln(0.55 + 0.45 * min(h, 1.0) + shift) * h;
    float glint = pow(max(dot(reflect(-light, n), toEye), 0.0), 40.0);
    float rim = pow(1.0 - max(dot(n, toEye), 0.0), 3.0);
    float3 c = body * (0.3 + 0.7 * diffuse) + burn * 0.9
             + float3(1.0, 0.95, 0.87) * glint * 0.9
             + kiln(0.7 + shift) * rim * 0.45;
    return float4(c, 1.0);
}

// Behind them: the dark of the kiln, with a glow round the cluster.
float4 background(float2 uv, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 q = (uv - float2(0.5, 0.5)) * float2(aspect, 1.0);
    float glow = exp(-length(q) * 2.6);
    return float4(kiln(0.03 + 0.17 * glow + (p.heat - 0.5) * 0.15), 1.0);
}

// The sound: each crystal chimes as it finishes growing, a struck glass
// note whose pitch falls the longer the crystal, from where it is on screen,
// over a low hum that swells while many are growing.

// Kaleidoscope's scale: D minor pentatonic over three octaves.
float scaleNote(int i) {
    const float notes[16] = {
        293.66, 349.23, 392.00, 440.00, 523.25, 587.33, 698.46, 783.99,
        880.00, 1046.50, 1174.66, 1396.91, 1567.98, 1760.00, 2093.00, 2349.32};
    return notes[clamp(i, 0, 15)];
}

// A struck glass note `since` seconds after it was struck: four partials,
// the higher ones dying sooner, as a glass's do, and softer on high notes.
float chime(float f, float since) {
    float on = 1.0 - exp(-since * 500.0);
    float soft = clamp(900.0 / f, 0.0, 1.0);
    float s = sin(TAU * f * since) * exp(-since * 2.4)
            + 0.5 * sin(TAU * f * 2.32 * since) * exp(-since * 4.0)
            + 0.25 * soft * sin(TAU * f * 4.25 * since) * exp(-since * 7.0)
            + 0.12 * soft * sin(TAU * f * 6.63 * since) * exp(-since * 11.0);
    return s * on;
}

float2 sound(float t, constant Params& p) {
    float T = t * p.speed;
    float rate = max(p.speed, 0.01);
    float3 eye = bloomEye(t);
    float3 f = normalize(float3(0.0, 0.25, 0.0) - eye);
    float3 r = normalize(cross(float3(0, 1, 0), f));

    float2 out = float2(0.0);
    float growing = 0.0;
    for (uint k = 0u; k < CRYSTALS; k++) {
        Crystal c = bloomCrystal(k, T, p);
        float2 life = bloomLife(k);
        float span = life.y;
        float local = T - life.x;
        if (local < 0.0) continue;
        float f0 = fract(local / span);
        growing += f0 < 0.3 ? 1.0 : 0.0;

        // When this crystal last finished growing, in seconds.
        float done = (local - 0.3 * span) / span;
        if (done < 0.0) continue;
        float since = fract(done) * span / rate;
        if (since > 3.0) continue;
        float h = hashU(k * 7u + uint(floor(done)) * 131u + 9u);
        int note = 13 - int(h * h * 9.0) - int(hashU(k + 1500u) * 3.0);
        float pan = clamp(dot(c.root - eye, r) * 0.9, -1.0, 1.0) * 0.5 + 0.5;
        float s = chime(scaleNote(note), since) * 0.1;
        out += s * float2(sqrt(1.0 - pan), sqrt(pan));
    }

    // The hum: a low D and the A above it, a little apart in each ear so
    // they beat slowly, louder while more crystals grow.
    float swell = 0.05 + 0.11 * min(growing / 25.0, 1.0);
    float hum = sin(TAU * 73.42 * t) + 0.6 * sin(TAU * 110.0 * t);
    float hum2 = sin(TAU * 73.62 * t) + 0.6 * sin(TAU * 110.3 * t);
    out += swell * float2(hum, hum2) * 0.8;

    out = tanh(out * 1.2) / 1.2;
    return out * p.volume * smoothstep(0.0, 0.5, t);
}

struct Uniforms {
    var resolution: SIMD2<Float>
    var time: Float
    var timeDelta: Float
    var mouse: SIMD2<Float>
    var frame: UInt32
    var pad: UInt32 = 0
}

let library = try device.makeLibrary(source: standaloneSource, options: nil)
let descriptor = MTLRenderPipelineDescriptor()
descriptor.vertexFunction = library.makeFunction(name: "rk_vertex")
descriptor.fragmentFunction = library.makeFunction(name: "rk_fragment")
descriptor.colorAttachments[0].pixelFormat = .bgra8Unorm
let pipeline = try device.makeRenderPipelineState(descriptor: descriptor)

// Each frame, inside a render pass:
var uniforms = Uniforms(resolution: size, time: time, timeDelta: dt, mouse: mouse, frame: frame)
var params: [Float] = [1, 1, 6, 0.75, 1, 0.5, 0.8]     // the @param defaults, in the order they are declared
encoder.setRenderPipelineState(pipeline)
encoder.setFragmentBytes(&uniforms, length: MemoryLayout<Uniforms>.stride, index: 0)
encoder.setFragmentBytes(&params, length: params.count * MemoryLayout<Float>.stride, index: 1)
encoder.drawPrimitives(type: .triangle, vertexStart: 0, vertexCount: 3)