Atmospheric Scattering Shader

Hello All!

So, i gave it a go at trying to adapt the Rayleigh Scattering vs. Mie Scattering :

http://http.developer.nvidia.com/GPUGems2/gpugems2_chapter16.html

I used bits and pieces of code present around the web and came up with the following:

shader OCC01(
    vector Eye = vector(0,0,1),
    vector Sun = vector(0,0,1),
    output closure color CL = holdout()
){
    // math const
    float PI = 3.14159265359;
    float DEG_TO_RAD = PI / 360.0;
    float MAX = 10000.0;
    // scatter const
    float K_R = 0.166;
    float K_M = 0.0025;
    float E = 20;       //14.3                         // light intensity
    vector C_R = vector( 0.3, 0.7, 1.0 );     // 1 / wavelength ^ 4
    float G_M = -0.85;    // -0.85                // Mie g
    float R = 1.0; // 1.0
    float R_INNER = 0.7;
    float SCALE_H = 4.0 / ( R - R_INNER );
    float SCALE_L = 1.0 / ( R - R_INNER );
    int NUM_OUT_SCATTER = 10;
    float FNUM_OUT_SCATTER = 10.0;
    int NUM_IN_SCATTER = 10;
    float FNUM_IN_SCATTER = 10.0;
    /* Begin */  
    // Mie
    // g : ( -0.75, -0.999 )
    //      3 * ( 1 - g^2 )               1 + c^2
    // F = ----------------- * -------------------------------
    //      2 * ( 2 + g^2 )     ( 1 + g^2 - 2 * g * c )^(3/2)
    float phase_mie( float g, float c, float cc ) {
        float gg = g * g;
        
        float a = ( 1.0 - gg ) * ( 1.0 + cc );

        float b = 1.0 + gg - 2.0 * g * c;
        b *= sqrt( b );
        b *= 2.0 + gg;    
        
        return 2.5 * a / b;
    }
    //
    // Reyleigh
    // g : 0
    // F = 3/4 * ( 1 + c^2 )
    float phase_reyleigh( float cc ) {
        return 0.75 * ( 1.0 + cc );
    }
    
    float density( vector p ){
       return exp( -( length( p ) - R_INNER ) * SCALE_H );
    }
    
    float optic( vector p, vector q ) {
       vector step = ( q - p ) / FNUM_OUT_SCATTER;
       vector v1 = p + step * 0.5;
       float sum = 0.0;
       for ( int i = 0; i < NUM_OUT_SCATTER; i++ ) {
          sum += density( v1 );
          v1 += step;
       }
       sum *= length( step ) * SCALE_L;
       return sum;
    }  
    //
    vector o = normalize(Eye); // eye
    vector dir = vector(N); // dir
    vector l = normalize(Sun); // sun light dir
    vector e = cross(o,dir);
    // 
    float len = ( e[1] - e[0] ) / FNUM_IN_SCATTER;
    vector step01 = dir * len;
    vector p = o + dir * e[0];
    vector v1 = p + dir * ( len * 0.5 );
    vector sum = vector( 0,0,0 );
    
    for ( int i = 0; i < NUM_IN_SCATTER; i++ ) {
        vector f = cross(I,-dir);
        vector u1 = v1 + l * f[1];    
        float n = ( optic( p, v1 ) + optic( v1, u1 ) ) * ( PI * 4.0 );    
        sum += density( v1 ) * exp( -n * ( K_R * C_R + K_M ) );
        v1 += step01;
       
    }   
    sum *= len * SCALE_L;    
    float c  = dot( dir, -l );
    float cc = c * c;
    vector Fact = sum * K_R * C_R * phase_reyleigh( cc ) + K_M * phase_mie( G_M, c, cc ) * E;
    CL =  emission()*1.5*color(Fact);
    return;    
}

This code produced this result in a sphere surface placed at the center of the scene:


I hope you can improve on it and have fun!

Attachments