Simulating Aurora Borealis Ribbons in GLSL

This article explains how to simulate the fluid, waving curtains of an aurora borealis in GLSL using sinusoidal field displacement. By feeding continuous coordinate space through multi-frequency sine and cosine functions—a technique known as domain warping—shaders perturb flat planes or raymarched volumes into intricate folds. The following sections detail the mathematical principles of sinusoidal displacement, how to layer waves for organic motion, the construction of the ribbon geometry, and the final color mapping.

Mathematical Basis of Sinusoidal Displacement

The fundamental visual characteristic of an aurora is its vertical, folded "curtain" shape drifting across the sky. In GLSL, this curtain begins as an implicit 2D or 3D coordinate field. Instead of manipulating individual mesh vertices, sinusoidal field displacement alters the coordinates before sampling a density function.

In its simplest 2D form, a point \((x, z)\) along the horizontal ground plane is displaced along one axis using a sine function of the other axis and time:

vec2 displace(vec2 p, float time) {
    p.x += sin(p.z * 1.5 + time * 0.8) * 0.4;
    return p;
}

By substituting the displaced coordinate into a distance function, straight lines or boundaries warp into smooth, wave-like curves.

Layering Frequencies and Domain Warping

A single sine wave produces an unnatural, uniform oscillation. Realistic auroral motion requires chaotic, organic variation created by layering multiple sine and cosine waves at different frequencies, amplitudes, and speeds, mirroring Fractional Brownian Motion (fBM):

vec2 auroraDisplacement(vec2 p, float time) {
    float wave1 = sin(p.y * 2.0 + time * 0.5) * 0.3;
    float wave2 = cos(p.x * 3.5 - time * 0.7) * 0.15;
    float wave3 = sin((p.x + p.y) * 5.0 + time * 1.2) * 0.05;
    
    return p + vec2(wave1 + wave3, wave2 + wave3);
}

To create the tight folds visible in high-intensity displays, shaders apply domain warping: displacing coordinates with a function that itself samples displaced coordinates. This feedback loop folds the coordinate space over itself, producing sharp crests and stretched troughs reminiscent of magnetic flux lines.

Defining Ribbon Density and Vertical Strands

Once the domain is displaced, the shader calculates ribbon density. The core of a ribbon is typically defined by calculating the distance of the displaced coordinates to an arbitrary centerline:

float ribbonPath = abs(displacedPos.x);
float ribbonMask = exp(-ribbonPath * 8.0); // Sharp falloff from the center

Aurora borealis curtains also feature distinctive vertical striations (rays) caused by charged particles aligning with the Earth's magnetic field lines. This is achieved by modulating the ribbon mask with a high-frequency horizontal wave applied along the vertical axis (\(y\)):

float striations = sin(displacedPos.x * 40.0 + displacedPos.y * 10.0) * 0.5 + 0.5;
float finalDensity = ribbonMask * mix(0.6, 1.0, striations);

Altitude Shaping and Volumetric Raymarching

Auroras exist within specific altitude bands, fading out both below and above their peak excitation layer. A vertical envelope restricts the ribbon's intensity:

float altitudeFade = smoothstep(0.0, 0.3, uv.y) * smoothstep(1.0, 0.4, uv.y);
float auroraIntensity = finalDensity * altitudeFade;

In 3D implementations, this density function is sampled along a camera ray using raymarching. The ray steps through the bounding volume, evaluating the displaced coordinates at each step, accumulating density and applying lighting via additive blending:

vec4 colorAccum = vec4(0.0);
for (int i = 0; i < 32; i++) {
    vec3 p = rayOrigin + rayDir * stepSize * float(i);
    vec2 displacedXZ = auroraDisplacement(p.xz, u_time);
    float density = calculateDensity(vec3(displacedXZ.x, p.y, displacedXZ.y));
    
    vec4 stepColor = mapAuroraColor(p.y) * density * stepSize;
    colorAccum += stepColor;
}

Color Mapping

The emission spectrum of an aurora depends on altitude and gas composition. At lower altitudes (100–150 km), molecular nitrogen and atomic oxygen produce bright green emissions. At higher altitudes (above 200 km), low-density atomic oxygen produces deep reds and purples.

In the shader, the altitude coordinate (or the vertical texture coordinate) maps directly to this gradient:

vec4 mapAuroraColor(float altitude) {
    vec4 baseGreen = vec4(0.1, 0.9, 0.4, 1.0);
    vec4 topPurple = vec4(0.6, 0.1, 0.8, 1.0);
    return mix(baseGreen, topPurple, smoothstep(0.2, 0.8, altitude));
}

Multiplying this gradient by the displaced density field completes the simulation, delivering continuous, undulating ribbon oscillations entirely computed within the GPU's fragment pipeline.