#version 330 core
out vec4 FragColor;

in vec3 TexCoords;

uniform sampler2D u_CloudTexture;
uniform vec3 u_SunPos;

const float PI = 3.14159265359;

/* Atmosphere parameters. */
const float ATMOSPHERE_RADIUS     = 6520275.0;
const float PLANET_RADIUS         = 6371000.0;
const float RAYLEIGH_SCALE_HEIGHT = 8000.0;
const float MIE_SCALE_HEIGHT      = 1200.0;

/* Scattering coefficients. */
const vec3 RAYLEIGH_SCATTERING_COEFF  = vec3(0.0000058, 0.0000135, 0.0000331);
const vec3 MIE_SCATTERING_COEFF       = vec3(0.00002);

/* Mie scattering parameters. */
const float MIE_G = 0.76;

/* Number of samples for ray marching. */
const int SAMPLES = 16;
const int LIGHT_SAMPLES = 8;

float ray_sphere_intersect(vec3 origin, vec3 dir, float radius) {
  float a = dot(dir, dir);
  float b = 2.0 * dot(origin, dir);
  float c = dot(origin, origin) - radius * radius;
  float d = b * b - 4.0 * a * c;
  if (d < 0.0) {
    return -1.0;
  }
  return (-b + sqrt(d)) / (2.0 * a);
}

vec3 get_sky_color(vec3 dir) {
  vec3 eye = vec3(0.0, PLANET_RADIUS + 1.0, 0.0);
  float dist = ray_sphere_intersect(eye, dir, ATMOSPHERE_RADIUS);
  vec3 totalRayleigh = vec3(0.0);
  vec3 totalMie = vec3(0.0);
  float stepSize = dist / float(SAMPLES);
  for (int i = 0; i < SAMPLES; i++) {
    vec3 p = eye + dir * (float(i) + 0.5) * stepSize;
    float h = length(p) - PLANET_RADIUS;
    float rayleighDensity = exp(-h / RAYLEIGH_SCALE_HEIGHT);
    float mieDensity = exp(-h / MIE_SCALE_HEIGHT);
    float lightDist = ray_sphere_intersect(p, normalize(u_SunPos), ATMOSPHERE_RADIUS);
    float lightStepSize = lightDist / float(LIGHT_SAMPLES);
    float opticalDepthLight = 0.0;
    for (int j = 0; j < LIGHT_SAMPLES; j++) {
      vec3 lp = p + normalize(u_SunPos) * (float(j) + 0.5) * lightStepSize;
      float lh = length(lp) - PLANET_RADIUS;
      opticalDepthLight += exp(-lh / RAYLEIGH_SCALE_HEIGHT) * lightStepSize;
    }
    vec3 transmittance = exp(-(opticalDepthLight * RAYLEIGH_SCATTERING_COEFF
          + opticalDepthLight * MIE_SCATTERING_COEFF));
    totalRayleigh += rayleighDensity * stepSize * transmittance;
    totalMie += mieDensity * stepSize * transmittance;
  }

  float mu = dot(dir, normalize(u_SunPos));
  float rayleighPhase = 3.0 / (16.0 * PI) * (1.0 + mu * mu);
  float miePhase = 3.0 / (8.0 * PI) * ((1.0 - MIE_G * MIE_G)
      * (1.0 + mu * mu)) / ((2.0 + MIE_G * MIE_G) * pow(1.0 + MIE_G
      * MIE_G - 2.0 * MIE_G * mu, 1.5));
  return 20.0 * (totalRayleigh * RAYLEIGH_SCATTERING_COEFF * rayleighPhase
      + totalMie * MIE_SCATTERING_COEFF * miePhase);
}

void main() {
  float u = atan(TexCoords.z, TexCoords.x) / (2.0 * PI) + 0.5;
  float v = asin(TexCoords.y) / PI + 0.5;
  float cloudAlpha = texture(u_CloudTexture, vec2(u,v)).r;

  vec3 skyColor = get_sky_color(normalize(TexCoords));
  vec3 cloudColor = mix(skyColor, vec3(1.0), smoothstep(0.5, 0.8, cloudAlpha));

  float poleFade = 1.0 - abs(TexCoords.y);
  float finalAlpha = smoothstep(0.0, 0.2, poleFade) * smoothstep(0.5, 0.8, cloudAlpha);

  FragColor = vec4(pow(cloudColor, vec3(1.0/2.2)), finalAlpha);
}