Curl-Noise Flow

// Stage 1 · Downhill on noise
// Noise makes a landscape of hills and hollows. The obvious way to make a flow from it is to let every spark roll
// downhill. It moves for a moment, then the sparks pile up in the valleys and the bottoms of the hollows, and the
// rest of the picture empties out.
//
// @particles 150000

// @param speed 1.0 0.0 3.0 How fast the sparks roll
// @param scale 1.6 0.6 4.0 How small the hills are
// @param trails 0.9 0.6 0.97 How much of the last image is left at each step

// A spark is just where it is.
struct Particle {
    float2 position;
};

// Episode 5's integer hash: a whole number in, a random number from 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);
}

// Episode 8's corner hash, keeping all its bits: a whole number for each whole-numbered point in 3D.
uint hashCorner(float3 p) {
    uint3 v = uint3(int3(floor(p))) * 1664525u + 1013904223u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    v ^= v >> 16u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    return v.x;
}

// A corner's slope: a random direction in 3D from its hash (a height, and an angle
// around), dotted with the way from the corner to the point.
float cornerSlope(uint h, float3 d) {
    float z = float(h & 0xffffu) / 32767.5 - 1.0;
    float a = float(h >> 16u) * (6.2831853 / 65536.0);
    float s = sqrt(max(1.0 - z * z, 0.0));
    return 1.5 * dot(float3(s * cos(a), s * sin(a), z), d);
}

// One corner's part of the noise at the point: its slope, from the corner to there.
float cornerPart(float3 i, float3 f, float3 c) {
    return cornerSlope(hashCorner(i + c), f - c);
}

// Gradient noise: a random slope at each corner instead of a random height, so the
// hills don't line up with the grid. Roughly -1 to 1.
float gnoise(float3 p) {
    float3 i = floor(p), f = p - i;
    float3 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);
    float n000 = cornerPart(i, f, float3(0, 0, 0));
    float n100 = cornerPart(i, f, float3(1, 0, 0));
    float n010 = cornerPart(i, f, float3(0, 1, 0));
    float n110 = cornerPart(i, f, float3(1, 1, 0));
    float n001 = cornerPart(i, f, float3(0, 0, 1));
    float n101 = cornerPart(i, f, float3(1, 0, 1));
    float n011 = cornerPart(i, f, float3(0, 1, 1));
    float n111 = cornerPart(i, f, float3(1, 1, 1));
    float lo = mix(mix(n000, n100, w.x), mix(n010, n110, w.x), w.y);
    float hi = mix(mix(n001, n101, w.x), mix(n011, n111, w.x), w.y);
    return mix(lo, hi, w.z);
}

// The landscape: one octave of gradient noise, standing still.
float potential(float2 x, constant Params& p) {
    return gnoise(float3(x * p.scale, 0.0)) / p.scale;
}

// The slope of the landscape, measured with four samples a small step apart: which way is uphill, and how steep.
float2 slope(float2 x, constant Params& p) {
    float2 a = float2(0.002, 0.0), b = float2(0.0, 0.002);
    float dx = potential(x + a, p) - potential(x - a, p);
    float dy = potential(x + b, p) - potential(x - b, p);
    return float2(dx, dy) / 0.004;
}

// Somewhere on the picture or a way past its edges: a spark off the picture isn't
// drawn, but it can come back on.
float2 scatter(uint k) {
    float2 r = float2(hashU(k), hashU(k + 1u)) * 2.0 - 1.0;
    return r * float2(16.0 / 9.0 + 0.4, 1.4);
}

// Each spark starts anywhere.
Particle spawn(uint id, constant Params& p) {
    Particle q;
    q.position = scatter(id * 4u);
    return q;
}

// One step: the spark rolls against the slope, downhill, and stops where the
// ground is flat, with every spark that rolled the same way.
void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    // Every three seconds or so, a spark is dropped somewhere new.
    if (hashU(id * 7u + u.frame * 7919u) < u.timeDelta / 3.0) {
        q.position = scatter(id * 5u + u.frame * 7919u);
    }
    q.position -= slope(q.position, p) * p.speed * 0.3 * u.timeDelta;
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    return float3(q.position / float2(aspect, 1.0) * 0.5 + 0.5, 0.003);
}

float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float edge = smoothstep(1.0, 0.0, length(at));
    float bright = edge * 0.3 * (1.0 - p.trails) / 0.1;
    return float4(float3(1.0, 0.59, 0.20) * bright, 1.0);
}

// Episode 10's trails: each step keeps most of the last one's image, so every spark draws its path.
float4 fade(float2 uv, Trails previous,
            constant Uniforms& u, constant Params& p) {
    return previous.at(uv) * p.trails;
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float3 c = trails.at(uv).rgb;
    return float4(c / (1.0 + c), 1.0);
}

// Stage 2 · Turned a quarter: the curl
// The same landscape, but each spark goes sideways to the slope instead of down it: along the contour lines, the
// way a path keeps to one height around a hill. Nothing can gather anywhere, so the sparks stream round and round.
//
// @particles 150000

