Atmospheric Scattering

                        
void DrawAtmosphere(const glm::mat4 & P)
{
	// Draw the sky
	progAtmos->bind();

	glm::mat4 V = SetAtmosphereView();
	//This is positioning the quad at the near plane of the perspective view frustum 
	glm::mat4 M = CalcualteFrustumFarBounds(mycam.getFOV(), mycam.getFarDist(), mycam.getAspect());

	//send the matrices to the shaders 
	glUniformMatrix4fv(progAtmos->getUniform("M"), 1, GL_FALSE, &M[0][0]);
	glUniformMatrix4fv(progAtmos->getUniform("V"), 1, GL_FALSE, &V[0][0]);
	glUniformMatrix4fv(progAtmos->getUniform("P"), 1, GL_FALSE, &P[0][0]);
	glUniform3fv(progAtmos->getUniform("campos"), 1, &mycam.getLocation()[0]);
	glUniform3fv(progAtmos->getUniform("uSunPos"), 1, &sun.position[0]);

	glDisable(GL_DEPTH_TEST);
	atmosQuad->draw(progAtmos);

	glEnable(GL_DEPTH_TEST);

	progAtmos->unbind();
}

// This takes a unit quad and places it at the far plane 
//      of the view frustum with the same dimensions
// @fovy The field of view in the y direction of the camera.
// @far The distance to the far plane from the camera.
// @aspect The aspect ratio of the camera.
glm::mat4 CalcualteFrustumFarBounds(float fovy, float far, float aspect) {
	// The length of the hypotenuse of the triangle defined by the camera, the middle of the atmosphere quad, and the top of the atmosphere quad.
	float hypot = far / sin(fovy * 2 * pi / 360.0f);
	// The angle between the hypotenuse and atmosphere quad.
	float angle = 90.0f - fovy;
	// The length of the quad from center to top.
	float side = sin(angle * 2 * pi / 360.0f) * hypot;

	auto S = glm::scale(glm::mat4(1.0f), glm::vec3(side * aspect, side, 1.0f));
	auto T = glm::translate(glm::mat4(1.0f), glm::vec3(0.0f, 0.0f, far));
	return S * T;
}

// Returns the rotation matrix to match the camera's rotation for the atmosphere plane
// Assumes two degrees of freedom for the camera, pitch and yaw. Roll is not currently supported.
glm::mat4 SetAtmosphereView()
{
	auto R = glm::rotate(glm::mat4(1.0f), -mycam.getTheta(), glm::vec3(0.0f, 1.0f, 0.0f));
	R *= glm::rotate(glm::mat4(1.0f), mycam.getPhi(), glm::vec3(1.0f, 0.0f, 0.0f));
	return R;
}
                        
                    

                        
#version 410 core

out vec4 color;
in vec4 fragPos;
uniform vec3 campos;
uniform vec3 uSunPos;

#define PI 3.141592
#define iSteps 16
#define jSteps 8
#define planetRadius 6372e3

vec2 rsi(vec3 r0, vec3 rd, float sr) {
    // ray-sphere intersection that assumes
    // the sphere is centered at the origin.
    // No intersection when result.x > result.y
    float a = dot(rd, rd);
    float b = 2.0 * dot(rd, r0);
    float c = dot(r0, r0) - (sr * sr);
    float d = (b*b) - 4.0*a*c;
    if (d < 0.0) return vec2(1e5,-1e5);
    return vec2(
        (-b - sqrt(d))/(2.0*a),
        (-b + sqrt(d))/(2.0*a)
    );
}

