#version 400 compatibility

/*
====================================================================================================

    Copyright (C) 2020 RRe36

    All Rights Reserved unless otherwise explicitly stated.


    By downloading this you have agreed to the license and terms of use.
    These can be found inside the included license-file
    or here: https://rre36.com/copyright-license

    Violating these terms may be penalized with actions according to the Digital Millennium
    Copyright Act (DMCA), the Information Society Directive and/or similar laws
    depending on your country.

====================================================================================================
*/

#include "/lib/head.glsl"

#if (defined lightHqFilterEnabled && defined lightmapSmoothingEnabled)
/*DRAWBUFFERS:356*/
layout(location = 0) out vec4 lightmapData;
layout(location = 1) out vec4 skyboxData;
layout(location = 2) out vec4 auxData;
#else
/*DRAWBUFFERS:56*/
layout(location = 0) out vec4 skyboxData;
layout(location = 1) out vec4 auxData;
#endif

#define ATROUS_SIZE 32

#include "/lib/util/colorspace.glsl"

in vec2 coord;

flat in mat4x3 lightColor;

uniform sampler2D colortex3;
uniform sampler2D colortex5;
uniform sampler2D colortex6;

uniform sampler2D depthtex0;

uniform sampler2D noisetex;
uniform sampler3D depthtex1;

uniform sampler3D colortex1;
uniform sampler3D colortex2;

uniform int frameCounter;
uniform int worldTime;

uniform float eyeAltitude;
uniform float far, near;
uniform float frameTimeCounter;
uniform float wetness;
uniform float worldAnimTime;

uniform vec2 pixelSize;
uniform vec2 viewSize;
uniform vec2 taaOffset;
uniform vec2 skyCaptureResolution;

uniform vec3 cameraPosition;
uniform vec3 upvec, upvecView;
uniform vec3 sunvec, sunvecView;
uniform vec3 moonvec, moonvecView, lightvec;
uniform vec3 cloudLightDir, cloudLightDirView;

uniform vec4 daytime;

uniform mat4 gbufferModelView, gbufferModelViewInverse;
uniform mat4 gbufferProjection, gbufferProjectionInverse;


/* ------ includes ------*/
#define FUTIL_LINDEPTH
#define FUTIL_D3X3
#define FUTIL_ROT2
#include "/lib/fUtil.glsl"

#include "/lib/util/kernel.glsl"
#include "/lib/light/atrousLightmap.glsl"


#include "/lib/util/transforms.glsl"
#include "/lib/atmos/phase.glsl"
#include "/lib/frag/bluenoise.glsl"
#include "/lib/frag/gradnoise.glsl"

#define airmassStepBias 0.33
#include "/lib/atmos/air/const.glsl"
#include "/lib/atmos/air/density.glsl"

#include "/lib/atmos/project.glsl"
#include "/lib/frag/noise.glsl"
#include "/lib/util/bicubic.glsl"


/* ------ functions ------ */
#define cloudVolumeSamples 50       //[30 40 50 60 70 80 90 100]
#define cloudTransmittanceThreshold 0.05
#define cloudVolumeClip 6.5e4
#define cloudCumulusAlt 800.0       //[300.0 400.0 500.0 600.0 800.0 1000.0 1200.0 1400.0]
#define cloudCumulusDepth 2000.0    //[1000.0 1200.0 1400.0 1600.0 1800.0 2000.0 2200.0 2400.0 2600.0 2800.0 3000.0]
#define cloudCumulusScaleMult 1.0   //[0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5]
#define cloudCumulusCoverageBias 0.0    //[-0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0.0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2]

const float cloudCumulusMaxY    = cloudCumulusAlt + cloudCumulusDepth;
const float cloudCumulusMidY    = cloudCumulusAlt + cloudCumulusDepth * 0.5;
const float cloudCumulusScale   = 15e-5 * cloudCumulusScaleMult * (cloudCumulusAlt / 800.0);

#ifdef freezeAtmosAnim
    const float cloudTime   = float(atmosAnimOffset) * 0.003;
#else
    #ifdef volumeWorldTimeAnim
        float cloudTime     = worldAnimTime * 1.8;
    #else
        float cloudTime     = frameTimeCounter * 0.003;
    #endif
#endif

uniform vec3 volumeCloudData;

const float phaseConst  = 1.0 / (tau);

float mieCloud(float cosTheta, float g) {
    float sqrG  = sqr(g);
    float a     = (1.0 - sqrG) * rcp(2.0 + sqrG);
    float b     = (1.0 + sqr(cosTheta)) * rcp((-2.0 * (g * cosTheta)) + 1.0 + sqrG);

    return max((1.5 * (a * b)) + (g * cosTheta), 0.0) * phaseConst;
}

float cloudPhase(float cosTheta, float g, vec3 gMult) {
    float x = mieCloud(cosTheta, gMult.x * g);
    float y = mieCloud(cosTheta, -gMult.y * g) * volumeCloudData.z;
    float z = mieHG(cosTheta, gMult.z * g);

    return mix(mix(x, y, 0.25), z, 0.15) + 0.016 * g;    //i assume this is more energy conserving than summing them
}

