float cloud_height = 1500.0;
float maxHeight = 3000.0;

int maxIT_clouds = 4;
int maxIT = 4;
#define cloudDensity 0.01 // Cloud Density, 0.02-0.04 is around irl values	[0.0010 0.0011 0.0013 0.0015 0.0017 0.0020 0.0023 0.0026 0.0030 0.0034 0.0039 0.0045 0.0051 0.0058 0.0067 0.0077 0.0088 0.0101 0.0115 0.0132 0.0151 0.0173 0.0199 0.0228 0.0261 0.0299 0.0342 0.0392 0.0449 0.0514 0.0589 0.0675 0.0773 0.0885 0.1014 0.1162 0.1331 0.1524 0.1746 0.2000]

float hash(vec3 p)
{
    p = fract(p * 0.3183099 + 0.1) * 17.0;
    return fract(dot(p, p.yzx + vec3(1.0)));
}
float customMix(vec3 i, vec3 f)
{
    float n000 = hash(i + vec3(0.0));
    float n100 = hash(i + vec3(1.0, 0.0, 0.0));
    float n010 = hash(i + vec3(0.0, 1.0, 0.0));
    float n110 = hash(i + vec3(1.0, 1.0, 0.0));
    float n001 = hash(i + vec3(0.0, 0.0, 1.0));
    float n101 = hash(i + vec3(1.0, 0.0, 1.0));
    float n011 = hash(i + vec3(0.0, 1.0, 1.0));
    float n111 = hash(i + vec3(1.0, 1.0, 1.0));

    float nx00 = mix(n000, n100, f.x);
    float nx10 = mix(n010, n110, f.x);
    float nx01 = mix(n001, n101, f.x);
    float nx11 = mix(n011, n111, f.x);

    float nxy0 = mix(nx00, nx10, f.y);
    float nxy1 = mix(nx01, nx11, f.y);

    return mix(nxy0, nxy1, f.z);
}

float noisem(vec3 x)
{
    vec3 i = floor(x);
    vec3 f = fract(x);
    f = f * f * (3.0 - 2.0 * f);

    return customMix(i, f);
}
float noise(vec3 x)
{

    vec3 p = floor(x);
    vec3 f = fract(x);

    f = f * f * (3.0 - 2.0 * f);
    vec2 uv = p.xz + f.xz + p.y;
    vec2 coord = uv / 16.0;

    vec2 noiseSample = texture2D(noisetex, coord).yx;
    return mix(noiseSample.x, noiseSample.y, f.y) * 2.0;
}

float cloud_density(in vec3 x)
{
    float mult = max(x.y - 2000.0, 0.0) / 2000.0;
    float mult2 = max(-x.y + 2000.0, 0.0) / 500.0;
    // Scale the input coordinates
    x *= 0.5;
    float bias_noise = noise((x * vec3(1, 1.2, 1)) / 1024.0);

    // Cloud octaves
    float result = 0.0;
    float w = 1.0;
    float w_sum = 0.0;
    float wx = 0.6;
    vec3 x_oct = x / 256.0;

    float fc2 = frameTimeCounter;
    for (int i = 0; i < 4; i++)
    {
        result += 1.0 * w * noisem(x_oct + fc2 * 0.01);
     
        x_oct *= 2.0;
        w_sum += w;
        w *= wx;
    }

    result /= w_sum;

    // Bias noise modification
    bias_noise = max(0.0, bias_noise - (-0.35) * exp(-(x.y - 400.0) / 8192.0));
    float bias = 1.0 - ((0.0) + (pow(bias_noise, 3.14)) * ((1.0) - (0.0)));
    result = max(0.0, result - bias) / (1.0 - bias);

    // Final shaping with gamma
    result = pow(result, 3.14);
    float cloud = result * result * 3.0 - mult * mult * mult * 3.0 - mult2 * mult2 * 0.9;

    return max(cloud, 0.0);
}

vec4 renderClouds(float dither, vec3 view, vec3 fogC)
{
    vec3 color = vec3(0);

    vec3 worldV = normalize(view.rgb);

    float VdotU = worldV.y;
    float maxIT_clouds = int(clamp(maxIT_clouds / sqrt(VdotU), 0.0, maxIT * 1.0));

    vec3 dV_view = worldV;
    worldV = normalize(worldV) * 100000.0 + cameraPosition;
    if (worldV.y < cloud_height)
    {
        return vec4(fogC, 0.0);
    }

    vec3 progress_view = dV_view * dither + cameraPosition;
    float total_extinction = 1.0;

    dV_view = normalize(dV_view);

    dV_view *= max(maxHeight - cloud_height, 0.0) / dV_view.y / maxIT_clouds;
    vec3 startOffset = dV_view * dither;
    progress_view = startOffset + cameraPosition + dV_view * (cloud_height - cameraPosition.y) / dV_view.y;

    float ftimeTemp = frameTimeCounter * 35.0;
    progress_view.x += ftimeTemp;

    float shadowStep = 400.0;
    vec3 dV_Sun = normalize(mat3(gbufferModelViewInverse) * sunVec) * shadowStep;
    float mult = length(dV_view);
    float SdotV = dot(sunVec, mat3(gbufferModelView) * view);

    float mieDay = max(phaseg(SdotV, 0.3), phaseg(SdotV, 0.12) * 0.7);
    vec3 sunContribution = mieDay * sunLight;
    float ambientMult = exp(-(0.24 + 0.8 * rainStrength) * cloudDensity * 750.0);
    vec3 skyCol0 = atmosphereUp * ambientMult;
#ifdef END
    skyCol0 = fogColor.rgb * fogColor.rgb;
    sunContribution = mieDay * skyCol0 * PI;

#endif
    const float threshold = 0.001;
    const float expThreshold = 1e-5;

    for (int i = 0; i < maxIT_clouds; i++)
    {
        float cloud = cloud_density(progress_view);

        if (cloud > threshold)
        {
            float muS = cloud * cloudDensity;
            float muE = muS;
            float muEshD = 0.0;

            vec3 shadowSamplePos = progress_view + dV_Sun;
            if (shadowSamplePos.y < maxHeight)
            {
                float cloudS = cloud_density(shadowSamplePos);
                muEshD += cloudS * cloudDensity * shadowStep;
            }

            float h = 0.5 - 0.5 * clamp((progress_view.y - 1500.0) / 4000.0, 0.0, 1.0);
            float powder = 1.0 - exp(-muE * mult);

            vec3 S = sunContribution * exp(-muEshD) + skyCol0 * mix(1.0, powder, h * ambientMult);
            vec3 Sint = (S - S * exp(-mult * muE)) / muE;
            color += muS * Sint * total_extinction;
            total_extinction *= exp(-muE * mult);

            if (total_extinction < expThreshold)
            {
                break;
            }
        }

        progress_view += dV_view;
    }

    float dist = (cloud_height - cameraPosition.y) / normalize(dV_view).y;
    float fog = clamp(exp(-dist / 6000.0 * (1.0 + rainStrength * 2.0)), 0.0, 1.0);

    return mix(vec4(fogC, 0.0), vec4(color, clamp(total_extinction, 0.0, 1.0)), fog);
}
