Atmo Scattering

/*
  
  Thanks to whomever wrote the shader at http://glslsandbox.com/e#36146.5. My code is based off of their's.

*/

#define INCLUDED_ATMOSCATTERING

#define ATMOSPHERIC_SCATTERING // Swaps out the sky gradient for a proper atmospheric scattering model. Atmospheric scattering models how light would interact with an actual atmosphere, producing realistic colours as the sun nears the horizon.

//#define ATMOSCATTERING_TONEMAP // Should atmospheric scattering use a tonemap? By default this is disabled as a tonemap is redundant due to the shader having a tonemap of its own, but this is left for personal preference.

#define ATMOSCATTERING_RAIN_DESATURATION -0.8
#define ATMOSCATTERING_RAIN_STRENGTH 0.5

const float turbidity = 1.5;
float rayleighCoefficient = mix(2.5, 5.0, rain);

const float mieCoefficient = 0.005;
const float mieDirectionalG = 0.76;

// Constants for atmo scattering
const float n = 1.00029;
const float N = 2.54743E25; // scattering width. lower = wider. [1.54743E25 2.54743E25]

const float pn = 0.02;

// Wavelength of used primaries, according to praetham
const vec3 primaryWavelengths = vec3(680, 550, 450) * 1.0E-9;

// Mie stuff
// K coefficient for the primaries
const vec3 K = vec3(0.686, 0.678, 0.666);
const float v = 4.0;

// Optical length at zenith for molecules
const float rayleighZenithLength = 8.4E3;
const float mieZenithLength = 1.25E3;

const float sunIntensity = 1000.0;
const float sunAngularDiameterCos = 0.99883194915 * 1.0;
const float moonAngularDiameterCos = 0.99883194915 * 0.999;

const float cutoffAngle = pi * 0.5128205128205128;
const float steepness = 1.0;

float rayleighPhase(in float cosViewSunAngle) {
  return (3.0 / (16.0 * pi)) * (1.0 + pow2(cosViewSunAngle));
}

vec3 totalRayleigh(in vec3 lambda, in float n, in float N, in float pn){
	return (24.0 * pow3(pi) * pow(pow2(n) - 1.0, 2.0) * (6.0 + 3.0 * pn)) / (N * pow(lambda, vec3(4.0)) * pow(pow2(n) + 2.0, 2.0) * (6.0 - 7.0 * pn));
}

vec3 totalMie(in vec3 primaryWavelengths, in vec3 K, in float T) {
  return 0.434 * ((0.2 * T) * 10E-18) * pi * pow((2.0 * pi) / primaryWavelengths, vec3(v - 2.0)) * K;
}

float hgPhase(in float cosViewSunAngle, in float g) {
  return (1.0 / (4.0 * pi)) * ((1.0 - pow2(g)) / pow(1.0 - 2.0 * g * cosViewSunAngle + pow2(g), 1.5));
}

float SunIntensity(in float zenithAngleCos, in float intensity) {
  return intensity * max(0.0, 1.0 - exp(-((cutoffAngle - acos(zenithAngleCos)) / steepness)));
}

#ifdef ATMOSCATTERING_TONEMAP
const float A = 0.10;
const float B = 0.50;
const float C = 0.10;
const float D = 0.20;
const float E = 0.03;
const float F = 0.30;
const float W = 1000.0;

vec3 Uncharted2Tonemap(in vec3 x) {
   return ((x * (A * x + C * B) + D * E) / (x * (A * x + B) + D * F)) - E / F;
}

vec3 ToneMap(in vec3 color) {
    vec3 toneMappedColor;
    
    toneMappedColor = color * 0.04;
    toneMappedColor = Uncharted2Tonemap(toneMappedColor);
    
    float sunfade = 1.0 - clamp(1.0 - exp(-(sunVector.z / 500.0)), 0.0, 1.0);
    toneMappedColor = pow(toneMappedColor, vec3(1.0 / (1.2 + (1.2 * sunfade))));
    
    return toneMappedColor;
}
#endif

