GLSL port of https://github.com/Fewes/MinimalAtmosphere

atmosphere.glsl

//MIT License
//Copyright (c) 2021 Felix Westin
//Source: https://github.com/Fewes/MinimalAtmosphere

//Ported to GLSL by Marcin Ignac

#ifndef ATMOSPHERE_INCLUDED
#define ATMOSPHERE_INCLUDED

// -------------------------------------
// Defines
#define EPS					1e-6
#define PI					3.14159265359
#define INFINITY 			1.0 / 0.0
#define PLANET_RADIUS		6371000.0
#define PLANET_CENTER		vec3(0.0, -PLANET_RADIUS, 0.0)
#define ATMOSPHERE_HEIGHT	1.0*100000.0
#define RAYLEIGH_HEIGHT		(ATMOSPHERE_HEIGHT * 0.08)
#define MIE_HEIGHT			(ATMOSPHERE_HEIGHT * 0.012)

// -------------------------------------
// Coefficients
#define C_RAYLEIGH			(vec3(5.802, 13.558, 33.100) * 1e-6)
#define C_MIE				(vec3(3.996,  3.996,  3.996) * 1e-6)
#define C_OZONE				(vec3(0.650,  1.881,  0.085) * 1e-6)

#define ATMOSPHERE_DENSITY	1.0
#define EXPOSURE			20.0

float saturate(float v) {
  return clamp(v, 0.0, 1.0);
}

// -------------------------------------
// Math
vec2 SphereIntersection (vec3 rayStart, vec3 rayDir, vec3 sphereCenter, float sphereRadius)
{
  rayStart -= sphereCenter;
  float a = dot(rayDir, rayDir);
  float b = 2.0 * dot(rayStart, rayDir);
  float c = dot(rayStart, rayStart) - (sphereRadius * sphereRadius);
  float d = b * b - 4.0 * a * c;
  if (d < 0.0)
  {
    return vec2(-1.0);
  }
  else
  {        
    d = sqrt(d);
    return vec2(-b - d, -b + d) / (2.0 * a);
  }
}

vec2 PlanetIntersection (vec3 rayStart, vec3 rayDir)
{
  return SphereIntersection(rayStart, rayDir, PLANET_CENTER, PLANET_RADIUS);
}

vec2 AtmosphereIntersection (vec3 rayStart, vec3 rayDir)
{
  return SphereIntersection(rayStart, rayDir, PLANET_CENTER, PLANET_RADIUS + ATMOSPHERE_HEIGHT);
}

// -------------------------------------
// Phase functions
float PhaseRayleigh (float costh)
{
  return 3.0 * (1.0 + costh*costh) / (16.0 * PI);
}
float PhaseMie (float costh, float g)
{      
  g = min(g, 0.9381);
  float k = 1.55 * g - 0.55*g*g*g;
  float kcosth = k*costh;
  return (1.0 - k*k) / ((4.0 * PI) * (1.0-kcosth) * (1.0-kcosth));
}
float PhaseMie (float costh)
{
  return PhaseMie(costh, 0.85);
}

// -------------------------------------
// Atmosphere
float AtmosphereHeight (vec3 positionWS)
{
  return distance(positionWS, PLANET_CENTER) - PLANET_RADIUS;
}
float DensityRayleigh (float h)
{
  return exp(-max(0.0, h / RAYLEIGH_HEIGHT));
}
float DensityMie (float h)
{
  return exp(-max(0.0, h / MIE_HEIGHT));
}
float DensityOzone (float h)
{
  // Ozone is represented as a tent function with a width of 30km, centered around an altitude of 25km.
  return max(0.0, 1.0 - abs(h - 25000.0) / 15000.0);
}
vec3 AtmosphereDensity (float h)
{
  return vec3(DensityRayleigh(h), DensityMie(h), DensityOzone(h));
}

