Shader

Spiral Galaxy

A spiral galaxy seen at a slant, from 160,000 stars on orbits. Each orbit is an oval turned a little further than the one inside it, so the ovals crowd into two arms that turn slowly while the stars drive through them. Young stars and ember gas light up where an arm squeezes them, a lane of dust runs inside each arm, and the core is a ball of old stars under a warm haze.

· Particles & Simulation · Compute shader

Teaches Stars on orbits, density waves, dust lanes

Builds on Flowing Color Field (value noise), Floating 3D Scene (rotation), Cracked Clay (integer hash, sound noise), Particle Swarm, Curl-Noise Flow, Murmuration

Watch

The idea

Particle Swarm, Curl-Noise Flow and Murmuration moved their particles a step at a time: each step started from where the particle was and added a little. A star in a galaxy needs none of that. Its orbit is known, so where it is at any moment is a function of its number and the clock. update here throws the last step away and works every star out again from scratch, which is why the galaxy can be scrubbed to any moment and still look right.

The arms are the puzzle. Stars near the middle go round faster than stars far out, so any arm made of the same stars would wind up tighter and tighter, and be gone in a few turns. Real arms are density waves: places where stars crowd for a while, as cars crowd in a traffic jam, then drive on. This galaxy makes its jam the way Kalnajs drew it in 1973. Every orbit is an oval, and each oval is turned a little further than the one inside it. Where neighbouring ovals swing close together, the stars on them crowd, and those crowded places line up along two spirals.

It is built in four steps, and each step is a shader you can run on its own.

Stage 1 · Stars on orbits

Every star gets a distance a from the middle, most of them near it (an exponential spread, 0.3 × −log(1 − h)) and none past 1.3, and an angle to start from. orbitRate is the rotation curve: how many radians a second a star goes round at that distance. A star's speed along its orbit rises through the core (CORE, 0.12) and then stays at 0.11 a second, so an outer star, with further to go, takes longer: once in about 13 seconds at 0.2, about 57 seconds at 1 and 75 at the edge, 1.3 (at a speed of 1). The disk is thicker in the middle, and viewOf turns it a little and leans it back by tilt degrees.

Two spokes of stars are marked in ember: the ones that started within a few degrees of a line through the middle. Run it and watch them. Their inner ends race ahead, the spokes bend into spirals, and the spirals wind tighter until they are a set of rings. That is the winding problem, and why the arms can't be made of the same stars.

stage-1-orbits.metal
// Stage 1 · Stars on orbits
// Every star goes round the middle on a circle of its own, the inner ones faster. Two spokes of stars are marked in
// ember: they wind up into tighter and tighter spirals, so arms can't be made of the same stars for long.
//
// @particles 160000

// @param speed 1.0 0.0 3.0 How fast the stars go round: the outer ones once a minute at 1
// @param tilt 58.0 0.0 80.0 How far the galaxy leans away from face-on, in degrees
// @param trails 0.6 0.0 0.9 How much of the last image is left at each step

constant float PI = 3.1415927;
constant float TAU = 6.2831853;
constant float CORE = 0.12;
constant float ZOOM = 1.1;

// A star: where it is (x and y in the disk, height out of it), and if it is marked.
struct Particle {
    float2 position;
    float height;
    float marked;
};

// Cracked Clay'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);
}

// Kaleidoscope's rotation: a point turned by angle a about the middle.
float2 rot(float2 p, float a) {
    float c = cos(a), s = sin(a);
    return float2(c * p.x - s * p.y, s * p.x + c * p.y);
}

// How many radians a second a star goes round at distance a: its speed rises through the core, then stays the same,
// so an outer star, with further to go, takes longer.
float orbitRate(float a, constant Params& p) {
    float v = 0.11 * a / sqrt(a * a + CORE * CORE);
    return p.speed * v / max(a, 0.01);
}

// Star number id at time t: nothing is carried between steps, a star is wherever its orbit has it at that time.
// Most stars are near the middle; the disk is thicker there.
Particle star(uint id, float t, constant Params& p) {
    float h1 = hashU(id * 5u + 1u), h2 = hashU(id * 5u + 2u);
    float h3 = hashU(id * 5u + 3u);
    Particle q;
    float a = min(0.012 + 0.3 * -log(1.0 - h1 * 0.995), 1.3);
    float m = h2 * TAU + orbitRate(a, p) * t;
    q.position = float2(cos(m), sin(m)) * a;
    q.height = (h3 - 0.5) * 0.05 * (1.0 + 2.0 * exp(-a / 0.15));
    // Two spokes: the stars that started within a few degrees of a line through the middle.
    q.marked = step(abs(fract(h2 * 2.0 + 0.5) - 0.5), 0.012);
    return q;
}

Particle spawn(uint id, constant Params& p) {
    return star(id, 0.0, p);
}

void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    q = star(id, u.time, p);
}

// Where a point of the galaxy is on the picture, one unit the half height: turned a little, then leaned by tilt.
float2 viewOf(float2 position, float height, constant Params& p) {
    float lean = p.tilt * PI / 180.0;
    float2 r = rot(position, 0.4);
    return float2(r.x, r.y * cos(lean) + height * sin(lean));
}

// Its radius a little more than the half height, or short of the half width in a picture taller than wide.
float zoomOf(constant Uniforms& u) {
    return min(ZOOM, 0.95 * u.resolution.x / u.resolution.y);
}

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

// Gold light, ember for the marked stars. The dots add up; longer trails keep more of each step, so each gives less.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float3 c = mix(float3(1.0, 0.8, 0.55) * 0.75, float3(1.0, 0.3, 0.08) * 3.0,
                   q.marked);
    float edge = smoothstep(1.0, 0.0, length(at));
    return float4(c * edge * (1.0 - p.trails) * 0.5, 0.0);
}

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

// The light the stars left, rolled off so it never clips, on a dark ground.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float3 col = 1.0 - exp(-trails.at(uv).rgb * 1.6);
    return float4(col + float3(0.035, 0.022, 0.024), 1.0);
}

Stage 2 · Ovals that crowd

Now each orbit is an oval: ovalness squashes its short axis by up to oval (0.3), round in the core and rounder again at the edge, and orbitTurn turns each oval by twist × a, so each one is turned a little further than the one inside it. The whole set of ovals also turns slowly, 0.035 radians a second, a turn in about three minutes; each star goes round its own oval at its own rate, less that turning.

Where is the jam? Going round the oval, the next oval out, turned a little further, swings in towards this one an eighth of a turn back from each end of the long axis. That is where armWeight puts it: −sin(2m), eased, where m is how far round its oval the star is. Stars glow ember as they pass through, and two arms appear, made of stars that are only visiting. Set oval to 0 and the arms vanish: round orbits never crowd.

stage-2-ovals.metal23 new lines since stage 1
// Stage 2 · Ovals that crowd
// Each orbit becomes an oval, and each oval is turned a little further than the one inside it. The ovals crowd
// along two spiral lines: the arms. Stars drive into them and out again, glowing ember while they are squeezed.
//
// @particles 160000

// @param speed 1.0 0.0 3.0 How fast the stars go round: the outer ones once a minute at 1
// @param twist 3.2 0.5 6.0 How far each orbit is turned from the one inside it: loose arms or tightly wound
// @param oval 0.3 0.0 0.45 How oval the orbits are: round orbits make no arms
// @param tilt 58.0 0.0 80.0 How far the galaxy leans away from face-on, in degrees
// @param trails 0.6 0.0 0.9 How much of the last image is left at each step

