scatter.frag

vec3 sundir;
vec3 I_R, I_M;
vec2 totalDepthRM;


    // consts 
    const vec3 rayScatterCoeff = vec3(58e-7, 135e-7, 331e-7);
    const vec3 rayEffectiveCoeff = rayScatterCoeff; // Rayleight doesn't absorb light 
    
    const vec3 mieScatterCoeff = vec3(2e-5);
    const vec3 mieEffectiveCoeff = mieScatterCoeff * 1.1; // Approximate absorption as a factor of scattering
    
    const vec3 sunIntensity = vec3(7); // sun is modelled as infinitely far away
    
    const float earthRadius = 6360e3;
    const float atmosphereRadius = 6380e3;
	const float hRay = 8e3;
	const float hMie = 12e2;

    const vec3 center = vec3(0., -earthRadius, 0.); // earth center point
	

// Basically a ray-sphere intersection. Find distance to where rays escapes a sphere with given radius.
// Used to calculate length at which ray escapes atmosphere
float escape(vec3 p, vec3 d, float R) {
	vec3 v = p - center;
	float b = dot(v, d);
	float det = b * b - dot(v, v) + R*R;
	if (det < 0.) return -1.;
	det = sqrt(det);
	float t1 = -b - det, t2 = -b + det;
	return (t1 >= 0.) ? t1 : t2;
}

vec2 densitiesRM (vec3 p) {
	float h = max(0.,length(p - center) - earthRadius);
    return vec2(exp(-h/hRay), exp(-h/hMie));
}

vec2 scatterDepthInt (vec3 o, vec3 d, float L, float steps) {
    // Approximate sampling the middle
    // return densitiesRM (o + d * (L / 2.)) * L;
    
    // Approximate by combining 2 samples
   	return (densitiesRM (o) * (L / 2.) + densitiesRM (o + d * L) * (L / 2.));
    
    // Integrate
    
    // Accumulator
	vec2 depthRMs = vec2(0.);

	// Set L to be step distance and pre-multiply d with it
	L /= steps; d *= L;
	
	// Go from point P to A
	for (float i = 0.; i < steps; ++i)
		// Simply accumulate densities
		depthRMs += densitiesRM(o + d * i);

	return depthRMs * L;
}

void scatterIn (vec3 origin, vec3 direction, float depth, float steps) {

    
    depth = depth / steps;
    
    for (float i = 0.; i < steps; i++) {
    	vec3 p = origin + direction * (depth * i);
    	vec2 dRM = densitiesRM(p) * depth;
    	totalDepthRM += dRM;
        
        // Calculate optical depth
        vec2 depthRMsum = totalDepthRM + scatterDepthInt(p, sundir, escape(p, sundir, atmosphereRadius), 4.);
    	
        // Calculate exponent part of both integrals
		vec3 A = exp(-rayEffectiveCoeff * depthRMsum.x - mieEffectiveCoeff * depthRMsum.y);
        
        I_R += A * dRM.x;
        I_M += A * dRM.y;
    }
}

// Calculate the complete scattering effect for a ray going from origin
// along direction for a distance of depth. 
// luminance represents the light coming from that point e.g. in a scene
// If luminance is zero the effective return is in-scattering from the sun
vec3 scatter (vec3 origin, vec3 direction, float depth, vec3 luminance) {
   
	I_R = I_M = vec3(0.);
    
    scatterIn (origin, direction, depth, 12.);
    
    float mu = dot(direction, sundir);
    
    vec3 extinction = luminance * exp(-rayEffectiveCoeff * totalDepthRM.x - mieEffectiveCoeff * totalDepthRM.y);
    
    return extinction
        	// Add in-scattering
		+ sunIntensity * (1. + mu * mu) * (
            // 3/16pi = 0.597
			I_R * rayEffectiveCoeff * .0597 +
            // mie phase function is complicated [Nisita et al. 93, Riley et al. 04]
			I_M * mieScatterCoeff * .0196 / pow(1.58 - 1.52 * mu, 1.5));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from -1 to 1)
    vec2 uv = fragCoord/iResolution.xy * 2. - 1.;
    
    // Fix aspect
    uv.x *= iResolution.x / iResolution.y;
    
    
	// sundir could be set as a uniform or calculated dynamically, this is for demo
    sundir = normalize (vec3(.5, .4 * (1. + sin(.5 * iTime)), -1.));
    
    
    vec3 origin = vec3(0., 0., 0.);
    vec3 direction = normalize(vec3(uv, -2.));
    vec3 col = vec3(0.);
    
    // Here you would do scene-intersection
    float L = escape (origin, direction, atmosphereRadius);
    col = scatter (origin, direction, L, col);
    
    // not sure what's with the sqrt tbh.
    fragColor = vec4(sqrt(col), 1.);
}