Space Travel · Interactive explainer · Scene 13

A Galaxy of Ellipses

Eight hundred thousand suns, glowing gas and dark dust — and not a single line of gravity. How a spiral galaxy is computed from nothing but each star's index and the clock, with the real shaders running right here in your browser.

The scene

Through the wormhole we arrive above another galaxy. In the demo the scene lasts sixty-four beats — 30.72 seconds at 125 BPM, from 3:42.7 to 4:13.4 — and fades in from the wormhole over a single beat. We start 260 units from the center, 55° above the disc, circling slowly while we close in. Two-thirds of the way through, the camera dives: down to 32 units from the core and only 6° above the galactic plane, skimming through the stars and the dust toward the bright center, where the next scene's black hole waits.

The galaxy has a warm, puffy bulge, a thin disc of white, blue and orange stars, pink and blue gas clouds, dark dust and a glowing core. Its two arms look like the result of gravity — but nothing is simulated. Every star is one GPU vertex whose position is computed from scratch each frame from two numbers: its index and the time. Each star travels on its own ellipse; the ellipses are slightly rotated relative to each other, and where neighbouring ellipses crowd together, the stars crowd together too. Those crowds are the spiral arms. This page builds the idea up piece by piece.

Live: the real shaders

This is the actual galaxy scene: the same background, star and cloud shaders as the demo, ported to WebGL2, driven by a JavaScript copy of the scene's C++ render() (camera path, uniforms and draw order are identical). Drag the time slider to scrub: because every frame is a pure function of time, scrubbing is just "render the frame for this time".

Background

Nine small ideas make up the whole scene. Each one is simple; the spiral appears only when they are combined.

1 · A galaxy from an index

A GPU usually draws triangles from a list of vertices stored in memory. Here there is no list. The C++ code asks the GPU to draw 800,000 points (GL_POINTS, one vertex each) from an empty vertex array, and the vertex shader receives just one number: gl_VertexID, the index of the vertex being processed, 0 to 799,999. That number is the star. Everything about it — bulge or disc, distance from the center, height, position on its orbit, color, brightness — must be derived from that index.

The tool is a hash function: a formula that scrambles an integer into a number that looks random but is perfectly repeatable. The galaxy uses hashU(id, channel), which returns a value in [0, 1]; each property reads a different channel (1 to 7), so the seven numbers of one star are unrelated to each other.

// common/noise.glsl — integer hash for large ids
uint pcgHash(uint v) {
    uint state = v * 747796405u + 2891336453u;                    // one step of a PCG generator
    uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u; // data-dependent shift: mix the bits
    return (word >> 22u) ^ word;
}
float hashU(uint id, uint channel) { return float(pcgHash(id ^ pcgHash(channel))) / 4294967295.0; }

