Ok! Here are the shaders.

Vertex shader:

```
/*
Shader for rendering a planet atmosphere.
*/
// Interpolated position in eye coordinates
varying vec4 eyePosition;
void main()
{
// Compute position in eye coordinates
eyePosition = gl_ModelViewMatrix * gl_Vertex;
// Obtain projected position
gl_Position = ftransform();
}
```

The fragment shader:

```
/*
Shader for rendering a planet atmosphere.
All lights are assumed to be positional. Maximum number of lights taken into account is 2
Lights are treated as directional because they are assumed to be far away, so the light position
is treated as a direction since they are in camera space. ( But the position must be normalized)
Uniforms:
numLights - The number of (positional) lights activated in OpenGL. GL_LIGHT0 through GL_LIGHT(numLights-1)
atmosphereAltitude - Altitude at wich atmosphere density becomes very small
atmosphericPressure - Atmosphere pressure at 0 altitude
atmosphericScale - Atmospheric pressure constant
atmosphereColor - atmosphere diffused light color
planetRadius - Planet radius
planetEyePosition - Planet position in eye coordinates
*/
const int NUM_SAMPLES_ATMOSPHERE = 5;
const float E = 2.718281828;
//const float ATMOSPHERE_DENSITY = 0.002;
const float ATMOSPHERE_DENSITY = 0.01;
const vec3 REDDISH_ATMOSPHERE = vec3( 0.7, 0.2, 0.1 );
// Number of active lights (set by application)
uniform int numLights;
// Atmosphere altitude
uniform float atmosphereAltitude;
// Atmospheric pressure at surface level (units = terrestrial pressures)
uniform float atmosphericPressure;
// atmosphericScale is -g / H (g: surface gravity in earths gravities, H = atmosphere height scale)
uniform float atmosphericScale;
// Atmosphere color
uniform vec3 atmosphereColor;
// Planet radius
uniform float planetRadius;
// Planet position in eye coordinates
uniform vec3 planetEyePosition;
// Interpolated position of vertex in eye coordinates
varying vec4 eyePosition;
// Computes ray-sphere intersection.
// rayDir must be unitary.
// Returns the ray parameters for first and second intersection in that order.
// Only positive intersections (parameter > 0) are computed.
// If the first is < 0, there is no intersection.
// Else, if the second is < 0, there is only one intersection.
vec2 raySphereIntersection( vec3 rayOrigin, vec3 rayDir, vec3 spherePos, float sphereRadius ) {
vec3 centerToOrigin = rayOrigin - spherePos;
float b = 2.0 * ( dot( rayDir, centerToOrigin ) );
float c = dot( centerToOrigin, centerToOrigin ) - sphereRadius * sphereRadius;
float d = b * b - 4.0 * c;
if ( d < 0.0 ) {
// No interesections
return vec2( -1.0, -1.0 );
}
d = sqrt( d );
float t1 = 0.5 * ( -b - d );
float t2 = 0.5 * ( -b + d );
if ( t1 < 0.0 ) {
return vec2( t2, -1.0 );
}
return vec2( t1, t2 );
}
void main()
{
vec3 rayDir = normalize( eyePosition.xyz );
vec2 tAtm = raySphereIntersection( vec3( 0.0 ), rayDir, planetEyePosition, planetRadius + atmosphereAltitude );
if ( tAtm.x < 0.0 ) {
// Ray lies outside of planet and atmosphere
discard;
}
vec2 tPlanet = raySphereIntersection( vec3( 0.0 ), rayDir, planetEyePosition, planetRadius );
// ti and t2 are the ray parameters for initial and end points of the ray inside the atmosphere
float t1 = 0.0;
float t2 = 0.0;
float rayTouchesOnlySky = 0.0;
if ( tPlanet.x < 0.0 ) {
// The ray doesn't intersect the planet.
if ( tAtm.y > 0.0 ) {
// The ray enters and exits the atmosphere without touching the planet
t1 = tAtm.x;
t2 = tAtm.y;
}
else {
// The camera is inside the atmosphere and the ray exits it without touching the planet.
t1 = 0.0;
t2 = tAtm.x;
rayTouchesOnlySky = 1.0;
}
}
else {
// The ray touches the planet.
if ( tAtm.y > 0.0 ) {
// The ray enters the atmosphere and then touches the planet.
t1 = tAtm.x;
t2 = tPlanet.x;
}
else {
// The camera is inside the atmosphere and the ray touches the planet.
t1 = 0.0;
t2 = tPlanet.x;
}
}
float factorHeight = 1.0 / ( atmosphereAltitude );
float planetEyePositionLength = length( planetEyePosition );
float cameraHeight = planetEyePositionLength - planetRadius;
//float density;
float densityIlluminated;
float alpha;
// Sampling of atmosphere
densityIlluminated = 0.0;
float dt = ( t2 - t1 ) / float( NUM_SAMPLES_ATMOSPHERE );
vec3 sampleIncrem = rayDir * dt;
vec3 samplePos = rayDir * t1 + sampleIncrem * 0.5;
for ( int i = 0; i < NUM_SAMPLES_ATMOSPHERE; i++ ) {
vec3 planet2sample = samplePos - planetEyePosition;
// Compute air density
float sampleHeight = max( 0.0, length( planet2sample ) - planetRadius );
// Compute illuminated air density
vec3 normal = normalize( planet2sample );
float dotLightNormal = dot( normal, normalize( gl_LightSource[ 0 ].position.xyz /*- samplePos*/ ) );
if ( dotLightNormal > -0.0 ) {
float d = exp( atmosphericScale * sampleHeight ) * dt;
densityIlluminated += d * dotLightNormal;
}
// Increment sample position
samplePos += sampleIncrem;
}
densityIlluminated *= ATMOSPHERE_DENSITY * atmosphericPressure;
alpha = densityIlluminated;
vec3 lightDir = normalize( gl_LightSource[ 0 ].position.xyz );
float dotLightNormal = max( 0.0, - dot( planetEyePosition, lightDir ) / planetEyePositionLength );
if ( cameraHeight < atmosphereAltitude ) {
// Camera is inside the atmosphere, add light diffusion
float d = exp( atmosphericScale * cameraHeight );
densityIlluminated = min( 1.0, max( densityIlluminated, rayTouchesOnlySky * sqrt( dotLightNormal ) * d ) );
alpha = densityIlluminated;
}
// Computes 'redding' of color due to angle to the sun and planet normal
float dotLightRayRaw = dot( rayDir, lightDir );
float dotLightRay = max( 0.0, dotLightRayRaw );
float dotRayNormal = max( 0.0, - dot( planetEyePosition, rayDir ) / planetEyePositionLength );
float redding = pow( dotLightRay, 50.0 ) + ( 0.5 + 0.5 * dotLightRayRaw ) * pow( 1.0 - dotRayNormal, 2.0 );
// Computes sun color
vec3 sunColor = mix( gl_LightSource[ 0 ].diffuse.xyz, REDDISH_ATMOSPHERE, 1.0 - dotLightNormal );
if ( numLights > 1 ) {
// Only two suns contribute diffused light
lightDir = normalize( gl_LightSource[ 1 ].position.xyz );
dotLightNormal = max( 0.0, - dot( planetEyePosition, lightDir ) / planetEyePositionLength );
dotLightRayRaw = dot( rayDir, lightDir );
dotLightRay = max( 0.0, dotLightRayRaw );
dotRayNormal = max( 0.0, - dot( planetEyePosition, rayDir ) / planetEyePositionLength );
float redding1 = pow( dotLightRay, 50.0 ) + ( 0.5 + 0.5 * dotLightRayRaw ) * pow( 1.0 - dotRayNormal, 2.0 );
vec3 sunColor1 = mix( gl_LightSource[ 1 ].diffuse.xyz, REDDISH_ATMOSPHERE, 1.0 - dotLightNormal );
redding = 0.5 * ( redding + redding1);
sunColor = mix( sunColor, sunColor1, redding );
}
// Store color with alpha
gl_FragColor = vec4( mix( atmosphereColor, sunColor, redding ), alpha );
}
```

Comments:

I use GLSL 1.2, but I’m planning to convert the engine to GL 3.3/GLSL 1.5

I use only 2 lights for performance, though i can have 4 suns in a solar system, only 2 of them (close binary) affect a planet.

I use a camera-centered world due to numeric precision. So camera is at origin (0,0,0) and planet is at planetEyePosition. So rays go from origin to eyePosition, wich is the fragment position on the quad.

The units in planet drawing are kilometers, the units for object drawing are meters.

Have fun!

Btw you can visit my blog and see some pictures and a video (sorry, no pictures of the atmosphere. I’m going to make a post in the blog with pictures for you)

The blog url is below.