// @param speed 1.0 0.0 3.0 How fast the sparks roll
// @param scale 1.6 0.6 4.0 How small the hills are
// @param trails 0.9 0.6 0.97 How much of the last image is left at each step

// A spark is just where it is.
struct Particle {
    float2 position;
};

// Episode 5's integer hash: a whole number in, a random number from 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);
}

// Episode 8's corner hash, keeping all its bits: a whole number for each whole-numbered point in 3D.
uint hashCorner(float3 p) {
    uint3 v = uint3(int3(floor(p))) * 1664525u + 1013904223u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    v ^= v >> 16u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    return v.x;
}

// A corner's slope: a random direction in 3D from its hash (a height, and an angle
// around), dotted with the way from the corner to the point.
float cornerSlope(uint h, float3 d) {
    float z = float(h & 0xffffu) / 32767.5 - 1.0;
    float a = float(h >> 16u) * (6.2831853 / 65536.0);
    float s = sqrt(max(1.0 - z * z, 0.0));
    return 1.5 * dot(float3(s * cos(a), s * sin(a), z), d);
}

// One corner's part of the noise at the point: its slope, from the corner to there.
float cornerPart(float3 i, float3 f, float3 c) {
    return cornerSlope(hashCorner(i + c), f - c);
}

// Gradient noise: a random slope at each corner instead of a random height, so the
// hills don't line up with the grid. Roughly -1 to 1.
float gnoise(float3 p) {
    float3 i = floor(p), f = p - i;
    float3 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);
    float n000 = cornerPart(i, f, float3(0, 0, 0));
    float n100 = cornerPart(i, f, float3(1, 0, 0));
    float n010 = cornerPart(i, f, float3(0, 1, 0));
    float n110 = cornerPart(i, f, float3(1, 1, 0));
    float n001 = cornerPart(i, f, float3(0, 0, 1));
    float n101 = cornerPart(i, f, float3(1, 0, 1));
    float n011 = cornerPart(i, f, float3(0, 1, 1));
    float n111 = cornerPart(i, f, float3(1, 1, 1));
    float lo = mix(mix(n000, n100, w.x), mix(n010, n110, w.x), w.y);
    float hi = mix(mix(n001, n101, w.x), mix(n011, n111, w.x), w.y);
    return mix(lo, hi, w.z);
}

// The landscape: one octave of gradient noise, standing still.
float potential(float2 x, constant Params& p) {
    return gnoise(float3(x * p.scale, 0.0)) / p.scale;
}

// The curl: the slope (dx, dy) turned a quarter, to (dy, -dx). It is as long as the slope is steep, but it points
// along the hill, never up or down it. What flows into any small patch flows out of it again.
float2 curl(float2 x, constant Params& p) {
    float2 a = float2(0.002, 0.0), b = float2(0.0, 0.002);
    float dx = potential(x + a, p) - potential(x - a, p);
    float dy = potential(x + b, p) - potential(x - b, p);
    return float2(dy, -dx) / 0.004;
}

// Somewhere on the picture or a way past its edges: a spark off the picture isn't
// drawn, but it can come back on.
float2 scatter(uint k) {
    float2 r = float2(hashU(k), hashU(k + 1u)) * 2.0 - 1.0;
    return r * float2(16.0 / 9.0 + 0.4, 1.4);
}

// Each spark starts anywhere.
Particle spawn(uint id, constant Params& p) {
    Particle q;
    q.position = scatter(id * 4u);
    return q;
}

// One step: the spark goes along its contour line. Since the landscape stands still, each spark goes round its own
// loop until it is dropped somewhere new, faster where the hill is steep.
void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    // Every three seconds or so, a spark is dropped somewhere new.
    if (hashU(id * 7u + u.frame * 7919u) < u.timeDelta / 3.0) {
        q.position = scatter(id * 5u + u.frame * 7919u);
    }
    q.position += curl(q.position, p) * p.speed * 0.3 * u.timeDelta;
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    return float3(q.position / float2(aspect, 1.0) * 0.5 + 0.5, 0.003);
}

float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float edge = smoothstep(1.0, 0.0, length(at));
    float bright = edge * 0.3 * (1.0 - p.trails) / 0.1;
    return float4(float3(1.0, 0.59, 0.20) * bright, 1.0);
}

// Episode 10's trails: each step keeps most of the last one's image, so every spark draws its path.
float4 fade(float2 uv, Trails previous,
            constant Uniforms& u, constant Params& p) {
    return previous.at(uv) * p.trails;
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float3 c = trails.at(uv).rgb;
    return float4(c / (1.0 + c), 1.0);
}

// Stage 3 · A current that changes
// The landscape now drifts through time, with a second, finer layer of noise on top, so the loops open up and the
// sparks travel. Each spark lives a few seconds and is born again somewhere else, fading in and out.
//
// @particles 90000

