GLSL Ice Crystal Halo Refraction Algorithm

This article provides an overview of the algorithms used to evaluate frost crystal and ice halo refraction rings around light sources in GLSL shaders. Simulating these atmospheric optics requires modeling how light interacts with hexagonal ice crystals, specifically focusing on the minimum deviation angle, chromatic dispersion, and angular intensity falloff. The following sections break down the underlying mathematical model and provide a direct GLSL implementation suitable for real-time rendering.

The Minimum Deviation Algorithm

Atmospheric halos, such as the common 22° and 46° halos, are produced by light refracting through hexagonal prism ice crystals suspended in the air. In real-time computer graphics, performing brute-force Monte Carlo ray tracing through millions of randomly oriented 3D crystals is computationally prohibitive. Instead, GLSL implementations utilize an analytical phase function based on the Prism Minimum Deviation Law.

When light passes through a crystal prism with an apex angle \(\alpha\) and a refractive index \(n\), the angle of minimum deviation \(\theta_{min}\) represents the threshold below which no light can be refracted. It is determined by:

\[\theta_{min} = 2 \arcsin\left(n \sin\left(\frac{\alpha}{2}\right)\right) - \alpha\]

  • 22° Halo: Light passes through alternating prism faces with an apex angle of \(\alpha = 60^\circ\) (\(\pi/3\) radians). For ice (\(n \approx 1.309\)), \(\theta_{min} \approx 21.84^\circ\).
  • 46° Halo: Light passes between a prism face and a basal face with an apex angle of \(\alpha = 90^\circ\) (\(\pi/2\) radians), yielding \(\theta_{min} \approx 45.7^\circ\).

Because light rays pile up at this minimum deviation threshold before dropping off gradually at wider angles, the intensity profile features a sharp inner boundary and a diffused outer tail.

Chromatic Dispersion Modeling

The characteristic red inner edge and blue/white outer fringe of an ice halo occur because the refractive index of ice varies with wavelength. This dispersion can be modeled using Cauchy's equation:

\[n(\lambda) = A + \frac{B}{\lambda^2}\]

In a fragment shader, evaluating three distinct wavelengths (red: ~650 nm, green: ~532 nm, blue: ~450 nm) allows the calculation of separate minimum deviation angles for each color channel:

  • \(n_{red} \approx 1.307 \rightarrow \theta_{min} \approx 21.6^\circ\)
  • \(n_{green} \approx 1.310 \rightarrow \theta_{min} \approx 21.9^\circ\)
  • \(n_{blue} \approx 1.317 \rightarrow \theta_{min} \approx 22.5^\circ\)

GLSL Implementation

The following fragment shader function calculates the scattering intensity of an ice halo around a point or directional light source. It calculates the angular separation between the view vector and the light vector, evaluates the deviation profile across the RGB spectrum, and produces the ring effect.

// Computes halo intensity for a given scattering angle and refractive index
float evaluateHaloComponent(float cosTheta, float n, float apexAngle) {
    // Calculate minimum deviation angle for the given refractive index
    float sinHalfAlpha = sin(apexAngle * 0.5);
    float sinVal = n * sinHalfAlpha;
    
    // Total internal reflection check
    if (sinVal > 1.0) return 0.0;
    
    float thetaMin = 2.0 * asin(sinVal) - apexAngle;
    float currentAngle = acos(clamp(cosTheta, -1.0, 1.0));
    
    // Light cannot refract below the minimum deviation angle
    if (currentAngle < thetaMin) {
        return 0.0;
    }
    
    // Angular distance from minimum deviation peak
    float deltaTheta = currentAngle - thetaMin;
    
    // Intensity spikes at thetaMin and decays asymptotically (1 / sqrt(deltaTheta))
    // Softened with an epsilon to prevent division by zero
    float peakIntensity = 1.0 / sqrt(deltaTheta + 0.001);
    
    // Exponential falloff for larger angles
    float falloff = exp(-deltaTheta * 12.0);
    
    return peakIntensity * falloff;
}

// Main evaluation function for the 22-degree ice halo
vec3 evaluateIceHalo(vec3 viewDir, vec3 lightDir) {
    // Cosine of the angle between view and light directions
    float cosTheta = dot(normalize(viewDir), normalize(lightDir));
    
    // Apex angle for 22-degree halo (60 degrees in radians)
    const float alpha = 1.04719755; 
    
    // Wavelength-dependent refractive indices for ice (RGB)
    const vec3 nIce = vec3(1.307, 1.310, 1.317);
    
    vec3 haloColor;
    haloColor.r = evaluateHaloComponent(cosTheta, nIce.r, alpha);
    haloColor.g = evaluateHaloComponent(cosTheta, nIce.g, alpha);
    haloColor.b = evaluateHaloComponent(cosTheta, nIce.b, alpha);
    
    // Scale intensity
    return haloColor * 0.15;
}

Integration into Atmospheric Scattering

To integrate this algorithm into a sky or lens shader:

  1. Normalize the interpolated view ray and the vector pointing toward the light source.
  2. Call evaluateIceHalo(viewDir, lightDir).
  3. Modulate the returned color by the atmospheric optical depth (crystal concentration) and light source radiance.
  4. Add the resulting halo value directly to the ambient sky color and forward Mie scattering contributions before tone mapping.