Venus Twilight Belt Formula in GLSL Shaders

This article explains how to calculate the twilight belt illumination of Venus in atmospheric GLSL shaders. Due to its thick carbon dioxide atmosphere and reflective sulfuric acid cloud layers, Venus exhibits an extended twilight zone that wraps around its dark limb beyond the geometric terminator. Real-time rendering engines simulate this phenomenon by modifying standard diffuse lighting through extended terminator offsets, exponential twilight falloff functions, and Henyey-Greenstein forward-scattering phase approximations.

The Physical Mechanism

Standard Lambertian shading evaluates light using \(I = \max(\mathbf{N} \cdot \mathbf{L}, 0.0)\), cutting off abruptly at the geometric terminator where the solar zenith angle (\(\theta\)) equals \(90^\circ\) (\(\mathbf{N} \cdot \mathbf{L} = 0.0\)).

On Venus, intense atmospheric refraction and multiple scattering within the upper aerosol deck extend diffuse light into the nightside up to a solar zenith angle of approximately \(98^\circ\) to \(102^\circ\). This creates the visible twilight belt and crescent cusp extensions.

The Mathematical Formula

To compute this twilight illumination analytically in a fragment shader, the angle of solar depression past the terminator is mapped to an exponential scattering falloff modulated by directional forward scattering.

1. Solar Zenith Angle and Twilight Depression

First, evaluate the cosine of the solar zenith angle:

\[\mu = \mathbf{N} \cdot \mathbf{L}\]

The angle \(\theta\) in radians is:

\[\theta = \arccos(\text{clamp}(\mu, -1.0, 1.0))\]

The twilight depression angle past the geometric terminator (\(\theta_{\text{term}} = \frac{\pi}{2}\)) is:

\[\Delta\theta = \max\left(0.0, \theta - \frac{\pi}{2}\right)\]

2. Twilight Transmittance Function

Atmospheric optical depth causes the twilight belt to decay exponentially as the sun dips lower below the local horizon:

\[T_{\text{twilight}} = \exp\left(-\frac{\Delta\theta}{\sigma}\right)\]

Where:

• \(\Delta\theta\) is the depression angle past the terminator.
• \(\sigma\) is the twilight scale factor (calibrated to Venus's atmospheric scale height and cloud optical depth, typically between \(0.03\) and \(0.08\) radians).

3. Forward Phase Function (Henyey-Greenstein)

Because Venusian twilight is heavily dominated by forward scattering of sunlight through the upper hazes toward the observer, the intensity is scaled by a phase function:

\[P(\phi, g) = \frac{1 - g^2}{4\pi (1 + g^2 - 2g \cos\phi)^{1.5}}\]

Where:

• \(\cos\phi = \mathbf{V} \cdot \mathbf{L}\) (the cosine of the angle between the view vector and the light vector).
• \(g\) is the asymmetry parameter (\(g \approx 0.7\) to \(0.85\) for Venusian sulfuric acid aerosols).

4. Combined Twilight Intensity

The total twilight belt illumination \(I_{\text{twilight}}\) is calculated as:

\[I_{\text{twilight}} = I_{\text{sun}} \cdot T_{\text{twilight}} \cdot P(\mathbf{V} \cdot \mathbf{L}, g) \cdot W(\mu)\]

Where \(W(\mu)\) is a smooth transition weight that blends the daytime diffuse illumination into the twilight belt without producing hard edges.

GLSL Shader Implementation

Below is the GLSL function implementing this formulation:

// Calculates the Venus twilight belt intensity
// N: Normalized surface normal
// L: Normalized light direction (towards Sun)
// V: Normalized view direction (towards camera)
vec3 calculateVenusTwilight(vec3 N, vec3 L, vec3 V, vec3 sunColor) {
    const float PI = 3.14159265359;
    const float HALF_PI = 1.57079632679;
    
    // Venus atmospheric constants
    const float twilightScale = 0.052; // Angular falloff rate (~3 degrees)
    const float maxTwilightAngle = 0.174; // Max twilight depth (~10 degrees)
    const float g = 0.80; // Aerosol forward-scattering asymmetry factor
    
    float NdotL = dot(N, L);
    float VdotL = dot(V, L);
    
    // Compute solar zenith angle
    float zenithAngle = acos(clamp(NdotL, -1.0, 1.0));
    
    // Angular distance past the terminator
    float depression = max(0.0, zenithAngle - HALF_PI);
    
    // Cut off evaluation beyond the active twilight belt limit
    if (depression > maxTwilightAngle) {
        return vec3(0.0);
    }
    
    // Exponential attenuation through atmospheric path
    float attenuation = exp(-depression / twilightScale);
    
    // Henyey-Greenstein phase function for forward haze scattering
    float g2 = g * g;
    float phase = (1.0 - g2) / (4.0 * PI * pow(1.0 + g2 - 2.0 * g * VdotL, 1.5));
    
    // Smooth transition near the terminator boundary
    float blend = smoothstep(HALF_PI + maxTwilightAngle, HALF_PI, zenithAngle);
    
    return sunColor * attenuation * phase * blend;
}

This formula isolates the illumination occurring strictly in the twilight region and can be added directly to the daytime multi-scattering terms to produce a physically representative illumination model of the Venusian atmosphere.