// @param speed 1.0 0.0 3.0 How fast the current carries the sparks
// @param scale 1.6 0.6 4.0 How small the eddies are
// @param evolve 0.12 0.0 0.5 How fast the current itself changes
// @param life 6.0 2.0 12.0 How many seconds a spark lives
// @param trails 0.96 0.8 0.985 How much of the last image is left at each step

// A spark now has a life, from 1 down to 0, and a random number of its own.
struct Particle {
    float2 position;
    float life;
    float seed;
};

// Episode 5's integer hash: a whole number in, a random number from 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);
}

// Episode 8's corner hash, keeping all its bits: a whole number for each whole-numbered point in 3D.
uint hashCorner(float3 p) {
    uint3 v = uint3(int3(floor(p))) * 1664525u + 1013904223u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    v ^= v >> 16u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    return v.x;
}

// A corner's slope: a random direction in 3D from its hash (a height, and an angle
// around), dotted with the way from the corner to the point.
float cornerSlope(uint h, float3 d) {
    float z = float(h & 0xffffu) / 32767.5 - 1.0;
    float a = float(h >> 16u) * (6.2831853 / 65536.0);
    float s = sqrt(max(1.0 - z * z, 0.0));
    return 1.5 * dot(float3(s * cos(a), s * sin(a), z), d);
}

// One corner's part of the noise at the point: its slope, from the corner to there.
float cornerPart(float3 i, float3 f, float3 c) {
    return cornerSlope(hashCorner(i + c), f - c);
}

// Gradient noise: a random slope at each corner instead of a random height, so the
// hills don't line up with the grid. Roughly -1 to 1.
float gnoise(float3 p) {
    float3 i = floor(p), f = p - i;
    float3 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);
    float n000 = cornerPart(i, f, float3(0, 0, 0));
    float n100 = cornerPart(i, f, float3(1, 0, 0));
    float n010 = cornerPart(i, f, float3(0, 1, 0));
    float n110 = cornerPart(i, f, float3(1, 1, 0));
    float n001 = cornerPart(i, f, float3(0, 0, 1));
    float n101 = cornerPart(i, f, float3(1, 0, 1));
    float n011 = cornerPart(i, f, float3(0, 1, 1));
    float n111 = cornerPart(i, f, float3(1, 1, 1));
    float lo = mix(mix(n000, n100, w.x), mix(n010, n110, w.x), w.y);
    float hi = mix(mix(n001, n101, w.x), mix(n011, n111, w.x), w.y);
    return mix(lo, hi, w.z);
}

// Two octaves of gradient noise, the second turned and finer. The third coordinate
// is time, so the hills slowly change shape.
float potential(float2 x, float t, constant Params& p) {
    float2 q = x * p.scale;
    float2 r = float2(0.8 * q.x - 0.6 * q.y, 0.6 * q.x + 0.8 * q.y) * 2.1 + 7.3;
    float z = t * p.evolve;
    float n = gnoise(float3(q, z)) + 0.45 * gnoise(float3(r, z * 1.4 + 3.1));
    return n / p.scale;
}

// The curl, as before, at a moment in time.
float2 curl(float2 x, float t, constant Params& p) {
    float2 a = float2(0.002, 0.0), b = float2(0.0, 0.002);
    float dx = potential(x + a, t, p) - potential(x - a, t, p);
    float dy = potential(x + b, t, p) - potential(x - b, t, p);
    return float2(dy, -dx) / 0.004;
}

// Anywhere on the picture, or a little past its edges.
float2 anywhere(uint k, float aspect) {
    float2 r = float2(hashU(k), hashU(k + 1u)) * 2.0 - 1.0;
    return r * float2(aspect + 0.05, 1.05);
}

// Each spark starts partway through its life, so they don't all end together.
Particle spawn(uint id, constant Params& p) {
    Particle q;
    q.position = anywhere(id * 4u, 16.0 / 9.0);
    q.life = hashU(id * 4u + 2u);
    q.seed = hashU(id * 4u + 3u);
    return q;
}

// One step along the current, measured halfway through the step (the midpoint method), so a spark keeps to its
// stream through tight bends. Then it burns down; when it burns out, or drifts off the picture, it is born again.
void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    float dt = u.timeDelta;
    float aspect = u.resolution.x / u.resolution.y;
    float2 v = curl(q.position, u.time, p) * p.speed * 0.3;
    float2 mid = q.position + v * dt * 0.5;
    q.position += curl(mid, u.time, p) * p.speed * 0.3 * dt;

    q.life -= dt / (p.life * (0.5 + q.seed));
    bool gone = abs(q.position.x) > aspect + 0.1 || abs(q.position.y) > 1.1;
    if (q.life <= 0.0 || gone) {
        q.position = anywhere(id * 3u + u.frame * 7919u, aspect);
        q.life = 1.0;
    }
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    return float3(q.position / float2(aspect, 1.0) * 0.5 + 0.5, 0.003);
}

// A spark fades in when it is born and out as it ends, so none pops in or out.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float alive = smoothstep(0.0, 0.2, q.life) * smoothstep(1.0, 0.9, q.life);
    float edge = smoothstep(1.0, 0.0, length(at));
    float bright = edge * alive * 0.25 * (1.0 - p.trails) / 0.04;
    return float4(float3(1.0, 0.59, 0.20) * bright, 1.0);
}

