          + vec3( 0.8, 0.9, 1.0 ) * 0.40 * pow( sun_body, ( 1.0 - sun_flare_size ) * 100.0 );
 }
 
+/*
+ * Henyey-Greenstein phase function for scattering
+ *
+ * q: the cosine (dot product) between incident (light to cloud pos) and scattered (cloud to eye pos) ray
+ *
+ */ 
+float hg( float q )
+{
+    float g = 0.5;
+    float t1 = 1.0 / ( 4.0 * PI );
+    float gsq = g * g;
+    float num = 1.0 - gsq;
+    float denom = pow( ( 1.0 + gsq - 2.0 * g * q ), 1.5 );
+
+    return t1 * ( num / denom );
+}
+
 vec4 clouds( vec3 v )
 {
     vec2 ofs = vec2( uTime * cloud_speed * 80.0, uTime * cloud_speed * 60.0 );
-    vec4 acc = vec4( 0, 0, 0, 0 );
+    vec4 acc = vec4( 0, 0, 0, 0 ); // this is the final color value we return
 
     // early exit if we're beneath a certain threshold
     // this doesn't seem to save any frames, though

// Worley noise implementation source:
// https://github.com/Erkaman/glsl-worley
// modified a touch; we only care about F1, so just return that in our 3D worley impl

vec3 permute( vec3 x ) {
    return mod( ( 34.0 * x + 1.0 ) * x, 289.0 );
}

        float height = ( float( i ) * 12.0 + 200.0 - cPosition_World.y ) / v.y;
        vec3 cloudPos = cPosition_World.xyz + height*v + vec3( 831.0, 321.0 + float( i ) * 0.15 - 0.2*ofs.x, 1330.0 + 0.3*ofs.y );
        float density = cloudiness * smoothstep( 0.4, 1.0, fbm( cloudPos * 0.0015 ) );
        /*
         * Beer's law: Energy = e^(-density)
         * Powdered sugar scattering approximation: Energy = 1.0 - e^-(2.0 * density)
         */

        float T = exp( -density ) * ( 1.0 - exp( -2.0 * density ) ) * hg( dot( normalize( cloudPos - cPosition_World.xyz ), sun_light_dir ) );

        vec3 light = mix( vec3( 1.0, 1.0, 1.0 ), cloud_base_color, T );
        T = ( 1.0 - acc.a ) * T;
        acc += vec4( light * T, T );
    }

    acc.rgb /= acc.a + 0.0001;

    // everything from here on out doesn't contribute much to the end result
    float alpha = smoothstep( 0.7, 1.0, acc.a );

    /*
     * Beer's law: transmittance = e^(-thickness)
     *
     * transmittance is the resulting alpha
     * thickness is the density
     */ 

    float T = exp( -alpha );

    acc.rgb -= 0.6 * vec3( 0.8, 0.75, 0.7 ) * alpha * pow( normalize( vec3( 1.0 ) ), vec3( 13.0 ) );
    acc.rgb += 0.2 * vec3( 1.3, 1.2, 1.0 ) * ( 1.0 - alpha ) * pow( normalize( vec3( 1.0 ) ), vec3( 5.0 ) );