constant float PI = 3.1415927;
constant float TAU = 6.2831853;
constant float CORE = 0.12;
// How fast the arms turn, in radians a second: far slower than the stars.
constant float PATTERN = 0.035;
constant float ZOOM = 1.1;

// A star: where it is (x and y in the disk, height out of it), and how deep in an arm.
struct Particle {
    float2 position;
    float height;
    float glow;
};

// Cracked Clay'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);
}

// Kaleidoscope's rotation: a point turned by angle a about the middle.
float2 rot(float2 p, float a) {
    float c = cos(a), s = sin(a);
    return float2(c * p.x - s * p.y, s * p.x + c * p.y);
}

// How many radians a second a star goes round at distance a: its speed rises through the core, then stays the same,
// so an outer star, with further to go, takes longer.
float orbitRate(float a, constant Params& p) {
    float v = 0.11 * a / sqrt(a * a + CORE * CORE);
    return p.speed * v / max(a, 0.01);
}

// How oval the orbit of size a is: round in the core, rounder at the edge.
float ovalness(float a, constant Params& p) {
    float edge = 1.0 - 0.6 * smoothstep(0.6, 1.3, a);
    return p.oval * smoothstep(0.03, 0.25, a) * edge;
}

// Each oval turned a little further than the one inside it: the arms.
float orbitTurn(float a, float t, constant Params& p) {
    return p.twist * a + PATTERN * p.speed * t;
}

// The point at angle m round the oval of size a, in the galaxy's frame at time t.
float2 onOrbit(float a, float m, float t, constant Params& p) {
    float2 q = float2(a * cos(m), a * (1.0 - ovalness(a, p)) * sin(m));
    return rot(q, orbitTurn(a, t, p));
}

// How deep in an arm a star at angle m round its oval is, 0 to 1: the ovals crowd an eighth of a turn back from each
// end of the long axis, where the next oval out swings in. None in the core, fading out at the edge.
float armWeight(float m, float a, constant Params& p) {
    float squeeze = smoothstep(0.2, 1.0, -sin(2.0 * m));
    float disk = smoothstep(0.08, 0.3, a) * (1.0 - smoothstep(0.75, 1.15, a));
    return squeeze * disk * smoothstep(0.0, 0.1, p.oval);
}

// Star number id at time t: nothing is carried between steps, a star is wherever its orbit has it at that time.
// Most stars are near the middle; the disk is thicker there.
Particle star(uint id, float t, constant Params& p) {
    float h1 = hashU(id * 5u + 1u), h2 = hashU(id * 5u + 2u);
    float h3 = hashU(id * 5u + 3u);
    Particle q;
    float a = min(0.012 + 0.3 * -log(1.0 - h1 * 0.995), 1.3);
    // Round its oval at the rate of its size, less the turning of the ovals.
    float m = h2 * TAU + (orbitRate(a, p) - PATTERN * p.speed) * t;
    q.position = onOrbit(a, m, t, p);
    q.height = (h3 - 0.5) * 0.05 * (1.0 + 2.0 * exp(-a / 0.15));
    q.glow = armWeight(m, a, p);
    return q;
}

Particle spawn(uint id, constant Params& p) {
    return star(id, 0.0, p);
}

void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    q = star(id, u.time, p);
}

// Where a point of the galaxy is on the picture, one unit the half height: turned a little, then leaned by tilt.
float2 viewOf(float2 position, float height, constant Params& p) {
    float lean = p.tilt * PI / 180.0;
    float2 r = rot(position, 0.4);
    return float2(r.x, r.y * cos(lean) + height * sin(lean));
}

// Its radius a little more than the half height, or short of the half width in a picture taller than wide.
float zoomOf(constant Uniforms& u) {
    return min(ZOOM, 0.95 * u.resolution.x / u.resolution.y);
}

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

// Gold light, ember as the arm squeezes a star. The dots add up; longer trails keep more of each step.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float3 c = mix(float3(1.0, 0.8, 0.55) * 0.75, float3(1.0, 0.3, 0.08) * 1.6,
                   q.glow);
    float edge = smoothstep(1.0, 0.0, length(at));
    return float4(c * edge * (1.0 - p.trails) * 0.5, 0.0);
}

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

// The light the stars left, rolled off so it never clips, on a dark ground.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float3 col = 1.0 - exp(-trails.at(uv).rgb * 1.6);
    return float4(col + float3(0.035, 0.022, 0.024), 1.0);
}
6 lines from stage 1 that are gone
    float marked;
⋯
    float m = h2 * TAU + orbitRate(a, p) * t;
    q.position = float2(cos(m), sin(m)) * a;
⋯
    q.marked = step(abs(fract(h2 * 2.0 + 0.5) - 0.5), 0.012);
⋯
    float3 c = mix(float3(1.0, 0.8, 0.55) * 0.75, float3(1.0, 0.3, 0.08) * 3.0,
                   q.marked);

Stage 3 · Young stars and the core

Squeezed gas makes new stars, so stars now come in kinds. About one in sixteen (at young 0.6) is a young star, thin in the disk, nearly white and bright, and lit only while the arm holds it (armWeight squared). A few more are clouds of glowing gas, drawn as soft ember dots about ten times wider than a star, lit only in the thick of the squeeze (armWeight cubed). Most are old disk stars, everywhere, gold and dim, fading far out so the disk has no edge.

One star in six is in the core: small round orbits, each on a plane of its own leaned by up to 1.3 radians and turned to any direction, so the core is a ball rather than part of the disk. shade adds a haze for the stars too many to see one by one: a warm glow round the middle, round on the picture as a ball is, and a faint one over the disk, worked out back on the disk itself.

stage-3-young-stars.metal48 new lines since stage 2
// Stage 3 · Young stars and the core
// The arm's squeeze makes new stars: young stars, nearly white, and clouds of glowing gas light up only while the arm
// holds them. In the middle, a ball of old stars on round orbits tilted every way, under a warm haze.
//
// @particles 160000

// @param speed 1.0 0.0 3.0 How fast the stars go round: the outer ones once a minute at 1
// @param twist 3.2 0.5 6.0 How far each orbit is turned from the one inside it: loose arms or tightly wound
// @param oval 0.3 0.0 0.45 How oval the orbits are: round orbits make no arms
// @param young 0.6 0.0 1.0 How many young stars and glowing clouds light up in the arms
// @param core 0.5 0.0 1.0 How bright the old stars of the core are
// @param tilt 58.0 0.0 80.0 How far the galaxy leans away from face-on, in degrees
// @param trails 0.6 0.0 0.9 How much of the last image is left at each step

constant float PI = 3.1415927;
constant float TAU = 6.2831853;
// The core's size in the rotation curve: inside it the speed rises from nothing, outside it stays the same.
constant float CORE = 0.12;
// How fast the arms themselves turn, in radians a second at speed 1: far slower than the stars.
constant float PATTERN = 0.035;
// How much of the picture's half height the galaxy's radius fills.
constant float ZOOM = 1.1;

// What one star is, where it is now: its place in the galaxy's own frame (x and y in the disk, height out of it),
// what kind of star it is (0 an old disk star, 1 a star of the core, 2 a young star, 3 a cloud of glowing gas), and
// how lit it is now.
struct Particle {
    float2 position;
    float height;
    float kind;
    float glow;
};

// Cracked Clay'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);
}

// Kaleidoscope's rotation: a point turned by angle a about the middle.
float2 rot(float2 p, float a) {
    float c = cos(a), s = sin(a);
    return float2(c * p.x - s * p.y, s * p.x + c * p.y);
}

