Rayleigh Scattering Phase Function in GLSL

This article explains how the Rayleigh scattering angular phase function is evaluated and implemented within a GLSL sky shader. Rayleigh scattering governs the interaction of light with air molecules smaller than the light's wavelength, dictating the characteristic blue of the daytime sky and the warm colors of sunrise and sunset. Evaluating its phase function in real-time shaders involves calculating the angle between the viewer's gaze and the incident sunlight, then applying the mathematically normalized Rayleigh distribution formula to weight the scattered light intensity.

The Rayleigh Phase Function Formula

The phase function describes the angular distribution of scattered light relative to the incident light direction. For Rayleigh scattering, the phase function \(P_R(\theta)\) is defined analytically as:

\[P_R(\theta) = \frac{3}{16\pi} (1 + \cos^2\theta)\]

Where:

• \(\theta\) is the scattering angle between the view ray direction and the sun direction.
• \(\cos\theta\) represents the directional alignment.
• \(\frac{3}{16\pi} \approx 0.0596831\) is the normalization factor ensuring the integral of the phase function over a unit sphere equals 1.

Because \(\cos^2\theta\) is symmetric, Rayleigh scattering produces equal intensity in both the forward (\(\theta = 0\)) and backward (\(\theta = \pi\)) directions, with minimum scattering occurring at \(90^\circ\) (\(\theta = \pi/2\)).

Evaluating the Cosine of the Scattering Angle

In a GLSL shader, you compute \(\cos\theta\) directly using the dot product of two normalized unit vectors: the ray direction vector (rayDir) pointing away from the camera into the sky, and the light direction vector (sunDir) pointing toward the sun.

float cosTheta = dot(rayDir, sunDir);

Using the dot product eliminates the need for expensive inverse trigonometric calls like acos().

GLSL Implementation

The calculation can be encapsulated into a standalone GLSL function:

#define PI 3.14159265359

float rayleighPhase(float cosTheta) {
    const float k = 3.0 / (16.0 * PI);
    return k * (1.0 + cosTheta * cosTheta);
}

Integration into the Atmospheric Scattering Pipeline

During raymarching or analytical sky evaluation, total in-scattered light is split into Rayleigh (molecular) and Mie (aerosol) components. The phase function is applied outside the primary transmittance integration loop because it depends solely on the scattering angle, which remains constant along a straight ray.

// Directions must be normalized
vec3 V = normalize(vWorldRayDirection);
vec3 L = normalize(uSunDirection);

// Compute cosine of the scattering angle
float cosTheta = dot(V, L);

// Evaluate phase functions
float phaseR = rayleighPhase(cosTheta);
float phaseM = miePhase(cosTheta, 0.76); // Using Henyey-Greenstein for Mie

// Accumulate optical depth along the ray (via raymarching or lookup table)
vec3 rayleighScattering = ...; 
vec3 mieScattering = ...;

// Apply phase functions to modulate accumulated in-scattering
vec3 finalSkyColor = rayleighScattering * phaseR + mieScattering * phaseM;

Precomputing the scalar factor 3.0 / (16.0 * PI) into a const float avoids redundant operations, maximizing per-pixel execution speed.