float4 fade(float2 uv, Trails previous,
            constant Uniforms& u, constant Params& p) {
    return previous.at(uv) * p.trails;
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float3 c = trails.at(uv).rgb;
    return float4(c / (1.0 + c), 1.0);
}

// The kiln palette from episode 1.
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.78, t));
    c = mix(c, c4, smoothstep(0.78, 1.00, t));
    return c;
}

// The picture, from the light the sparks left.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    // How much light piled up here, rolled off so it never clips, then through the
    // kiln palette: a stream's faint edges are ember, its thick middle gold.
    float light = trails.at(uv).r;
    float level = 1.0 - exp(-light * 1.4);
    float2 d = uv - 0.5;
    float3 ground = float3(0.028, 0.017, 0.019) * (1.0 - dot(d, d) * 1.6);
    float t = level * (0.7 + 0.8 * p.heat) + (p.heat - 0.5) * 0.2;
    float3 c = kiln(t) * smoothstep(0.0, 0.08, level);
    return float4(ground + c, 1.0);
}

// Ray Kiln · Episode 14 · Curl-Noise Flow
//
// @episode 14 Curl-Noise Flow
// @length 60
// @variation 30
// @variations 1,3,4,5,6,8,9,11
// @still 20.0
// @teaches Compute particles following a flow field, with trails
// @category Compute
// @tags compute shaders, particles, curl noise, gradient noise, flow fields, trails, procedural sound
// @final Curl-noise flow | A slow whirl around the middle, eddies of curl noise on top, and every spark coloured by how fast it goes.
//
// Ninety thousand sparks drift on a current that never stops changing. The current is made from noise, but
// not in the obvious way: the sparks don't run downhill on the noise, they run along its contour lines, sideways to
// the slope. That one turn (the curl) makes a flow that never piles sparks up or leaves holes, so the picture stays
// evenly full of streams that swirl, part and join.
// 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
//
// @particles 90000

// @param speed 1.0 0.0 3.0 How fast the current carries the sparks
// @param scale 1.6 0.6 4.0 How small the eddies are: a few wide swirls or many tight ones
// @param evolve 0.12 0.0 0.5 How fast the current itself changes
// @param swirl 0.5 -1.5 1.5 A slow whirl around the middle, one way or the other
// @param life 6.0 2.0 12.0 How many seconds a spark lives before it is born again somewhere else
// @param trails 0.96 0.8 0.985 How much of the last image is left at each step: short sparks or long streaks
// @param heat 0.5 0.0 1.0 Slides the palette from deep ember to white-hot
// @param volume 0.8 0.0 1.0 Loudness of the sound

constant float TAU = 6.2831853;

// What one spark is: where it is (in a space one unit tall from the middle, as wide as the picture) and where it was
// a step ago, how much of its life is left (1 to 0), a random number of its own, and how fast it went at its last
// step, for its brightness.
struct Particle {
    float2 position;
    // Where it was a step ago.
    float2 last;
    float life;
    float seed;
    // How fast it went at its last step.
    float speed;
};

// Episode 5's integer hash: a whole number in, a random number from 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);
}

// Episode 8's corner hash, keeping all its bits: a whole number for each whole-numbered point in 3D.
uint hashCorner(float3 p) {
    uint3 v = uint3(int3(floor(p))) * 1664525u + 1013904223u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    v ^= v >> 16u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    return v.x;
}

// A corner's slope: a random direction in 3D from its hash (a height, and an angle
// around), dotted with the way from the corner to the point.
float cornerSlope(uint h, float3 d) {
    float z = float(h & 0xffffu) / 32767.5 - 1.0;
    float a = float(h >> 16u) * (6.2831853 / 65536.0);
    float s = sqrt(max(1.0 - z * z, 0.0));
    return 1.5 * dot(float3(s * cos(a), s * sin(a), z), d);
}

// One corner's part of the noise at the point: its slope, from the corner to there.
float cornerPart(float3 i, float3 f, float3 c) {
    return cornerSlope(hashCorner(i + c), f - c);
}

// Gradient noise: where value noise gives each corner a height, this gives each corner a slope, so the noise is 0
// at every corner and its slope there is random. Value noise is flat at its corners, and its contour lines would all
// cross the grid square on; gradient noise hides its grid far better. Roughly -1 to 1. z is time, so the field can change.
float gnoise(float3 p) {
    float3 i = floor(p), f = p - i;
    float3 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);
    float n000 = cornerPart(i, f, float3(0, 0, 0));
    float n100 = cornerPart(i, f, float3(1, 0, 0));
    float n010 = cornerPart(i, f, float3(0, 1, 0));
    float n110 = cornerPart(i, f, float3(1, 1, 0));
    float n001 = cornerPart(i, f, float3(0, 0, 1));
    float n101 = cornerPart(i, f, float3(1, 0, 1));
    float n011 = cornerPart(i, f, float3(0, 1, 1));
    float n111 = cornerPart(i, f, float3(1, 1, 1));
    float lo = mix(mix(n000, n100, w.x), mix(n010, n110, w.x), w.y);
    float hi = mix(mix(n001, n101, w.x), mix(n011, n111, w.x), w.y);
    return mix(lo, hi, w.z);
}