vec3 curl3D(vec3 pos) {
    return texture(depthtex1, fract(pos)).xyz;
}
float erosion3D(vec3 pos) {
    return texture(colortex2, fract(pos)).x;
}
float shape3D(vec3 pos) {
    return texture(colortex1, fract(pos)).x;
}

float noisePerlinWorleyCubic(sampler2D tex, vec2 pos) {
        pos        *= 256.0;
    ivec2 location  = ivec2(floor(pos));

    float samples[4]    = float[4](
        texelFetch(tex, location                 & 255, 0).z, texelFetch(tex, (location + ivec2(1, 0)) & 255, 0).z,
        texelFetch(tex, (location + ivec2(0, 1)) & 255, 0).z, texelFetch(tex, (location + ivec2(1, 1)) & 255, 0).z
    );

    vec2 weights    = cubeSmooth(fract(pos));


    return mix(
        mix(samples[0], samples[1], weights.x),
        mix(samples[2], samples[3], weights.x), weights.y
    );
}

float cloudCumulusShape(vec3 pos, out float powderBias) {
    powderBias = 0.0;
    float altitude  = pos.y;
    float altRemapped = saturate((altitude - cloudCumulusAlt) * rcp(cloudCumulusDepth));

    float erodeLow  = 1.0 - sstep(altRemapped, 0.0, 0.22);
    float erodeHigh = sstep(altRemapped, 0.19, 1.0);
    float fadeLow   = sstep(altRemapped, 0.0, 0.2);
    float fadeHigh  = 1.0 - sstep(altRemapped, 0.6, 1.0);

    vec3 wind       = vec3(cloudTime, 0.0, cloudTime * 0.4);

        pos         = (pos * cloudCumulusScale) + wind;

    #ifdef cloudLocalCoverage
    float localCoverage = noisePerlinWorleyCubic(noisetex, pos.xz * 0.055 + wind.xz * 0.03);
        localCoverage = linStep(localCoverage, 0.35, 0.65);

        //return localCoverage;
        localCoverage = mix(0.47, 0.66, sqr(localCoverage));

    float coverageBias = localCoverage + cloudCumulusCoverageBias + volumeCloudData.x;
    #else
    float coverageBias = 0.47 + cloudCumulusCoverageBias + volumeCloudData.x;
    #endif

    float coverage  = noise2D(pos.xz * 0.24).z;
        coverage    = max0(coverage - coverageBias) * rcp(1.0 - saturate(coverageBias) * 0.999);

        coverage    = (coverage * fadeLow * fadeHigh) - erodeLow * 0.26 - erodeHigh * 0.65;

        coverage    = saturate(coverage * 1.05);

    if (coverage <= 1e-8) return 0.0;

    float wf1       = sstep(altRemapped, 0.0, 0.33);
    float wfade     = 1.0 - wf1 * 0.93;

    float dfade     = 0.001 + sstep(altRemapped, 0.0, 0.2) * 0.1;
        dfade      += sstep(altRemapped, 0.1, 0.45) * 0.4;
        dfade      += sstep(altRemapped, 0.2, 0.60) * 0.8;
        dfade      += sstep(altRemapped, 0.3, 0.85) * 0.9;
        dfade      /= 0.001 + 0.1 + 0.4 + 0.8 + 0.9;

    float shape     = coverage;

    float hardness  = sqrt(sstep(altRemapped, 0.26, 0.64));

    powderBias      = max(1.0 - hardness * 0.85, 0.36 * rcp(volumeCloudData.y));

    vec3 curl       = curl3D(pos * 4.0) * 2.0 - 1.0;
        pos        += wind * 0.3;

    vec3 posCurl    = pos + (curl * 0.12 * wfade);

        shape      -= (1.0 - shape3D(pos * 3.0)) * 0.66;

    if (shape <= 0.0) return 0.0;

        shape      -= (1.0 - erosion3D(posCurl * 20.0)) * 0.17;

        shape       = saturate(shape);
        shape       = 1.0 - pow(1.0 - shape, 3.0 + hardness * 5.0);

    return max(shape * dfade * volumeCloudData.y, 0.0);
}

float cloudCumulusShape(vec3 pos) {
    float t = 0.0;
    return cloudCumulusShape(pos, t);
}

float cloudVolumeLightOD(vec3 pos, const uint steps) {
    float basestep = 22.2;
    float exponent = 2.0;

    float stepsize  = basestep;
    float prevStep  = stepsize;

    float od = 0.0;

    for(uint i = 0; i < steps; ++i, pos += cloudLightDir * stepsize) {

        if(pos.y > cloudCumulusMaxY || pos.y < cloudCumulusAlt) continue;

            prevStep  = stepsize;
            stepsize *= exponent;
        
        float density = cloudCumulusShape(pos);
        if (density <= 0.0) continue;

            od += density * prevStep;
    }

    return od;
}
float cloudVolumeSkyOD(vec3 pos, const uint steps) {
    const vec3 dir = vec3(0.0, 1.0, 0.0);

    float stepsize = (cloudCumulusDepth / float(steps));
        stepsize  *= 1.0-linStep(pos.y, cloudCumulusAlt, cloudCumulusMaxY) * 0.9;

    float od = 0.0;

    for(uint i = 0; i < steps; ++i, pos += dir * stepsize) {

        if(pos.y > cloudCumulusMaxY || pos.y < cloudCumulusAlt) continue;
        
        float density = cloudCumulusShape(pos);
        if (density <= 0.0) continue;

            od += density * stepsize;
    }

    return od;
}