vec3 atmosphere(vec3 r, vec3 r0, vec3 pSun, float iSun, float rPlanet, float rAtmos, vec3 kRlh, float kMie, float shRlh, float shMie, float g) {
    // Normalize the sun and view directions.
    pSun = normalize(pSun);
    r = normalize(r);

    // Calculate the step size of the primary ray.
    vec2 p = rsi(r0, r, rAtmos);
    if (p.x > p.y) return vec3(0,0,0);
    p.y = min(p.y, rsi(r0, r, rPlanet).x);
    float iStepSize = (p.y - p.x) / float(iSteps);

    // Initialize the primary ray time.
    float iTime = 0.0;

    // Initialize accumulators for Rayleigh and Mie scattering.
    vec3 totalRlh = vec3(0,0,0);
    vec3 totalMie = vec3(0,0,0);

    // Initialize optical depth accumulators for the primary ray.
    float iOdRlh = 0.0;
    float iOdMie = 0.0;

    // Calculate the Rayleigh and Mie phases.
    float mu = dot(r, pSun);
    float mumu = mu * mu;
    float gg = g * g;
    float pRlh = 3.0 / (16.0 * PI) * (1.0 + mumu);
    float pMie = 3.0 / (8.0 * PI) * ((1.0 - gg) * (mumu + 1.0)) / (pow(1.0 + gg - 2.0 * mu * g, 1.5) * (2.0 + gg));

    // Sample the primary ray.
    for (int i = 0; i < iSteps; i++) {

        // Calculate the primary ray sample position.
        vec3 iPos = r0 + r * (iTime + iStepSize * 0.5);

        // Calculate the height of the sample.
        float iHeight = length(iPos) - rPlanet;

        // Calculate the optical depth of the Rayleigh and Mie scattering for this step.
        float odStepRlh = exp(-iHeight / shRlh) * iStepSize;
        float odStepMie = exp(-iHeight / shMie) * iStepSize;

        // Accumulate optical depth.
        iOdRlh += odStepRlh;
        iOdMie += odStepMie;

        // Calculate the step size of the secondary ray.
        float jStepSize = rsi(iPos, pSun, rAtmos).y / float(jSteps);

        // Initialize the secondary ray time.
        float jTime = 0.0;

        // Initialize optical depth accumulators for the secondary ray.
        float jOdRlh = 0.0;
        float jOdMie = 0.0;

        // Sample the secondary ray.
        for (int j = 0; j < jSteps; j++) {

            // Calculate the secondary ray sample position.
            vec3 jPos = iPos + pSun * (jTime + jStepSize * 0.5);

            // Calculate the height of the sample.
            float jHeight = length(jPos) - rPlanet;

            // Accumulate the optical depth.
            jOdRlh += exp(-jHeight / shRlh) * jStepSize;
            jOdMie += exp(-jHeight / shMie) * jStepSize;

            // Increment the secondary ray time.
            jTime += jStepSize;
        }

        // Calculate attenuation.
        vec3 attn = exp(-(kMie * (iOdMie + jOdMie) + kRlh * (iOdRlh + jOdRlh)));

        // Accumulate scattering.
        totalRlh += odStepRlh * attn;
        totalMie += odStepMie * attn;

        // Increment the primary ray time.
        iTime += iStepSize;

    }

    // Calculate and return the final color.
    return iSun * (pRlh * kRlh * totalRlh + pMie * kMie * totalMie);
}



void main() {
    vec3 col = atmosphere(
        normalize(vec3(fragPos)),               // normalized ray direction
        vec3(0,planetRadius,0) + campos,        // ray origin
        uSunPos,                                // position of the sun
        22.0,                                   // intensity of the sun
        planetRadius,                           // radius of the planet in meters
        planetRadius + 100e3,                   // radius of the atmosphere in meters
        vec3(5.5e-6, 13.0e-6, 22.4e-6),         // Rayleigh scattering coefficient
        21e-6,                                  // Mie scattering coefficient
        8e3,                                    // Rayleigh scale height
        1.2e3,                                  // Mie scale height
        0.758                                   // Mie preferred scattering direction
    );

    // Apply exposure.
    col = 1.0 - exp(-1.0 * col);

    color = vec4(col, 1);
}