// The landscape the flow comes from: two octaves of gradient noise, the second
// turned and finer, drifting through time. Divided by scale, so smaller eddies
// don't also make a faster current.
float potential(float2 x, float t, constant Params& p) {
    float2 q = x * p.scale;
    float2 r = float2(0.8 * q.x - 0.6 * q.y, 0.6 * q.x + 0.8 * q.y) * 2.1 + 7.3;
    float z = t * p.evolve;
    float n = gnoise(float3(q, z)) + 0.45 * gnoise(float3(r, z * 1.4 + 3.1));
    // Plus a smooth hill in the middle. Its contour lines are circles, so its curl
    // is a slow whirl around the middle, under the eddies; swirl sets which way.
    return n / p.scale + p.swirl * 1.2 * exp(-dot(x, x) / 0.6);
}

// The curl: the slope of the landscape, turned a quarter. Measured with four samples a small step apart. The water
// of a flow made this way never gathers anywhere or runs dry: what flows into any small patch flows out of it.
float2 curl(float2 x, float t, constant Params& p) {
    float2 a = float2(0.002, 0.0), b = float2(0.0, 0.002);
    float dx = potential(x + a, t, p) - potential(x - a, t, p);
    float dy = potential(x + b, t, p) - potential(x - b, t, p);
    return float2(dy, -dx) / 0.004;
}

// Anywhere on the picture, a little past its edges, so sparks drift in rather than appear at the border.
float2 anywhere(uint k, float aspect) {
    float2 r = float2(hashU(k), hashU(k + 1u)) * 2.0 - 1.0;
    return r * float2(aspect + 0.05, 1.05);
}

// Each spark starts somewhere on the picture, partway through its life, so they don't all end together.
Particle spawn(uint id, constant Params& p) {
    Particle q;
    q.position = anywhere(id * 4u, 16.0 / 9.0);
    q.last = q.position;
    q.life = hashU(id * 4u + 2u);
    q.seed = hashU(id * 4u + 3u);
    q.speed = 0.0;
    return q;
}

// One step: the spark goes where the current at its place carries it. No velocity is kept: a spark has no weight,
// it is wherever the water takes it. Two half steps (the midpoint method) keep it on its stream through tight bends.
void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    float dt = u.timeDelta;
    float aspect = u.resolution.x / u.resolution.y;
    // Where it was, before it moves, so place can draw it anywhere along the step.
    q.last = q.position;
    float2 v = curl(q.position, u.time, p) * p.speed * 0.3;
    float2 mid = q.position + v * dt * 0.5;
    v = curl(mid, u.time, p) * p.speed * 0.3;
    q.position += v * dt;
    // How fast it went, for its brightness.
    q.speed = length(v);

    // Burning down; a spark that burns out, or drifts off the picture, is born again somewhere else.
    q.life -= dt / (p.life * (0.5 + q.seed));
    bool gone = abs(q.position.x) > aspect + 0.1 || abs(q.position.y) > 1.1;
    if (q.life <= 0.0 || gone) {
        uint k = id * 3u + u.frame * 7919u;
        q.position = anywhere(k, aspect);
        q.last = q.position;
        q.life = 1.0;
        q.speed = 0.0;
    }
}

// Where to draw a spark (the picture's uv) and how big.
float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    // A fast spark jumps a few pixels a step, and its dots would line up like
    // stitches. Drawn at a random point along the step it took, its trail fills in.
    float along = hashU(uint(q.seed * 16777216.0) + u.frame * 7919u);
    float2 spot = mix(q.last, q.position, along);
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, 0.003);
}

// The kiln palette from episode 1.
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.78, t));
    c = mix(c, c4, smoothstep(0.78, 1.00, t));
    return c;
}

// A spark's light. It fades in when it is born and out as it ends, so no spark
// pops, and longer trails dim each dot by as much. A narrower picture holds the
// same sparks closer together, so each dot is dimmer by as much as it is narrower.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    // Brighter the faster the current carries it: the fast streams glow and the
    // slow pools smoulder. Only light is kept here; shade gives it its colour.
    float fast = q.speed / (0.3 * max(p.speed, 0.05));
    float glow = 0.15 + smoothstep(0.1, 1.6, fast);
    float alive = smoothstep(0.0, 0.2, q.life) * smoothstep(1.0, 0.9, q.life);
    float edge = smoothstep(1.0, 0.0, length(at));
    float keep = (1.0 - p.trails) / 0.04;
    float narrow = min(u.resolution.x / u.resolution.y / (16.0 / 9.0), 1.0);
    return float4(float3(edge * alive * glow * 0.16 * keep * narrow), 1.0);
}