float cloudVolumeLightOD(vec3 pos, const uint steps, float noise) {
    float basestep = 22.2;
    float exponent = 2.0;

    float stepsize  = basestep;
    float prevStep  = stepsize;

    float od = 0.0;

    pos += cloudLightDir * noise * stepsize;

    for(uint i = 0; i < steps; ++i, pos += cloudLightDir * stepsize) {

        pos += cloudLightDir * noise * (stepsize - prevStep);

        if(pos.y > cloudCumulusMaxY || pos.y < cloudCumulusAlt) continue;

            prevStep  = stepsize;
            stepsize *= exponent;
        
        float density = cloudCumulusShape(pos);
        if (density <= 0.0) continue;

            od += density * prevStep;
    }

    return od;
}
float cloudVolumeSkyOD(vec3 pos, const uint steps, float noise) {
    const vec3 dir = vec3(0.0, 1.0, 0.0);

    float stepsize = (cloudCumulusDepth / float(steps));
        stepsize  *= 1.0-linStep(pos.y, cloudCumulusAlt, cloudCumulusMaxY) * 0.9;

        pos += dir * stepsize * noise;

    float od = 0.0;

    for(uint i = 0; i < steps; ++i, pos += dir * stepsize) {

        if(pos.y > cloudCumulusMaxY || pos.y < cloudCumulusAlt) continue;
        
        float density = cloudCumulusShape(pos);
        if (density <= 0.0) continue;

            od += density * stepsize;
    }

    return od;
}


/* ------ cirrus clouds ------ */
#define cloudCirrusClip 2e5
#define cloudCirrusAlt 8000.0   //[5000.0 6000.0 7000.0 8000.0 9000.0 10000.0 12000.0]
#define cloudCirrusDepth 3000.0 //[2000.0 3000.0 4000.0 5000.0 6000.0]
#define cloudCirrusScale 3e-5
#define cloudCirrusCoverageBias 0.0 //[-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5]

const float cloudCirrusMaxY     = cloudCirrusAlt + cloudCirrusDepth;
const float cloudCirrusPlaneY   = cloudCirrusAlt + cloudCirrusDepth * 0.2;

uniform vec2 volumeCirrusData;

vec2 noise2DCubic(sampler2D tex, vec2 pos) {
        pos        *= 256.0;
    ivec2 location  = ivec2(floor(pos));

    vec2 samples[4]    = vec2[4](
        texelFetch(tex, location                 & 255, 0).xy, texelFetch(tex, (location + ivec2(1, 0)) & 255, 0).xy,
        texelFetch(tex, (location + ivec2(0, 1)) & 255, 0).xy, texelFetch(tex, (location + ivec2(1, 1)) & 255, 0).xy
    );

    vec2 weights    = cubeSmooth(fract(pos));


    return mix(
        mix(samples[0], samples[1], weights.x),
        mix(samples[2], samples[3], weights.x), weights.y
    );
}

float cloudCirrusShape(vec3 pos) {
    float altitude  = pos.y;
    float altRemapped = saturate((altitude - cloudCirrusAlt) * rcp(cloudCirrusDepth));

    float erodeLow  = 1.0 - sstep(altRemapped, 0.0, 0.19);
    float erodeHigh = sstep(altRemapped, 0.22, 1.0);
    float fadeLow   = sstep(altRemapped, 0.0, 0.2);
    float fadeHigh  = 1.0 - sstep(altRemapped, 0.6, 1.0);

    vec3 wind       = vec3(cloudTime, 0.0, -cloudTime * 0.2) * 0.5;

        pos        *= rcp(1.0 + distance(cameraPosition, pos) * rcp(cloudCirrusClip));

        pos.xz     += noise2DCubic(noisetex, pos.xz * 0.004) * 1.66;

        pos.x      *= 0.6;

        pos         = (pos * cloudCirrusScale) + wind;

    float coverageBias = 0.43 + cloudCirrusCoverageBias + volumeCirrusData.x;

    float coverage  = noise2D(pos.xz * 0.23 * vec2(0.75, 1.0)).z;
        coverage    = (coverage - coverageBias) * rcp(1.0 - saturate(coverageBias) * 0.999);

        coverage    = (coverage * fadeLow * fadeHigh) - erodeLow * 0.3 - erodeHigh * 0.6;

        coverage    = saturate(coverage * 1.05);

    if (coverage <= 1e-8) return 0.0;

    float shape     = coverage;
    float slope     = sqrt(1.0 - saturate(shape));

    if (slope <= 1e-12) return 0.0;

        pos.x      += shape * 0.1;
        pos.y      -= sqr(shape) * cloudCirrusScale * 250.0;

        shape      -= value3D(pos * 48.0) * 0.14;

        pos.xz     += shape * 0.13;

        slope     = sqrt(1.0 - saturate(shape));

        shape      -= value3D(pos * 128.0 + wind) * 0.07 * slope;

    if (shape <= 0.0) return 0.0;

        slope    = (1.0 - saturate(shape));

        pos.x   *= pi;

        shape   -= value3D(pos * 256.0 + wind) * 0.04 * slope;

        shape    = max0(shape);
        shape    = sqr(shape) * sqrt(shape);

    return max(shape * volumeCirrusData.y, 0.0);
}

