Revision control
++ // Copyright (C) 2023 Fernando García Liñán
++ // Copyright (C) 2026 The Qt Company Ltd.
++ // SPDX-License-Identifier: MIT
++
++ // Configurable parameters
++
++ // Ray marching steps. More steps mean better accuracy but worse performance
++ const int TRANSMITTANCE_STEPS = 32;
++ const int IN_SCATTERING_STEPS = 16;
++
++ // Debug
++ #define ENABLE_SPECTRAL 1
++ #define ENABLE_MULTIPLE_SCATTERING 1
++ #define ENABLE_AEROSOLS 1
++
++ //-----------------------------------------------------------------------------
++ // Constants
++
++ // All parameters that depend on wavelength (vec4) are sampled at
++ // 630, 560, 490, 430 nanometers
++
++ const float PI = 3.14159265358979323846;
++ const float INV_PI = 0.31830988618379067154;
++ const float INV_4PI = 0.25 * INV_PI;
++ const float PHASE_ISOTROPIC = INV_4PI;
++ const float RAYLEIGH_PHASE_SCALE = (3.0 / 16.0) * INV_PI;
++ const float g = 0.8;
++ const float gg = g*g;
++
++ const float EARTH_RADIUS = 6371.0; // km
++ const float ATMOSPHERE_THICKNESS = 100.0; // km
++ const float ATMOSPHERE_RADIUS = EARTH_RADIUS + ATMOSPHERE_THICKNESS;
++
++ #if ENABLE_SPECTRAL == 1
++ // Extraterrestial Solar Irradiance Spectra, units W * m^-2 * nm^-1
++ const vec4 sunSpectralIrradiance = vec4(1.679, 1.828, 1.986, 1.307);
++ // Rayleigh scattering coefficient at sea level, units km^-1
++ // "Rayleigh-scattering calculations for the terrestrial atmosphere"
++ // by Anthony Bucholtz (1995).
++ const vec4 molecularScatteringCoefficientBase = vec4(6.605e-3, 1.067e-2, 1.842e-2, 3.156e-2);
++ // Ozone absorption cross section, units m^2 / molecules
++ // "High spectral resolution ozone absorption cross-sections"
++ // by V. Gorshelev et al. (2014).
++ const vec4 ozoneAbsorptionCrossSection = vec4(3.472e-21, 3.914e-21, 1.349e-21, 11.03e-23) * 1e-4f;
++ #else
++ // Same as above but for the following "RGB" wavelengths: 680, 550, 440 nm
++ // The Sun spectral irradiance is also multiplied by a constant factor to
++ // compensate for the fact that we use the spectral samples directly as RGB,
++ // which is incorrect.
++ const vec4 sunSpectralIrradiance = vec4(1.500, 1.864, 1.715, 0.0) * 150.0;
++ const vec4 molecularScatteringCoefficientBase = vec4(4.847e-3, 1.149e-2, 2.870e-2, 0.0);
++ const vec4 ozoneAbsorptionCrossSection = vec4(3.36e-21f, 3.08e-21f, 20.6e-23f, 0.0) * 1e-4f;
++ #endif
++
++ /*
++ * Every aerosol type expects 5 parameters:
++ * - Scattering cross section
++ * - Absorption cross section
++ * - Base density (km^-3)
++ * - Background density (km^-3)
++ * - Height scaling parameter
++ * These parameters can be sent as uniforms.
++ *
++ * This model for aerosols and their corresponding parameters come from
++ * "A Physically-Based Spatio-Temporal Sky Model"
++ * by Guimera et al. (2018).
++ */
++
++ vec3 getSunDirection(float time)
++ {
++ float a = sin(time*0.5 - 1.5) * 0.55 + 0.45;
++ return vec3(-sqrt(1.0 - a*a), 0.0, a);
++ }
++
++ /*
++ * Helper function to obtain the transmittance to the top of the atmosphere
++ * from Buffer A.
++ */
++ vec4 transmittanceFromLut(sampler2D lut, float cosTheta, float normalizedAltitude)
++ {
++ float u = clamp(cosTheta * 0.5 + 0.5, 0.0, 1.0);
++ float v = clamp(normalizedAltitude, 0.0, 1.0);
++ return texture(lut, vec2(u, v));
++ }
++
++ /*
++ * Returns the distance between ro and the first intersection with the sphere
++ * or -1.0 if there is no intersection. The sphere's origin is (0,0,0).
++ * -1.0 is also returned if the ray is pointing away from the sphere.
++ */
++ float raySphereIntersection(vec3 ro, vec3 rd, float radius)
++ {
++ float b = dot(ro, rd);
++ float c = dot(ro, ro) - radius*radius;
++ if (c > 0.0 && b > 0.0) return -1.0;
++ float d = b*b - c;
++ if (d < 0.0) return -1.0;
++ if (d > b*b) return (-b+sqrt(d));
++ return (-b-sqrt(d));
++ }
++
++ /*
++ * Rayleigh phase function.
++ */
++ float molecularPhaseFunction(float cosTheta)
++ {
++ return RAYLEIGH_PHASE_SCALE * (1.0 + cosTheta*cosTheta);
++ }
++
++ /*
++ * Henyey-Greenstrein phase function.
++ */
++ float aerosolPhaseFunction(float cosTheta)
++ {
++ float den = 1.0 + gg + 2.0 * g * cosTheta;
++ return INV_4PI * (1.0 - gg) / (den * sqrt(den));
++ }
++
++ vec4 getMultipleScattering(sampler2D transmittanceLut, float cosTheta, float normalizedHeight, float d, vec4 groundAlbedo)
++ {
++ #if ENABLE_MULTIPLE_SCATTERING == 1
++ // Solid angle subtended by the planet from a point at d distance
++ // from the planet center.
++ float omega = 2.0 * PI * (1.0 - sqrt(d*d - EARTH_RADIUS*EARTH_RADIUS) / d);
++
++ vec4 tToGround = transmittanceFromLut(transmittanceLut, cosTheta, 0.0);
++
++ vec4 tGroundToSample =
++ transmittanceFromLut(transmittanceLut, 1.0, 0.0) /
++ transmittanceFromLut(transmittanceLut, 1.0, normalizedHeight);
++
++ // 2nd order scattering from the ground
++ vec4 lGround = PHASE_ISOTROPIC * omega * (groundAlbedo / PI) * tToGround * tGroundToSample * cosTheta;
++
++ // Fit of Earth's multiple scattering coming from other points in the atmosphere
++ vec4 lMs = 0.02 * vec4(0.217, 0.347, 0.594, 1.0) * (1.0 / (1.0 + 5.0 * exp(-17.92 * cosTheta)));
++
++ return lMs + lGround;
++ #else
++ return vec4(0.0);
++ #endif
++ }
++
++ /*
++ * Return the molecular volume scattering coefficient (km^-1) for a given altitude
++ * in kilometers.
++ */
++ vec4 getMolecularScatteringCoefficient(float h)
++ {
++ return molecularScatteringCoefficientBase * exp(-0.07771971 * pow(h, 1.16364243));
++ }
++
++ /*
++ * Return the molecular volume absorption coefficient (km^-1) for a given altitude
++ * in kilometers.
++ */
++ vec4 getMolecularAbsorptionCoefficient(float h, float ozoneMeanMonthlyDobson)
++ {
++ h += 1e-4; // Avoid division by 0
++ float t = log(h) - 3.22261;
++ float density = 3.78547397e20 * (1.0 / h) * exp(-t * t * 5.55555555);
++ return ozoneAbsorptionCrossSection * ozoneMeanMonthlyDobson * density;
++ }
++
++ float getAerosolDensity(float h,
++ int aerosolType,
++ float aerosolBaseDensity,
++ float aerosolHeightScale,
++ float aerosolBackgroundDividedByBaseDensity
++ )
++ {
++ if (aerosolType == 0) { // Only for the Background aerosol type, no dependency on height
++ return aerosolBaseDensity * (1.0 + aerosolBackgroundDividedByBaseDensity);
++ } else {
++ return aerosolBaseDensity * (exp(-h / aerosolHeightScale)
++ + aerosolBackgroundDividedByBaseDensity);
++ }
++ }
++
++ /*
++ * Get the collision coefficients (scattering and absorption) of the
++ * atmospheric medium for a given point at an altitude h.
++ */
++ void getAtmosphereCollisionCoefficients(in float h,
++ in vec4 aerosolAbsorptionCrossSection,
++ in vec4 aerosolScatteringCrossSection,
++ in float aerosolBaseDensity,
++ in float aerosolHeightScale,
++ in float aerosolBackgroundDividedByBaseDensity,
++ in int aerosolType,
++ in float ozoneMeanMonthlyDobson,
++ in float aerosolTurbidity,
++ out vec4 aerosolAbsorption,
++ out vec4 aerosolScattering,
++ out vec4 molecularAbsorption,
++ out vec4 molecularScattering,
++ out vec4 extinction)
++ {
++ h = max(h, 0.0); // In case height is negative
++ #if ENABLE_AEROSOLS == 0
++ aerosolAbsorption = vec4(0.0);
++ aerosolScattering = vec4(0.0);
++ #else
++ float aerosolDensity = getAerosolDensity(h, aerosolType, aerosolBaseDensity, aerosolHeightScale, aerosolBackgroundDividedByBaseDensity);
++ aerosolAbsorption = aerosolAbsorptionCrossSection * aerosolDensity * aerosolTurbidity;
++ aerosolScattering = aerosolScatteringCrossSection * aerosolDensity * aerosolTurbidity;
++ #endif
++ molecularAbsorption = getMolecularAbsorptionCoefficient(h, ozoneMeanMonthlyDobson);
++ molecularScattering = getMolecularScatteringCoefficient(h);
++ extinction = aerosolAbsorption + aerosolScattering + molecularAbsorption + molecularScattering;
++ }
++
++ //-----------------------------------------------------------------------------
++ // Spectral rendering stuff
++
++ const mat4 M = mat4(
++ 137.672389239975, -8.632904716299537, -1.7181567391931372, 0,
++ 32.549094028629234, 91.29801417199785, -12.005406444382531, 0,
++ -38.91428392614275, 34.31665471469816, 29.89044807197628, 0,
++ 8.572844237945445, -11.103384660054624, 117.47585277566478, 0
++ );
++
++ vec3 linearSrgbFromSpectralSamples(vec4 L)
++ {
++ return (M * L).xyz;
++ }
++
++ /*
++ * Buffer B: Sky texture
++ *
++ * "A Scalable and Production Ready Sky and Atmosphere Rendering Technique"
++ * by Sébastien Hillaire (2020).
++ *
++ * We render the sky to a texture instead of raymarching on the entire screen.
++ * This is not very useful in Shadertoy, but very useful for someone looking
++ * to implement this on a real application.
++ *
++ * It is important to note that quality decreases significantly when rendering
++ * space views. To avoid this, the computeInscattering() function can be used
++ * directly when rendering to a fullscreen quad.
++ */
++
++ vec4 computeInscattering(vec3 rayOrigin, vec3 rayDir, float tD, out vec4 transmittance)
++ {
++
++ float cosTheta = dot(-rayDir, sunDirection);
++
++ float molecularPhase = molecularPhaseFunction(cosTheta);
++ float aerosolPhase = aerosolPhaseFunction(cosTheta);
++
++ float dt = tD / float(IN_SCATTERING_STEPS);
++
++ vec4 lInscattering = vec4(0.0);
++ transmittance = vec4(1.0);
++
++ for (int i = 0; i < IN_SCATTERING_STEPS; ++i) {
++ float t = (float(i) + 0.5) * dt;
++ vec3 xT = rayOrigin + rayDir * t;
++
++ float distanceToEarthCenter = length(xT);
++ vec3 zenithDir = xT / distanceToEarthCenter;
++ float altitude = distanceToEarthCenter - EARTH_RADIUS;
++ float normalizedAltitude = altitude / ATMOSPHERE_THICKNESS;
++
++ float sampleCosTheta = dot(zenithDir, sunDirection);
++
++ vec4 aerosolAbsorption, aerosolScattering;
++ vec4 molecularAbsorption, molecularScattering;
++ vec4 extinction;
++ getAtmosphereCollisionCoefficients(
++ altitude,
++ aerosolAbsorptionCrossSection,
++ aerosolScatteringCrossSection,
++ aerosolBaseDensity,
++ aerosolHeightScale,
++ aerosolBackgroundDividedByBaseDensity,
++ aerosolType,
++ ozoneMeanMonthlyDobson,
++ aerosolTurbidity,
++ aerosolAbsorption, aerosolScattering,
++ molecularAbsorption, molecularScattering,
++ extinction);
++
++ vec4 transmittanceToSun = transmittanceFromLut(
++ transmittanceBuffer, sampleCosTheta, normalizedAltitude);
++
++ vec4 ms = getMultipleScattering(
++ transmittanceBuffer, sampleCosTheta, normalizedAltitude,
++ distanceToEarthCenter,
++ groundAlbedo);
++
++ vec4 S = sunSpectralIrradiance *
++ (molecularScattering * (molecularPhase * transmittanceToSun + ms) +
++ aerosolScattering * (aerosolPhase * transmittanceToSun + ms));
++
++ vec4 stepTransmittance = exp(-dt * extinction);
++
++ // Energy-conserving analytical integration
++ // "Physically Based Sky, Atmosphere and Cloud Rendering in Frostbite"
++ // by Sébastien Hillaire
++ vec4 sInt = (S - S * stepTransmittance) / max(extinction, 1e-7);
++ lInscattering += transmittance * sInt;
++ transmittance *= stepTransmittance;
++ }
++
++ return lInscattering;
++ }
++
++ void MAIN()
++ {
++ vec2 uv = INPUT_UV;
++ const float EYE_DISTANCE_TO_EARTH_CENTER = EARTH_RADIUS + eyeAltitude;
++
++ float azimuth = 2.0 * PI * uv.x;
++
++ // Apply a non-linear transformation to the elevation to dedicate more
++ // texels to the horizon, where having more detail matters.
++ float l = uv.y * 2.0 - 1.0;
++ float elev = l*l * sign(l) * PI * 0.5; // [-pi/2, pi/2]
++
++ vec3 rayDir = vec3(cos(elev) * cos(azimuth),
++ cos(elev) * sin(azimuth),
++ sin(elev));
++
++ vec3 rayOrigin = vec3(0.0, 0.0, EYE_DISTANCE_TO_EARTH_CENTER);
++
++ float atmosDist = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS);
++ float groundDist = raySphereIntersection(rayOrigin, rayDir, EARTH_RADIUS);
++ float tD;
++ if (eyeAltitude < ATMOSPHERE_THICKNESS) {
++ // We are inside the atmosphere
++ if (groundDist < 0.0) {
++ // No ground collision, use the distance to the outer atmosphere
++ tD = atmosDist;
++ } else {
++ // We have a collision with the ground, use the distance to it
++ tD = groundDist;
++ }
++ } else {
++ // We are in outer space
++ if (atmosDist < 0.0) {
++ // No collision with the atmosphere, just return black
++ FRAGCOLOR = vec4(0.0, 0.0, 0.0, 1.0);
++ return;
++ } else {
++ // Move the ray origin to the atmosphere intersection
++ rayOrigin = rayOrigin + rayDir * (atmosDist + 1e-3);
++ if (groundDist < 0.0) {
++ // No collision with the ground, so the ray is exiting through
++ // the atmosphere.
++ float secondAtmosDist = raySphereIntersection(
++ rayOrigin, rayDir, ATMOSPHERE_RADIUS);
++ tD = secondAtmosDist;
++ } else {
++ tD = groundDist - atmosDist;
++ }
++ }
++ }
++
++ vec4 transmittance;
++ vec4 L = computeInscattering(rayOrigin, rayDir, tD, transmittance);
++
++ #if ENABLE_SPECTRAL == 1
++ FRAGCOLOR = vec4(linearSrgbFromSpectralSamples(L), 1.0);
++ #else
++ FRAGCOLOR = vec4(L.rgb, 1.0);
++ #endif
++
++
++ }