vec3 atmosphericScattering(in vec3 fragpos, in bool reflection) {
  // Cos angles
  vec4 dotVector = vec4(0.0);

  #define cosViewSunAngle dotVector.x
  #define cosSunUpAngle dotVector.y
  #define cosUpViewAngle dotVector.z
  #define cosViewMoonAngle dotVector.w

  cosViewSunAngle = max(0.0, dot(fragpos, sunVector));
  cosSunUpAngle = dot(sunVector, upVector) * 0.95 + 0.05;
  cosUpViewAngle = dot(upVector, fragpos);
  cosViewMoonAngle = max(0.0, dot(fragpos, -sunVector));

  float sunE = SunIntensity(cosSunUpAngle, sunIntensity);

  // Extinction (absorption + out scattering)
  mat2x3 extinctionMatrix;

  #define rayleighAtX extinctionMatrix[0]
  #define mieAtX extinctionMatrix[1]

  // Rayleigh coefficients
  rayleighAtX = totalRayleigh(primaryWavelengths, n, N, pn) * rayleighCoefficient;

  // Mie coefficients
  mieAtX = totalMie(primaryWavelengths, K, turbidity) * mieCoefficient;

  // Optical length
  // Cutoff angle at 90 to avoid singularity in next formula
  vec3 opticalLengthVector;

  #define zenithAngle opticalLengthVector.x
  #define rayleighOpticalLength opticalLengthVector.y
  #define mieOpticalLength opticalLengthVector.z

  zenithAngle = max(0.0, cosUpViewAngle);

  rayleighOpticalLength = rayleighZenithLength / zenithAngle;
  mieOpticalLength = mieZenithLength / zenithAngle;

  // Combined extinction factor
  vec3 Fex = exp(-(rayleighAtX * rayleighOpticalLength + mieAtX * mieOpticalLength));

  #undef zenithAngle
  #undef rayleighOpticalLength
  #undef mieOpticalLength

  // In scattering
  mat4x3 inscatteringMatrix;

  #define rayleighXtoEye inscatteringMatrix[0]
  #define mieXtoEye inscatteringMatrix[1]
  #define totalLightAtX inscatteringMatrix[2]
  #define lightFromXtoEye inscatteringMatrix[3]

  rayleighXtoEye = rayleighAtX * rayleighPhase(cosViewSunAngle);
  mieXtoEye = mieAtX * hgPhase(cosViewSunAngle, mieDirectionalG);

  totalLightAtX = rayleighAtX + mieAtX;
  lightFromXtoEye = rayleighXtoEye + mieXtoEye;

  vec3 somethingElse = sunE * (lightFromXtoEye / totalLightAtX);

  #undef rayleighAtX
  #undef mieAtX

  #undef rayleighXtoEye
  #undef mieXtoEye
  #undef totalLightAtX
  #undef lightFromXtoEye

  mat3 componentMatrix;

  #define sky componentMatrix[0]
  #define sun componentMatrix[1]
  #define moon componentMatrix[2]

  sky = somethingElse * (1.0 - Fex);
  sky = colourSaturate(sky, mix(0.0, ATMOSCATTERING_RAIN_DESATURATION, rain)) * mix(1.0, ATMOSCATTERING_RAIN_STRENGTH, rain);
  sky *= mix(vec3(1.0), pow(somethingElse * Fex, vec3(0.5)), clamp(pow(1.0 - dot(upVector, sunVector), 5.0), 0.0, 1.0));
  sky = max(sky, 0.0);

  // Composition + solar disc
  float sundisk = smoothstep(sunAngularDiameterCos, sunAngularDiameterCos + 0.00002, cosViewSunAngle);
  if(!reflection) sun = (sunE * 8.0 * Fex) * sundisk;

  float moondisk = smoothstep(moonAngularDiameterCos, moonAngularDiameterCos + 0.00002, cosViewMoonAngle);
  if(!reflection) moon = vec3(dlMoonlight * 1000.0) * moondisk;

  #ifdef ATMOSCATTERING_TONEMAP
    return ToneMap(sky + sun);
  #else
    return (sky + sun + moon) * 0.01;
  #endif

  #undef sky
  #undef sun
  #undef moon

  #undef cosViewSunAngle
  #undef cosSunUpAngle
  #undef cosUpViewAngle
}