#ifndef CLOUD_GLSL
    #define CLOUD_GLSL a

    const float cloudHighRealHeight = EarthRadius + CLOUD_HIGH_ALTITUDE;
    const float cloudHighEpsilon = 0.001;

    #define CloudScatteringCoefficient 50.0
    #define CloudAbsorptionCoefficient 20.0
    #define CloudMultiScatteringD 0.5

    float sampleCloudDensity(vec3 pos, float cloudBottomHeight, float cloudThickness) {
        float frequency = mix(CLOUD_CLEAR_SCALE, CLOUD_RAIN_SCALE, rainStrength);
        float weight = 1.0;
        float time = 0.02 * frameTimeCounter;

        vec3 worldPos = pos;
        worldPos.xz += cameraPosition.xz * 0.001;

        float n = 0.0;
        float c = 0.0;
        for(int i = 0; i < 4; i++) {
            n += perlin2D(worldPos.xz * frequency + time) * weight;
            c += weight;

            frequency *= 2.0;
            weight *= 0.5;
            time *= 1.1;
        }
        n /= c;

        float cloudmap = saturate((sample2DNoise(worldPos.xz * 0.2 + frameTimeCounter * 0.01) * 3.0 - 0.5) * 0.5 + 0.5);
        n *= cloudmap;

        float coverage = mix(CLOUD_CLEAR_COVERAGE, CLOUD_RAIN_COVERAGE, rainStrength);
        float bias = 1.0 - coverage;
        n = max0(n - bias) / max(1.0 - bias, 1e-4);
        n = saturate(n);

        float frequency1 = 5.0;
        float weight1 = 1.0;
        float time1 = 0.1 * frameTimeCounter;

        float n1 = 0.0;
        float c1 = 0.0;
        for(int i = 0; i < 4; i++) {
            n1 += sample3DNoise(worldPos * frequency1 + time1) * weight1;
            c1 += weight1;

            frequency1 *= 2.5;
            weight1 *= 0.5;
            time1 *= 1.1;
        }
        n1 /= c1;

        if(n <= 0.0) return 0.0;

        float hRaw = (length(pos) - cloudBottomHeight) / cloudThickness;

        float heightJitter = (n1 - 0.5) * 0.6;
        float h = saturate(hRaw + heightJitter);

        float s = smoothstep(0.0, 1.0, h);
        s = s * 0.8 + 0.25;

        n = remapSaturate(n, n1 * s * 0.6, 1.0, 0.0, 1.0);

        float bottomFade = pow(saturate(hRaw + heightJitter * 1.5), 0.35);
        float topFade = pow(saturate(1.0 - hRaw + heightJitter * 0.5), 0.6);
        n *= bottomFade * topFade;

        float density = mix(CLOUD_CLEAR_DENSITY, CLOUD_RAIN_DENSITY, rainStrength);
        return saturate(n) * density;
    }

    float CloudOpticalDepth(vec3 pos, vec3 dir, float len, float currDensity, int N_SAMPLE, float cloudBottomHeight, float cloudThickness) {
        if(len <= 0.0) return 0.0;

        float ds = len / float(N_SAMPLE);
        vec3 samplePos = pos;
        vec3 stepVec = dir * ds;

        float density0 = currDensity;
        float opticalDepth = 0.0;
        for(int i = 0; i < N_SAMPLE; i++) {
            samplePos += stepVec;

            float density = sampleCloudDensity(samplePos, cloudBottomHeight, cloudThickness);

            opticalDepth += density0 + density;
            density0 = density;
        }

        return opticalDepth * 0.5 * ds;
    }

    float CloudTransmittance(vec3 pos, vec3 dir, float currDensity, int N_SAMPLE, float cloudTopHeight, float cloudBottomHeight, float cloudThickness) {
        float lenToCloudTop = rayIntersectSphere(pos, dir, cloudTopHeight * cloudTopHeight);

        float len = min(lenToCloudTop, cloudThickness * 0.5);
        float optic = CloudOpticalDepth(pos, dir, len, currDensity, N_SAMPLE, cloudBottomHeight, cloudThickness);
        return exp(-(CloudScatteringCoefficient + CloudAbsorptionCoefficient) * optic);
    }

    bool CameraInsideCloudLayer() {
        float cloudAltitude = mix(CLOUD_CLEAR_ALTITUDE, CLOUD_RAIN_ALTITUDE, rainStrength);
        float cloudThickness = mix(CLOUD_CLEAR_THICKNESS, CLOUD_RAIN_THICKNESS, rainStrength);
        float cloudBottomHeight = EarthRadius + cloudAltitude;
        float cloudTopHeight = cloudBottomHeight + cloudThickness;
        float margin = cloudThickness * 0.25;
        return cameraHeight > cloudBottomHeight - margin && cameraHeight < cloudTopHeight + margin;
    }

    float GetAttenuationProbability(float opticalDepth) {
        return exp(-opticalDepth);
    }
    float GetAttenuationProbability(float opticalDepth, float secondSpread, float secondIntensity) {
        return max(exp(-opticalDepth), exp(-opticalDepth * secondSpread) * secondIntensity);
    }

    float GetInScatterProbability(float heightFraction, float opticalDepth) {
        float heightFactor = remapSaturate(heightFraction, 0.0, 1.0, 0.7, 1.0);
        float depthFactor = remapSaturate(opticalDepth, 0.0, 1.0, 0.25, 1.0);
        return 0.05 + pow(saturate(opticalDepth), depthFactor * heightFactor);
    }

    float GetSilverPhase(float cosTheta, float baseEcc, float silverIntensity, float silverSpread) {
        return max(getPhase(cosTheta, baseEcc), silverIntensity * getPhase(cosTheta, 0.99 - silverSpread));
    }

    vec4 RenderCloud(vec3 pos, vec3 dir, vec3 lightDir, vec3 lightLuminance, vec3 skylight) {
        float cloudAltitude = mix(CLOUD_CLEAR_ALTITUDE,  CLOUD_RAIN_ALTITUDE, rainStrength);
        float cloudThickness = mix(CLOUD_CLEAR_THICKNESS, CLOUD_RAIN_THICKNESS, rainStrength);
        float cloudBottomHeight = EarthRadius + cloudAltitude;
        float cloudTopHeight = cloudBottomHeight + cloudThickness;

        float lenToCloudBottom = rayIntersectSphere(pos, dir, cloudBottomHeight * cloudBottomHeight);
        float lenToCloudTop = rayIntersectSphere(pos, dir, cloudTopHeight * cloudTopHeight);

        float lenToCloud = 0.0;
        float lenToCloudEnd = 0.0;
        if(pos.y < cloudBottomHeight) {
            float lenToEarth = rayIntersectSphere(pos, dir, pow2(EarthRadius + ReferenceHeight));
            if(lenToEarth > 0.0) {
                return vec4(vec3(0.0), 1.0);
            } else {
                lenToCloud = max0(lenToCloudBottom);
                lenToCloudEnd = max0(lenToCloudTop);
            }
        } else if(pos.y < cloudTopHeight) {
            lenToCloud = 0.0;
            if(lenToCloudBottom > 0.0) {
                lenToCloudEnd = max0(lenToCloudBottom);
            } else {
                lenToCloudEnd = max0(lenToCloudTop);
            }
        } else {
            lenToCloud = max0(lenToCloudTop);
            lenToCloudEnd = max0(lenToCloudBottom);
        }
        float len = max0(lenToCloudEnd - lenToCloud);

        vec3 cloudTopPos = pos + dir * lenToCloudTop;
        lightLuminance *= remapSaturate(lightDir.y, 0.02, 0.1, 0.0, 1.0);
        vec3 lightColor = lightLuminance * TransToAtmos(cloudTopPos, lightDir);
        vec3 ambientColor = saturation(skylight, 0.33) * 3.0;
        ambientColor *= remapSaturate(lightDir.y, 0.0, 0.25, 0.05, 1.0);

        float cosTheta = dot(lightDir, dir);

        float phase = GetSilverPhase(cosTheta, 0.5, 0.5, 0.25);
        float rearPhase = getPhase(cosTheta, -0.5);
        phase = mix(phase, rearPhase, 0.15);
        float u_phase = 1.0 / (4.0 * PI);

        int sampleCount = max(1, int(mix(72.0, 18.0, abs(dir.y)) * CLOUD_QUALITY));

        float ds = len / float(sampleCount);
        vec3 stepVec = dir * ds;
        vec3 samplePos = pos + dir * lenToCloud;

        vec3 jitter = stepVec * blueNoise;
        samplePos += jitter;

        float transmittance = 1.0;
        vec3 scattering = vec3(0.0);
        for(int i = 0; i < sampleCount; i++) {
            float stepCloudDensity = sampleCloudDensity(samplePos, cloudBottomHeight, cloudThickness);
            float stepTransmittance = 1.0;

            if(stepCloudDensity > 1e-5) {
                float sigmaS = stepCloudDensity * CloudScatteringCoefficient;
                float sigmaA = stepCloudDensity * CloudAbsorptionCoefficient;
                float sigmaE = sigmaS + sigmaA;
                stepTransmittance = exp(-sigmaE * ds);

                float lightOpticalDepth = CloudOpticalDepth(samplePos, lightDir, min(rayIntersectSphere(samplePos, lightDir, cloudTopHeight * cloudTopHeight), cloudThickness * 0.5), stepCloudDensity, 3, cloudBottomHeight, cloudThickness) * (CloudScatteringCoefficient + CloudAbsorptionCoefficient);
                float lightVisibility = GetAttenuationProbability(lightOpticalDepth, 0.35, 0.6);

                vec3 lightEnergy = lightColor * lightVisibility * sigmaS;

                float D = CloudMultiScatteringD;
                float f_ms = (sigmaS / sigmaE) * (1.0 - exp(-D * sigmaE));
                vec3 multiScattering = lightEnergy * f_ms / max(1.0 - f_ms, 1e-6) * u_phase;

                float ambientVisibility = CloudTransmittance(samplePos, vec3(0.0, 1.0, 0.0), stepCloudDensity, 2, cloudTopHeight, cloudBottomHeight, cloudThickness) * 0.5 + 0.5;
                vec3 ambientlight = ambientColor * ambientVisibility * sigmaS;

                float h = saturate((length(samplePos) - cloudBottomHeight) / cloudThickness);
                float inScatterBoost = GetInScatterProbability(h, lightOpticalDepth);

                vec3 stepScattering = (lightEnergy * phase * inScatterBoost) + multiScattering + ambientlight;
                stepScattering = stepScattering * (1.0 - stepTransmittance) / max(sigmaE, 1e-6);

                scattering += stepScattering * transmittance;
            }

            samplePos += stepVec;
            transmittance *= stepTransmittance;
        }

        mat2x3 atmosphere = AtmosphereScattering(pos, dir, lenToCloud, lightDir, lightLuminance, 4);
        vec3 insca = atmosphere[0];
        vec3 outsca = atmosphere[1];
        scattering = scattering * outsca + insca * (1.0 - transmittance);

        return vec4(scattering, transmittance);
    }

    #define Cloud2DHeight 8.0

    const float Cloud2DReallyHeight = EarthRadius + Cloud2DHeight;

    vec4 RenderCloud2D(vec3 pos, vec3 dir, vec3 lightDir, vec3 lightLuminance) {
        float lenToEarth = rayIntersectSphere(pos, dir, pow2(EarthRadius + ReferenceHeight));
        if(pos.y < Cloud2DReallyHeight && lenToEarth > 0.0) return vec4(0.0, 0.0, 0.0, 1.0);

        float lenToCloud = rayIntersectSphere(pos, dir, Cloud2DReallyHeight * Cloud2DReallyHeight);
        if(lenToCloud < 0.0) return vec4(0.0, 0.0, 0.0, 1.0);

        vec3 cloudpos = pos + dir * lenToCloud;
        lightLuminance *= remapSaturate(lightDir.y, 0.02, 0.1, 0.0, 1.0);
        vec3 lightColor = lightLuminance * TransToAtmos(cloudpos, lightDir);

        float frequency = 700.0;
        float weight = 1.0;
        float time = 0.01 * frameTimeCounter;

        vec2 np = cloudpos.xz / Cloud2DReallyHeight;
        np -= (1.0 - worley2DCell(np * frequency * 1.0 + time)) * 0.001;
        float n = 0.0;
        float c = 0.0;
        for(int i = 0; i < 6; i++) {
            n += perlin2D(np * frequency + time) * weight;
            c += weight;

            frequency *= 2.5;
            weight *= 0.5;
            time *= 1.1;
        }
        n /= c;

        float cloudRange = 0.35;
        float cloudDensity = remapSaturate(n - 1.0 + cloudRange, 0.0, cloudRange, 0.0, 1.0) * 0.035 * mix(CLOUD_CLEAR_DENSITY, CLOUD_RAIN_DENSITY, rainStrength);
        if(cloudDensity < 1e-5) return vec4(0.0, 0.0, 0.0, 1.0);

        float cosTheta = dot(lightDir, dir);
        float forwardPhase = getPhase(cosTheta,  0.5);
        float rearPhase = getPhase(cosTheta, -0.5);
        float phase = mix(forwardPhase, rearPhase, 0.2);

        float sigmaS = cloudDensity * CloudScatteringCoefficient;
        float sigmaA = cloudDensity * CloudAbsorptionCoefficient;
        float sigmaE = sigmaS + sigmaA;
        float transmittance = exp(-sigmaE);
        vec3 scattering = lightColor * sigmaS * phase;

        mat2x3 atmosphere = AtmosphereScattering(pos, dir, lenToCloud, lightDir, lightLuminance, 4);
        vec3 insca = atmosphere[0];
        vec3 outsca = atmosphere[1];
        scattering = scattering * outsca + insca * (1.0 - transmittance);

        return vec4(scattering, transmittance);
    }

    // from Revelation by HaringPro
    // Apache 2.0
    float CloudHighDensity(vec2 rayPos) {
        const float windAngle = radians(CLOUD_HIGH_WIND_ANGLE);
        const vec2 windVelocity = vec2(cos(windAngle), sin(windAngle)) * CLOUD_HIGH_WIND_SPEED;
        vec2 windOffset = windVelocity * worldTimeCounter;

        rayPos -= windOffset;
        rayPos += cameraPosition.xz;

        vec2 curlN = texture(curlNoise2D, rayPos * 5e-5).xy * 0.25;
        vec2 position = rayPos * 2e-4 + curlN;

        float density = 0.0;

        #ifdef CLOUD_CIRRUS
            float ciCoverage = CLOUD_CI_COVERAGE - 0.5 + texture(noisetex, position * 0.01).z;
            ciCoverage = saturate(ciCoverage - texture(cloudMapTex, position * 0.01).y);

            if (ciCoverage > 0.25) {
                vec2 p = position + ciCoverage * 0.5 - windOffset * 1e-4;
                float cirrus = texture(cirroLutTex, p * 0.25).y;
                cirrus *= smoothstep(0.25, 1.0, ciCoverage);
                density += cirrus * cirrus;
            }
        #endif

        #ifdef CLOUD_CIRROCUMULUS
            float ccCoverage = CLOUD_CC_COVERAGE - saturate(texture(noisetex, position * 0.01).z * 1.5);
            ccCoverage = saturate(texture(cloudMapTex, position * 0.01).x * 0.75 + ccCoverage);

            if (ccCoverage > 0.25) {
                vec2 p = position + ccCoverage * 0.5 - windOffset * 1e-4;
                float cirrocumulus = pow2(texture(cirroLutTex, p * 0.25).x);
                cirrocumulus *= smoothstep(0.25, 1.0, ccCoverage);
                density += cirrocumulus;
            }
        #endif

        return density;
    }

    vec4 RenderCloudHigh(vec3 pos, vec3 dir, vec3 lightDir, vec3 lightLuminance) {
        float lenToEarth = rayIntersectSphere(pos, dir, pow2(EarthRadius + ReferenceHeight));
        if (pos.y < cloudHighRealHeight && lenToEarth > 0.0) return vec4(0.0, 0.0, 0.0, 1.0);

        float lenToCloud = rayIntersectSphere(pos, dir, cloudHighRealHeight * cloudHighRealHeight);
        if (lenToCloud < 0.0) return vec4(0.0, 0.0, 0.0, 1.0);

        vec3 cloudPos = pos + dir * lenToCloud;

        float density = CloudHighDensity(cloudPos.xz * 1000.0);
        if (density < cloudHighEpsilon) return vec4(0.0, 0.0, 0.0, 1.0);

        lightLuminance *= remapSaturate(lightDir.y, 0.02, 0.1, 0.0, 1.0);
        vec3 lightColor = lightLuminance * TransToAtmos(cloudPos, lightDir);

        float cosTheta = dot(lightDir, dir);
        float phase = mix(getPhase(cosTheta, 0.5), getPhase(cosTheta, -0.5), 0.3);

        float sigmaS = density * CloudScatteringCoefficient * 0.08;
        float sigmaA = density * CloudAbsorptionCoefficient * 0.02;
        float sigmaE = sigmaS + sigmaA;
    
        float transmittance = exp(-sigmaE);
        vec3 scattering = lightColor * sigmaS * phase * (1.0 - transmittance) / max(sigmaE, 1e-6);

        mat2x3 atmo = AtmosphereScattering(pos, dir, lenToCloud, lightDir, lightLuminance, 4);
        scattering = scattering * atmo[1] + atmo[0] * (1.0 - transmittance);

        return vec4(scattering, transmittance);
    }

#endif