float cloudCirrusLightOD(vec3 pos, const uint steps, vec3 dir) {
    float stepsize = (cloudCirrusDepth / float(steps));
        stepsize  *= 1.0 - linStep(pos.y, cloudCirrusAlt, cloudCirrusMaxY) * 0.9;

    float od = 0.0;

    for(uint i = 0; i < steps; ++i, pos += dir * stepsize) {

        if(pos.y > cloudCirrusMaxY || pos.y < cloudCirrusAlt) continue;
        
        float density = cloudCirrusShape(pos);
        if (density <= 0.0) continue;

            od += density * stepsize;
    }

    return od;
}

/* ------ cloud system function ------ */

vec4 cloudSystem(vec3 worldDir, float vDotL, float dither, vec3 skyColor) {
    vec3 totalScattering    = vec3(0.0);
    float totalTransmittance = 1.0;

    vec3 sunlight       = (worldTime>23000 || worldTime<12900) ? lightColor[0] : lightColor[2];
    vec3 skylight       = lightColor[1];

    float pFade         = saturate(mieCS(vDotL, 0.5));
        pFade           = mix(pFade, vDotL * 0.5 + 0.5, 1.0 / sqrt2);

    // --- volumetric clouds --- //
    #ifdef cloudVolumeEnabled
    float within    = sstep(eyeAltitude, cloudCumulusAlt - 75.0, cloudCumulusAlt) * (1.0 - sstep(eyeAltitude, cloudCumulusMaxY, cloudCumulusMaxY + 75.0));
    bool isBelowVol = eyeAltitude < cloudCumulusMidY;
    bool visibleVol = worldDir.y > 0.0 && isBelowVol || worldDir.y < 0.0 && !isBelowVol;

    float lightNoise = ditherGradNoiseTemporal();

    if (visibleVol || within > 0.0) {

        vec3 bottom     = worldDir * ((cloudCumulusAlt - eyeAltitude) * rcp(worldDir.y));
        vec3 top        = worldDir * ((cloudCumulusMaxY - eyeAltitude) * rcp(worldDir.y));

        float airDist   = length(worldDir.y < 0.0 ? bottom : top);
            airDist     = min(airDist, cloudVolumeClip);

        vec3 airmass    = getAirmass(vec3(0.0, planetRad + eyeAltitude, 0.0), worldDir, airDist, 0.15, 3) * rcp(airDist);

        if (worldDir.y < 0.0 && isBelowVol || worldDir.y > 0.0 && !isBelowVol) {
            bottom      = vec3(0.0);
            top         = vec3(0.0);
        }

        vec3 start      = isBelowVol ? bottom : top;
        vec3 end        = isBelowVol ? top : bottom;
            start       = mix(start, gbufferModelViewInverse[3].xyz, within);
            end         = mix(end, worldDir * 6e4, within);

        const float baseLength  = cloudCumulusDepth / cloudVolumeSamples;
        float stepLength    = length(end - start) * rcp(float(cloudVolumeSamples));
        float stepCoeff     = stepLength * rcp(baseLength);
            stepCoeff       = 0.45 + clamp(stepCoeff - 1.1, 0.0, 5.0) * 0.5;
            stepCoeff       = mix(stepCoeff, 10.0, sqr(within));

        uint steps          = uint(cloudVolumeSamples * stepCoeff);

        vec3 rStep          = (end - start) * rcp(float(steps));
        vec3 rPos           = rStep * dither + start + cameraPosition;
        float rLength       = length(rStep);

        vec3 scattering     = vec3(0.0);
        float transmittance = 1.0;

        vec3 bouncelight    = vec3(0.6, 1.0, 0.8) * sunlight * rcp(pi * 14.0 * sqrt2) * max0(dot(cloudLightDir, upvec));

        const float sigmaA  = 1.0;
        const float sigmaT  = 0.066;

        for (uint i = 0; i < steps; ++i, rPos += rStep) {
            if (transmittance < cloudTransmittanceThreshold) break;
            if (rPos.y < cloudCumulusAlt || rPos.y > cloudCumulusMaxY) continue;

            float dist  = distance(rPos, cameraPosition);
            if (dist > cloudVolumeClip) continue;

            float powderBias;
            float density = cloudCumulusShape(rPos, powderBias);
            if (density <= 0.0) continue;

            float extinction    = density * sigmaT;
            float stepT         = exp(-extinction * rLength);
            float integral      = (1.0 - stepT) * rcp(sigmaT);

            //float airmassMult = 1.0 + sstep(dist / cloudVolumeClip, 0.5, 1.0) * pi; 
            vec3 atmosFade = expf(-max(airScatterMat * airmass.xy * (dist), 0.0));
            //float atmosFade = exp(-max0(dist * 2e-5));
                atmosFade   = mix(atmosFade, vec3(0.0), sqr(linStep(dist / cloudVolumeClip, 0.5, 1.0)));

            if (max3(atmosFade) < 1e-4) {
                scattering     += (skyColor * sigmaT) * (integral * transmittance);
                transmittance  *= stepT;

                continue;
            }

            vec3 stepScatter    = vec3(0.0);

            #ifdef cloudVolumeSunlightDither
                float lightOD       = cloudVolumeLightOD(rPos, 5, lightNoise);
            #else
                float lightOD       = cloudVolumeLightOD(rPos, 5);
            #endif
            
            #ifdef cloudVolumeSkylightDither
                float skyOD         = cloudVolumeSkyOD(rPos, 4, lightNoise);
            #else
                float skyOD         = cloudVolumeSkyOD(rPos, 4);
            #endif

            float bounceOD      = cube(1.0 - linStep(rPos.y, cloudCumulusAlt, cloudCumulusAlt + cloudCumulusDepth * 0.3));
            /*
            float powder        = 8.0 * (1.0 - 0.97 * expf(-extinction * 16.0 * powderBias));     //this is loosely based on spectrum
            float anisoPowder   = powder * (1.0 - expf(-max0(lightOD - sigmaT) * powderBias * (1.0 + sigmaT)) * rpi);
                anisoPowder     = mix(anisoPowder, 1.0, pFade);
            */
            float powderCoeff   = 16.0 * powderBias;

            float powderLo1 = 1.0 - exp(-extinction *  2.0 * powderCoeff);
            float powderLo2 = 1.0 - exp(-extinction *  4.0 * powderCoeff);
            float powderLo  = (powderLo1 + powderLo2) / 4.0;

            float powderHi1 = 1.0 - exp(-extinction *  5.0 * powderCoeff);
            float powderHi2 = 1.0 - exp(-extinction * 20.0 * powderCoeff);
            float powderHi  = (powderHi1 + powderHi2) / 6.0;

            float powder    = (powderLo + powderHi) * 2.2;
            float anisoPowder = mix(powder, 1.0, pFade);

            vec3 phaseG         = pow(vec3(0.45, 0.45, 0.95), vec3(1.0 + lightOD));

            for (uint j = 0; j < 6; ++j) {
                float n     = float(j);

                float td    = sigmaT * pow(0.5, n);
                float phase = cloudPhase(vDotL, pow(0.5, n), phaseG);

                stepScatter.x  += expf(-lightOD * td) * phase * td;
                stepScatter.y  += expf(-skyOD * td) * td;
            }

            stepScatter.xy *= vec2(anisoPowder, powder);
            stepScatter.z  += bounceOD * powder;

            stepScatter     = (sunlight * stepScatter.x) + (skylight * stepScatter.y) + (bouncelight * stepScatter.z);
            stepScatter     = mix(skyColor * sigmaT, stepScatter, atmosFade);
            scattering     += stepScatter * (integral * transmittance);

            transmittance  *= stepT;
        }

        transmittance       = linStep(transmittance, cloudTransmittanceThreshold, 1.0);

        totalScattering    += scattering;
        totalTransmittance *= transmittance;
    }
    #endif

    #ifdef cloudCirrusEnabled
    // --- planar cirrus clouds --- //
    bool isBelowPlane   = eyeAltitude < cloudCirrusPlaneY;
    bool visiblePlane   = worldDir.y > 0.0 && isBelowPlane || worldDir.y < 0.0 && !isBelowPlane;

    if (visiblePlane && (!isBelowPlane || totalTransmittance > 1e-10)) {
        vec3 plane      = worldDir * ((cloudCirrusPlaneY - eyeAltitude) * rcp(worldDir.y));

        float airDist   = length(plane);
            //airDist     = mix(min(airDist, cloudVolumeClip), 9e4, within);

        vec3 airmass    = getAirmass(vec3(0.0, planetRad + eyeAltitude, 0.0), worldDir, airDist, 0.25, 3) * rcp(airDist);

        vec3 rPos       = plane + cameraPosition;

        const float sigmaT  = 0.01;

        float dist      = distance(rPos, cameraPosition);

        vec3 scattering = vec3(0.0);
        float transmittance = 1.0;

        mat2x3 planarBounds = mat2x3(worldDir*((cloudCirrusMaxY-eyeAltitude)/worldDir.y),
                                     worldDir*((cloudCirrusAlt-eyeAltitude)/worldDir.y));

        float rLength   = distance(planarBounds[0], planarBounds[1]);

        if (dist < cloudCirrusClip) {
            float density   = cloudCirrusShape(rPos);

            if (density > 0.0) {
                float extinction    = density * sigmaT;
                float stepT         = exp(-extinction * rLength);
                float integral      = (1.0 - stepT) * rcp(sigmaT);

                vec3 atmosFade  = expf(-max(airScatterMat * airmass.xy * dist, 0.0));
                    atmosFade   = mix(atmosFade, vec3(0.0), sqr(linStep(dist / cloudCirrusClip, 0.5, 1.0)));

                if (max3(atmosFade) < 1e-4) {
                    scattering     += (skyColor * sigmaT) * (integral * transmittance);
                    transmittance  *= stepT;
                } else {
                    float lightOD       = cloudCirrusLightOD(rPos, 5, cloudLightDir);
                    float skyOD         = cloudCirrusLightOD(rPos, 4, vec3(0.0, 1.0, 0.0));

                    float powder        = 8.0 * (1.0 - 0.97 * expf(-extinction * 18.0));     //this is loosely based on spectrum
                    float anisoPowder   = mix(powder, 1.0, pFade);

                    vec3 phaseG         = pow(vec3(0.45, 0.45, 0.95), vec3(1.0 + lightOD));

                    for (uint j = 0; j < 6; ++j) {
                        float n     = float(j);

                        float td    = sigmaT * pow(0.5, n);
                        float phase = cloudPhase(vDotL, pow(0.5, n), phaseG);

                        scattering.x   += expf(-lightOD * td) * anisoPowder * phase * td;
                        scattering.y   += expf(-skyOD * td) * powder * td;
                    }

                    scattering      = (sunlight * scattering.x) + (skylight * scattering.y);
                    scattering      = mix(skyColor * sigmaT, scattering, atmosFade);
                    scattering      = scattering * (integral * transmittance);

                    transmittance  *= stepT;
                }
            }
        }

        if (isBelowPlane) {
            totalScattering    += scattering * totalTransmittance;
            totalTransmittance *= transmittance;
        } else {
            totalScattering     = totalScattering * transmittance + scattering;
            totalTransmittance *= transmittance;
        }
    }
    #endif

    return vec4(totalScattering, totalTransmittance);
}

