Shader

Black Hole

A black hole with a disk of hot gas around it, seen from close by. Every pixel's ray is walked outwards and bent towards the hole at each step, so the far side of the disk shows over and under the shadow and the stars behind are drawn into a ring. The gas orbits near the speed of light, so the side coming towards you is beamed white-hot and the side going away falls to ember.

· 3D Scenes · Fragment shader

Teaches Rays bent by gravity, an accretion disk, Doppler colour

Builds on Floating 3D Scene, Cracked Clay (integer hash, sound noise), Kiln Fire (3D fractal noise), Glass Orb, Spinning Planet (3D hash, 3D value noise), Spiral Galaxy (orbit rate)

Watch

The idea

A black hole is the one thing in space you can't draw by finding its surface, because what you see of it is not a surface at all. It is a shadow: the directions from which no light reaches you, because any ray that came from there would have had to pass too close to the hole and fell in. Around the shadow, light that passed near the hole arrives bent, so the disk of gas feeding it, which lies flat in one plane, appears wrapped over and under the shadow, and the stars behind it are smeared into a ring.

So this episode doesn't march to a surface, as Glass Orb and Floating 3D Scene do. It walks every pixel's ray outwards from the camera in small steps and, at each step, bends the ray's direction towards the hole. Everything else, the shadow, the wrapped disk, the Einstein ring, falls out of that one bend. The units are the hole's own: the horizon is at a radius of 1, light can circle the hole at 1.5, and the disk runs from the innermost stable orbit at 3 out to 9. The camera sits 19 radii away.

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

Stage 1 · Straight rays

The walk first, with no gravity. trace starts at the camera with the pixel's ray direction and steps along it, a step of about a twentieth of the distance to the hole, never less than 0.03 and never more than 0.5. If the ray's distance from the centre drops below the horizon's radius, it has fallen in and the pixel is black. If it gets past 28 radii and is heading outwards, it has escaped, and stars gives it the sky: a grid of cells over the sky's directions, a star in about one cell in six, as in Spinning Planet.

The disk is a plane, y = 0. When a step would take the ray from one side of that plane to the other, the crossing point is found by interpolation, and diskGlow gives the light there: nothing inside the innermost orbit or beyond the outer edge, and between them a profile that is brightest just outside the inner edge and falls as the square of the radius, which is roughly how a real disk burns. With straight rays, the picture is what you would expect: a black ball with a flat ring of fire around it, the ring hiding the ball's lower half.

stage-1-straight.metal
// Stage 1 · Straight rays
// A camera looks at the hole from nineteen radii away. Each pixel's ray is
// walked outwards in small steps: if it runs into the horizon it goes black,
// where it crosses the disk's plane it takes the ring's light, and if it
// gets away it sees the stars. No gravity yet: the rays are straight.

// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;

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

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

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    return tint * star * (0.3 + 2.0 * bright);
}

// The disk's glow where a ray crosses its plane at the point q: a flat ring
// from the innermost orbit outwards, brightest just outside its inner edge
// and fading as the square of r.
float diskGlow(float3 q) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return 0.7 * profile * edge * 24.0;
}

// Walks a ray from the camera, straight, and gathers what it sees: the
// disk where it crosses the disk's plane, black if it runs into the hole,
// the stars if it escapes.
float3 trace(float3 ro, float3 rd, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    float3 light = float3(0.0);
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where there is more to see.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q);
            if (glow > 0.0) { light += glow * p.disk; break; }
        }
        pos = next;
    }
    if (escaped) light += stars(vel);
    return light;
}

// The camera: at rest, looking at the hole, tilted up out of the disk's
// plane as far as the slider says.
float3 cameraPos(constant Params& p) {
    float a = 0.35;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// The hum of the hole, for now.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

float2 sound(float t, constant Params& p) {
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    float2 v = hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * fadeIn * p.volume;
}

Stage 2 · Gravity bends the ray

One line of physics. Light in Schwarzschild's geometry feels a pull towards the hole of 1.5 h² / r⁵ times the position, where h is the ray's angular momentum about the hole, cross(pos, vel), fixed when the ray starts and conserved after. Each step adds that pull times the step to the ray's direction, which is then normalised again, so the ray only turns. The inverse fifth power is why nothing much happens far away and everything happens near the hole: a ray passing at 1.5 radii can be turned right round.

The gravity slider scales h². At 0 the rays are stage 1's; at 1 they are Einstein's; at 1.6 the hole bends harder than any real one, and the shadow grows. At 1, rays passing above the shadow are bent down behind it and cross the disk's plane on its far side, so the far side of the disk appears above the shadow; rays passing below are bent up and show it again, smaller, beneath. The flat ring has become the shape every picture of a black hole has. The steps are smaller near the hole, where the bend is strong, so the ring at the shadow's edge, where light has circled the hole, stays sharp.

stage-2-gravity.metal4 new lines since stage 1
// Stage 2 · Gravity bends the ray
// One line of physics. At each step the ray's direction is pulled towards
// the hole by Schwarzschild's bend, which grows as the inverse fifth power
// of the distance, times the ray's angular momentum. Rays passing behind
// the hole are bent over and under it, so the far side of the disk shows
// above and below the shadow, and the stars are drawn into a ring.

// @param gravity 1 0 1.6 How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;

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

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

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    return tint * star * (0.3 + 2.0 * bright);
}

// The disk's glow where a ray crosses its plane at the point q: a flat ring
// from the innermost orbit outwards, brightest just outside its inner edge
// and fading as the square of r.
float diskGlow(float3 q) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return 0.7 * profile * edge * 24.0;
}

// Walks a ray from the camera, bent by the hole, and gathers what it sees:
// the disk where it crosses the disk's plane, black if it falls into the
// hole, the stars if it escapes.
float3 trace(float3 ro, float3 rd, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    // Angular momentum about the hole, which the bending conserves.
    float3 L = cross(pos, vel);
    float h2 = dot(L, L) * p.gravity;
    float3 light = float3(0.0);
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where the bending is strong.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q);
            if (glow > 0.0) { light += glow * p.disk; break; }
        }
        // Gravity: the pull on a ray of light in Schwarzschild's geometry,
        // felt as a bend towards the hole that grows fast as r shrinks.
        vel += -1.5 * h2 * pos / (r2 * r2 * sqrt(r2)) * dt;
        vel = normalize(vel);
        pos = next;
    }
    if (escaped) light += stars(vel);
    return light;
}