// What is left of the last step's image: most of it, so each spark draws its stream behind it.
float4 fade(float2 uv, Trails previous,
            constant Uniforms& u, constant Params& p) {
    return previous.at(uv) * p.trails;
}

// The picture, from the light the sparks left.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    // How much light piled up here, rolled off so it never clips, then through the
    // kiln palette: a stream's faint edges are ember, its thick middle gold.
    float light = trails.at(uv).r;
    float level = 1.0 - exp(-light * 1.4);
    float2 d = uv - 0.5;
    float3 ground = float3(0.028, 0.017, 0.019) * (1.0 - dot(d, d) * 1.6);
    float t = level * (0.7 + 0.8 * p.heat) + (p.heat - 0.5) * 0.2;
    float3 c = kiln(t) * smoothstep(0.0, 0.08, level);
    return float4(ground + c, 1.0);
}

// A smooth random value from 0 to 1 that changes `rate` times a second: two hashed values eased into each other.
float noise1(float t, float rate, uint salt) {
    float x = t * rate;
    float i = floor(x);
    float f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    uint n = uint(int(i) + 1000000);
    return mix(hashU(n * 31u + salt), hashU((n + 1u) * 31u + salt), f);
}

// The sound listens to the current itself, which is a function of place and time alone: three fixed points on the
// picture, each with a voice of its own (D, A and F, a D minor chord) that swells as the water there runs faster; and
// wind, a rumble of smooth noise that follows how fast the current runs on average.
float2 sound(float t, constant Params& p) {
    float2 spots[3] = {
        float2(-0.9, 0.35), float2(0.15, -0.45), float2(0.95, 0.3)
    };
    float notes[3] = { 146.83, 220.0, 349.23 };
    float2 v = float2(0.0);
    float total = 0.0;
    for (int i = 0; i < 3; i++) {
        float run = length(curl(spots[i], t, p)) * 0.3 * p.speed;
        total += run;
        float swell = 0.2 + 0.8 * smoothstep(0.02, 0.5, run);
        float x = notes[i] * t;
        float tone = sin(TAU * x) + 0.3 * sin(2.0 * TAU * x + 0.7)
                   + 0.12 * sin(3.0 * TAU * x + 1.3);
        float sway = 0.8 + 0.2 * sin(TAU * 0.23 * t + float(i) * 2.0);
        tone *= 0.09 * swell * sway;
        float pan = 0.5 + 0.4 * spots[i].x / 1.2;
        v += float2(1.0 - pan, pan) * tone;
    }
    float wind = (noise1(t, 260.0, 3u) - 0.5)
               + 0.5 * (noise1(t, 900.0, 9u) - 0.5);
    float gust = 0.3 + 0.7 * smoothstep(0.0, 1.2, total / 3.0);
    gust *= 0.6 + 0.4 * noise1(t, 0.7, 17u);
    float lean = noise1(t, 0.3, 23u);
    v += float2(1.0 - lean * 0.6, 0.4 + lean * 0.6) * wind * 0.22 * gust;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 3.6) / 1.5 * fadeIn * p.volume;
}

// Ray Kiln · Episode 14 · Curl-Noise Flow
//
// @episode 14 Curl-Noise Flow
// @length 60
// @variation 30
// @variations 1,3,4,5,6,8,9,11
// @still 20.0
// @teaches Compute particles following a flow field, with trails
// @category Compute
// @tags compute shaders, particles, curl noise, gradient noise, flow fields, trails, procedural sound
// @final Curl-noise flow | A slow whirl around the middle, eddies of curl noise on top, and every spark coloured by how fast it goes.
//
// Ninety thousand sparks drift on a current that never stops changing. The current is made from noise, but
// not in the obvious way: the sparks don't run downhill on the noise, they run along its contour lines, sideways to
// the slope. That one turn (the curl) makes a flow that never piles sparks up or leaves holes, so the picture stays
// evenly full of streams that swirl, part and join.
// 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
//
// @particles 90000

// @param speed 1.0 0.0 3.0 How fast the current carries the sparks
// @param scale 1.6 0.6 4.0 How small the eddies are: a few wide swirls or many tight ones
// @param evolve 0.12 0.0 0.5 How fast the current itself changes
// @param swirl 0.5 -1.5 1.5 A slow whirl around the middle, one way or the other
// @param life 6.0 2.0 12.0 How many seconds a spark lives before it is born again somewhere else
// @param trails 0.96 0.8 0.985 How much of the last image is left at each step: short sparks or long streaks
// @param heat 0.5 0.0 1.0 Slides the palette from deep ember to white-hot
// @param volume 0.8 0.0 1.0 Loudness of the sound

constant float TAU = 6.2831853;

// What one spark is: where it is (in a space one unit tall from the middle, as wide as the picture) and where it was
// a step ago, how much of its life is left (1 to 0), a random number of its own, and how fast it went at its last
// step, for its brightness.
struct Particle {
    float2 position;
    // Where it was a step ago.
    float2 last;
    float life;
    float seed;
    // How fast it went at its last step.
    float speed;
};

