Earth Shadow Terminator Formulation in GLSL

This article details the mathematical formulation and GLSL implementation of the Earth's shadow terminator across atmospheric layers. It covers the analytical geometry of the planetary shadow cone, the derivation of the dip angle at varying altitudes, atmospheric transmission extinction, and an optimized GLSL routine for physically based volumetric atmospheric rendering.

Geometric Shadow Formulation

Consider the Earth as a sphere of radius \(R_E\) centered at the coordinate origin \(\mathbf{O} = (0, 0, 0)\). An atmospheric shell surrounds the Earth up to a maximum radius \(R_A = R_E + H_{atm}\), where \(H_{atm}\) is the atmospheric thickness.

Let:

  • \(\mathbf{p}\) be a sample point within the atmosphere, with radial distance \(r = \|\mathbf{p}\|\).
  • \(\mathbf{s}\) be the normalized unit vector pointing directly toward the Sun.

Ray-Sphere Occlusion Test

A point \(\mathbf{p}\) is occluded by the planetary core if a ray cast from \(\mathbf{p}\) along the light direction \(\mathbf{s}\), parameterized as \(\mathbf{r}(t) = \mathbf{p} + t\mathbf{s}\) for \(t > 0\), intersects the sphere of radius \(R_E\).

Substituting \(\mathbf{r}(t)\) into the sphere equation \(\|\mathbf{r}(t)\|^2 = R_E^2\) yields the quadratic equation:

\[t^2 + 2(\mathbf{p} \cdot \mathbf{s})t + (\|\mathbf{p}\|^2 - R_E^2) = 0\]

Let \(b = \mathbf{p} \cdot \mathbf{s}\) and \(c = \|\mathbf{p}\|^2 - R_E^2\). The discriminant is:

\[\Delta = b^2 - c = (\mathbf{p} \cdot \mathbf{s})^2 - \|\mathbf{p}\|^2 + R_E^2\]

The point \(\mathbf{p}\) lies in the geometric umbra of the planet if:

  1. \(\Delta \ge 0\) (the infinite line intersects the sphere).
  2. The closest approach along the ray occurs behind the observer relative to the Sun: \(b = \mathbf{p} \cdot \mathbf{s} < 0\).
  3. The point is outside the Earth: \(c > 0\) (i.e., \(\|\mathbf{p}\| \ge R_E\)).

The distance to the closest point on the ray to the Earth's center is given by:

\[d_\perp = \sqrt{\|\mathbf{p}\|^2 - (\mathbf{p} \cdot \mathbf{s})^2}\]

Occlusion occurs when \(d_\perp \le R_E\) and \(\mathbf{p} \cdot \mathbf{s} < 0\).

Altitude-Dependent Terminator Zenith Angle

At altitude \(h = r - R_E\), the geometric horizon dips. The critical solar zenith angle \(\theta_{term}\) at which the Sun touches the geometric horizon is given by:

\[\cos \theta_{term} = -\frac{\sqrt{r^2 - R_E^2}}{r}\]

Points with \(\cos \theta_s = \frac{\mathbf{p} \cdot \mathbf{s}}{r} < \cos \theta_{term}\) lie inside the geometric shadow.

Atmospheric Penumbra and Optical Extinction

In a real atmosphere, the transition across the terminator is softened by two factors: the non-zero angular radius of the Sun (\(\alpha_\odot \approx 0.00465\text{ rad}\)) and the exponential density decay of atmospheric constituents (Rayleigh scattering and ozone absorption).

Along a ray passing at a minimum tangent height \(h_{min}\) above the surface, the density profile follows:

\[\rho(h) = \rho_0 \exp\left(-\frac{h}{H_0}\right)\]

where \(H_0\) is the scale height (approximately \(8.0\text{ km}\) for Rayleigh scattering and \(1.2\text{ km}\) for Mie aerosols).

The effective optical depth \(\tau(\mathbf{p}, \mathbf{s})\) along the solar ray can be approximated near the terminator by the grazing transmittance function:

\[T(\mathbf{p}, \mathbf{s}) = \exp\left(-\beta_{ext} \int_0^\infty \rho(\mathbf{p} + t\mathbf{s}) \, dt\right)\]

For real-time evaluation, the shadow visibility factor \(V(\mathbf{p}, \mathbf{s}) \in [0, 1]\) accounts for both the geometric occultation of the solar disk and density-based transmission:

\[V(\mathbf{p}, \mathbf{s}) = \text{clamp}\left(\frac{d_\perp - R_E}{R_E \tan \alpha_\odot}, 0.0, 1.0\right) \times T(\mathbf{p}, \mathbf{s})\]

GLSL Implementation

The following GLSL function computes the planetary shadow factor across atmospheric layers, returning a smooth visibility factor between 0.0 (full shadow) and 1.0 (full sunlight), including twilight penumbra.

// Earth parameters (in meters)
const float EARTH_RADIUS = 6371000.0;
const float ATMOS_RADIUS = 6471000.0;
const float SUN_ANGULAR_RADIUS = 0.00465; // ~0.267 degrees in radians

/**
 * Calculates the solar visibility factor at point 'p' considering Earth's shadow.
 * @param p World-space position relative to Earth center (meters)
 * @param s Normalized unit vector pointing to the Sun
 * @return Scalar visibility factor in range [0.0, 1.0]
 */
float getEarthShadowTerminator(vec3 p, vec3 s) {
    float r2 = dot(p, p);
    float pDotS = dot(p, s);
    
    // If the point is facing the Sun, it cannot be in the planet's shadow
    if (pDotS >= 0.0) {
        return 1.0;
    }
    
    // Perpendicular distance squared from Earth center to the light ray
    float d2 = r2 - (pDotS * pDotS);
    float rEarth2 = EARTH_RADIUS * EARTH_RADIUS;
    
    // Ray misses the Earth entirely
    if (d2 > rEarth2) {
        return 1.0;
    }
    
    float d = sqrt(max(0.0, d2));
    
    // Penumbra width at the distance from the occulting limb to point p
    float limbDist = -pDotS;
    float penumbraRadius = limbDist * tan(SUN_ANGULAR_RADIUS);
    
    // Distance from the ray to the Earth's surface at the closest approach
    float distToSurface = d - EARTH_RADIUS;
    
    // Smooth transition between total shadow (umbra) and full light
    float shadowFactor = smoothstep(-penumbraRadius, penumbraRadius, distToSurface);
    
    // Atmospheric optical extinction near the horizon limb
    // Adds exponential absorption for rays passing through dense lower air
    float tangentAltitude = max(0.0, d - EARTH_RADIUS);
    const float RAYLEIGH_SCALE_HEIGHT = 8000.0;
    float airMassFactor = exp(-tangentAltitude / RAYLEIGH_SCALE_HEIGHT);
    float atmosphericAttenuation = exp(-airMassFactor * 4.0);
    
    return shadowFactor * atmosphericAttenuation;
}

This formulation integrates directly into single-scattering and multiple-scattering raymarch loops within atmospheric pixel and compute shaders.