// The camera: at rest, looking at the hole, tilted up out of the disk's
// plane as far as the slider says.
float3 cameraPos(constant Params& p) {
    float a = 0.35;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// The hum of the hole, still.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

float2 sound(float t, constant Params& p) {
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    float2 v = hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * fadeIn * p.volume;
}

Stage 3 · The disk

Gas, not a ring. orbitRate gives how fast the gas goes round at each radius, Kepler's law with the inner edge going round once every 10 s at spin 1 and the outer gas five times slower. diskGlow takes each crossing point's angle round the hole, turns it back by the orbit's rate times the time, and reads fractal noise there. The noise is sampled on a cylinder, (cos φ, sin φ, r), so there is no seam where the angle wraps, and it is stretched round the orbit by sampling the angle coarsely and the radius finely, so it reads as streaks. The inner streaks shear past the outer ones as the picture runs.

The gas is now not quite opaque. A crossing adds its light times what the ray still carries, and the ray carries on with 35% of it, so where the bent far side lies over the near side, a little of what is behind shows through, and the disk looks like something with depth rather than a painted ring.

stage-3-disk.metal60 new lines since stage 2
// Stage 3 · The disk
// Gas, not a ring. Streaks of noise are carried round at Kepler's rate for
// their radius, faster nearer in, so the inner edge shears past the outer
// gas; the noise is read on a cylinder so there is no seam. The gas is not
// quite opaque, and a little of what is behind it shows through.

// @param gravity 1 0 1.6 How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param spin 1 0 3 How fast the gas goes round: at 1 the inner edge orbits once in 10 s
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;

// At spin 1 the disk's inner edge goes round once every ten seconds.
constant float ORBIT = 3.27;

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

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

// A random 0 to 1 for each whole-numbered point in space.
float hash3(float3 c) {
    uint x = uint(int(c.x) + 4096);
    uint y = uint(int(c.y) + 4096);
    uint z = uint(int(c.z) + 4096);
    return hashU(x * 73856093u ^ y * 19349663u ^ z * 83492791u);
}

// Value noise in three dimensions: random values at the whole-numbered
// points, blended across the cube between eight of them.
float valueNoise3(float3 x) {
    float3 i = floor(x);
    float3 f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    float a = mix(hash3(i), hash3(i + float3(1, 0, 0)), f.x);
    float b = mix(hash3(i + float3(0, 1, 0)), hash3(i + float3(1, 1, 0)), f.x);
    float c = mix(hash3(i + float3(0, 0, 1)), hash3(i + float3(1, 0, 1)), f.x);
    float d = mix(hash3(i + float3(0, 1, 1)), hash3(i + float3(1, 1, 1)), f.x);
    return mix(mix(a, b, f.y), mix(c, d, f.y), f.z);
}

// Fractal noise: `octaves` layers of value noise, each twice as fine and half
// as strong as the one before, shifted so their grids don't line up. 0 to 1.
float fbm3(float3 x, int octaves) {
    float sum = 0.0, amp = 0.5, total = 0.0;
    for (int k = 0; k < octaves; k++) {
        sum += amp * valueNoise3(x);
        total += amp;
        x = x * 2.03 + float3(1.7, 9.2, 4.1);
        amp *= 0.5;
    }
    return sum / total;
}

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    return tint * star * (0.3 + 2.0 * bright);
}

// How fast the gas goes round at radius r, in turns a second at spin 1:
// faster nearer in, as the planets are.
float orbitRate(float r) {
    return ORBIT / (r * sqrt(r));
}

// The disk's glow where a ray crosses its plane at q. Streaks of noise are
// carried round at the orbit's rate for that radius, so the inner edge
// shears past the outer gas, and it burns brightest just outside the edge.
float diskGlow(float3 q, float t, constant Params& p) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float phi = atan2(q.z, q.x) - orbitRate(r) * p.spin * t;
    // Noise on a cylinder, so there is no seam where phi wraps round.
    float3 at = float3(cos(phi) * 1.6, sin(phi) * 1.6, r * 2.2);
    float n = fbm3(at, 4);
    float streak = 0.35 + 0.65 * smoothstep(0.3, 0.8, n);
    // Brightest just outside the inner edge, fading as the square of r.
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return streak * profile * edge * 24.0;
}

// Walks a ray from the camera, bent by the hole, and gathers what it sees:
// the disk where it crosses the disk's plane, black if it falls into the
// hole, the stars if it escapes.
float3 trace(float3 ro, float3 rd, float t, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    // Angular momentum about the hole, which the bending conserves.
    float3 L = cross(pos, vel);
    float h2 = dot(L, L) * p.gravity;
    float3 light = float3(0.0);
    float carry = 1.0;
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where the bending is strong.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q, t, p);
            if (glow > 0.0) {
                light += carry * glow * p.disk;
                // The gas is not quite opaque: a little of what is behind
                // shows through.
                carry *= 0.35;
            }
        }
        // Gravity: the pull on a ray of light in Schwarzschild's geometry,
        // felt as a bend towards the hole that grows fast as r shrinks.
        vel += -1.5 * h2 * pos / (r2 * r2 * sqrt(r2)) * dt;
        vel = normalize(vel);
        pos = next;
    }
    if (escaped) light += carry * stars(vel);
    return light;
}

