#version 400 compatibility

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

    Copyright (C) 2021 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"

/* RENDERTARGETS: 3,5 */
layout(location = 0) out vec4 cloudCurrentFrame;
layout(location = 1) out vec4 skyboxData;

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

in vec2 coord;

flat in mat4x3 lightColor;

uniform sampler2D colortex5;

uniform sampler2D depthtex0;
uniform sampler2D depthtex2;

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 upDir, upDirView;
uniform vec3 sunDir, sunDirView;
uniform vec3 moonDir, moonDirView;
uniform vec3 lightDir;
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/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"

/*
float bluenoiseLookup1D() {     //surprisingly this works less good with taa
    ivec2 xy    = ivec2(gl_FragCoord.xy) & 255;
    uint z      = frameCounter & 7;

    return texelFetch(colortex2, ivec3(xy, z), 0).x;
}*/

/* ------ functions ------ */
#define cloudVolumeSamples 40       //[20 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
}
float cloudSkyPhase(float cosTheta, float g, vec3 gMult) {
    float x = mieHG(cosTheta, gMult.x * g);
    float y = mieCloud(cosTheta, -gMult.y * g) * volumeCloudData.z;
    //float z = mieCloud(cosTheta, gMult.z * g);

    return mix(x, y, 0.19);    //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;
}

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
    );
}

vec2 noise2DCubic(vec2 pos) {
    return noise2DCubic(noisetex, pos);
}

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 + vec2(0.17, 0.34));
        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;
        //shape       = thresholdStep(shape, 0.07 * hardness);

    if (shape <= 0.0) return 0.0;

        shape      -= (1.0 - erosion3D(posCurl * 20.0)) * 0.17;
        //shape       = thresholdStep(shape, 0.03 * hardness);

        shape      -= (1.0 - erosion3D((pos + curl * 0.2) * 6.0)) * 0.1 * sqr(1.0 - hardness);

        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;