// The rotation curve: how many radians a second a star goes round at distance a. Its speed along the orbit rises
// through the core and then stays the same, so an outer star, with further to go, takes longer to go round.
float orbitRate(float a, constant Params& p) {
    float v = 0.11 * a / sqrt(a * a + CORE * CORE);
    return p.speed * v / max(a, 0.01);
}

// How oval the orbit of size a is: round in the core, most oval in the disk, rounder again at the edge.
float ovalness(float a, constant Params& p) {
    float edge = 1.0 - 0.6 * smoothstep(0.6, 1.3, a);
    return p.oval * smoothstep(0.03, 0.25, a) * edge;
}

// How far the oval of size a is turned: twist turns each oval a little further than the one inside it, and the
// whole set turns slowly with time. This one line is where the arms come from.
float orbitTurn(float a, float t, constant Params& p) {
    return p.twist * a + PATTERN * p.speed * t;
}

// The point at angle m round the oval of size a, in the galaxy's frame at time t.
float2 onOrbit(float a, float m, float t, constant Params& p) {
    float2 q = float2(a * cos(m), a * (1.0 - ovalness(a, p)) * sin(m));
    return rot(q, orbitTurn(a, t, p));
}

// How deep in an arm a star at angle m round its oval is, 0 to 1. The ovals crowd where the next one out, turned a
// little further, swings in towards this one: an eighth of a turn back from each end of the long axis. None in the
// core, where the orbits are round, and fading out at the edge of the disk.
float armWeight(float m, float a, constant Params& p) {
    float squeeze = smoothstep(0.2, 1.0, -sin(2.0 * m));
    float disk = smoothstep(0.08, 0.3, a) * (1.0 - smoothstep(0.75, 1.15, a));
    return squeeze * disk * smoothstep(0.0, 0.1, p.oval);
}

// Star number id at time t. Nothing is carried from one step to the next: a star on an orbit is wherever its orbit
// has it at that time, so every step works out each star from its number and the clock.
Particle star(uint id, float t, constant Params& p) {
    float pick = hashU(id * 5u);
    float h1 = hashU(id * 5u + 1u), h2 = hashU(id * 5u + 2u);
    float h3 = hashU(id * 5u + 3u), h4 = hashU(id * 5u + 4u);
    Particle q;
    if (pick < 0.16) {
        // The core: small round orbits, each on a plane of its own, leaned by lean and turned by node, so the core
        // is a ball, not a disk.
        float a = 0.004 + 0.09 * -log(1.0 - h1 * 0.97);
        float m = h2 * TAU + orbitRate(a, p) * t;
        float lean = (h3 * 2.0 - 1.0) * 1.3;
        float2 flat = float2(cos(m), sin(m) * cos(lean)) * a;
        q.position = rot(flat, h4 * TAU);
        q.height = sin(m) * sin(lean) * a * 0.8;
        q.kind = 1.0;
        q.glow = 1.0;
        return q;
    }
    // The disk: most stars near the middle, fewer further out. Each goes round its oval at the rate of its size,
    // less the turning of the ovals themselves.
    float a = min(0.012 + 0.3 * -log(1.0 - h1 * 0.995), 1.3);
    float m = h2 * TAU + (orbitRate(a, p) - PATTERN * p.speed) * t;
    q.position = onOrbit(a, m, t, p);
    float arm = armWeight(m, a, p);
    if (pick < 0.16 + 0.1 * p.young) {
        // Young stars: thin in the disk, and lit only while the arm squeezes them.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 2.0;
        q.glow = arm * arm;
    } else if (pick < 0.27 && h4 < p.young) {
        // Glowing gas, fewer, lit only in the thick of the squeeze.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 3.0;
        q.glow = arm * arm * arm;
    } else {
        // Old stars: everywhere in the disk, which is thicker in the middle, and dimmer far out so it has no edge.
        q.height = (h3 - 0.5) * 0.05 * (1.0 + 2.0 * exp(-a / 0.15));
        q.kind = 0.0;
        q.glow = 1.0 - 0.8 * smoothstep(0.6, 1.3, a);
    }
    return q;
}

Particle spawn(uint id, constant Params& p) {
    return star(id, 0.0, p);
}

void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    q = star(id, u.time, p);
}

// Where a point of the galaxy is on the picture, one unit the half height: turned a little, then leaned by tilt.
float2 viewOf(float2 position, float height, constant Params& p) {
    float lean = p.tilt * PI / 180.0;
    float2 r = rot(position, 0.4);
    return float2(r.x, r.y * cos(lean) + height * sin(lean));
}

// How big the galaxy is drawn: its radius a little more than the half height (it leans, so it is less tall than
// wide), or a little short of the half width in a picture taller than wide.
float zoomOf(constant Uniforms& u) {
    return min(ZOOM, 0.95 * u.resolution.x / u.resolution.y);
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 spot = viewOf(q.position, q.height, p) * zoomOf(u);
    float size = q.kind == 3.0 ? 0.03 : q.kind == 2.0 ? 0.0034 : 0.0024;
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, size);
}

// A star's light, by kind: the old stars gold and dim, the core warmer, the young stars nearly white and bright,
// the gas ember and soft. Only light: the dots add up, and shade rolls the sum off so it never clips. Longer
// trails keep more of each step, so each dot gives less.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float d = length(at);
    float3 c = float3(1.0, 0.8, 0.55) * 0.75;
    float edge = smoothstep(1.0, 0.0, d);
    if (q.kind == 1.0) c = float3(1.0, 0.62, 0.3) * (0.11 + 0.64 * p.core);
    if (q.kind == 2.0) c = float3(1.0, 0.95, 0.88) * 2.6;
    if (q.kind == 3.0) {
        c = float3(1.0, 0.28, 0.1) * 0.1;
        edge = exp(-d * d * 4.0) * step(d, 1.0);
    }
    return float4(c * edge * q.glow * (1.0 - p.trails) * 0.5, 0.0);
}

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

// The picture: the starlight, and a smooth haze for the stars too many to see one by one.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float aspect = u.resolution.x / u.resolution.y;
    float3 light = trails.at(uv).rgb;
    // Back from the picture to the disk: undo the zoom, the lean and the turn.
    float lean = p.tilt * PI / 180.0;
    float2 s = (uv - 0.5) * 2.0 * float2(aspect, 1.0) / zoomOf(u);
    float2 x = rot(float2(s.x, s.y / max(cos(lean), 0.15)), -0.4);
    // The core's glow, round on the picture since the core is a ball, and the disk's, flat with it.
    float bulge = exp(-length(s) / 0.07) * (0.2 + 0.8 * p.core) * 0.9;
    light += float3(1.0, 0.6, 0.28) * bulge;
    light += float3(1.0, 0.7, 0.45) * exp(-length(x) / 0.3) * 0.06;
    float3 col = 1.0 - exp(-light * 1.6);
    return float4(col + float3(0.035, 0.022, 0.024), 1.0);
}
8 lines from stage 2 that are gone
    float h3 = hashU(id * 5u + 3u);
⋯
    q.glow = armWeight(m, a, p);
⋯
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, 0.0024);
⋯
    float3 c = mix(float3(1.0, 0.8, 0.55) * 0.75, float3(1.0, 0.3, 0.08) * 1.6,
                   q.glow);
    float edge = smoothstep(1.0, 0.0, length(at));
    return float4(c * edge * (1.0 - p.trails) * 0.5, 0.0);
⋯
    float3 col = 1.0 - exp(-trails.at(uv).rgb * 1.6);