// The camera: at rest, looking at the hole, tilted up out of the disk's
// plane as far as the slider says.
float3 cameraPos(constant Params& p) {
    float a = 0.35;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float t = u.time;
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, t, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// Three voices go round the hole as the gas does, each at its own Kepler
// rate, over the hum.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

// One orbiting voice: a tone at `f0` for gas at radius r, going round at
// that radius's Kepler rate, from one ear to the other and back.
float2 orbitVoice(float t, float f0, float r, float phase0,
                  constant Params& p) {
    float w = orbitRate(r) * p.spin;
    float ang = w * t + phase0;
    float ph = TAU * f0 * t;
    float tone = sin(ph) + 0.5 * sin(2.0 * ph + 0.3) + 0.2 * sin(3.0 * ph);
    // Which ear it is nearer to, as it goes round.
    float pan = 0.5 + 0.4 * sin(ang);
    return float2(1.0 - pan, pan) * tone;
}

float2 sound(float t, constant Params& p) {
    float2 v = float2(0.0);
    v += orbitVoice(t, 82.41, DISK_IN, 0.0, p) * 0.22;
    v += orbitVoice(t, 123.47, 5.0, 2.1, p) * 0.14;
    v += orbitVoice(t, 164.81, 7.5, 4.0, p) * 0.09;
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    v += hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * fadeIn * p.volume;
}
8 lines from stage 2 that are gone
float diskGlow(float3 q) {
⋯
    return 0.7 * profile * edge * 24.0;
⋯
float3 trace(float3 ro, float3 rd, constant Params& p) {
⋯
            float glow = diskGlow(q);
            if (glow > 0.0) { light += glow * p.disk; break; }
⋯
    if (escaped) light += stars(vel);
⋯
    float3 seen = trace(ro, rd, p);
⋯
    float2 v = hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;

Stage 4 · Doppler

The gas at the inner edge orbits at four tenths of the speed of light, and at that speed light itself changes. dopplerAt gives the Doppler factor for the gas at each crossing: its speed sqrt(0.5 / r) along the orbit, and the factor 1 / (γ (1 + β · d)), where d is the direction the ray is walking, which is backwards along the light. Gas coming towards you has a factor above 1: its light arrives brighter and hotter, and the brightness goes as the cube of the factor, which is relativistic beaming. Gas going away has a factor below 1 and fades. A second factor, sqrt(1 − 1 / r), is the gravitational redshift: light climbing out of the hole's well loses energy, so the inner edge is dimmed a little more.

The doppler slider mixes the factor between 1 and its true value. At 0 the disk is symmetric; at 1 the right side, where the gas swings towards the camera, burns white, and the left falls to ember. That lopsidedness is the one asymmetry in the picture, and it is the thing real images of a black hole show most clearly.

stage-4-doppler.metal18 new lines since stage 3
// Stage 4 · Doppler
// The gas orbits at close to the speed of light, so its light is beamed:
// the side coming towards you arrives brighter and hotter, the side going
// away dimmer and redder, by a factor that depends on the gas's speed and
// which way the ray leaves it. Climbing out of the hole's gravity costs the
// light a little more, so the inner edge is dimmed too.

// @param gravity 1 0 1.6 How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param spin 1 0 3 How fast the gas goes round: at 1 the inner edge orbits once in 10 s
// @param doppler 1 0 1 How much the gas's speed beams the near side and dims the far side
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;

// At spin 1 the disk's inner edge goes round once every ten seconds.
constant float ORBIT = 3.27;

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

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

// A random 0 to 1 for each whole-numbered point in space.
float hash3(float3 c) {
    uint x = uint(int(c.x) + 4096);
    uint y = uint(int(c.y) + 4096);
    uint z = uint(int(c.z) + 4096);
    return hashU(x * 73856093u ^ y * 19349663u ^ z * 83492791u);
}

// Value noise in three dimensions: random values at the whole-numbered
// points, blended across the cube between eight of them.
float valueNoise3(float3 x) {
    float3 i = floor(x);
    float3 f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    float a = mix(hash3(i), hash3(i + float3(1, 0, 0)), f.x);
    float b = mix(hash3(i + float3(0, 1, 0)), hash3(i + float3(1, 1, 0)), f.x);
    float c = mix(hash3(i + float3(0, 0, 1)), hash3(i + float3(1, 0, 1)), f.x);
    float d = mix(hash3(i + float3(0, 1, 1)), hash3(i + float3(1, 1, 1)), f.x);
    return mix(mix(a, b, f.y), mix(c, d, f.y), f.z);
}

// Fractal noise: `octaves` layers of value noise, each twice as fine and half
// as strong as the one before, shifted so their grids don't line up. 0 to 1.
float fbm3(float3 x, int octaves) {
    float sum = 0.0, amp = 0.5, total = 0.0;
    for (int k = 0; k < octaves; k++) {
        sum += amp * valueNoise3(x);
        total += amp;
        x = x * 2.03 + float3(1.7, 9.2, 4.1);
        amp *= 0.5;
    }
    return sum / total;
}

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    return tint * star * (0.3 + 2.0 * bright);
}

// How fast the gas goes round at radius r, in turns a second at spin 1:
// faster nearer in, as the planets are.
float orbitRate(float r) {
    return ORBIT / (r * sqrt(r));
}

// The disk's glow where a ray crosses its plane at q. Streaks of noise are
// carried round at the orbit's rate for that radius, so the inner edge
// shears past the outer gas, and it burns brightest just outside the edge.
float diskGlow(float3 q, float t, constant Params& p) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float phi = atan2(q.z, q.x) - orbitRate(r) * p.spin * t;
    // Noise on a cylinder, so there is no seam where phi wraps round.
    float3 at = float3(cos(phi) * 1.6, sin(phi) * 1.6, r * 2.2);
    float n = fbm3(at, 4);
    float streak = 0.35 + 0.65 * smoothstep(0.3, 0.8, n);
    // Brightest just outside the inner edge, fading as the square of r.
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return streak * profile * edge * 24.0;
}

// The Doppler factor for gas at q seen along the ray direction d: above 1
// where the gas comes towards you, below 1 where it goes away. The gas
// orbits at a fraction beta of the speed of light; the light travels
// back along -d.
float dopplerAt(float3 q, float3 d, constant Params& p) {
    float r = length(q.xz);
    float beta = sqrt(0.5 / r);
    float3 along = normalize(float3(-q.z, 0.0, q.x));
    float gamma = 1.0 / sqrt(1.0 - beta * beta);
    float shift = 1.0 / (gamma * (1.0 + beta * dot(along, d)));
    // Climbing out of the hole's gravity costs the light some energy too.
    float well = sqrt(1.0 - HORIZON / r);
    return mix(1.0, shift, p.doppler) * well;
}

// Walks a ray from the camera, bent by the hole, and gathers what it sees:
// the disk where it crosses the disk's plane, black if it falls into the
// hole, the stars if it escapes.
float3 trace(float3 ro, float3 rd, float t, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    // Angular momentum about the hole, which the bending conserves.
    float3 L = cross(pos, vel);
    float h2 = dot(L, L) * p.gravity;
    float3 light = float3(0.0);
    float carry = 1.0;
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where the bending is strong.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q, t, p);
            if (glow > 0.0) {
                float dop = dopplerAt(q, vel, p);
                float seen = glow * dop * dop * dop * p.disk;
                light += carry * seen;
                // The gas is not quite opaque: a little of what is behind
                // shows through.
                carry *= 0.35;
            }
        }
        // Gravity: the pull on a ray of light in Schwarzschild's geometry,
        // felt as a bend towards the hole that grows fast as r shrinks.
        vel += -1.5 * h2 * pos / (r2 * r2 * sqrt(r2)) * dt;
        vel = normalize(vel);
        pos = next;
    }
    if (escaped) light += carry * stars(vel);
    return light;
}

