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));
++ }
++
++ //-----------------------------------------------------------------------------
++ // Volumetric clouds
++ // Density model: Schneider, "Real-time Volumetric Cloudscapes of Horizon
++ // Zero Dawn", GPU Pro 7 (2015).
++ // Lighting model: Wrenninge & Zafar, "Oz: The Great and Volumetric"
++ // (SIGGRAPH 2013 course on Production Volume Rendering).
++
++ float remap01(float v, float lo, float hi)
++ {
++ return clamp((v - lo) / max(hi - lo, 1e-5), 0.0, 1.0);
++ }
++
++ float cloudDensityAt(vec3 pos, float heightFraction)
++ {
++ // pos is in earth-centered coordinates (km).
++ // Detail UVW — one full tile of the noise volume = cloudScale km. The XZ
++ // offset shifts the camera column off (0,0,0) in noise space so a noise
++ // valley at the origin doesn't empty the zenith view. cloudTimeOffset
++ // animates the noise slowly through its 3D volume so clouds evolve.
++ vec3 local = vec3(pos.x, pos.y - EARTH_RADIUS, pos.z);
++ vec3 baseUvw = (local + vec3(cloudWindOffset.x + 73.5, 0.0, cloudWindOffset.y + 41.2))
++ / max(cloudScale, 1e-3);
++ baseUvw.y += cloudTimeOffset;
++
++ // Domain warp breaks the periodic repetition of the noise tile.
++ vec3 warpSample = textureLod(noiseVolume, baseUvw * 0.35 + vec3(7.2, 1.1, 13.5), 0.0).gba;
++ vec3 uvw = baseUvw + (warpSample - 0.5) * 0.35;
++
++ // Coverage UVW at 5× larger period with independent offset.
++ vec3 covUvw = uvw * 0.2 + vec3(1.7, 0.31, 2.3);
++ float coverageR = textureLod(noiseVolume, covUvw, 0.0).r;
++
++ float threshold = mix(0.85, 0.0, cloudCoverage);
++ float baseShape = remap01(coverageR, threshold, threshold + 0.10);
++ if (baseShape < 0.001)
++ return 0.0;
++
++ float bottomShape = smoothstep(0.0, 0.05, heightFraction);
++ float topShape = smoothstep(1.0, 0.4, heightFraction);
++ float vshape = bottomShape * topShape;
++
++ // Worley FBM channels carve cauliflower edges into the base shape.
++ vec4 nDetail = textureLod(noiseVolume, uvw, 0.0);
++ float erosion = nDetail.g * 0.5 + nDetail.b * 0.3 + nDetail.a * 0.2;
++ float erodeStrength = mix(0.1, 0.5, smoothstep(1.0, 0.0, heightFraction));
++
++ float density = baseShape * vshape - erosion * erodeStrength;
++ density = clamp(density, 0.0, 1.0);
++
++ return density * cloudDensityScale * 1.6;
++ }
++
++ vec3 raymarchClouds(vec3 ro, vec3 rd, vec3 sunDirWorld, vec3 sunRadiance,
++ vec3 skyZenithRGB, vec3 skyGroundRGB, vec3 skySunSideRGB,
++ out float transmittance)
++ {
++ transmittance = 1.0;
++ if (cloudCoverage <= 0.0 || cloudDensityScale <= 0.0)
++ return vec3(0.0);
++ if (rd.y <= 0.0)
++ return vec3(0.0);
++
++ float cloudBottomR = EARTH_RADIUS + cloudBottomKm;
++ float cloudTopR = EARTH_RADIUS + cloudTopKm;
++ float tStart = raySphereIntersection(ro, rd, cloudBottomR);
++ float tEnd = raySphereIntersection(ro, rd, cloudTopR);
++ if (tStart < 0.0 || tEnd < 0.0 || tEnd <= tStart)
++ return vec3(0.0);
++ tEnd = min(tEnd, tStart + cloudMaxDistanceKm);
++
++ int steps = max(cloudPrimarySteps, 1);
++ int lightSteps = max(cloudLightSteps, 1);
++ float stepLen = (tEnd - tStart) / float(steps);
++ float lStepLen = (cloudTopKm - cloudBottomKm) / float(lightSteps);
++
++ // Two-octave Wrenninge multi-scattering approximation (see paper above).
++ const int MS_COUNT = 2;
++ const float MS_A = 0.1;
++ const float MS_B = 1.0;
++ const float MS_C = 0.4;
++
++ // Dual-lobe Henyey-Greenstein — backward lobe brightens front-lit clouds.
++ const float gBack = -0.3;
++ const float backLobeBlend = 0.2;
++
++ float cosThetaView = dot(rd, sunDirWorld);
++ vec3 scattering = vec3(0.0);
++
++ for (int i = 0; i < steps; ++i) {
++ float t = tStart + (float(i) + 0.5) * stepLen;
++ vec3 pos = ro + rd * t;
++ float h = length(pos) - EARTH_RADIUS;
++ float hf = (h - cloudBottomKm) / max(cloudTopKm - cloudBottomKm, 1e-3);
++
++ float distanceFade = 1.0 - smoothstep(0.5 * cloudMaxDistanceKm, cloudMaxDistanceKm, t);
++ float density = cloudDensityAt(pos, hf) * distanceFade;
++
++ if (density > 1e-4) {
++ float ext = density * cloudExtinction;
++ float stepT = exp(-ext * stepLen);
++
++ float odSun = 0.0;
++ for (int j = 0; j < lightSteps; ++j) {
++ vec3 lpos = pos + sunDirWorld * lStepLen * (float(j) + 0.5);
++ float lh = length(lpos) - EARTH_RADIUS;
++ float lhf = (lh - cloudBottomKm) / max(cloudTopKm - cloudBottomKm, 1e-3);
++ odSun += cloudDensityAt(lpos, lhf) * cloudExtinction * lStepLen;
++ }
++
++ vec3 msLight = vec3(0.0);
++ float a = 1.0, b = 1.0, c = 1.0;
++ for (int n = 0; n < MS_COUNT; ++n) {
++ float gFwd = c * cloudPhaseG;
++ float gBk = c * gBack;
++ float ggF = gFwd * gFwd;
++ float ggB = gBk * gBk;
++ float denF = 1.0 + ggF - 2.0 * gFwd * cosThetaView;
++ float denB = 1.0 + ggB - 2.0 * gBk * cosThetaView;
++ float pFwd = INV_4PI * (1.0 - ggF) / (denF * sqrt(denF));
++ float pBk = INV_4PI * (1.0 - ggB) / (denB * sqrt(denB));
++ float phaseN = mix(pFwd, pBk, backLobeBlend);
++ float sunTN = exp(-odSun * a);
++ msLight += b * phaseN * sunTN * sunRadiance;
++ a *= MS_A; b *= MS_B; c *= MS_C;
++ }
++
++ // Powder brightening for dense interiors.
++ float powderBoost = max(1.0, pow(density * 20.0, 0.5));
++ msLight *= powderBoost;
++
++ // Additive silver lining toward the sun.
++ float silver = pow(clamp(cosThetaView, 0.0, 1.0), 4.0) * (1.0 - exp(-odSun * 2.0));
++ msLight += silver * sunRadiance;
++
++ vec3 ambient = mix(skyGroundRGB, skyZenithRGB, clamp(hf, 0.0, 1.0)) * 0.6;
++ float sunFacing = clamp(cosThetaView * 0.5 + 0.5, 0.0, 1.0);
++ ambient += skySunSideRGB * sunFacing * 0.6;
++
++ vec3 stepScatter = (msLight + ambient) * (1.0 - stepT);
++ scattering += stepScatter * transmittance;
++ transmittance *= stepT;
++
++ if (transmittance < 0.01) {
++ transmittance = 0.0;
++ break;
++ }
++ }
++ }
++
++ return scattering;
++ }
++
++ /*
++ * 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;
++ }
++
++
++
++ float ease(float x, float c)
++ {
++ x = clamp(x, 0.0, 1.0);
++ if (c > 0.0) {
++ if (c < 1.0)
++ return 1.0 - pow(1.0 - x, 1.0 / c);
++ return pow(x, c);
++ }
++ if (c < 0.0) {
++ if (x < 0.5)
++ return pow(x * 2.0, -c) * 0.5;
++ return (1.0 - pow(1.0 - (x - 0.5) * 2.0, -c)) * 0.5 + 0.5;
++ }
++ return 0.0;
++ }
++
++
++ vec2 sampleSphericalMap(vec3 v)
++ {
++ // Linear azimuth/elevation
++ float azimuth = atan(v.z, v.x); // [-pi, pi]
++ float theta = asin(v.y); // [-pi/2, pi/2]
++
++ // Non-linear LUT warp: undo the warp applied during LUT generation
++ float elevWarp = sqrt(abs(theta) / (PI * 0.5)) * sign(theta) * 0.5 + 0.5;
++
++ // Map azimuth to [0,1]
++ float azimuthUv = azimuth / (2.0 * PI) + 0.5;
++
++ return vec2(azimuthUv, elevWarp);
++ }
++
++ /*
++ * ACES tonemapping fit for the sRGB color space
++ */
++ // sRGB => XYZ => D65_2_D60 => AP1 => RRT_SAT
++ const mat3 acesInputMat = mat3(
++ 0.59719, 0.07600, 0.02840,
++ 0.35458, 0.90834, 0.13383,
++ 0.04823, 0.01566, 0.83777
++ );
++
++ // ODT_SAT => XYZ => D60_2_D65 => sRGB
++ const mat3 acesOutputMat = mat3(
++ 1.60475, -0.10208, -0.00327,
++ -0.53108, 1.10813, -0.07276,
++ -0.07367, -0.00605, 1.07602
++ );
++
++ vec3 rrtAndOdtFit(vec3 v)
++ {
++ vec3 a = v * (v + 0.0245786) - 0.000090537;
++ vec3 b = v * (0.983729 * v + 0.4329510) + 0.238081;
++ return a / b;
++ }
++
++ vec3 acesFitted(vec3 color)
++ {
++ color = acesInputMat * color;
++ color = rrtAndOdtFit(color);
++ color = acesOutputMat * color;
++ return clamp(color, 0.0, 1.0);
++ }
++
++ //-----------------------------------------------------------------------------
++
++ vec3 gammaCorrect(vec3 linearSrgb)
++ {
++ vec3 a = 12.92 * linearSrgb;
++ vec3 b = 1.055 * pow(linearSrgb, vec3(1.0 / 2.4)) - 0.055;
++ vec3 c = step(vec3(0.0031308), linearSrgb);
++ return mix(a, b, c);
++ }
++
++ void MAIN()
++ {
++ vec3 direction = normalize(qt_eyeDir);
++ vec2 uv = sampleSphericalMap(direction);
++
++ vec3 col = texture(skytextureBuffer, uv).rgb;
++
++ vec3 sunDirWorld = normalize(vec3(-sunDirection.x, sunDirection.z, -sunDirection.y));
++ vec3 sunTint = sunColor.rgb * sunEnergy;
++
++ vec3 sunCol = vec3(0.0);
++ if (direction.y > 0.0) {
++ float sunDot = clamp(dot(sunDirWorld, direction), -1.0, 1.0);
++ float sunAngle = degrees(acos(sunDot));
++
++ if (sunAngle < sunDiskInnerAngle) {
++ sunCol = sunTint * sunAlpha;
++ } else if (sunAngle < sunDiskOuterAngle) {
++ float t = (sunAngle - sunDiskInnerAngle) / (sunDiskOuterAngle - sunDiskInnerAngle);
++ t = ease(t, sunDiskFalloff);
++ sunCol = sunTint * sunAlpha * (1.0 - t);
++ }
++ }
++
++ // Volumetric clouds: raymarch a cumulus layer that attenuates the sky and
++ // adds its own scattered light. See cloudDensityAt / raymarchClouds above.
++ float cloudT = 1.0;
++ vec3 cloudScatter = vec3(0.0);
++ if (cloudsEnabled) {
++ // Sun radiance feeding the cloud raymarch. Kept modest because the
++ // example has no DirectionalLight — all PBR shading on the scene
++ // comes from the SkyMaterial's IBL cubemap, and very bright cloud
++ // HDR there blows out the metallic spheres' reflections. Paired
++ // with the SceneEnvironment's probeExposure < 1 to keep the IBL
++ // energy comparable to a directional-light-anchored scene. (The
++ // linear-clamp tonemap means the visible sky still saturates to
++ // white at this scale, so this knob is effectively an IBL-only
++ // brightness control.)
++ const float SUN_RADIANCE_SCALE = 100.0;
++
++ // Cheap analytical atmospheric reddening for low sun.
++ float sunCosZenith = max(sunDirWorld.y, 0.0);
++ float slant = 1.0 / max(sunCosZenith + 0.025, 0.025);
++ vec3 sunAtten = exp(-slant * vec3(0.05, 0.15, 0.40));
++ vec3 sunRadianceLinear = sunTint * sunAtten * SUN_RADIANCE_SCALE;
++
++ // Ambient sky samples used for cloud underside / top / sun-side fill.
++ vec3 skyZenithRGB = texture(skytextureBuffer, sampleSphericalMap(vec3(0.0, 1.0, 0.0))).rgb;
++ vec3 skyGroundRGB = texture(skytextureBuffer, sampleSphericalMap(normalize(vec3(0.0, -0.25, 0.0)))).rgb;
++ vec3 sunSideDir = normalize(sunDirWorld + vec3(0.0, 0.1, 0.0));
++ vec3 skySunSideRGB = texture(skytextureBuffer, sampleSphericalMap(sunSideDir)).rgb;
++
++ vec3 ro = vec3(0.0, EARTH_RADIUS + eyeAltitude, 0.0);
++ cloudScatter = raymarchClouds(ro, direction, sunDirWorld, sunRadianceLinear,
++ skyZenithRGB, skyGroundRGB, skySunSideRGB, cloudT);
++ }
++
++ col = col * cloudT + sunCol * cloudT + cloudScatter;
++ col = max(col, vec3(0.0));
++
++ // Output linear HDR values — skybox.frag applies exposure + tonemapping.
++ // The 0.1 factor normalizes physical radiance (W/m²/sr) to display units.
++ // Matches the manual test's calibration so cloud lighting math produces
++ // the same visual result.
++ FRAGCOLOR = vec4(col * 0.1, 1.0);
++ }