vec4 cloudSystemReflected(vec3 worldDir, float vDotL, float dither, vec3 skyColor) {
    vec3 totalScattering    = vec3(0.0);
    float totalTransmittance = 1.0;

    vec3 sunlight       = (worldTime>23000 || worldTime<12900) ? lightColor[0] : lightColor[2];
    vec3 skylight       = lightColor[1];

    float pFade         = saturate(mieCS(vDotL, 0.5));
        pFade           = mix(pFade, vDotL * 0.5 + 0.5, 1.0 / sqrt2);

    const float eyeAltitude = 64.0;
    vec3 camPos     = vec3(cameraPosition.x, eyeAltitude, cameraPosition.z);
    const uint volumeSamples = 40;

    // --- volumetric clouds --- //
    #ifdef cloudVolumeEnabled
    bool visibleVol = worldDir.y > 0.0;

    if (visibleVol) {

        vec3 bottom     = worldDir * ((cloudCumulusAlt - eyeAltitude) * rcp(worldDir.y));
        vec3 top        = worldDir * ((cloudCumulusMaxY - eyeAltitude) * rcp(worldDir.y));

        float airDist   = length(worldDir.y < 0.0 ? bottom : top);
            airDist     = min(airDist, cloudVolumeClip);

        vec3 airmass    = getAirmass(vec3(0.0, planetRad + eyeAltitude, 0.0), worldDir, airDist, 0.15, 2) * rcp(40000.0) * sqrt2;

        vec3 start      = bottom;
        vec3 end        = top;

        const float baseLength  = cloudCumulusDepth / float(volumeSamples);
        float stepLength    = length(end - start) * rcp(float(volumeSamples));
        float stepCoeff     = stepLength * rcp(baseLength);
            stepCoeff       = 0.45 + clamp(stepCoeff - 1.1, 0.0, 2.5) * 0.5;

        uint steps          = uint(volumeSamples * stepCoeff);

        vec3 rStep          = (end - start) * rcp(float(steps));
        vec3 rPos           = rStep * dither + start + camPos;
        float rLength       = length(rStep);

        vec3 scattering     = vec3(0.0);
        float transmittance = 1.0;

        const float sigmaA  = 1.0;
        const float sigmaT  = 0.066;

        for (uint i = 0; i < steps; ++i, rPos += rStep) {
            if (transmittance < 0.1) break;
            if (rPos.y < cloudCumulusAlt || rPos.y > cloudCumulusMaxY) continue;

            float dist  = distance(rPos, camPos);
            if (dist > cloudVolumeClip) continue;

            float powderBias;
            float density = cloudCumulusShape(rPos, powderBias);
            if (density <= 0.0) continue;

            float extinction    = density * sigmaT;
            float stepT         = exp(-extinction * rLength);
            float integral      = (1.0 - stepT) * rcp(sigmaT);

            float airmassMult = 1.0 + sstep(dist / cloudVolumeClip, 0.5, 1.0) * pi; 
            vec3 atmosFade = expf(-max(airScatterMat * airmass.xy * (dist * airmassMult), 0.0));

            if (max3(atmosFade) < 1e-4) {
                scattering     += (skyColor * sigmaT) * (integral * transmittance);
                transmittance  *= stepT;

                continue;
            }

            vec3 stepScatter    = vec3(0.0);

            float lightOD       = cloudVolumeLightOD(rPos, 4);
            
            float skyOD         = cloudVolumeSkyOD(rPos, 3);

            /*
            float powder        = 8.0 * (1.0 - 0.97 * expf(-extinction * 16.0 * powderBias));     //this is loosely based on spectrum
            float anisoPowder   = powder * (1.0 - expf(-max0(lightOD - sigmaT) * powderBias * (1.0 + sigmaT)) * rpi);
                anisoPowder     = mix(anisoPowder, 1.0, pFade);
            */
            float powderCoeff   = 16.0 * powderBias;

            float powderLo1 = 1.0 - exp(-extinction *  2.0 * powderCoeff);
            float powderLo2 = 1.0 - exp(-extinction *  4.0 * powderCoeff);
            float powderLo  = (powderLo1 + powderLo2) / 4.0;

            float powderHi1 = 1.0 - exp(-extinction *  5.0 * powderCoeff);
            float powderHi2 = 1.0 - exp(-extinction * 20.0 * powderCoeff);
            float powderHi  = (powderHi1 + powderHi2) / 6.0;

            float powder    = (powderLo + powderHi) * 2.2;
            float anisoPowder = mix(powder, 1.0, pFade);

            vec3 phaseG         = pow(vec3(0.45, 0.45, 0.95), vec3(1.0 + lightOD / pi));

            for (uint j = 0; j < 4; ++j) {
                float n     = float(j);

                float td    = sigmaT * pow(0.5, n);
                float phase = cloudPhase(vDotL, pow(0.5, n), phaseG);

                stepScatter.x  += expf(-lightOD * td) * phase * td;
                stepScatter.y  += expf(-skyOD * td) * td;
            }

            stepScatter.xy *= vec2(anisoPowder, powder);

            stepScatter     = (sunlight * stepScatter.x) + (skylight * stepScatter.y);
            stepScatter     = mix(skyColor * sigmaT, stepScatter, atmosFade);
            scattering     += stepScatter * (integral * transmittance);

            transmittance  *= stepT;
        }

        transmittance       = linStep(transmittance, 0.1, 1.0);

        totalScattering    += scattering;
        totalTransmittance *= transmittance;
    }
    #endif

    #ifdef cloudCirrusEnabled
    // --- planar cirrus clouds --- //
    bool visiblePlane   = worldDir.y > 0.0;

    if (visiblePlane && totalTransmittance > 1e-10) {
        vec3 plane      = worldDir * ((cloudCirrusPlaneY - eyeAltitude) * rcp(worldDir.y));

        float airDist   = length(plane);
            //airDist     = mix(min(airDist, cloudVolumeClip), 9e4, within);

        vec3 airmass    = getAirmass(vec3(0.0, planetRad + eyeAltitude, 0.0), worldDir, airDist, 0.25, 3) * rcp(60000.0);

        vec3 rPos       = plane + camPos;

        const float sigmaT  = 0.016;

        float dist      = distance(rPos, camPos);

        vec3 scattering = vec3(0.0);
        float transmittance = 1.0;
        const float rLength = cloudCirrusDepth / tau;

        if (dist < cloudCirrusClip) {
            float density   = cloudCirrusShape(rPos);

            if (density > 0.0) {
                float extinction    = density * sigmaT;
                float stepT         = exp(-extinction * rLength);
                float integral      = (1.0 - stepT) * rcp(sigmaT);

                vec3 atmosFade  = expf(-max(airScatterMat * airmass.xy * dist, 0.0));

                if (max3(atmosFade) < 1e-4) {
                    scattering     += (skyColor * sigmaT) * (integral * transmittance);
                    transmittance  *= stepT;
                } else {
                    float lightOD       = cloudCirrusLightOD(rPos, 4, cloudLightDir);
                    float skyOD         = cloudCirrusLightOD(rPos, 3, vec3(0.0, 1.0, 0.0));

                    float powder        = 8.0 * (1.0 - 0.97 * expf(-extinction * 18.0));     //this is loosely based on spectrum
                        powder          = powder * 0.5 + 0.5;
                    float anisoPowder   = mix(powder, 1.0, pFade);

                    vec3 phaseG         = pow(vec3(0.45, 0.45, 0.95), vec3(1.0 + lightOD));

                    for (uint j = 0; j < 5; ++j) {
                        float n     = float(j);

                        float td    = sigmaT * pow(0.5, n);
                        float phase = cloudPhase(vDotL, pow(0.5, n), phaseG);

                        scattering.x   += expf(-lightOD * td) * anisoPowder * phase * td;
                        scattering.y   += expf(-skyOD * td) * powder * td;
                    }

                    scattering      = (sunlight * scattering.x) + (skylight * scattering.y);
                    scattering      = mix(skyColor * sigmaT, scattering, atmosFade);
                    scattering      = scattering * (integral * transmittance);

                    transmittance  *= stepT;
                }
            }
        }

        totalScattering    += scattering * totalTransmittance;
        totalTransmittance *= transmittance;
    }
    #endif

    return vec4(totalScattering, totalTransmittance);
}