// The camera: at rest, looking at the hole, tilted up out of the disk's
// plane as far as the slider says.
float3 cameraPos(constant Params& p) {
    float a = 0.35;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float t = u.time;
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, t, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// The voices are Doppler shifted as they come and go, up in pitch and
// louder as each swings towards you, by the factor that beams the picture.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

// One orbiting voice: a tone at `f0` whose pitch rides the Doppler shift of
// gas at radius r. The shift is 1 + beta cos(angle) to a first order, and
// the phase is the integral of the frequency, so the pitch glides smoothly
// and the sound stays a pure function of time.
float2 orbitVoice(float t, float f0, float r, float phase0,
                  constant Params& p) {
    float beta = sqrt(0.5 / r) * p.doppler;
    float w = TAU * orbitRate(r) * p.spin / TAU;
    float ang = w * t + phase0;
    float ph = TAU * f0 * (t + (w > 1e-4 ? beta * sin(ang) / w : 0.0));
    // Louder coming towards you, as the picture is, by the cube of the shift.
    float shift = 1.0 + beta * cos(ang);
    float loud = shift * shift * shift;
    float tone = sin(ph) + 0.5 * sin(2.0 * ph + 0.3) + 0.2 * sin(3.0 * ph);
    // Which ear it is nearer to, as it goes round.
    float pan = 0.5 + 0.4 * sin(ang);
    return float2(1.0 - pan, pan) * tone * loud;
}

float2 sound(float t, constant Params& p) {
    float2 v = float2(0.0);
    v += orbitVoice(t, 82.41, DISK_IN, 0.0, p) * 0.22;
    v += orbitVoice(t, 123.47, 5.0, 2.1, p) * 0.14;
    v += orbitVoice(t, 164.81, 7.5, 4.0, p) * 0.09;
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    v += hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * fadeIn * p.volume;
}
4 lines from stage 3 that are gone
                light += carry * glow * p.disk;
⋯
    float w = orbitRate(r) * p.spin;
⋯
    float ph = TAU * f0 * t;
⋯
    return float2(1.0 - pan, pan) * tone;

Stage 5 · Black hole

The finished shader adds the camera's drift, a slow circuit of the hole that takes the shadow and the beamed side slowly round, and a faint band of dust across the sky, so the lensing has something to bend besides stars. The palette is the kiln's: the light a pixel gathered sets its heat through 1 − exp(−lum), so the beamed gas saturates to white-hot while the far side stays ember, and the heat slider scales that curve instead of shifting the palette, so the sky stays black at any setting. A little of the light's own colour tints the result, so a blue star stays faintly blue.

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

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float t = u.time;
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(t, p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, t, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}
black-hole.metal, the whole file19 new lines since stage 4
// Ray Kiln · Black Hole
//
// @episode Black Hole
// @length 60
// @variation 30
// @variations 1,2,3,7,8,9,11,12
// @still 22.0
// @teaches Rays bent by gravity, an accretion disk, Doppler colour
// @category Fragment
// @tags raymarching, black holes, gravitational lensing, doppler, accretion disk, value noise, stars, procedural sound
// @sweeps gravity, tilt, doppler
// @final Black hole | A disk of fire wraps over and under the shadow it feeds, the near side rushing towards you white-hot, the far side falling away in ember, the stars behind it bent into a ring.
//
// A black hole with a disk of hot gas around it, seen from close by. Every
// pixel's ray is walked outwards step by step, and at every step gravity
// bends it towards the hole: a ray that comes near enough circles it, and one
// that comes nearer falls in, which is the black shadow in the middle. The
// disk is in one plane; where the ray crosses that plane it picks up the
// disk's light, and because the rays behind the hole are bent up and over it,
// the far side of the disk shows above and below the shadow. The gas orbits
// at close to the speed of light, so the side coming towards you is beamed
// bright and white-hot and the side going away is dim and ember: the Doppler
// effect, and the one asymmetry in an otherwise symmetric picture.
// 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

// @param gravity 1 0 1.6 How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param spin 1 0 3 How fast the gas goes round: at 1 the inner edge orbits once in 10 s
// @param doppler 1 0 1 How much the gas's speed beams the near side and dims the far side
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;
// At spin 1 the disk's inner edge goes round once every ten seconds.
constant float ORBIT = 3.27;

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

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

// A random 0 to 1 for each whole-numbered point in space.
float hash3(float3 c) {
    uint x = uint(int(c.x) + 4096);
    uint y = uint(int(c.y) + 4096);
    uint z = uint(int(c.z) + 4096);
    return hashU(x * 73856093u ^ y * 19349663u ^ z * 83492791u);
}

// Value noise in three dimensions: random values at the whole-numbered
// points, blended across the cube between eight of them.
float valueNoise3(float3 x) {
    float3 i = floor(x);
    float3 f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    float a = mix(hash3(i), hash3(i + float3(1, 0, 0)), f.x);
    float b = mix(hash3(i + float3(0, 1, 0)), hash3(i + float3(1, 1, 0)), f.x);
    float c = mix(hash3(i + float3(0, 0, 1)), hash3(i + float3(1, 0, 1)), f.x);
    float d = mix(hash3(i + float3(0, 1, 1)), hash3(i + float3(1, 1, 1)), f.x);
    return mix(mix(a, b, f.y), mix(c, d, f.y), f.z);
}

// Fractal noise: `octaves` layers of value noise, each twice as fine and half
// as strong as the one before, shifted so their grids don't line up. 0 to 1.
float fbm3(float3 x, int octaves) {
    float sum = 0.0, amp = 0.5, total = 0.0;
    for (int k = 0; k < octaves; k++) {
        sum += amp * valueNoise3(x);
        total += amp;
        x = x * 2.03 + float3(1.7, 9.2, 4.1);
        amp *= 0.5;
    }
    return sum / total;
}

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    // A faint band of dust across the sky, so the lensing has something
    // to bend.
    float band = exp(-dir.y * dir.y * 18.0) * 0.15;
    band *= 0.4 + 0.6 * fbm3(dir * 6.0, 3);
    return kiln(0.45) * band + tint * star * (0.3 + 2.0 * bright);
}