/*
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 cloudCirrusShape(vec3 pos) {  //altocumulus ish
    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         = (pos * 2e-4) + wind;

        pos.xz     += noise2DCubic(noisetex, pos.xz * 0.002) * 1.0;

    float localCoverage = noisePerlinWorleyCubic(noisetex, pos.xz * 0.025 + wind.xz * 0.03);
        localCoverage = linStep(localCoverage, 0.25, 0.65);

        localCoverage = mix(0.25, 0.7, sqr(localCoverage));

    float coverageBias = localCoverage + cloudCumulusCoverageBias;

    float coverage  = noise2D(pos.xz * 0.16).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.xz += shape * 0.15;
        pos.y  += sqr(shape) * 0.2;

        shape  -= value3D(pos * 24.0) * 0.3;

        //pos.xz     += shape * 0.1;

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

        shape  -= value3D(pos * 72.0 + wind) * 0.13;

    if (shape <= 0.0) return 0.0;

        slope   = (1.0 - saturate(shape));

        //pos.xz -= shape * 0.1;

        shape  -= value2D((pos * 216.0 + wind).xz) * 0.065;

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

    return max(shape * 2.0, 0.0);
}
*/

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         = (pos * 6e-5) + wind;

    float localCoverage = noisePerlinWorleyCubic(noisetex, pos.xz * 0.03 + wind.xz * 0.03);
        localCoverage = linStep(localCoverage, 0.25, 0.65);

        localCoverage = mix(0.25, 0.66, sqr(localCoverage));

    float coverageBias = 0.4 + cloudCirrusCoverageBias + volumeCirrusData.x;

    float posBias   = noise2DCubic(noisetex, pos.xz * 0.01).x;
        posBias    += noise2DCubic(pos.xz * 0.02).x * 0.5;
        posBias    += noise2DCubic(pos.xz * 0.04).x * 0.25;
        posBias    /= 1.75;
        posBias    *= sqrt(posBias);

        pos        += posBias * vec3(0.3, 0.6, 0.5);

    float coverage  = noise2DCubic(pos.xz * 0.01).x;
        coverage   += noise2DCubic(pos.xz * 0.02).x * 0.5;
        coverage   += noise2DCubic(pos.xz * 0.04).x * 0.25;

        coverage   /= 1.75;

        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.xz += shape * 0.1;
        pos.y  -= sqr(shape) * cloudCirrusScale * 250.0;

        shape  -= value3D(pos * 32.0) * 0.15;

        //pos.xz     += shape * 0.1;

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

        shape  -= value3D(pos * 96.0 + wind) * 0.075;

    if (shape <= 0.0) return 0.0;

        slope   = (1.0 - saturate(shape));

        pos.xz -= shape * 0.1;

        shape  -= value2D((pos * 288.0 + wind).xz) * 0.035;

        shape   = max0(shape);
        shape   = cube(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, float lightNoise, 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;

    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, upDir));

        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 (maxOf(atmosFade) < 1e-4) {
                scattering     += (skyColor * sigmaT) * (integral * transmittance);
                transmittance  *= stepT;

                continue;
            }

            vec3 stepScatter    = vec3(0.0);

            #ifdef cloudVolumeSunlightDither
                float lightOD       = cloudVolumeLightOD(rPos, 6, lightNoise);
            #else
                float lightOD       = cloudVolumeLightOD(rPos, 6);
            #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 < 8; ++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 (maxOf(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 * 12.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 = 45;

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

    float lightNoise = ditherGradNoise();

    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;

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

        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, 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 (maxOf(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 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 < 5; ++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 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.03;

        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 (maxOf(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 * 38.0));     //this is loosely based on spectrum
                    float anisoPowder   = mix(powder, 1.0, pFade);

                    vec3 phaseG         = pow(vec3(0.45, 0.25, 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);
}

/* ------ sky reflection ------ */
#include "/lib/frag/capture.glsl"

/* ------ 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;
}

/* --- TEMPORAL CHECKERBOARD --- */

#define checkerboardDivider 9
#define ditherPass
#include "/lib/frag/checkerboard.glsl"

void main() {
    cloudCurrentFrame   = vec4(0.0, 0.0, 0.0, 1.0);

    #if (defined cloudVolumeEnabled || defined cloudCirrusEnabled)

    #ifdef cloudTemporalUpscaleEnabled
        const float cLOD    = sqrt(CLOUD_RENDER_LOD)*3.0;
    #else
        const float cLOD    = sqrt(CLOUD_RENDER_LOD);
    #endif

    vec2 cloudCoord = coord * cLOD;

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

        float depth = depthMax3x3(depthtex2, cloudCoord, pixelSize * cLOD);

        if (!landMask(depth)) {
            #ifdef cloudTemporalUpscaleEnabled
            int frame       = (frameCounter) % 9;
            ivec2 offset    = temporalOffset9[frame];
            ivec2 pixel     = ivec2(floor(cloudCoord * viewSize * rcp(cLOD/3.0) + offset));
            cloudCoord      = vec2(pixel) / 3.0 / viewSize * cLOD;

            float dither    = ditherBluenoiseCheckerboard(vec2(offset));
            float noise     = ditherGradNoiseCheckerboard(vec2(offset));
            #else
            float dither    = ditherBluenoise();
            float noise     = ditherGradNoiseTemporal();
            #endif

            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;

            cloudCurrentFrame = cloudSystem(worldDir, dot(viewDir, cloudLightDirView), dither, noise, skyColor);
        }
    }
    #endif

    cloudCurrentFrame   = clamp16F(cloudCurrentFrame);

    skyboxData      = texture(colortex5, coord);

    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;

        #if (defined cloudVolumeEnabled || defined cloudCirrusEnabled)

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

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

        skyboxData  = vec4(skyColor, 1.0);
    }

    skyboxData.a = 1.0;

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

    if (saturate(shadowmapCoord) == shadowmapCoord) {
        #ifdef cloudVolumeEnabled
        skyboxData.a    = generateCloudShadowmap(shadowmapCoord, lightDir, ditherBluenoiseStatic());
        #endif
    }
    #endif

    skyboxData  = clamp16F(skyboxData);
}