#ifdef lightHqFilterEnabled
 /* - */
#endif

/* ------ sky reflection ------ */
vec2 sincos(float x) {
    return vec2(sin(x), cos(x));
}

vec3 projectSphere(vec2 coord) {
    coord  *= vec2(tau, pi);
    vec2 lon = sincos(coord.x) * sin(coord.y);
    return vec3(lon.x, cos(coord.y), lon.y);
}

/* ------ cloud shadows ------ */

uniform mat4 shadowModelViewInverse;

float generateCloudShadowmap(vec2 coord, vec3 lightDir, float dither) {
    vec3 position   = vec3(coord, 0.0) * 2.0 - 1.0;
        position.xy *= cloudShadowmapRenderDistance;
        position    = mat3(shadowModelViewInverse) * position;
        position.xz += cameraPosition.xz;
        //position    += lightDir * (256.0 - position.y) * rcp(lightDir.y);

    float lightFade = cubeSmooth(sqr(linStep(lightDir.y, 0.1, 0.15)));

    float transmittance = 1.0;

    if (lightFade > 0.0) {
        vec3 bottom     = lightDir * ((cloudCumulusAlt - position.y) * rcp(lightDir.y));
        vec3 top        = lightDir * ((cloudCumulusMaxY - position.y) * rcp(lightDir.y));

        vec3 start      = bottom;
        vec3 end        = top;

        uint steps          = 40;

        vec3 rStep          = (end - start) * rcp(float(steps));
        vec3 rPos           = rStep * dither + start + position;
        float rLength       = length(rStep);

        const float sigmaT  = 0.066;

        for (uint i = 0; i < steps; ++i, rPos += rStep) {
            if (transmittance < 0.05) break;

            float fade      = 1.0 - sstep(distance(rPos.xz, cameraPosition.xz), cloudShadowmapRenderDistance * 0.5, cloudShadowmapRenderDistance);
            if (fade < 1e-3) {
                transmittance = 1.0 - wetness * 0.90;
                continue;
            }

            float density   = cloudCumulusShape(rPos);

            float extinction = density * sigmaT;
            float stepT     = exp(-extinction * rLength);
                stepT       = mix(1.0 - wetness * 0.90, stepT, fade);

            transmittance  *= stepT;
        }
        transmittance       = linStep(transmittance, 0.05, 1.0);
        transmittance       = mix(1.0 - wetness * 0.95, transmittance, lightFade);
    } else {
        transmittance   = 1.0 - wetness * 0.95;
    }

    return transmittance;
}