// How fast the gas goes round at radius r, in turns a second at spin 1:
// faster nearer in, as the planets are.
float orbitRate(float r) {
    return ORBIT / (r * sqrt(r));
}

// The disk's glow where a ray crosses its plane at q. Streaks of noise are
// carried round at the orbit's rate for that radius, so the inner edge
// shears past the outer gas, and it burns brightest just outside the edge.
float diskGlow(float3 q, float t, constant Params& p) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float phi = atan2(q.z, q.x) - orbitRate(r) * p.spin * t;
    // Noise on a cylinder, so there is no seam where phi wraps round.
    float3 at = float3(cos(phi) * 1.6, sin(phi) * 1.6, r * 2.2);
    float n = fbm3(at, 4);
    float streak = 0.35 + 0.65 * smoothstep(0.3, 0.8, n);
    // Brightest just outside the inner edge, fading as the square of r.
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return streak * profile * edge * 24.0;
}

// The Doppler factor for gas at q seen along the ray direction d: above 1
// where the gas comes towards you, below 1 where it goes away. The gas
// orbits at a fraction beta of the speed of light; the light travels
// back along -d.
float dopplerAt(float3 q, float3 d, constant Params& p) {
    float r = length(q.xz);
    float beta = sqrt(0.5 / r);
    float3 along = normalize(float3(-q.z, 0.0, q.x));
    float gamma = 1.0 / sqrt(1.0 - beta * beta);
    float shift = 1.0 / (gamma * (1.0 + beta * dot(along, d)));
    // Climbing out of the hole's gravity costs the light some energy too.
    float well = sqrt(1.0 - HORIZON / r);
    return mix(1.0, shift, p.doppler) * well;
}

// Walks a ray from the camera, bent by the hole, and gathers what it sees:
// the disk where it crosses the disk's plane, black if it falls into the
// hole, the stars if it escapes.
float3 trace(float3 ro, float3 rd, float t, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    // Angular momentum about the hole, which the bending conserves.
    float3 L = cross(pos, vel);
    float h2 = dot(L, L) * p.gravity;
    float3 light = float3(0.0);
    float carry = 1.0;
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where the bending is strong.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q, t, p);
            if (glow > 0.0) {
                float dop = dopplerAt(q, vel, p);
                float seen = glow * dop * dop * dop * p.disk;
                light += carry * seen;
                // The gas is not quite opaque: a little of what is behind
                // shows through.
                carry *= 0.35;
            }
        }
        // Gravity: the pull on a ray of light in Schwarzschild's geometry,
        // felt as a bend towards the hole that grows fast as r shrinks.
        vel += -1.5 * h2 * pos / (r2 * r2 * sqrt(r2)) * dt;
        vel = normalize(vel);
        pos = next;
    }
    if (escaped) light += carry * stars(vel);
    return light;
}