// Episode 5's integer hash: a whole number in, a random number from 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);
}

// Episode 8's corner hash, keeping all its bits: a whole number for each whole-numbered point in 3D.
uint hashCorner(float3 p) {
    uint3 v = uint3(int3(floor(p))) * 1664525u + 1013904223u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    v ^= v >> 16u;
    v.x += v.y * v.z; v.y += v.z * v.x; v.z += v.x * v.y;
    return v.x;
}

// A corner's slope: a random direction in 3D from its hash (a height, and an angle
// around), dotted with the way from the corner to the point.
float cornerSlope(uint h, float3 d) {
    float z = float(h & 0xffffu) / 32767.5 - 1.0;
    float a = float(h >> 16u) * (6.2831853 / 65536.0);
    float s = sqrt(max(1.0 - z * z, 0.0));
    return 1.5 * dot(float3(s * cos(a), s * sin(a), z), d);
}

// One corner's part of the noise at the point: its slope, from the corner to there.
float cornerPart(float3 i, float3 f, float3 c) {
    return cornerSlope(hashCorner(i + c), f - c);
}

// Gradient noise: where value noise gives each corner a height, this gives each corner a slope, so the noise is 0
// at every corner and its slope there is random. Value noise is flat at its corners, and its contour lines would all
// cross the grid square on; gradient noise hides its grid far better. Roughly -1 to 1. z is time, so the field can change.
float gnoise(float3 p) {
    float3 i = floor(p), f = p - i;
    float3 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);
    float n000 = cornerPart(i, f, float3(0, 0, 0));
    float n100 = cornerPart(i, f, float3(1, 0, 0));
    float n010 = cornerPart(i, f, float3(0, 1, 0));
    float n110 = cornerPart(i, f, float3(1, 1, 0));
    float n001 = cornerPart(i, f, float3(0, 0, 1));
    float n101 = cornerPart(i, f, float3(1, 0, 1));
    float n011 = cornerPart(i, f, float3(0, 1, 1));
    float n111 = cornerPart(i, f, float3(1, 1, 1));
    float lo = mix(mix(n000, n100, w.x), mix(n010, n110, w.x), w.y);
    float hi = mix(mix(n001, n101, w.x), mix(n011, n111, w.x), w.y);
    return mix(lo, hi, w.z);
}

// The landscape the flow comes from: two octaves of gradient noise, the second
// turned and finer, drifting through time. Divided by scale, so smaller eddies
// don't also make a faster current.
float potential(float2 x, float t, constant Params& p) {
    float2 q = x * p.scale;
    float2 r = float2(0.8 * q.x - 0.6 * q.y, 0.6 * q.x + 0.8 * q.y) * 2.1 + 7.3;
    float z = t * p.evolve;
    float n = gnoise(float3(q, z)) + 0.45 * gnoise(float3(r, z * 1.4 + 3.1));
    // Plus a smooth hill in the middle. Its contour lines are circles, so its curl
    // is a slow whirl around the middle, under the eddies; swirl sets which way.
    return n / p.scale + p.swirl * 1.2 * exp(-dot(x, x) / 0.6);
}

// The curl: the slope of the landscape, turned a quarter. Measured with four samples a small step apart. The water
// of a flow made this way never gathers anywhere or runs dry: what flows into any small patch flows out of it.
float2 curl(float2 x, float t, constant Params& p) {
    float2 a = float2(0.002, 0.0), b = float2(0.0, 0.002);
    float dx = potential(x + a, t, p) - potential(x - a, t, p);
    float dy = potential(x + b, t, p) - potential(x - b, t, p);
    return float2(dy, -dx) / 0.004;
}

// Anywhere on the picture, a little past its edges, so sparks drift in rather than appear at the border.
float2 anywhere(uint k, float aspect) {
    float2 r = float2(hashU(k), hashU(k + 1u)) * 2.0 - 1.0;
    return r * float2(aspect + 0.05, 1.05);
}

// Each spark starts somewhere on the picture, partway through its life, so they don't all end together.
Particle spawn(uint id, constant Params& p) {
    Particle q;
    q.position = anywhere(id * 4u, 16.0 / 9.0);
    q.last = q.position;
    q.life = hashU(id * 4u + 2u);
    q.seed = hashU(id * 4u + 3u);
    q.speed = 0.0;
    return q;
}

// One step: the spark goes where the current at its place carries it. No velocity is kept: a spark has no weight,
// it is wherever the water takes it. Two half steps (the midpoint method) keep it on its stream through tight bends.
void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    float dt = u.timeDelta;
    float aspect = u.resolution.x / u.resolution.y;
    // Where it was, before it moves, so place can draw it anywhere along the step.
    q.last = q.position;
    float2 v = curl(q.position, u.time, p) * p.speed * 0.3;
    float2 mid = q.position + v * dt * 0.5;
    v = curl(mid, u.time, p) * p.speed * 0.3;
    q.position += v * dt;
    // How fast it went, for its brightness.
    q.speed = length(v);

    // Burning down; a spark that burns out, or drifts off the picture, is born again somewhere else.
    q.life -= dt / (p.life * (0.5 + q.seed));
    bool gone = abs(q.position.x) > aspect + 0.1 || abs(q.position.y) > 1.1;
    if (q.life <= 0.0 || gone) {
        uint k = id * 3u + u.frame * 7919u;
        q.position = anywhere(k, aspect);
        q.last = q.position;
        q.life = 1.0;
        q.speed = 0.0;
    }
}

