Just found this rather nice Atmospheric Scattering rendering C++ code over at scratchapixel.com and thought I’d do a quick conversion to a GLSL shader as a test for the Timeline software I’m working on. Works rather nicely…

My (none optimised) fragment shader conversion is:

#version 150
#define M_PI 3.1415926535897932384626433832795
uniform float TimeOfDay; // range 0.0 -> 1.0 (0.0 = Midnight, 0.5 = Midday, etc)
const float RADIUS_EARTH = 6360e3;
const float RADIUS_ATMOSPHERE = 6420e3;
const float RAYLEIGH_SCALE_HEIGHT = 7994;
const float MIE_SCALE_HEIGHT = 1200;
const float SUN_INTENSITY = 20;
const float g = 0.76;
const vec3 betaR = vec3( 5.5e-6, 13.0e-6, 22.4e-6 ); // Rayleigh scattering coefficients at sea level
const vec3 betaM = vec3( 21e-6 ); // Mie scattering coefficients at sea level
vec3 sunDirection = vec3( 0, 1, 0 );
const int numSamples = 16;
const int numSamplesLight = 8;
struct Ray
{
vec3 o; //origin
vec3 d; //direction (should always be normalized)
};
struct Sphere
{
vec3 pos; //center of sphere position
float rad; //radius
};
const Sphere SPHERE_EARTH = Sphere( vec3( 0 ), RADIUS_EARTH );
const Sphere SPHERE_ATMOSPHERE = Sphere( vec3( 0 ), RADIUS_ATMOSPHERE );
bool intersect( in Ray ray, in Sphere sphere, out float t0, out float t1 )
{
vec3 oc = ray.o - sphere.pos;
float b = 2.0 * dot(ray.d, oc);
float c = dot(oc, oc) - sphere.rad*sphere.rad;
float disc = b * b - 4.0 * c;
if (disc < 0.0)
return false;
float q;
if (b < 0.0) q = (-b - sqrt(disc))/2.0; else q = (-b + sqrt(disc))/2.0; t0 = q; t1 = c / q; // make sure t0 is smaller than t1 if (t0 > t1) {
// if t0 is bigger than t1 swap them around
float temp = t0;
t0 = t1;
t1 = temp;
}
// if t1 is less than zero, the object is in the ray's negative direction
// and consequently the ray misses the sphere
if (t1 < 0.0)
return false;
if( t0 < 0.0 )
{
t0 = 0;
}
return( true );
}
vec3 computeIncidentLight( in Ray r )
{
float t0, t1;
if( !intersect( r, SPHERE_ATMOSPHERE, t0, t1 ) )
{
return vec3( 1 );
}
float segmentLength = ( t1 - t0 ) / numSamples;
float tCurrent = t0;
vec3 sumR = vec3( 0 );
vec3 sumM = vec3( 0 );
float opticalDepthR = 0;
float opticalDepthM = 0;
float mu = dot( r.d, sunDirection );
float phaseR = 3 / ( 16 * M_PI ) * ( 1 + mu * mu );
float phaseM = 3 / ( 8 * M_PI ) * ( ( 1 - g * g ) * ( 1 + mu * mu ) ) / ( ( 2 + g * g ) * pow( 1 + g * g - 2 * g * mu, 1.5 ) );
for( int i = 0; i < numSamples ; i++ )
{
vec3 samplePosition = r.o + r.d * ( tCurrent + 0.5 * segmentLength );
float height = length( samplePosition ) - RADIUS_EARTH;
// compute optical depth for light
float hr = exp( -height / RAYLEIGH_SCALE_HEIGHT ) * segmentLength;
float hm = exp( -height / MIE_SCALE_HEIGHT ) * segmentLength;
opticalDepthR += hr;
opticalDepthM += hm;
// light optical depth
Ray lightRay = Ray( samplePosition, sunDirection );
float lmin, lmax;
intersect( lightRay, SPHERE_ATMOSPHERE, lmin, lmax );
float segmentLengthLight = lmax / numSamplesLight;
float tCurrentLight = 0;
float opticalDepthLightR = 0;
float opticalDepthLightM = 0;
int j = 0;
for( ; j < numSamplesLight ; j++ )
{
vec3 samplePositionLight = lightRay.o + lightRay.d * ( tCurrentLight + 0.5 * segmentLengthLight );
float heightLight = length( samplePositionLight ) - RADIUS_EARTH;
if( heightLight < 0 )
{
break;
}
opticalDepthLightR += exp( -heightLight / RAYLEIGH_SCALE_HEIGHT ) * segmentLengthLight;
opticalDepthLightM += exp( -heightLight / MIE_SCALE_HEIGHT ) * segmentLengthLight;
tCurrentLight += segmentLengthLight;
}
if( j == numSamplesLight )
{
vec3 tau = betaR * ( opticalDepthR + opticalDepthLightR ) + betaM * 1.1 * ( opticalDepthM + opticalDepthLightM );
vec3 attenuation = exp( -tau );
sumR += hr * attenuation;
sumM += hm * attenuation;
}
tCurrent += segmentLength;
}
return( SUN_INTENSITY * ( sumR * phaseR * betaR + sumM * phaseM * betaM ) );
}
void main()
{
const int width = 512;
const int height = 512;
float a = mod( TimeOfDay - 0.5, 1 ) * 2.0 * M_PI;
sunDirection = normalize( vec3( 0, cos( a ), sin( a ) ) );
float x = 2 * ( gl_FragCoord.x + 0.5 ) / ( width - 1 ) - 1;
float y = 2 * ( gl_FragCoord.y + 0.5 ) / ( height - 1 ) - 1;
float z2 = x * x + y * y;
if( z2 <= 1 )
{
float phi = atan( y, x );
float theta = acos( 1 - z2 );
vec3 dir = vec3( sin( theta ) * cos( phi ), cos( theta ), sin( theta ) * sin( phi ) );
vec3 pos = vec3( 0, RADIUS_EARTH + 1, 0 );
gl_FragColor = vec4( computeIncidentLight( Ray( pos, normalize( dir ) ) ), 1 );
}
else
{
gl_FragColor = vec4( 0 );
}
}