Stage 4 · Spiral Galaxy

The last step adds what lies in front and behind.

  • Dust lanes. shade undoes the view to find the point of the disk under each pixel, and dustAt measures where it is against the turning ovals. A lane runs a little inside each arm (pow(max(−sin(2θ − 1.5), 0), 12)), torn into wisps by two layers of value noise that turn with the arms, from 0.08 out to about 0.95. The starlight there is dimmed by exp(−dust × 2.2); dust sets how dark.
  • Far stars. The stars of our own sky, behind the galaxy: one in about four hundred cells of 2.5 pixels holds one.
  • The pointer turns the galaxy about its axis as it moves across, and leans it more or less as it moves up and down.
shade(), from the finished shader
// The picture: the starlight, dimmed where dust lies in front, a smooth haze for the stars too many to see one by
// one, and the stars of our own sky far behind.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float aspect = u.resolution.x / u.resolution.y;
    float3 light = trails.at(uv).rgb;
    // Back from the picture to the disk: undo the zoom, the lean and the turn, for the point of the disk here.
    float2 v = viewAngles(u, p);
    float2 s = (uv - 0.5) * 2.0 * float2(aspect, 1.0) / zoomOf(u);
    float2 x = rot(float2(s.x, s.y / max(cos(v.y), 0.15)), -v.x);
    // The haze: the core's glow, round on the picture since the core is a ball, and the disk's, flat with it.
    float bulge = exp(-length(s) / 0.07) * (0.2 + 0.8 * p.core) * 0.9;
    light += float3(1.0, 0.6, 0.28) * bulge;
    light += float3(1.0, 0.7, 0.45) * exp(-length(x) / 0.3) * 0.06;
    light *= exp(-dustAt(x, u.time, p) * p.dust * 2.2);
    float3 col = 1.0 - exp(-light * 1.6);
    // Far stars: one in a few hundred little cells.
    float2 cell = floor(uv * u.resolution / 2.5);
    uint n = uint(cell.x + 8192.0) * 73856093u
           ^ uint(cell.y + 8192.0) * 19349663u;
    col += float3(1.0, 0.92, 0.8) * smoothstep(0.9975, 1.0, hashU(n)) * 0.45;
    col += float3(0.035, 0.022, 0.024);
    return float4(col, 1.0);
}
spiral-galaxy.metal, the whole file85 new lines since stage 3
// Ray Kiln · Spiral Galaxy
//
// @episode Spiral Galaxy
// @length 60
// @variation 30
// @variations 1,2,3,4,7,10,11,12
// @still 24.0
// @teaches Stars on orbits, density waves, dust lanes
// @category Compute
// @tags compute shaders, particles, orbits, density waves, spiral galaxy, value noise, dust, stars, procedural sound
// @final Spiral galaxy | Stars on turning ovals crowd into two arms, young stars light where they pass through, and dust darkens the arms' inner edges.
//
// A spiral galaxy, seen at a slant. Its arms look like solid things, but no star stays in them: every star is on an
// orbit of its own, the inner ones going round faster than the outer. Each orbit is an oval, and each oval is turned
// a little further than the one inside it, so the ovals crowd together along two spiral lines. That crowding is the
// arm: a traffic jam that stars drive into and out of, while the jam itself turns slowly. Young stars and glowing gas
// light up where the jam squeezes them, and a lane of dust runs along each arm's inner edge.
// 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 160000

// @param speed 1.0 0.0 3.0 How fast the stars go round: the outer ones once a minute at 1
// @param twist 3.2 0.5 6.0 How far each orbit is turned from the one inside it: loose arms or tightly wound
// @param oval 0.3 0.0 0.45 How oval the orbits are: round orbits make no arms
// @param young 0.6 0.0 1.0 How many young stars and glowing clouds light up in the arms
// @param dust 0.7 0.0 1.5 How dark the dust lanes are
// @param core 0.5 0.0 1.0 How bright the old stars of the core are
// @param tilt 58.0 0.0 80.0 How far the galaxy leans away from face-on, in degrees
// @param trails 0.6 0.0 0.9 How much of the last image is left at each step: crisp stars or a soft glow
// @param volume 0.8 0.0 1.0 Loudness of the sound

constant float PI = 3.1415927;
constant float TAU = 6.2831853;
// The core's size in the rotation curve: inside it the speed rises from nothing, outside it stays the same.
constant float CORE = 0.12;
// How fast the arms themselves turn, in radians a second at speed 1: far slower than the stars.
constant float PATTERN = 0.035;
// How much of the picture's half height the galaxy's radius fills.
constant float ZOOM = 1.1;

// What one star is, where it is now: its place in the galaxy's own frame (x and y in the disk, height out of it),
// what kind of star it is (0 an old disk star, 1 a star of the core, 2 a young star, 3 a cloud of glowing gas), and
// how lit it is now.
struct Particle {
    float2 position;
    float height;
    float kind;
    float glow;
};

// Cracked Clay'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);
}

// Kaleidoscope's rotation: a point turned by angle a about the middle.
float2 rot(float2 p, float a) {
    float c = cos(a), s = sin(a);
    return float2(c * p.x - s * p.y, s * p.x + c * p.y);
}

// The rotation curve: how many radians a second a star goes round at distance a. Its speed along the orbit rises
// through the core and then stays the same, so an outer star, with further to go, takes longer to go round.
float orbitRate(float a, constant Params& p) {
    float v = 0.11 * a / sqrt(a * a + CORE * CORE);
    return p.speed * v / max(a, 0.01);
}

// How oval the orbit of size a is: round in the core, most oval in the disk, rounder again at the edge.
float ovalness(float a, constant Params& p) {
    float edge = 1.0 - 0.6 * smoothstep(0.6, 1.3, a);
    return p.oval * smoothstep(0.03, 0.25, a) * edge;
}

// How far the oval of size a is turned: twist turns each oval a little further than the one inside it, and the
// whole set turns slowly with time. This one line is where the arms come from.
float orbitTurn(float a, float t, constant Params& p) {
    return p.twist * a + PATTERN * p.speed * t;
}

// The point at angle m round the oval of size a, in the galaxy's frame at time t.
float2 onOrbit(float a, float m, float t, constant Params& p) {
    float2 q = float2(a * cos(m), a * (1.0 - ovalness(a, p)) * sin(m));
    return rot(q, orbitTurn(a, t, p));
}

// How deep in an arm a star at angle m round its oval is, 0 to 1. The ovals crowd where the next one out, turned a
// little further, swings in towards this one: an eighth of a turn back from each end of the long axis. None in the
// core, where the orbits are round, and fading out at the edge of the disk.
float armWeight(float m, float a, constant Params& p) {
    float squeeze = smoothstep(0.2, 1.0, -sin(2.0 * m));
    float disk = smoothstep(0.08, 0.3, a) * (1.0 - smoothstep(0.75, 1.15, a));
    return squeeze * disk * smoothstep(0.0, 0.1, p.oval);
}