// Optical depth is a unitless measurement of the amount of absorption of a participating medium (such as the atmosphere).
// This function calculates just that for our three atmospheric elements:
// R: Rayleigh
// G: Mie
// B: Ozone
// If you find the term "optical depth" confusing, you can think of it as "how much density was found along the ray in total".
vec3 IntegrateOpticalDepth (vec3 rayStart, vec3 rayDir)
{
  vec2 intersection = AtmosphereIntersection(rayStart, rayDir);
  float  rayLength    = intersection.y;

  int    sampleCount  = 8;
  float  stepSize     = rayLength / float(sampleCount);
  vec3 opticalDepth = vec3(0.0);

  //for (int i = 0; i < sampleCount; i++)
  for (int i = 0; i < 8; i++)
  {
    vec3 localPosition = rayStart + rayDir * (float(i) + 0.5) * stepSize;
    float  localHeight   = AtmosphereHeight(localPosition);
    vec3 localDensity  = AtmosphereDensity(localHeight) * stepSize;

    opticalDepth += localDensity;
  }

  return opticalDepth;
}

// Calculate a luminance transmittance value from optical depth.
vec3 Absorb (vec3 opticalDepth)
{
  // Note that Mie results in slightly more light absorption than scattering, hence * 1.1
  return exp(-(opticalDepth.x * C_RAYLEIGH + opticalDepth.y * C_MIE * 1.1 + opticalDepth.z * C_OZONE) * ATMOSPHERE_DENSITY);
}

// Integrate scattering over a ray for a single directional light source.
// Also return the transmittance for the same ray as we are already calculating the optical depth anyway.
vec3 IntegrateScattering (vec3 rayStart, vec3 rayDir, float rayLength, vec3 lightDir, vec3 lightColor, out vec3 transmittance)
{
  // We can reduce the number of atmospheric samples required to converge by spacing them exponentially closer to the camera.
  // This breaks space view however, so let's compensate for that with an exponent that "fades" to 1 as we leave the atmosphere.
  float  rayHeight = AtmosphereHeight(rayStart);
  float  sampleDistributionExponent = 1.0 + saturate(1.0 - rayHeight / ATMOSPHERE_HEIGHT) * 8.0; // Slightly arbitrary max exponent of 9

  vec2 intersection = AtmosphereIntersection(rayStart, rayDir);
  rayLength = min(rayLength, intersection.y);
  if (intersection.x > 0.0)
  {
    // Advance ray to the atmosphere entry point
    rayStart += rayDir * intersection.x;
    rayLength -= intersection.x;
  }

  float  costh     = dot(rayDir, lightDir);
  float  phaseRayleigh = PhaseRayleigh(costh);
  float  phaseMie = PhaseMie(costh);
  vec3 rayleigh = vec3(0.0);
  vec3 mie      = vec3(0.0);

  int  sampleCount = 64;
  vec3 opticalDepth = vec3(0.0);

  float prevRayTime = 0.0;
  // for (int i = 0; i < sampleCount; i++)
  for (int i = 0; i < 64; i++)
  {
    float  rayTime = pow(float(i) / float(sampleCount), sampleDistributionExponent) * rayLength;
    // Because we are distributing the samples exponentially, we have to calculate the step size per sample.
    float  stepSize = (rayTime - prevRayTime);

    vec3 localPosition = rayStart + rayDir * rayTime;
    float  localHeight   = AtmosphereHeight(localPosition);
    vec3 localDensity  = AtmosphereDensity(localHeight) * stepSize;

    opticalDepth += localDensity;

    vec3 opticalDepthlight  = IntegrateOpticalDepth(localPosition, lightDir);
    vec3 lightTransmittance = Absorb(opticalDepth + opticalDepthlight);

    rayleigh += lightTransmittance * phaseRayleigh * localDensity.x;
    mie      += lightTransmittance * phaseMie      * localDensity.y;

    prevRayTime = rayTime;
  }

  transmittance = Absorb(opticalDepth);


  return (rayleigh * C_RAYLEIGH + mie * C_MIE) * lightColor * EXPOSURE;
}

#endif // ATMOSPHERE_INCLUDED