// Where to draw a spark (the picture's uv) and how big.
float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    // A fast spark jumps a few pixels a step, and its dots would line up like
    // stitches. Drawn at a random point along the step it took, its trail fills in.
    float along = hashU(uint(q.seed * 16777216.0) + u.frame * 7919u);
    float2 spot = mix(q.last, q.position, along);
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, 0.003);
}

// The kiln palette from episode 1.
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.78, t));
    c = mix(c, c4, smoothstep(0.78, 1.00, t));
    return c;
}

// A spark's light. It fades in when it is born and out as it ends, so no spark
// pops, and longer trails dim each dot by as much. A narrower picture holds the
// same sparks closer together, so each dot is dimmer by as much as it is narrower.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    // Brighter the faster the current carries it: the fast streams glow and the
    // slow pools smoulder. Only light is kept here; shade gives it its colour.
    float fast = q.speed / (0.3 * max(p.speed, 0.05));
    float glow = 0.15 + smoothstep(0.1, 1.6, fast);
    float alive = smoothstep(0.0, 0.2, q.life) * smoothstep(1.0, 0.9, q.life);
    float edge = smoothstep(1.0, 0.0, length(at));
    float keep = (1.0 - p.trails) / 0.04;
    float narrow = min(u.resolution.x / u.resolution.y / (16.0 / 9.0), 1.0);
    return float4(float3(edge * alive * glow * 0.16 * keep * narrow), 1.0);
}

// What is left of the last step's image: most of it, so each spark draws its stream behind it.
float4 fade(float2 uv, Trails previous,
            constant Uniforms& u, constant Params& p) {
    return previous.at(uv) * p.trails;
}

// The picture, from the light the sparks left.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    // How much light piled up here, rolled off so it never clips, then through the
    // kiln palette: a stream's faint edges are ember, its thick middle gold.
    float light = trails.at(uv).r;
    float level = 1.0 - exp(-light * 1.4);
    float2 d = uv - 0.5;
    float3 ground = float3(0.028, 0.017, 0.019) * (1.0 - dot(d, d) * 1.6);
    float t = level * (0.7 + 0.8 * p.heat) + (p.heat - 0.5) * 0.2;
    float3 c = kiln(t) * smoothstep(0.0, 0.08, level);
    return float4(ground + c, 1.0);
}

// A smooth random value from 0 to 1 that changes `rate` times a second: two hashed values eased into each other.
float noise1(float t, float rate, uint salt) {
    float x = t * rate;
    float i = floor(x);
    float f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    uint n = uint(int(i) + 1000000);
    return mix(hashU(n * 31u + salt), hashU((n + 1u) * 31u + salt), f);
}

// The sound listens to the current itself, which is a function of place and time alone: three fixed points on the
// picture, each with a voice of its own (D, A and F, a D minor chord) that swells as the water there runs faster; and
// wind, a rumble of smooth noise that follows how fast the current runs on average.
float2 sound(float t, constant Params& p) {
    float2 spots[3] = {
        float2(-0.9, 0.35), float2(0.15, -0.45), float2(0.95, 0.3)
    };
    float notes[3] = { 146.83, 220.0, 349.23 };
    float2 v = float2(0.0);
    float total = 0.0;
    for (int i = 0; i < 3; i++) {
        float run = length(curl(spots[i], t, p)) * 0.3 * p.speed;
        total += run;
        float swell = 0.2 + 0.8 * smoothstep(0.02, 0.5, run);
        float x = notes[i] * t;
        float tone = sin(TAU * x) + 0.3 * sin(2.0 * TAU * x + 0.7)
                   + 0.12 * sin(3.0 * TAU * x + 1.3);
        float sway = 0.8 + 0.2 * sin(TAU * 0.23 * t + float(i) * 2.0);
        tone *= 0.09 * swell * sway;
        float pan = 0.5 + 0.4 * spots[i].x / 1.2;
        v += float2(1.0 - pan, pan) * tone;
    }
    float wind = (noise1(t, 260.0, 3u) - 0.5)
               + 0.5 * (noise1(t, 900.0, 9u) - 0.5);
    float gust = 0.3 + 0.7 * smoothstep(0.0, 1.2, total / 3.0);
    gust *= 0.6 + 0.4 * noise1(t, 0.7, 17u);
    float lean = noise1(t, 0.3, 23u);
    v += float2(1.0 - lean * 0.6, 0.4 + lean * 0.6) * wind * 0.22 * gust;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 3.6) / 1.5 * fadeIn * p.volume;
}

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.0, 1.6, 0.12, 0.5, 6.0, 0.96, 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)