// Star number id at time t. Nothing is carried from one step to the next: a star on an orbit is wherever its orbit
// has it at that time, so every step works out each star from its number and the clock.
Particle star(uint id, float t, constant Params& p) {
    float pick = hashU(id * 5u);
    float h1 = hashU(id * 5u + 1u), h2 = hashU(id * 5u + 2u);
    float h3 = hashU(id * 5u + 3u), h4 = hashU(id * 5u + 4u);
    Particle q;
    if (pick < 0.16) {
        // The core: small round orbits, each on a plane of its own, leaned by lean and turned by node, so the core
        // is a ball, not a disk.
        float a = 0.004 + 0.09 * -log(1.0 - h1 * 0.97);
        float m = h2 * TAU + orbitRate(a, p) * t;
        float lean = (h3 * 2.0 - 1.0) * 1.3;
        float2 flat = float2(cos(m), sin(m) * cos(lean)) * a;
        q.position = rot(flat, h4 * TAU);
        q.height = sin(m) * sin(lean) * a * 0.8;
        q.kind = 1.0;
        q.glow = 1.0;
        return q;
    }
    // The disk: most stars near the middle, fewer further out. Each goes round its oval at the rate of its size,
    // less the turning of the ovals themselves.
    float a = min(0.012 + 0.3 * -log(1.0 - h1 * 0.995), 1.3);
    float m = h2 * TAU + (orbitRate(a, p) - PATTERN * p.speed) * t;
    q.position = onOrbit(a, m, t, p);
    float arm = armWeight(m, a, p);
    if (pick < 0.16 + 0.1 * p.young) {
        // Young stars: thin in the disk, and lit only while the arm squeezes them.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 2.0;
        q.glow = arm * arm;
    } else if (pick < 0.27 && h4 < p.young) {
        // Glowing gas, fewer, lit only in the thick of the squeeze.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 3.0;
        q.glow = arm * arm * arm;
    } else {
        // Old stars: everywhere in the disk, which is thicker in the middle, and dimmer far out so it has no edge.
        q.height = (h3 - 0.5) * 0.05 * (1.0 + 2.0 * exp(-a / 0.15));
        q.kind = 0.0;
        q.glow = 1.0 - 0.8 * smoothstep(0.6, 1.3, a);
    }
    return q;
}

Particle spawn(uint id, constant Params& p) {
    return star(id, 0.0, p);
}

void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    q = star(id, u.time, p);
}

// How the galaxy is seen: the pointer across turns it about its axis, the pointer up and down leans it more or less.
// Gives the turn and the lean.
float2 viewAngles(constant Uniforms& u, constant Params& p) {
    float2 m = u.mouse / max(u.resolution, float2(1.0)) - 0.5;
    float lean = clamp(p.tilt * PI / 180.0 + m.y * 1.2, 0.0, 1.45);
    return float2(0.4 - m.x * 2.0, lean);
}

// Where a point of the galaxy is on the picture, one unit the half height.
float2 viewOf(float2 position, float height,
              constant Uniforms& u, constant Params& p) {
    float2 v = viewAngles(u, p);
    float2 r = rot(position, v.x);
    return float2(r.x, r.y * cos(v.y) + height * sin(v.y));
}

// How big the galaxy is drawn: its radius a little more than the half height (it leans, so it is less tall than
// wide), or a little short of the half width in a picture taller than wide.
float zoomOf(constant Uniforms& u) {
    return min(ZOOM, 0.95 * u.resolution.x / u.resolution.y);
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 spot = viewOf(q.position, q.height, u, p) * zoomOf(u);
    float size = q.kind == 3.0 ? 0.03 : q.kind == 2.0 ? 0.0034 : 0.0024;
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, size);
}

// A star's light, by kind: the old stars gold and dim, the core warmer, the young stars nearly white and bright,
// the gas ember and soft. Only light: the dots add up, and shade rolls the sum off so it never clips. Longer
// trails keep more of each step, so each dot gives less.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float d = length(at);
    float3 c = float3(1.0, 0.8, 0.55) * 0.75;
    float edge = smoothstep(1.0, 0.0, d);
    if (q.kind == 1.0) c = float3(1.0, 0.62, 0.3) * (0.11 + 0.64 * p.core);
    if (q.kind == 2.0) c = float3(1.0, 0.95, 0.88) * 2.6;
    if (q.kind == 3.0) {
        c = float3(1.0, 0.28, 0.1) * 0.1;
        edge = exp(-d * d * 4.0) * step(d, 1.0);
    }
    return float4(c * edge * q.glow * (1.0 - p.trails) * 0.5, 0.0);
}

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

// Value noise from Flowing Color Field: hashed corners eased into each other.
float vhash(float2 c) {
    uint x = uint(int(c.x) + 4096), y = uint(int(c.y) + 4096);
    return hashU(x * 7919u + y * 104729u);
}

float vnoise(float2 x) {
    float2 i = floor(x), f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    return mix(mix(vhash(i), vhash(i + float2(1, 0)), f.x),
               mix(vhash(i + float2(0, 1)), vhash(i + float2(1, 1)), f.x), f.y);
}

// The dust, at a point of the disk: a lane a little inside each arm, where the gas crowds first, torn into wisps by
// two layers of noise that turn with the arms.
float dustAt(float2 x, float t, constant Params& p) {
    float r = length(x);
    float angle = atan2(x.y, x.x) - orbitTurn(r, t, p);
    float lane = pow(max(-sin(2.0 * angle - 1.5), 0.0), 12.0);
    float2 w = rot(x, -PATTERN * p.speed * t);
    float wisps = vnoise(w * 9.0) * 0.65 + vnoise(w * 23.0 + 5.0) * 0.35;
    float reach = smoothstep(0.08, 0.2, r) * (1.0 - smoothstep(0.6, 0.95, r));
    float torn = smoothstep(0.25, 0.75, wisps);
    return lane * torn * reach * smoothstep(0.0, 0.1, p.oval);
}

// The picture: the starlight, dimmed where dust lies in front, a smooth haze for the stars too many to see one by
// one, and the stars of our own sky far behind.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float aspect = u.resolution.x / u.resolution.y;
    float3 light = trails.at(uv).rgb;
    // Back from the picture to the disk: undo the zoom, the lean and the turn, for the point of the disk here.
    float2 v = viewAngles(u, p);
    float2 s = (uv - 0.5) * 2.0 * float2(aspect, 1.0) / zoomOf(u);
    float2 x = rot(float2(s.x, s.y / max(cos(v.y), 0.15)), -v.x);
    // The haze: the core's glow, round on the picture since the core is a ball, and the disk's, flat with it.
    float bulge = exp(-length(s) / 0.07) * (0.2 + 0.8 * p.core) * 0.9;
    light += float3(1.0, 0.6, 0.28) * bulge;
    light += float3(1.0, 0.7, 0.45) * exp(-length(x) / 0.3) * 0.06;
    light *= exp(-dustAt(x, u.time, p) * p.dust * 2.2);
    float3 col = 1.0 - exp(-light * 1.6);
    // Far stars: one in a few hundred little cells.
    float2 cell = floor(uv * u.resolution / 2.5);
    uint n = uint(cell.x + 8192.0) * 73856093u
           ^ uint(cell.y + 8192.0) * 19349663u;
    col += float3(1.0, 0.92, 0.8) * smoothstep(0.9975, 1.0, hashU(n)) * 0.45;
    col += float3(0.035, 0.022, 0.024);
    return float4(col, 1.0);
}