// The camera: on a slow drift around the hole, looking at it, tilted up out
// of the disk's plane as far as the slider says.
float3 cameraPos(float t, constant Params& p) {
    float a = 0.35 + 0.04 * t;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float t = u.time;
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(t, p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, t, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// Two voices orbit the hole as the gas does, a low one at the inner edge of
// the disk and a higher one farther out, each going round at its own
// Kepler rate. Each is Doppler shifted as it comes and goes, up in pitch and
// louder as it swings towards you, down and quieter as it swings away, by the
// same factor that beams the picture. Under them the hole hums, a deep slow
// swell, and now and then a star glints.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

// One orbiting voice: a tone at `f0` whose pitch rides the Doppler shift of
// gas at radius r. The shift is 1 + beta cos(angle) to a first order, and
// the phase is the integral of the frequency, so the pitch glides smoothly
// and the sound stays a pure function of time.
float2 orbitVoice(float t, float f0, float r, float phase0,
                  constant Params& p) {
    float beta = sqrt(0.5 / r) * p.doppler;
    float w = TAU * orbitRate(r) * p.spin / TAU;
    float ang = w * t + phase0;
    float ph = TAU * f0 * (t + (w > 1e-4 ? beta * sin(ang) / w : 0.0));
    // Louder coming towards you, as the picture is, by the cube of the shift.
    float shift = 1.0 + beta * cos(ang);
    float loud = shift * shift * shift;
    float tone = sin(ph) + 0.5 * sin(2.0 * ph + 0.3) + 0.2 * sin(3.0 * ph);
    // Which ear it is nearer to, as it goes round.
    float pan = 0.5 + 0.4 * sin(ang);
    return float2(1.0 - pan, pan) * tone * loud;
}

float2 sound(float t, constant Params& p) {
    float2 v = float2(0.0);
    v += orbitVoice(t, 82.41, DISK_IN, 0.0, p) * 0.22;
    v += orbitVoice(t, 123.47, 5.0, 2.1, p) * 0.14;
    v += orbitVoice(t, 164.81, 7.5, 4.0, p) * 0.09;
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    v += hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    // Star glints: a short high ping every so often, somewhere in the stereo.
    float slot = floor(t / 1.7);
    for (int k = -1; k <= 0; k++) {
        float j = slot + float(k);
        if (j < 0.0) continue;
        uint id = uint(j) * 5u + 11u;
        if (hashU(id) > 0.45) continue;
        float s = t - j * 1.7 - 1.2 * hashU(id + 1u);
        if (s < 0.0) continue;
        float f = 1760.0 * exp2(floor(hashU(id + 2u) * 5.0) * 2.0 / 12.0);
        float ping = sin(TAU * f * s) * exp(-s * 9.0) * (1.0 - exp(-s * 400.0));
        float pan = hashU(id + 3u);
        v += float2(1.0 - pan, pan) * ping * 0.05;
    }
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * fadeIn * p.volume;
}
4 lines from stage 4 that are gone
    return tint * star * (0.3 + 2.0 * bright);
⋯
float3 cameraPos(constant Params& p) {
    float a = 0.35;
⋯
    float3 ro = cameraPos(p);

The sound

Three voices go round the hole as the gas does. Each is a tone, a low E, a B and an E above it, with two overtones, placed at a radius of the disk, 3, 5 and 7.5, and going round at that radius's Kepler rate from orbitRate, so the inner one laps the others. Each voice is Doppler shifted by the same first-order factor that beams the picture, 1 + β cos(angle): its pitch glides up as it swings towards you and down as it swings away, and it is louder coming than going, by the cube of the factor, as the gas is. The phase is the integral of the frequency, f₀ (t + β sin(angle) / ω), so the glide is smooth and the sound stays a pure function of time. Each voice is panned round the stereo with its orbit. Under them the hole hums, a low E at 41 Hz and its octave, a little detuned between the ears so it beats, swelling slowly, and now and then a star glints: a short high ping somewhere in the stereo. tanh keeps it under full scale and it fades in over two seconds.

At the defaults the peak is about −7.2 dBFS and the average (RMS) about −18.1 dBFS over the minute. With the gas going round fast and the disk bright (spin 3, disk 1.5) the peak is −6.3 dBFS and the RMS −16.3; with the disk still and no Doppler (spin 0, doppler 0) the voices stop gliding and it is −10.1 and −18.7. The doppler slider changes the sound as it changes the picture: at 0 the voices hold their pitch.

New in Ray Kiln

Sound while you edit, in the Mac app. The app now plays a shader's sound live, following the picture's clock and the sliders. The sound is the export's own kernel, run on the GPU a chunk ahead: each chunk of 43 ms is asked for the shader time at which it will be heard, and whatever changes between two chunks, an edit that compiles, a slider, a restart, a pause, is joined with a 10 ms crossfade between what was playing and what plays now, so nothing clicks. A slider's change is heard about a tenth of a second after you move it. The inspector's Sound switch and Shader > Mute Sound silence it, and the readout beside the switch gives the peak of the last chunk in dBFS, or says why nothing plays. A pause fades out, and a shader with no sound costs nothing. It came in Ray Kiln 0.12; see the app.

Black Hole is the shader to try it with. Open it in the app and drag doppler from 1 to 0 while it plays: the three voices stop gliding and the right side of the disk stops burning white, at the same moment. Drag spin and the voices go round faster with the gas.

Try this

  • Set gravity to 0 and back to 1: a flat ring and a ball, then the same ring wrapped over and under the shadow.
  • Set tilt to 0.5 to look down on the disk from above: a ring of streaks with the shadow in the middle, and the beamed side still brighter.
  • Set doppler to 0: the disk is symmetric, and the sound's voices hold their pitch. Then 1 again, and watch the right side.
  • Set gravity to 1.6: the shadow grows and the Einstein ring moves out; the far side of the disk is bent higher.
  • In trace, change carry *= 0.35 to carry = 0.0: the disk is opaque again, and the far side no longer shows through the near side.
  • In orbitRate, change r * sqrt(r) to r: the gas goes round at the wrong law, and the streaks shear less.
  • Set heat to 1 with disk at 1.5: the whole near side goes white; at 0 only the beamed gas gets past orange.

Parameters

ParameterDefaultRangeWhat it does
gravity10 to 1.6How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
tilt0.140 to 0.5How far the camera looks down on the disk: 0 is edge on, 0.5 from above
spin10 to 3How fast the gas goes round: at 1 the inner edge orbits once in 10 s
doppler10 to 1How much the gas's speed beams the near side and dims the far side
disk10 to 1.5How bright the disk burns
heat0.50 to 1Slides the palette from deep ember to white-hot
volume0.80 to 1Loudness 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.

black-hole.metal
// Ray Kiln · Black Hole
//
// @episode Black Hole
// @length 60
// @variation 30
// @variations 1,2,3,7,8,9,11,12
// @still 22.0
// @teaches Rays bent by gravity, an accretion disk, Doppler colour
// @category Fragment
// @tags raymarching, black holes, gravitational lensing, doppler, accretion disk, value noise, stars, procedural sound
// @sweeps gravity, tilt, doppler
// @final Black hole | A disk of fire wraps over and under the shadow it feeds, the near side rushing towards you white-hot, the far side falling away in ember, the stars behind it bent into a ring.
//
// A black hole with a disk of hot gas around it, seen from close by. Every
// pixel's ray is walked outwards step by step, and at every step gravity
// bends it towards the hole: a ray that comes near enough circles it, and one
// that comes nearer falls in, which is the black shadow in the middle. The
// disk is in one plane; where the ray crosses that plane it picks up the
// disk's light, and because the rays behind the hole are bent up and over it,
// the far side of the disk shows above and below the shadow. The gas orbits
// at close to the speed of light, so the side coming towards you is beamed
// bright and white-hot and the side going away is dim and ember: the Doppler
// effect, and the one asymmetry in an otherwise symmetric picture.
// 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

// @param gravity 1 0 1.6 How strongly the hole bends light: 0 is straight rays, 1 is Einstein's
// @param tilt 0.14 0 0.5 How far the camera looks down on the disk: 0 is edge on, 0.5 from above
// @param spin 1 0 3 How fast the gas goes round: at 1 the inner edge orbits once in 10 s
// @param doppler 1 0 1 How much the gas's speed beams the near side and dims the far side
// @param disk 1 0 1.5 How bright the disk burns
// @param heat 0.5 0 1 Slides the palette from deep ember to white-hot
// @param volume 0.8 0 1 Loudness of the sound

constant float TAU = 6.2831853;
// Distances are in Schwarzschild radii: the horizon is at r = 1 and the
// disk runs from r = 3 out to r = 9.
constant float HORIZON = 1.0;
constant float DISK_IN = 3.0;
constant float DISK_OUT = 9.0;
// How far the camera is from the hole.
constant float DISTANCE = 19.0;
// The most steps a ray is walked, and how far out counts as escaped.
constant int STEPS = 200;
constant float FAR = 28.0;
// At spin 1 the disk's inner edge goes round once every ten seconds.
constant float ORBIT = 3.27;

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

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

// A random 0 to 1 for each whole-numbered point in space.
float hash3(float3 c) {
    uint x = uint(int(c.x) + 4096);
    uint y = uint(int(c.y) + 4096);
    uint z = uint(int(c.z) + 4096);
    return hashU(x * 73856093u ^ y * 19349663u ^ z * 83492791u);
}

// Value noise in three dimensions: random values at the whole-numbered
// points, blended across the cube between eight of them.
float valueNoise3(float3 x) {
    float3 i = floor(x);
    float3 f = x - i;
    f = f * f * (3.0 - 2.0 * f);
    float a = mix(hash3(i), hash3(i + float3(1, 0, 0)), f.x);
    float b = mix(hash3(i + float3(0, 1, 0)), hash3(i + float3(1, 1, 0)), f.x);
    float c = mix(hash3(i + float3(0, 0, 1)), hash3(i + float3(1, 0, 1)), f.x);
    float d = mix(hash3(i + float3(0, 1, 1)), hash3(i + float3(1, 1, 1)), f.x);
    return mix(mix(a, b, f.y), mix(c, d, f.y), f.z);
}

// Fractal noise: `octaves` layers of value noise, each twice as fine and half
// as strong as the one before, shifted so their grids don't line up. 0 to 1.
float fbm3(float3 x, int octaves) {
    float sum = 0.0, amp = 0.5, total = 0.0;
    for (int k = 0; k < octaves; k++) {
        sum += amp * valueNoise3(x);
        total += amp;
        x = x * 2.03 + float3(1.7, 9.2, 4.1);
        amp *= 0.5;
    }
    return sum / total;
}

// The stars behind everything, from the direction a ray finally leaves in:
// a grid of cells over the sky, a star in some, a few of them bright.
float3 stars(float3 dir) {
    float2 sph = float2(atan2(dir.z, dir.x), asin(clamp(dir.y, -1.0, 1.0)));
    sph *= 90.0;
    float2 cell = floor(sph);
    uint k = uint(int(cell.x) + 1000) * 7919u;
    k += uint(int(cell.y) + 1000) * 104729u;
    float2 at = float2(hashU(k), hashU(k + 1u)) * 0.8 + 0.1;
    float star = step(0.84, hashU(k + 2u));
    star *= smoothstep(0.11, 0.0, length(fract(sph) - at));
    float bright = pow(hashU(k + 3u), 3.0);
    float3 tint = mix(float3(0.7, 0.8, 1.0), kiln(0.85), hashU(k + 4u));
    // A faint band of dust across the sky, so the lensing has something
    // to bend.
    float band = exp(-dir.y * dir.y * 18.0) * 0.15;
    band *= 0.4 + 0.6 * fbm3(dir * 6.0, 3);
    return kiln(0.45) * band + tint * star * (0.3 + 2.0 * bright);
}

// How fast the gas goes round at radius r, in turns a second at spin 1:
// faster nearer in, as the planets are.
float orbitRate(float r) {
    return ORBIT / (r * sqrt(r));
}

// The disk's glow where a ray crosses its plane at q. Streaks of noise are
// carried round at the orbit's rate for that radius, so the inner edge
// shears past the outer gas, and it burns brightest just outside the edge.
float diskGlow(float3 q, float t, constant Params& p) {
    float r = length(q.xz);
    if (r < DISK_IN || r > DISK_OUT) return 0.0;
    float phi = atan2(q.z, q.x) - orbitRate(r) * p.spin * t;
    // Noise on a cylinder, so there is no seam where phi wraps round.
    float3 at = float3(cos(phi) * 1.6, sin(phi) * 1.6, r * 2.2);
    float n = fbm3(at, 4);
    float streak = 0.35 + 0.65 * smoothstep(0.3, 0.8, n);
    // Brightest just outside the inner edge, fading as the square of r.
    float profile = 1.0 - sqrt(DISK_IN / r) * 0.9;
    profile *= (DISK_IN * DISK_IN) / (r * r);
    float edge = smoothstep(DISK_IN, DISK_IN + 0.4, r);
    edge *= smoothstep(DISK_OUT, DISK_OUT - 2.5, r);
    return streak * profile * edge * 24.0;
}

// The Doppler factor for gas at q seen along the ray direction d: above 1
// where the gas comes towards you, below 1 where it goes away. The gas
// orbits at a fraction beta of the speed of light; the light travels
// back along -d.
float dopplerAt(float3 q, float3 d, constant Params& p) {
    float r = length(q.xz);
    float beta = sqrt(0.5 / r);
    float3 along = normalize(float3(-q.z, 0.0, q.x));
    float gamma = 1.0 / sqrt(1.0 - beta * beta);
    float shift = 1.0 / (gamma * (1.0 + beta * dot(along, d)));
    // Climbing out of the hole's gravity costs the light some energy too.
    float well = sqrt(1.0 - HORIZON / r);
    return mix(1.0, shift, p.doppler) * well;
}

// Walks a ray from the camera, bent by the hole, and gathers what it sees:
// the disk where it crosses the disk's plane, black if it falls into the
// hole, the stars if it escapes.
float3 trace(float3 ro, float3 rd, float t, constant Params& p) {
    float3 pos = ro;
    float3 vel = rd;
    // Angular momentum about the hole, which the bending conserves.
    float3 L = cross(pos, vel);
    float h2 = dot(L, L) * p.gravity;
    float3 light = float3(0.0);
    float carry = 1.0;
    bool escaped = false;
    for (int i = 0; i < STEPS; i++) {
        float r2 = dot(pos, pos);
        if (r2 < HORIZON * HORIZON) break;
        if (r2 > FAR * FAR && dot(pos, vel) > 0.0) { escaped = true; break; }
        // Smaller steps near the hole, where the bending is strong.
        float dt = clamp(0.045 * sqrt(r2), 0.03, 0.5);
        float3 next = pos + vel * dt;
        // Crossing the disk's plane: find where, and take its light.
        if (pos.y * next.y < 0.0) {
            float f = pos.y / (pos.y - next.y);
            float3 q = mix(pos, next, f);
            float glow = diskGlow(q, t, p);
            if (glow > 0.0) {
                float dop = dopplerAt(q, vel, p);
                float seen = glow * dop * dop * dop * p.disk;
                light += carry * seen;
                // The gas is not quite opaque: a little of what is behind
                // shows through.
                carry *= 0.35;
            }
        }
        // Gravity: the pull on a ray of light in Schwarzschild's geometry,
        // felt as a bend towards the hole that grows fast as r shrinks.
        vel += -1.5 * h2 * pos / (r2 * r2 * sqrt(r2)) * dt;
        vel = normalize(vel);
        pos = next;
    }
    if (escaped) light += carry * stars(vel);
    return light;
}

// The camera: on a slow drift around the hole, looking at it, tilted up out
// of the disk's plane as far as the slider says.
float3 cameraPos(float t, constant Params& p) {
    float a = 0.35 + 0.04 * t;
    float up = sin(p.tilt * TAU * 0.25);
    float flat = cos(p.tilt * TAU * 0.25);
    return DISTANCE * float3(flat * sin(a), up, flat * cos(a));
}

float4 shade(float2 uv, constant Uniforms& u, constant Params& p) {
    float t = u.time;
    float side = min(u.resolution.x, u.resolution.y);
    float2 q = (uv - 0.5) * u.resolution / side * 2.0;
    float3 ro = cameraPos(t, p);
    float3 fw = normalize(-ro);
    float3 rt = normalize(cross(fw, float3(0.0, 1.0, 0.0)));
    float3 up = cross(rt, fw);
    float3 rd = normalize(q.x * rt + q.y * up + 1.7 * fw);
    float3 seen = trace(ro, rd, t, p);
    // The light through the palette: its amount picks the heat, and what
    // colour it had (a blue star, the white of beamed gas) tints it.
    float lum = dot(seen, float3(0.3, 0.5, 0.2));
    float heat = 1.0 - exp(-lum * (0.7 + 0.8 * p.heat));
    float3 col = kiln(heat);
    col *= mix(float3(1.0), seen / max(lum, 1e-3), 0.25);
    return float4(clamp(col, 0.0, 1.0), 1.0);
}

// ---- Sound -----------------------------------------------------------------------------------------------------
// The sound is part of the shader: a function of time, run on the GPU 48,000 times a second and put in the video.
// Two voices orbit the hole as the gas does, a low one at the inner edge of
// the disk and a higher one farther out, each going round at its own
// Kepler rate. Each is Doppler shifted as it comes and goes, up in pitch and
// louder as it swings towards you, down and quieter as it swings away, by the
// same factor that beams the picture. Under them the hole hums, a deep slow
// swell, and now and then a star glints.

// A value that drifts smoothly: random values `rate` times a second, blended.
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);
}

// One orbiting voice: a tone at `f0` whose pitch rides the Doppler shift of
// gas at radius r. The shift is 1 + beta cos(angle) to a first order, and
// the phase is the integral of the frequency, so the pitch glides smoothly
// and the sound stays a pure function of time.
float2 orbitVoice(float t, float f0, float r, float phase0,
                  constant Params& p) {
    float beta = sqrt(0.5 / r) * p.doppler;
    float w = TAU * orbitRate(r) * p.spin / TAU;
    float ang = w * t + phase0;
    float ph = TAU * f0 * (t + (w > 1e-4 ? beta * sin(ang) / w : 0.0));
    // Louder coming towards you, as the picture is, by the cube of the shift.
    float shift = 1.0 + beta * cos(ang);
    float loud = shift * shift * shift;
    float tone = sin(ph) + 0.5 * sin(2.0 * ph + 0.3) + 0.2 * sin(3.0 * ph);
    // Which ear it is nearer to, as it goes round.
    float pan = 0.5 + 0.4 * sin(ang);
    return float2(1.0 - pan, pan) * tone * loud;
}

float2 sound(float t, constant Params& p) {
    float2 v = float2(0.0);
    v += orbitVoice(t, 82.41, DISK_IN, 0.0, p) * 0.22;
    v += orbitVoice(t, 123.47, 5.0, 2.1, p) * 0.14;
    v += orbitVoice(t, 164.81, 7.5, 4.0, p) * 0.09;
    // The hum of the hole: a low E and its octave, swelling slowly, a little
    // apart in each ear.
    float swell = 0.5 + 0.5 * sin(TAU * t / 13.0 + 2.0 * noise1(t, 0.1, 7u));
    float2 hum = float2(sin(TAU * 41.2 * t), sin(TAU * 41.3 * t + 1.0));
    hum += 0.5 * float2(sin(TAU * 82.4 * t + 0.5), sin(TAU * 82.6 * t));
    v += hum * 0.13 * (0.4 + 0.6 * swell) * p.disk;
    // Star glints: a short high ping every so often, somewhere in the stereo.
    float slot = floor(t / 1.7);
    for (int k = -1; k <= 0; k++) {
        float j = slot + float(k);
        if (j < 0.0) continue;
        uint id = uint(j) * 5u + 11u;
        if (hashU(id) > 0.45) continue;
        float s = t - j * 1.7 - 1.2 * hashU(id + 1u);
        if (s < 0.0) continue;
        float f = 1760.0 * exp2(floor(hashU(id + 2u) * 5.0) * 2.0 / 12.0);
        float ping = sin(TAU * f * s) * exp(-s * 9.0) * (1.0 - exp(-s * 400.0));
        float pan = hashU(id + 3u);
        v += float2(1.0 - pan, pan) * ping * 0.05;
    }
    float fadeIn = smoothstep(0.0, 2.0, t);
    return tanh(v * 1.3) / 1.5 * 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.14, 1, 1, 1, 0.5, 0.8]     // the @param defaults, in the order they are declared
encoder.setRenderPipelineState(pipeline)
encoder.setFragmentBytes(&uniforms, length: MemoryLayout<Uniforms>.stride, index: 0)
encoder.setFragmentBytes(&params, length: params.count * MemoryLayout<Float>.stride, index: 1)
encoder.drawPrimitives(type: .triangle, vertexStart: 0, vertexCount: 3)

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.

More in 3D Scenes