#include <common>
uniform vec3 lightDirection; // The direction vector to the light source
uniform vec3 invWavelength;  // 1 / pow(wavelength, 4) for the red, green, and blue channels
uniform float outerRadius;   // The outer (atmosphere) radius
uniform float innerRadius;   // The inner (planetary) radius
uniform float ESun;          // Sun exposure
uniform float Km;            // Mie coefficient
uniform float Kr;            // Rayleigh coefficient
uniform float scale;         // 1 / (outerRadius - innerRadius)
uniform float scaleDepth;    // The scale depth (i.e. the altitude at which the atmosphere's average density is found)

const int nSamples = 2;
const float fSamples = 1.0;

varying vec3 v3Direction;
varying vec3 c0;
varying vec3 c1;

float dscale(float fCos) {
  float x = 1.0 - fCos;
  return scaleDepth * exp(-0.00287 + x * (0.459 + x * (3.83 + x * (-6.80 + x * 5.25))));
}

float calculateNearScatter(vec3 v3CameraPosition, vec3 v3Ray, float fCameraHeight, float fOuterRadius) {
  float B = 2.0 * dot(v3CameraPosition, v3Ray);
  float C = pow(fCameraHeight, 2.0) - pow(fOuterRadius, 2.0);
  float fDet = max(0.0, B * B - 4.0 * C);
  return 0.5 * (-B - sqrt(fDet));
}

void main(void) {
  vec3 mvPosition = (modelMatrix * vec4(position, 1.0)).xyz;

  // Initialize variables
  float cameraHeight = length(vec3(0.0, 0.0, 0.0) - cameraPosition);
  float KmESun = Km * ESun;
  float KrESun = Kr * ESun;
  float Kr4PI = Kr * 4.0 * PI;
  float Km4PI = Km * 4.0 * PI;
  float scaleOverScaleDepth = scale / scaleDepth;

  // Get the ray from the camera to the vertex and its length (which is the far point of the ray passing through the atmosphere)
  vec3 v3Ray = mvPosition - cameraPosition;
  float fFar = length(v3Ray);
  v3Ray /= fFar;

  vec3 v3Start;
  float fStartAngle;
  float fStartDepth;
  float fStartOffset;

  if (cameraHeight > outerRadius) {
    // Sky from space
    // Calculate the closest intersection of the ray with the outer atmosphere (which is the near point of the ray passing through the atmosphere)
    float fNear = calculateNearScatter(cameraPosition, v3Ray, cameraHeight, outerRadius);

    // Calculate the ray's starting position, then calculate its scattering offset
    v3Start = cameraPosition + v3Ray * fNear;
    fFar -= fNear;
    fStartAngle = dot(v3Ray, v3Start) / outerRadius;
    fStartDepth = exp(-1.0 / scaleDepth);
    fStartOffset = fStartDepth * dscale(fStartAngle);
  } else {
    // Sky from within the atmosphere
    v3Start = cameraPosition;
    fStartDepth = exp(scaleOverScaleDepth * (innerRadius - cameraHeight));
    fStartAngle = dot(v3Ray, v3Start) / length(v3Start);
    fStartOffset = fStartDepth * dscale(fStartAngle);
  }

  // Initialize the scattering loop variables
  float fSampleLength = fFar / fSamples;
  float scaledLength = fSampleLength * scale;
  vec3 v3SampleRay = v3Ray * fSampleLength;
  vec3 v3SamplePoint = v3Start + v3SampleRay * 0.5;

  // Now loop through the sample rays
  vec3 v3FrontColor = vec3(0.0, 0.0, 0.0);
  for(int i=0; i<nSamples; i++)
  {
    float fHeight = length(v3SamplePoint);
    float fDepth = exp(scaleOverScaleDepth * (innerRadius - fHeight));
    float fLightAngle = dot(lightDirection, v3SamplePoint) / fHeight;
    float fCameraAngle = dot(v3Ray, v3SamplePoint) / fHeight;
    float fScatter = (fStartOffset + fDepth * (dscale(fLightAngle) - dscale(fCameraAngle)));
    vec3 v3Attenuate = exp(-fScatter * (invWavelength * Kr4PI + Km4PI));
    v3FrontColor += v3Attenuate * (fDepth * scaledLength);
    v3SamplePoint += v3SampleRay;
  }

  // Finally, scale the Mie and Rayleigh colors and set up the varying variables for the pixel shader
  gl_Position = projectionMatrix * viewMatrix * modelMatrix * vec4(position, 1);
  c0 = v3FrontColor * (invWavelength * KrESun);
  c1 = v3FrontColor * KmESun;
  v3Direction = cameraPosition - mvPosition;
}