// A smooth random value from 0 to 1 that changes `rate` times a second.
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 follows two stars of the disk, one near the middle on the left and one further out on the right. Each
// has a drone (D and A) that swells as its star goes into an arm, and while it is squeezed, young stars light round
// it as chimes: short notes of a D pentatonic, a little above, at random moments. Under them, the core's low hum.
float2 sound(float t, constant Params& p) {
    float radii[2] = { 0.38, 0.72 };
    float drones[2] = { 146.83, 220.0 };
    float pent[5] = { 587.33, 659.25, 739.99, 880.0, 987.77 };
    float2 v = float2(0.0);
    for (int i = 0; i < 2; i++) {
        float a = radii[i];
        float rate = orbitRate(a, p) - PATTERN * p.speed;
        float start = 1.3 + float(i) * 2.1;
        float arm = armWeight(start + rate * t, a, p);
        float pan = i == 0 ? 0.3 : 0.7;
        float x = drones[i] * t;
        float tone = sin(TAU * x) + 0.35 * sin(2.0 * TAU * x + 0.4)
                   + 0.15 * sin(3.0 * TAU * x + 1.1);
        tone *= 0.07 * (0.25 + 0.75 * arm);
        v += float2(1.0 - pan, pan) * tone;
        // Chimes, eight slots a second: a slot rings when its hash is under how squeezed the star was as it began,
        // so a note never starts halfway through.
        for (int k = 0; k < 6; k++) {
            float slot = floor(t * 8.0) - float(k);
            uint n = uint(int(slot) + 100000) * 2u + uint(i);
            float then = armWeight(start + rate * slot / 8.0, a, p);
            float chance = then * then * 0.35 * (0.3 + p.young);
            float ring = step(hashU(n * 13u + 5u), chance);
            float age = t - slot / 8.0;
            float f = pent[uint(hashU(n * 13u + 7u) * 4.99)];
            f *= i == 0 ? 1.0 : 1.5;
            float env = ring * exp(-age * 5.0) * smoothstep(0.0, 0.004, age);
            float wide = (hashU(n * 13u + 9u) - 0.5) * 0.5;
            float2 side = float2(0.5 - wide, 0.5 + wide);
            v += side * sin(TAU * f * age) * env * 0.05;
        }
    }
    float hum = sin(TAU * 73.42 * t) * (0.6 + 0.4 * noise1(t, 0.4, 3u))
              + 0.4 * sin(TAU * 36.71 * t);
    v += hum * 0.06 * (0.3 + 0.7 * p.core);
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 5.0) / 1.3 * fadeIn * p.volume;
}
8 lines from stage 3 that are gone
float2 viewOf(float2 position, float height, constant Params& p) {
    float lean = p.tilt * PI / 180.0;
    float2 r = rot(position, 0.4);
    return float2(r.x, r.y * cos(lean) + height * sin(lean));
⋯
    float2 spot = viewOf(q.position, q.height, p) * zoomOf(u);
⋯
    float lean = p.tilt * PI / 180.0;
⋯
    float2 x = rot(float2(s.x, s.y / max(cos(lean), 0.15)), -0.4);
⋯
    return float4(col + float3(0.035, 0.022, 0.024), 1.0);

The sound

The sound rides with two stars of the disk: one at 0.38 from the middle, on the left, and one at 0.72, on the right. Each has a drone, D3 (146.83 Hz) for the inner one and A3 (220 Hz) for the outer, that swells from a quarter to full strength as its star drives into an arm, using the same armWeight as the picture. The inner star reaches an arm about every 13 seconds and the outer about every 27, so the two swells drift in and out of step.

While a star is squeezed, young stars light round it as chimes: eight times a second a slot may ring, more likely the deeper the star is in the arm (and the higher young), each a short note of the D pentatonic (D5, E5, F♯5, A5, B5), a fifth higher for the outer star, dying in about a fifth of a second and spread a little across the ears. Whether a slot rings is decided by how squeezed the star was as the slot began, so no note starts halfway through. Under them hums the core, D2 and D1 (73.42 and 36.71 Hz), louder with core. It all goes through tanh and fades in over the first two seconds. At the defaults the peak is about -8.0 dBFS and the average (RMS) about -19.9 dBFS.

New in Ray Kiln

Particle count live. The particle readout under a particle episode's picture in the playground now has a Particles slider for its // @particles N line, from 1 to 4,000,000 (or as many as one buffer on that GPU can hold), on round numbers. While you drag, the swarm on screen changes at once, with no compile: past the new count, particles are dropped; new ones are made by spawn and join the swarm at the step it is at. Beside it the readout says the steps a second the swarm actually gets and, where the browser can time the GPU, what a step costs, in all and per particle. When you let go, the count is written into the code's // @particles line, so a Share link, a save and an export all have it, and the swarm carries on rather than starting again.

Open Spiral Galaxy in the playground and drag Particles down to 5,000: the galaxy thins to a faint scatter round its core. Since each step works every star out from the clock, the new stars are in their places from the next step. Then drag it up to 4,000,000: the dots add up to a bright, smooth disk. The step itself still costs under a nanosecond a star; it is drawing four million dots that takes the time, about 20 ms a step on the Mac it was checked on, so the swarm gets 40 or so steps a second instead of 60 and the readout says how far it falls behind the clock. Bring it back down and it catches up.

Try this

  • Set oval to 0: the orbits are round, nothing crowds, and the arms, their young stars and the dust are gone.
  • Set twist to 6 for tightly wound arms, or 0.5 for a pair of loose ones like a bar.
  • Set tilt to 80: the galaxy is nearly edge-on, and the dust lane crosses in front of the core.
  • In armWeight, change -sin(2.0 * m) to sin(2.0 * m): the young stars light up between the arms instead, out of step with the dust.
  • Set speed to 0: everything stops, and the picture is one moment of the galaxy.
  • Press Inspect and click dots in an arm until you find one of kind 2, a young star: its glow rises as the arm squeezes it and falls to 0 as it drives out.

Parameters

ParameterDefaultRangeWhat it does
speed1.00.0 to 3.0How fast the stars go round: the outer ones once a minute at 1
twist3.20.5 to 6.0How far each orbit is turned from the one inside it: loose arms or tightly wound
oval0.30.0 to 0.45How oval the orbits are: round orbits make no arms
young0.60.0 to 1.0How many young stars and glowing clouds light up in the arms
dust0.70.0 to 1.5How dark the dust lanes are
core0.50.0 to 1.0How bright the old stars of the core are
tilt58.00.0 to 80.0How far the galaxy leans away from face-on, in degrees
trails0.60.0 to 0.9How much of the last image is left at each step: crisp stars or a soft glow
volume0.80.0 to 1.0Loudness of the sound

More from this shader

The same shader 8 more ways: each from a different moment, with different settings, one after another. Each chapter below says which settings moved most, then lists them all.

Wallpapers

A still from the shader, rendered at full size for a desktop or a phone screen. Free to use as your wallpaper.

Clips to post

Fifteen seconds of the shader with its sound, sized for posting: square for feeds, vertical for stories and short videos. Free to post, ideally with a link to this page.

The full source

This is the whole episode: the picture, the parameters and the sound.

spiral-galaxy.metal
// Ray Kiln · Spiral Galaxy
//
// @episode Spiral Galaxy
// @length 60
// @variation 30
// @variations 1,2,3,4,7,10,11,12
// @still 24.0
// @teaches Stars on orbits, density waves, dust lanes
// @category Compute
// @tags compute shaders, particles, orbits, density waves, spiral galaxy, value noise, dust, stars, procedural sound
// @final Spiral galaxy | Stars on turning ovals crowd into two arms, young stars light where they pass through, and dust darkens the arms' inner edges.
//
// A spiral galaxy, seen at a slant. Its arms look like solid things, but no star stays in them: every star is on an
// orbit of its own, the inner ones going round faster than the outer. Each orbit is an oval, and each oval is turned
// a little further than the one inside it, so the ovals crowd together along two spiral lines. That crowding is the
// arm: a traffic jam that stars drive into and out of, while the jam itself turns slowly. Young stars and glowing gas
// light up where the jam squeezes them, and a lane of dust runs along each arm's inner edge.
// 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 160000

// @param speed 1.0 0.0 3.0 How fast the stars go round: the outer ones once a minute at 1
// @param twist 3.2 0.5 6.0 How far each orbit is turned from the one inside it: loose arms or tightly wound
// @param oval 0.3 0.0 0.45 How oval the orbits are: round orbits make no arms
// @param young 0.6 0.0 1.0 How many young stars and glowing clouds light up in the arms
// @param dust 0.7 0.0 1.5 How dark the dust lanes are
// @param core 0.5 0.0 1.0 How bright the old stars of the core are
// @param tilt 58.0 0.0 80.0 How far the galaxy leans away from face-on, in degrees
// @param trails 0.6 0.0 0.9 How much of the last image is left at each step: crisp stars or a soft glow
// @param volume 0.8 0.0 1.0 Loudness of the sound

constant float PI = 3.1415927;
constant float TAU = 6.2831853;
// The core's size in the rotation curve: inside it the speed rises from nothing, outside it stays the same.
constant float CORE = 0.12;
// How fast the arms themselves turn, in radians a second at speed 1: far slower than the stars.
constant float PATTERN = 0.035;
// How much of the picture's half height the galaxy's radius fills.
constant float ZOOM = 1.1;

// What one star is, where it is now: its place in the galaxy's own frame (x and y in the disk, height out of it),
// what kind of star it is (0 an old disk star, 1 a star of the core, 2 a young star, 3 a cloud of glowing gas), and
// how lit it is now.
struct Particle {
    float2 position;
    float height;
    float kind;
    float glow;
};

// Cracked Clay'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);
}