pcgHash is the GPU hash recommended by Jarzynski and Olano ("Hash Functions for GPU Rendering", 2020), built from one step of Melissa O'Neill's PCG random number generator (2014). Why not the float hash used by other scenes (such as the warp's hash12)? Because of precision. A 32-bit float carries about 7 significant digits. hash12 begins by multiplying the id by 0.1031; for star 800,000 that is 82,480 — and at that size a float can only step in increments of 1/128. The fractional part, which is all the hash keeps, has only 128 possible values left. Integer arithmetic has no such problem: every one of the 32 bits takes part in the scramble.

2 · Sampling a distribution

A hash gives a uniform number u: every value between 0 and 1 is equally likely. But a galaxy is not uniform. Real disc galaxies get fainter roughly exponentially with radius (Ken Freeman described this in 1970): every few thousand light-years outward, the light drops by the same factor. Inspired by that profile, the demo samples an exponential distribution of orbital semi-major axes a. This is an artistic shortcut, not an exact exponential surface-density model.

The recipe is inverse transform sampling. Write down the cumulative distribution F(a): the fraction of stars with a semi-major axis smaller than a. It climbs from 0 to 1 at the maximum axis length. Pick u uniformly, and find the axis length where F(a) = u. Where F climbs steeply (many stars), a wide range of u maps into a small range of axis lengths, so those lengths come up often. For p(a) ∝ exp(−a/h) with scale length h = 28, cut off at E = 110, the equation F(a) = u can be solved exactly:

That is diskRadius(u) in galaxy.glsl. Half of the disc stars have semi-major axes below ≈ 18.9; their instantaneous distances also depend on their elliptical orbits. In a circular disc, an exponential surface density Σ(r) ∝ exp(−r/h) would instead require p(r) ∝ r · exp(−r/h), because outer annuli have more area. The bulge uses a simpler shortcut, a = 14 · u1.6: raising u to a power above 1 pushes values toward 0, so bulge stars pile up toward the center (half have semi-major axes below 4.6 units). There is another classic way to sample, rejection sampling — throw random candidates and discard some — but a vertex shader must produce exactly one star per vertex, with no retries, so a direct formula is the right tool.

3 · A bulge and a thin disc

The first hash channel decides the star's family: if hashU(id, 1) < 0.18 it belongs to the bulge, the round swarm of old stars at the center of many spiral galaxies; otherwise to the disc. (18 % is the uniform uBulgeFraction.) Each star then gets a height above or below the galactic plane, drawn from a Gaussian — the bell curve — whose width σ depends on the family and the radius:

• bulge: σ = 1 + 7 · (1 − a / 16) — 8 units at the very center, under 2 at the bulge edge;
• disc: σ = mix(2.5, 0.8, smoothstep(10, 60, a)) — 2.5 units near the bulge, thinning to 0.8 beyond radius 60.

Hashes give uniform numbers, so the Gaussian needs its own trick: the Box–Muller transform (George Box and Mervin Muller, 1958) turns two uniform numbers into one normally-distributed number:

The √(−2 ln u₁) part is a radius whose square is exponentially distributed, cos(2π u₂) projects it on a random direction; together they give an exact standard normal value (the shader guards u₁ ≥ 10⁻⁶ so the logarithm stays finite).

4 · Orbits and the winding problem

Every star orbits the center. How fast? Planets around the Sun obey Kepler: the farther out, the slower they move, speed falling as 1/√r, because almost all the mass is in the center. Galaxies surprised astronomers: in the 1970s Vera Rubin and Kent Ford measured that stars far from the center move about as fast as those closer in — a flat rotation curve (one of the main clues for dark matter). The demo copies that: each star's angle advances at

The + 3 keeps the center from spinning infinitely fast; near it the speed rises gently, like a solid body. The angle is computed directly, φ = phase + ω(a) · t, never accumulated, so the scene can be scrubbed.

A flat rotation curve means differential rotation: inner stars complete their circles much sooner than outer ones (over this scene, a star at radius 10 turns 135°, one at radius 100 only 17°). Now suppose the arms were made of particular stars — a line of stars that stay together. Differential rotation would wrap that line into an ever tighter spiral. Galaxies have turned dozens of times since they formed, so such "material" arms would be wound to invisibility by now. This is the winding problem, and it means spiral arms cannot be fixed groups of stars.

5 · Arms from twisted ellipses

The accepted answer to the winding problem is the density-wave theory of C. C. Lin and Frank Shu (1964): spiral arms are a pattern, like a traffic jam on a highway. Cars enter the jam, slow down, leave it — the jam persists even though the cars in it keep changing. In a galaxy, stars and gas drift into an arm, linger, and move on.

A beautiful kinematic picture of this, shown by Agris Kalnajs in 1973, is what the demo implements. Give every star an elliptical orbit instead of a circle, and rotate each ellipse a little more than the one inside it. Every ellipse on its own is perfectly smooth. But where neighbouring ellipses come close together, stars from many orbits are packed into a small space — and those crowded places line up into two spiral arms. No star belongs to an arm; each one passes through both arms on every orbit.

// scenes/galaxy/galaxy.glsl
vec3 orbitPosition(float a, float phase, float height) {
    float phi = phase + uOrbitSpeed / (a + 3.0) * uTime;       // where along the orbit (§4)
    vec2 e = vec2(a * cos(phi), a * axisRatio(a) * sin(phi)); // point on an ellipse, axes a and a·k
    float theta = a * uTwist + uPattern;                   // bigger orbits are rotated more
    float c = cos(theta), s = sin(theta);
    return vec3(c * e.x - s * e.y, height, s * e.x + c * e.y);  // rotate in the disc plane (xz)
}

Three numbers shape the galaxy. The axis ratio k(a) — minor axis over major axis — is 0.6 in the center (long, bar-like orbits, so the core looks like a short bar), 0.82 from radius 18 to 60, and rounds off to 0.95 at the rim, where the arms fade. The twist (0.045 radians per unit of radius) is how much each orbit is rotated relative to its inner neighbour: between radius 10 and 110 the orientation turns 4.5 radians, about 258°. And uPattern = 0.02 · t rotates the whole family of ellipses slowly — the pattern speed of the density wave — while the stars run along their ellipses at their own, faster pace.

Because an ellipse looks the same after half a turn, the crowding happens on two opposite sides: the model always makes two arms. And because the orientation grows linearly with radius, the arms are close to Archimedean spirals (angle proportional to radius), a little different from the roughly logarithmic spirals measured in real galaxies.

6 · Star sprites: size and light

Each star is drawn as a point sprite: the vertex shader sets gl_PointSize, the GPU covers a square of that many pixels around the projected position, and the fragment shader receives gl_PointCoord — where inside the square the pixel lies, from 0 to 1. The demo turns that into a soft round dot, exp(−3 · |d|²) with d going from −1 to 1 across the square.

Luminosity. 0.3 + hashU(id, 7)⁸ · 25: raising a uniform number to the 8th power pushes almost all values toward 0, so half the stars have luminosity under 0.4 and only about one in ten exceeds 10 — "mostly faint, a few giants", as in a real stellar population. Color comes from channel 6: bulge stars are warm orange to cream; disc stars are 25 % young blue, 45 % white-yellow and 30 % orange.

Size. size = clamp(pointScale / range · (0.7 + 0.15 · lum), 1, 6) pixels, where range is the distance to the camera and pointScale = 220 at 1080 lines (scaled with the resolution). Near and bright stars get bigger, but never more than 6 pixels, and never under 1 so distant stars don't flicker away.

Brightness divided by size. The color is divided by size "to keep total light independent of size", says the shader comment. Strictly, a sprite's area grows with size², so no division by size can make the total exactly constant. But the soft profile is concentrated — its effective area is only about size² / 4 pixels — so over the allowed range of 1 to 6 pixels the total light stays between about half and one and a half times that of a 1-pixel star (the figure below adds up the actual pixels). Without the division, a 6-pixel star would deliver about nine times the light of a 1-pixel one. Note also that there is no 1/distance² dimming: once a star is down to 1 pixel, its light no longer depends on distance at all.

Additive blending. All stars are drawn with glBlendFunc(GL_ONE, GL_ONE): each sprite's color is added to the pixel. Light behaves like that, the order doesn't matter (so 800,000 stars need no sorting), and thousands of faint stars overlapping in the core add up to values far above 1.0. The image is rendered in 16-bit floating point (HDR) and only compressed for the screen at the end by a tone-mapping curve, which also feeds a glow (bloom) around the brightest areas.

7 · Gas glows, dust absorbs

Gas and dust are drawn as billboards: flat squares that always face the camera. Each one is built from six vertices (two triangles) with no vertex data, as in the warp: gl_VertexID / 6 picks the cloud, gl_VertexID % 6 the corner, and the corner is placed at center + (camRight · x + camUp · y) · size — using the camera's own right and up vectors is what makes the square face the viewer. The fragment shader paints a soft puff, exp(−3 · |corner|²). The clouds' centers use the same orbitPosition as the stars, so they gather in the same arms and turn with them. Their ids start at 7,000,000 (gas) and 9,000,000 (dust), so they hash differently from any star.

Gas (6,000 puffs, radius 15 to 95, layer σ = 0.9): 70 % pink — the color of the hydrogen glow of star-forming regions — and 30 % blue. It is drawn with additive blending, like the stars: it adds light.

Dust (14,000 puffs, radius 10 to 95, a thinner layer σ = 0.45) must remove light. It is drawn with glBlendFunc(GL_ZERO, GL_ONE_MINUS_SRC_COLOR), which makes the GPU compute

so the dust's "color" is the fraction of light it absorbs. That color is (0.10, 0.13, 0.18) · 0.25: more blue absorbed than red. Real interstellar dust does the same — it dims and reddens the starlight behind it, which is why the dust lanes of real galaxies look brown. Many puffs overlapping multiply: k layers let through (1 − s)^k ≈ e^−k·s of the light, the Beer–Lambert law of absorption in a medium.

8 · Haze, core and sky

Before any star is drawn, a fullscreen fragment shader paints the backdrop. For each pixel it builds the camera ray — the line from the eye through that pixel (cameraRay, the pinhole camera of the warp explainer) — and adds three things:

• Distant stars: the shared starSky of the demo, at 35 % strength, plus a faint blue-black.
• Disc haze: the unresolved glow of billions of stars too faint to draw. The ray meets the disc plane y = 0 at distance t = −ro.y / rd.y (from ro.y + t · rd.y = 0); if t > 0, the crossing point's distance from the center r gives the glow e−r/30 · smoothstep(115, 70, r), warm white, times uHaze = 0.25. (The reversed edges 115 → 70 make the smoothstep fall from 1 to 0 at the rim — a common idiom, formally undefined in the GLSL specification but computed the same way by every GPU.)
• Core glow: the brightness depends on how close the ray passes to the center. The closest approach of a line to a point is found by projecting: tc = −dot(ro, rd), r = |ro + rd · tc|. Done as is, the glow would be a sphere; the shader first scales y by 2.5, measuring the distance in a squashed frame, so the glow is a flattened ellipsoid like a real bulge. Two exponentials make a hot center and a wide halo: 2 · e−r/5 + 0.35 · e−r/16. On the C++ side uCoreGlow is multiplied by 1 + 0.25 · beatPulse, a pulse that jumps on every beat and decays as e^−4·s: the core breathes with the music.

9 · The camera: orbit, then dive

The camera always looks at the center, from a position given in spherical coordinates: a distance, an elevation above the disc and an azimuth around the vertical axis. Two smoothstep ramps drive it — smoothstep goes from 0 to 1 between two edges with a gentle start and finish:

p        = t / 30.72                         // progress through the scene
approach = smoothstep(0.00, 0.65, p)          // first two thirds: close in
dive     = smoothstep(0.62, 1.00, p)          // last third: dive (overlaps slightly, so there's no pause)
distance  = mix(mix(260, 150, approach), 32, dive)
elevation = mix(mix(55°, 32°, approach), 6°, dive)
azimuth   = 0.6 + 0.12 · t                    // constant circling, 211° over the scene
position  = distance · (cos elev · sin azim,  sin elev,  cos elev · cos azim)

Nesting mix calls chains two moves: the inner mix is the approach, the outer one blends from wherever the approach is toward the dive target. Camera::lookAt turns the position into the forward, right and up vectors, with a focal length of 1.3 screen heights (a vertical field of view of 2 · atan(0.5 / 1.3) ≈ 42°); the stars and clouds use the matching perspective matrix, the background the matching rays.

The algorithm, step by step

Here is the complete scene, in order, with each step linked to the idea above that explains it.

// CPU, once per frame
p        = t / 30.72
approach = smoothstep(0, 0.65, p)
dive     = smoothstep(0.62, 1, p)
dist     = mix(mix(260, 150, approach), 32, dive)
elev     = mix(mix(55°, 32°, approach), 6°, dive)
azim     = 0.6 + 0.12·t
camera   = lookAt(dist · dir(elev, azim) → 0)
viewProj = perspective(2·atan(0.5/1.3)) · view
pattern  = 0.02·t

// draw 1: background, fullscreen, replaces
for each pixel:
  rd   = cameraRay(pixel)
  col  = starSky(rd)·0.35 + (0.001, 0.001, 0.002)
  t    = −ro.y / rd.y
  if t > 0:
    r    = |(ro + rd·t).xz|
    col += warm · e^(−r/30)
           · smoothstep(115, 70, r) · 0.25
  q, dq = (ro, rd) with y × 2.5
  r     = |q + dq · max(−q·dq, 0)|
  col  += gold · (2e^(−r/5) + 0.35e^(−r/16))
          · (1 + 0.25·beatPulse)

// draw 2: stars, GL_POINTS, additive
for each vertex id in 0 … 799,999:
  bulge = hashU(id,1) < 0.18
  u     = hashU(id,2)
  a     = bulge ? 14·u^1.6
                : −28·ln(1 − u·(1 − e^(−110/28)))
  σ     = bulge ? 1 + 7·(1 − a/16)
                : mix(2.5, 0.8, smoothstep(10, 60, a))
  y     = boxMuller(hashU(id,3), hashU(id,4)) · σ

  φ   = hashU(id,5)·2π + orbitSpeed/(a+3) · t
  k   = axisRatio(a)     // 0.6 → 0.82 → 0.95
  e   = (a·cos φ, a·k·sin φ)
  θ   = a · twist + pattern
  pos = (rotate(θ)·e).x, y, (rotate(θ)·e).y

  color = bulge ? mix(orange, cream, hashU(id,6))
                : blue 25% | white 45% | orange 30%
  lum   = 0.3 + hashU(id,7)^8 · 25
  size  = clamp(220·H/1080 / |pos − cam|
                · (0.7 + 0.15·lum), 1, 6)
  gl_PointSize = size
  vColor = color · lum · 0.09 / size
// fragment, per covered pixel
d   = 2·gl_PointCoord − 1
out = vColor · e^(−3|d|²)

// draw 3: gas, 6,000 quads, additive
for each vertex v:
  id     = v/6 + 7,000,000
  corner = quadCorner(v % 6)
  a      = mix(15, 95, hashU(id,2))
  center = orbitPosition(a, hashU(id,5)·2π,
                         boxMuller · 0.9)
  size   = mix(2.5, 6, hashU(id,6))
  pos    = center + (right·corner.x + up·corner.y)·size
  color  = (70% pink | 30% blue) · 0.06 · 0.6
           · smoothstep(2, 12, |center − cam|)
// fragment
out = color · e^(−3|corner|²)

// draw 4: dust, 14,000 quads
// blend (ZERO, ONE_MINUS_SRC_COLOR)
  id     = v/6 + 9,000,000
  a      = mix(10, 95, hashU(id,2))
  center = orbitPosition(a, phase − 0.12,
                         boxMuller · 0.45)
  size   = mix(3, 8, hashU(id,6))
  color  = (0.10, 0.13, 0.18) · 0.25 · nearFade
  dst    = dst · (1 − color · e^(−3|corner|²))

// engine post chain
bloom → chromatic aberration → ACES
→ vignette, grain → sRGB

The real source

These are the files as they are in the demo's repository. The browser version above uses the same code.

Glossary

Vertex shader

A small program the GPU runs once per vertex; it decides where the vertex lands on screen.

Fragment shader

A program run once per covered pixel (fragment); it decides the pixel's color.

Uniform

A value set by the CPU before a draw and shared by every vertex and pixel of that draw (time, twist, camera…).

gl_VertexID

The index of the vertex being processed. With no vertex data at all, it is enough to build geometry procedurally.

Point sprite

A point drawn as a small screen-aligned square of gl_PointSize pixels; gl_PointCoord tells the fragment shader where inside it a pixel lies.

Billboard

A flat quad turned to always face the camera, used to draw soft round things such as clouds.

Hash

A deterministic scrambling function: same input, same "random-looking" output. Used instead of stored random numbers.

PCG

Permuted congruential generator (O'Neill, 2014): a multiply-add step followed by a bit permutation. One step of it makes a fast, high-quality integer hash.

Floating-point precision

A 32-bit float keeps 24 significant bits; the larger the number, the coarser the steps between representable values.

Probability density / CDF

The density says how likely each value is; the cumulative distribution F(x) is the probability of a value below x.

Inverse transform sampling

Turning a uniform random number u into a value with a chosen distribution by solving F(x) = u.

Exponential disc

A disc whose density falls by the same factor for every step of the scale length outward (here 28 units).

Bulge

The thick, round concentration of mostly old stars at the center of a spiral galaxy.

Gaussian (normal) distribution

The bell curve; σ (sigma) is its width — about two thirds of values fall within ±σ.

Box–Muller transform

A formula turning two uniform random numbers into a normally distributed one (Box and Muller, 1958).

Axis ratio

An ellipse's minor axis divided by its major axis: 1 is a circle, small values are long and thin.

Rotation curve

Orbital speed as a function of distance from a galaxy's center; "flat" when the speed stays nearly constant outward.

Differential rotation

Inner parts turning faster (in angle) than outer parts, as in a galaxy or a hurricane.

Winding problem

Arms made of fixed groups of stars would be wound up tight by differential rotation within a few orbits.

Density wave

A pattern of higher density that moves through a medium while the material itself flows through it — like a traffic jam.

Additive blending

Adding a new color to the existing pixel instead of replacing it — how light combines.

Multiplicative blending

Multiplying the existing pixel by a factor, here 1 − source: how an absorbing layer dims what is behind it.

Extinction and reddening

Dust dims starlight, and dims blue light more than red, so stars seen through dust look redder.

Beer–Lambert law

Light through an absorbing medium falls off exponentially with the amount of material it crosses.

Smoothstep

smoothstep(e0, e1, x): 0 before e0, 1 after e1, an S-shaped ramp in between.

Ray–plane intersection

Where a line ro + t·rd meets a plane; for the plane y = 0, t = −ro.y / rd.y.

HDR

High dynamic range: keeping brightness values above 1.0 until a final tone-mapping step.

Tone mapping (ACES)

A curve compressing HDR brightness into the screen's range while keeping highlights pleasant.

Bloom

Letting very bright areas glow into their neighborhood, like light scattering in a lens or eye.

Part of the Space Travel documentation. Direction & design: blacklizard · Code: Claude. The figures on this page are computed live in your browser; nothing is pre-rendered except the header image.

← 12 · Wormhole  ·  All scenes  ·  14 · Black hole →