void main() {
    #if (defined lightHqFilterEnabled && defined lightmapSmoothingEnabled)
    lightmapData    = clamp16F(lightmap3x3(colortex3, depthtex0, colortex6, coord, ATROUS_SIZE));
    #endif
    skyboxData      = texture(colortex5, coord);

    #if (defined cloudVolumeEnabled || defined cloudCirrusEnabled)
    const float clod = sqrt(CLOUD_RENDER_LOD);

    vec2 cloudCoord = (coord - vec2(1.0 - rcp(clod), 0.0)) * clod;

    if (clamp(cloudCoord, -pixelSize * clod, 1.0 + pixelSize * clod) == cloudCoord) {
        cloudCoord  = saturate(cloudCoord);

        float depth = depthMax3x3(depthtex0, cloudCoord, pixelSize * clod);

        if (!landMask(depth)) {
            vec3 viewPos    = screenToViewSpace(vec3(cloudCoord, 1.0));
            vec3 viewDir    = normalize(viewPos);
            vec3 scenePos   = viewToSceneSpace(viewPos);
            vec3 worldDir   = normalize(scenePos);

            vec3 skyColor   = textureBicubic(colortex5, projectSky(worldDir)).rgb;

            skyboxData      = cloudSystem(worldDir, dot(viewDir, cloudLightDirView), ditherBluenoise(), skyColor);
        } else {
            skyboxData.a    = 1.0;
        }
    }
    #endif

    vec2 reflectionCoord    = (coord - vec2(0.0, exp2(-SKY_RENDER_LOD) * 2.0)) * skyCaptureResolution;

    if (isOnScreenDownscale(reflectionCoord, pixelSize, 4)) {
        vec2 coord          = saturate(reflectionCoord);
        vec3 worldDir       = projectSphere(coord) * vec3(1.0, 1.0, -1.0);

        vec3 skyColor       = texture(colortex5, projectSky(worldDir)).rgb;

        #ifdef cloudReflectionsToggle
        vec4 cloudData      = cloudSystemReflected(worldDir, dot(mat3(gbufferModelView) * worldDir, cloudLightDirView), ditherBluenoiseStatic(), skyColor);

        skyColor    = skyColor * cloudData.a + cloudData.rgb;
        #endif

        skyboxData  = vec4(skyColor, 1.0);
    }

    skyboxData  = clamp16F(skyboxData);

    auxData     = stex(colortex6);

    #ifdef cloudShadowsEnabled
    vec2 shadowmapCoord     = coord * max(viewSize / cloudShadowmapResolution, 1.0);

    if (saturate(shadowmapCoord) == shadowmapCoord) {
        auxData.a           = generateCloudShadowmap(shadowmapCoord, lightvec, ditherBluenoiseStatic());
    }
    #endif
}