// Kaleidoscope's rotation: a point turned by angle a about the middle.
float2 rot(float2 p, float a) {
    float c = cos(a), s = sin(a);
    return float2(c * p.x - s * p.y, s * p.x + c * p.y);
}

// The rotation curve: how many radians a second a star goes round at distance a. Its speed along the orbit rises
// through the core and then stays the same, so an outer star, with further to go, takes longer to go round.
float orbitRate(float a, constant Params& p) {
    float v = 0.11 * a / sqrt(a * a + CORE * CORE);
    return p.speed * v / max(a, 0.01);
}

// How oval the orbit of size a is: round in the core, most oval in the disk, rounder again at the edge.
float ovalness(float a, constant Params& p) {
    float edge = 1.0 - 0.6 * smoothstep(0.6, 1.3, a);
    return p.oval * smoothstep(0.03, 0.25, a) * edge;
}

// How far the oval of size a is turned: twist turns each oval a little further than the one inside it, and the
// whole set turns slowly with time. This one line is where the arms come from.
float orbitTurn(float a, float t, constant Params& p) {
    return p.twist * a + PATTERN * p.speed * t;
}

// The point at angle m round the oval of size a, in the galaxy's frame at time t.
float2 onOrbit(float a, float m, float t, constant Params& p) {
    float2 q = float2(a * cos(m), a * (1.0 - ovalness(a, p)) * sin(m));
    return rot(q, orbitTurn(a, t, p));
}

// How deep in an arm a star at angle m round its oval is, 0 to 1. The ovals crowd where the next one out, turned a
// little further, swings in towards this one: an eighth of a turn back from each end of the long axis. None in the
// core, where the orbits are round, and fading out at the edge of the disk.
float armWeight(float m, float a, constant Params& p) {
    float squeeze = smoothstep(0.2, 1.0, -sin(2.0 * m));
    float disk = smoothstep(0.08, 0.3, a) * (1.0 - smoothstep(0.75, 1.15, a));
    return squeeze * disk * smoothstep(0.0, 0.1, p.oval);
}

// Star number id at time t. Nothing is carried from one step to the next: a star on an orbit is wherever its orbit
// has it at that time, so every step works out each star from its number and the clock.
Particle star(uint id, float t, constant Params& p) {
    float pick = hashU(id * 5u);
    float h1 = hashU(id * 5u + 1u), h2 = hashU(id * 5u + 2u);
    float h3 = hashU(id * 5u + 3u), h4 = hashU(id * 5u + 4u);
    Particle q;
    if (pick < 0.16) {
        // The core: small round orbits, each on a plane of its own, leaned by lean and turned by node, so the core
        // is a ball, not a disk.
        float a = 0.004 + 0.09 * -log(1.0 - h1 * 0.97);
        float m = h2 * TAU + orbitRate(a, p) * t;
        float lean = (h3 * 2.0 - 1.0) * 1.3;
        float2 flat = float2(cos(m), sin(m) * cos(lean)) * a;
        q.position = rot(flat, h4 * TAU);
        q.height = sin(m) * sin(lean) * a * 0.8;
        q.kind = 1.0;
        q.glow = 1.0;
        return q;
    }
    // The disk: most stars near the middle, fewer further out. Each goes round its oval at the rate of its size,
    // less the turning of the ovals themselves.
    float a = min(0.012 + 0.3 * -log(1.0 - h1 * 0.995), 1.3);
    float m = h2 * TAU + (orbitRate(a, p) - PATTERN * p.speed) * t;
    q.position = onOrbit(a, m, t, p);
    float arm = armWeight(m, a, p);
    if (pick < 0.16 + 0.1 * p.young) {
        // Young stars: thin in the disk, and lit only while the arm squeezes them.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 2.0;
        q.glow = arm * arm;
    } else if (pick < 0.27 && h4 < p.young) {
        // Glowing gas, fewer, lit only in the thick of the squeeze.
        q.height = (h3 - 0.5) * 0.01;
        q.kind = 3.0;
        q.glow = arm * arm * arm;
    } else {
        // Old stars: everywhere in the disk, which is thicker in the middle, and dimmer far out so it has no edge.
        q.height = (h3 - 0.5) * 0.05 * (1.0 + 2.0 * exp(-a / 0.15));
        q.kind = 0.0;
        q.glow = 1.0 - 0.8 * smoothstep(0.6, 1.3, a);
    }
    return q;
}

Particle spawn(uint id, constant Params& p) {
    return star(id, 0.0, p);
}

void update(thread Particle& q, uint id,
            constant Uniforms& u, constant Params& p) {
    q = star(id, u.time, p);
}

// How the galaxy is seen: the pointer across turns it about its axis, the pointer up and down leans it more or less.
// Gives the turn and the lean.
float2 viewAngles(constant Uniforms& u, constant Params& p) {
    float2 m = u.mouse / max(u.resolution, float2(1.0)) - 0.5;
    float lean = clamp(p.tilt * PI / 180.0 + m.y * 1.2, 0.0, 1.45);
    return float2(0.4 - m.x * 2.0, lean);
}

// Where a point of the galaxy is on the picture, one unit the half height.
float2 viewOf(float2 position, float height,
              constant Uniforms& u, constant Params& p) {
    float2 v = viewAngles(u, p);
    float2 r = rot(position, v.x);
    return float2(r.x, r.y * cos(v.y) + height * sin(v.y));
}

// How big the galaxy is drawn: its radius a little more than the half height (it leans, so it is less tall than
// wide), or a little short of the half width in a picture taller than wide.
float zoomOf(constant Uniforms& u) {
    return min(ZOOM, 0.95 * u.resolution.x / u.resolution.y);
}

float3 place(Particle q, constant Uniforms& u, constant Params& p) {
    float aspect = u.resolution.x / u.resolution.y;
    float2 spot = viewOf(q.position, q.height, u, p) * zoomOf(u);
    float size = q.kind == 3.0 ? 0.03 : q.kind == 2.0 ? 0.0034 : 0.0024;
    return float3(spot / float2(aspect, 1.0) * 0.5 + 0.5, size);
}

// A star's light, by kind: the old stars gold and dim, the core warmer, the young stars nearly white and bright,
// the gas ember and soft. Only light: the dots add up, and shade rolls the sum off so it never clips. Longer
// trails keep more of each step, so each dot gives less.
float4 ink(Particle q, float2 at, constant Uniforms& u, constant Params& p) {
    float d = length(at);
    float3 c = float3(1.0, 0.8, 0.55) * 0.75;
    float edge = smoothstep(1.0, 0.0, d);
    if (q.kind == 1.0) c = float3(1.0, 0.62, 0.3) * (0.11 + 0.64 * p.core);
    if (q.kind == 2.0) c = float3(1.0, 0.95, 0.88) * 2.6;
    if (q.kind == 3.0) {
        c = float3(1.0, 0.28, 0.1) * 0.1;
        edge = exp(-d * d * 4.0) * step(d, 1.0);
    }
    return float4(c * edge * q.glow * (1.0 - p.trails) * 0.5, 0.0);
}

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

// Value noise from Flowing Color Field: hashed corners eased into each other.
float vhash(float2 c) {
    uint x = uint(int(c.x) + 4096), y = uint(int(c.y) + 4096);
    return hashU(x * 7919u + y * 104729u);
}

float vnoise(float2 x) {
    float2 i = floor(x), f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    return mix(mix(vhash(i), vhash(i + float2(1, 0)), f.x),
               mix(vhash(i + float2(0, 1)), vhash(i + float2(1, 1)), f.x), f.y);
}

// The dust, at a point of the disk: a lane a little inside each arm, where the gas crowds first, torn into wisps by
// two layers of noise that turn with the arms.
float dustAt(float2 x, float t, constant Params& p) {
    float r = length(x);
    float angle = atan2(x.y, x.x) - orbitTurn(r, t, p);
    float lane = pow(max(-sin(2.0 * angle - 1.5), 0.0), 12.0);
    float2 w = rot(x, -PATTERN * p.speed * t);
    float wisps = vnoise(w * 9.0) * 0.65 + vnoise(w * 23.0 + 5.0) * 0.35;
    float reach = smoothstep(0.08, 0.2, r) * (1.0 - smoothstep(0.6, 0.95, r));
    float torn = smoothstep(0.25, 0.75, wisps);
    return lane * torn * reach * smoothstep(0.0, 0.1, p.oval);
}

// The picture: the starlight, dimmed where dust lies in front, a smooth haze for the stars too many to see one by
// one, and the stars of our own sky far behind.
float4 shade(float2 uv, constant Uniforms& u, constant Params& p,
             Trails trails) {
    float aspect = u.resolution.x / u.resolution.y;
    float3 light = trails.at(uv).rgb;
    // Back from the picture to the disk: undo the zoom, the lean and the turn, for the point of the disk here.
    float2 v = viewAngles(u, p);
    float2 s = (uv - 0.5) * 2.0 * float2(aspect, 1.0) / zoomOf(u);
    float2 x = rot(float2(s.x, s.y / max(cos(v.y), 0.15)), -v.x);
    // The haze: the core's glow, round on the picture since the core is a ball, and the disk's, flat with it.
    float bulge = exp(-length(s) / 0.07) * (0.2 + 0.8 * p.core) * 0.9;
    light += float3(1.0, 0.6, 0.28) * bulge;
    light += float3(1.0, 0.7, 0.45) * exp(-length(x) / 0.3) * 0.06;
    light *= exp(-dustAt(x, u.time, p) * p.dust * 2.2);
    float3 col = 1.0 - exp(-light * 1.6);
    // Far stars: one in a few hundred little cells.
    float2 cell = floor(uv * u.resolution / 2.5);
    uint n = uint(cell.x + 8192.0) * 73856093u
           ^ uint(cell.y + 8192.0) * 19349663u;
    col += float3(1.0, 0.92, 0.8) * smoothstep(0.9975, 1.0, hashU(n)) * 0.45;
    col += float3(0.035, 0.022, 0.024);
    return float4(col, 1.0);
}

// A smooth random value from 0 to 1 that changes `rate` times a second.
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 follows two stars of the disk, one near the middle on the left and one further out on the right. Each
// has a drone (D and A) that swells as its star goes into an arm, and while it is squeezed, young stars light round
// it as chimes: short notes of a D pentatonic, a little above, at random moments. Under them, the core's low hum.
float2 sound(float t, constant Params& p) {
    float radii[2] = { 0.38, 0.72 };
    float drones[2] = { 146.83, 220.0 };
    float pent[5] = { 587.33, 659.25, 739.99, 880.0, 987.77 };
    float2 v = float2(0.0);
    for (int i = 0; i < 2; i++) {
        float a = radii[i];
        float rate = orbitRate(a, p) - PATTERN * p.speed;
        float start = 1.3 + float(i) * 2.1;
        float arm = armWeight(start + rate * t, a, p);
        float pan = i == 0 ? 0.3 : 0.7;
        float x = drones[i] * t;
        float tone = sin(TAU * x) + 0.35 * sin(2.0 * TAU * x + 0.4)
                   + 0.15 * sin(3.0 * TAU * x + 1.1);
        tone *= 0.07 * (0.25 + 0.75 * arm);
        v += float2(1.0 - pan, pan) * tone;
        // Chimes, eight slots a second: a slot rings when its hash is under how squeezed the star was as it began,
        // so a note never starts halfway through.
        for (int k = 0; k < 6; k++) {
            float slot = floor(t * 8.0) - float(k);
            uint n = uint(int(slot) + 100000) * 2u + uint(i);
            float then = armWeight(start + rate * slot / 8.0, a, p);
            float chance = then * then * 0.35 * (0.3 + p.young);
            float ring = step(hashU(n * 13u + 5u), chance);
            float age = t - slot / 8.0;
            float f = pent[uint(hashU(n * 13u + 7u) * 4.99)];
            f *= i == 0 ? 1.0 : 1.5;
            float env = ring * exp(-age * 5.0) * smoothstep(0.0, 0.004, age);
            float wide = (hashU(n * 13u + 9u) - 0.5) * 0.5;
            float2 side = float2(0.5 - wide, 0.5 + wide);
            v += side * sin(TAU * f * age) * env * 0.05;
        }
    }
    float hum = sin(TAU * 73.42 * t) * (0.6 + 0.4 * noise1(t, 0.4, 3u))
              + 0.4 * sin(TAU * 36.71 * t);
    v += hum * 0.06 * (0.3 + 0.7 * p.core);
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 5.0) / 1.3 * fadeIn * p.volume;
}

Running it in your own project

Every Ray Kiln shader is the same shape. The standalone download above is the shader with the small wrapper around it that the Ray Kiln host adds, so it compiles with the ordinary Metal compiler and runs in your own app. It defines three entry points: rk_vertex (a full-screen triangle), rk_fragment (calls the shader's shade) and, for shaders with sound, the compute kernel rk_sound.

To draw it, pass the uniforms at fragment buffer 0 and the parameter values at buffer 1, in the order they are declared:

Drawing it, in Swift
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, 3.2, 0.3, 0.6, 0.7, 0.5, 58.0, 0.6, 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)

Time is in seconds, mouse is in pixels from the bottom left, and uv in the shader runs 0 to 1 with (0, 0) at the bottom left. Write the output as sRGB: the shader's numbers go to the screen as they are.

Built on later

More in Particles & Simulation