uniform sampler3D depthtex1;

#define cloudDensity CLOUD_DENSITY
#define sstep(x, low, high) smoothstep(low, high, x)
#define remap01(x, out_lo, out_hi) ((out_lo) + x * ((out_hi) - (out_lo)))

#ifdef END
float cloud_height = 350; // Start of the cloud layer
#else
float cloud_height = heightLimit + 128; // Start of the cloud layer
#endif

float cloud_depth = CLOUD_LAYER_HEIGHT;     // Depth of the cloud layer
float cloud_depth_half = cloud_depth * 0.5; // Half Depth of the cloud layer

float cloud_height_max = cloud_height + cloud_depth;
float cloud_height_mid = cloud_height + cloud_depth * 0.5;
#define PRIMARY_STEP_COUNT CLOUD_MARCHING_STEPS
#define SHADOW_STEP_COUNT 6
#define SHADOW_MAX_DISTANCE 8192

#ifdef DISTANT_HORIZONS
    #define PRIMARY_MAX_DISTANCE 32768
    #define SECONDARY_MAX_DISTANCE 8192

#else
    #define PRIMARY_MAX_DISTANCE 32768
    #define SECONDARY_MAX_DISTANCE 8192
#endif
#define STEP_NONLINEARITY 1.5

///

// Utility function that maps a value from one range to another.
float remap2(const float originalValue, const float originalMin, const float originalMax, const float newMin, const float newMax)
{
    return newMin + (((originalValue - originalMin) / (originalMax - originalMin)) * (newMax - newMin));
}

// fractional value for sample position in the cloud layer
float GetHeightFractionForPoint(float inPosition)
{ // get global fractional position in cloud zone
    float height_fraction = (inPosition - cloud_height) / cloud_depth;
    return clamp(height_fraction, 0.0, 1.0);
}
vec3 tonemap(vec3 color)
{
    float white = PI;
    float luma = dot(color, vec3(0.2126, 0.7152, 0.0722));
    float toneMappedLuma = luma * (1. + luma / (white * white)) / (1. + luma);
    color *= toneMappedLuma / luma;

    return color;
}

vec4 mixGradients(const float cloudType)
{

    const vec4 STRATUS_GRADIENT = vec4(0.02, 0.05, 0.09, 0.11);
    const vec4 STRATOCUMULUS_GRADIENT = vec4(0.02, 0.2, 0.48, 0.625);
    const vec4 CUMULUS_GRADIENT = vec4(0.01, 0.0625, 0.78, 1.0); // these fractions would need to be altered if cumulonimbus are added to the same pass
    float stratus = 1.0 - clamp(cloudType * 2.0, 0.0, 1.0);
    float stratocumulus = 1.0 - abs(cloudType - 0.5) * 2.0;
    float cumulus = clamp(cloudType - 0.5, 0.0, 1.0) * 2.0;
    return STRATUS_GRADIENT * stratus + STRATOCUMULUS_GRADIENT * stratocumulus + CUMULUS_GRADIENT * cumulus;
}
vec4 getCloudGradient(float height_fraction)
{
    const float volumeHeight = 1028.0;
    float height = height_fraction * volumeHeight;

    float cloudType;

    if (height < 257.0)
    {
        // Stratus layer
        cloudType = height / 257.0 * 0.25;
    }
    else if (height < 771.0)
    {
        // Stratocumulus layer
        cloudType = 0.25 + (height - 257.0) / (771.0 - 257.0) * 0.5;
    }
    else
    {
        // Cumulus layer
        cloudType = 0.75 + (height - 771.0) / (1028.0 - 771.0) * 0.25;
    }

    return mixGradients(cloudType);
}
float densityHeightGradient(const float heightFrac, const float cloudType)
{
    vec4 cloudGradient = mixGradients(cloudType);
    return smoothstep(cloudGradient.x, cloudGradient.y, heightFrac) - smoothstep(cloudGradient.z, cloudGradient.w, heightFrac);
}
const vec4 turbulenceTable[17] = vec4[17](
    vec4(0.031000, -0.015000, 0.022500, 1.080000),
    vec4(0.041500, 0.062500, -0.004500, 0.895000),
    vec4(-0.073500, 0.061000, 0.067500, 0.575000),
    vec4(0.058000, 0.060500, -0.081000, 0.755000),
    vec4(0.026500, 0.037500, -0.069000, 0.545000),
    vec4(0.058500, -0.039000, -0.050500, 1.125000),
    vec4(0.074500, -0.035500, -0.052500, 0.940000),
    vec4(0.009000, 0.043500, 0.055500, 0.620000),
    vec4(-0.061000, 0.048000, 0.003000, 0.805000),
    vec4(0.082000, -0.021000, 0.070000, 1.020000),
    vec4(0.018000, 0.048000, -0.052000, 0.700000),
    vec4(0.013000, -0.052500, -0.032500, 0.530000),
    vec4(0.059500, -0.003500, -0.059500, 0.890000),
    vec4(-0.080500, 0.005500, -0.050000, 0.615000),
    vec4(0.022000, -0.066500, 0.085000, 0.825000),
    vec4(0.048500, 0.063000, -0.049000, 0.980000),
    vec4(-0.021500, -0.054000, 0.061500, 0.650000));

vec4 GetNoise(sampler3D noiseSampler, vec3 position)
{
    return texture(noiseSampler, fract(position));
}
vec4 GetNoiseMip(sampler3D noiseSampler, vec3 position, float mip)
{
    return texture(noiseSampler, fract(position));
}
float cloud_coverage = mix(CLOUD_COVER, 1.0, rainStrength) * 0.5;
float cloudtype = mix(CLOUD_TYPE, 1.0, rainStrength);
vec2 windGenerator(float t)
{
    float tx = t;

    // Calculate the raw values
    vec2 p = vec2(sin(2.2 * tx) + -cos(1.4 * tx), cos(1.3 * t) + sin(-1.9 * t));

    // Apply scaling to fit the desired range
    p.y *= 0.0015;
    p.x *= 0.001;

    // Normalize to [0, 1]
    p = (p + 1.0) * 0.5;

    // Clamp to ensure values are within [0, 1]
    p = clamp(p, 0.0, 1.0);

    return p;
}

float calculateCloudDensity(vec3 pointInSpace, vec3 position, float mipLevel, float currentTime, vec2 windDirection)
{
    float windSpeed = sin(windDirection.x + windDirection.y);
    vec3 weather = vec3(cloudtype, 1, 1);
    vec3 currentPosition = pointInSpace;
    float heightFraction = GetHeightFractionForPoint(position.y);
    float gradient = densityHeightGradient(heightFraction, weather.x);
    // Early exit if there are no clouds
    if (gradient <= 0.001)
        return 0.0;

    // Apply wind influence to the position
    currentPosition.xz += currentTime * 20.0 * normalize(windDirection) * windSpeed * 0.6;
    currentPosition.y *= 2.0;

    // Get noise values
    vec2 noise = GetNoiseMip(depthtex1, currentPosition.xyz * 0.00011, mipLevel - 2.0).xy;

    float baseCloudDensity = remap2(noise.x, -(1.0 - noise.y * noise.x), 1.0, 0.0, 1.0);

    // Adjust cloud coverage
    float weatherCoverage = cloud_coverage * weather.b;
    baseCloudDensity = remap2(baseCloudDensity * gradient, 1.0 - (weatherCoverage), 1.0, 0.0, 1.0);
    baseCloudDensity *= weatherCoverage;

    // Adjust position for time
    currentPosition.xz -= currentTime * normalize(windDirection) * 5.0;
    currentPosition.y -= currentTime * 5.0;

    // Initialize variables for fractal brownian motion (FBM)
    float fractalNoise = 0;
    float weight = 1.0;
    float weightSum = 0.0;
    currentPosition *= 0.003;
    currentPosition.y *= 0.5;

    // Calculate FBM
    for (int i = 0; i < mipLevel; i++)
    {
        vec4 turbulence = turbulenceTable[i];
        fractalNoise += dot(turbulence.w * weight * GetNoise(depthtex1, currentPosition + 1.8 * (currentTime * 0.05) * turbulence.xyz).xy, vec2(0.5, 2.0));

        currentPosition *= 2.0;
        weightSum += weight;
        weight *= 0.5;
    }
    fractalNoise /= weightSum;

    // Blend fractal noise based on height
    fractalNoise = mix(fractalNoise, 1.0 - fractalNoise, clamp(heightFraction * 4.0, 0.0, 1.0));
    baseCloudDensity = remap2(baseCloudDensity, fractalNoise * 0.4 * heightFraction, 1.0, 0.0, 1.0);

    // return gradient;
    return pow(clamp(baseCloudDensity, 0.0, 1.0), (1.0 - heightFraction) * 0.8 + 0.5);
}

vec2 hash(in vec2 x) // this hash is not production ready, please
{                    // replace this by something better
    const vec2 k = vec2(0.3183099, 0.3678794);
    x = x * k + k.yx;
    return -1.0 + 2.0 * fract(16.0 * k * fract(x.x * x.y * (x.x + x.y)));
}

float getCloudDensity(in vec3 pos, in int LoD)
{
    float time = cloudTime;

    vec2 wind = windGenerator(time * 0.01);
    return calculateCloudDensity(pos / 3.0, pos, LoD, time, wind);
}
float gbf10(in vec3 eyePlayerPos, in float far)
{
    float distXZ = length((vec3(eyePlayerPos.x, 0.0, eyePlayerPos.z)).xyz);
    float distY = length((vec3(0.0, eyePlayerPos.y, 0.0)).xyz);
    float eyeDist = max(distXZ, distY);
#ifdef DISTANT_HORIZONS
    float borderFogFactor = smoothstep(0.0, dhRenderDistance, eyeDist);
    borderFogFactor *= borderFogFactor * borderFogFactor;
    return borderFogFactor;
#else

    float borderFogFactor = smoothstep(far * 0.5, far, eyeDist);
    borderFogFactor *= borderFogFactor * borderFogFactor;

    return borderFogFactor;

#endif
}
vec3 calculateLighting(float stepPy, float muE, float muEshD, float mult, vec3 sunContribution, vec3 sunContributionMulti, vec3 skyCol0, float transmittance)
{

    float h = 0.35 - 0.35 * GetHeightFractionForPoint(stepPy);
    float powder = 1.0 - exp(-muE * mult);
    float sunShadowMulti = exp(-log(muEshD * 0.15 + 0.5)) * powder;
    float sunShadow = exp(-muEshD);

    float ambientPowder = mix(1.0, powder, h) * (1.0 - h);
    vec3 S = vec3(sunContribution * sunShadow + sunShadowMulti * sunContributionMulti + skyCol0 * ambientPowder);

    vec3 Sint = (S - S * exp(-mult * muE)) / muE;
    vec3 stepScatter = muE * Sint * transmittance;

    return stepScatter;
}
float calculateSunShadow(float sunContributionG, vec2 stepPosXZ, vec2 cameraPositionXZ, vec3 stepPos, vec3 dV_Sun, float lss)
{

    float shadow = 0.0;

    if (sunContributionG > 1e-5 && distance(stepPosXZ, cameraPositionXZ) < SHADOW_MAX_DISTANCE)
    {
        for (int j = 1; j < SHADOW_STEP_COUNT; j++)
        {
            vec3 shadowSamplePos = stepPos + dV_Sun * float(j);
            if (shadowSamplePos.y < cloud_height_max)
            {
                float cloudS = getCloudDensity(shadowSamplePos, 1);
                shadow += cloudS * cloudDensity * lss;
            }
        }
    }

    return shadow;
}
vec3 calculateStepPosition(float i, float noise, float expFactor, float stepCount,
                           float max_distance, vec3 rStep, vec3 rPos,
                           vec3 cameraPosition, vec3 worldDir, float within)
{
    float step_factor = (pow(expFactor, float(i + noise) / float(stepCount)) / expFactor - 1.0 / expFactor) / (1.0 - 1.0 / expFactor);
    float step_distance = mix(step_factor, step_factor * float(max_distance), step_factor);

    vec3 step_p_main = rStep * i + rPos;
    vec3 step_p_alt = cameraPosition + step_distance * worldDir;

    vec3 stepPos = mix(step_p_main, step_p_alt, within);

    return stepPos;
}

vec4 renderClouds(float dither, vec3 worldPos, vec3 fogC, bool terrain, float usecamera)
{
    float eyeAlt;
    vec3 cameraPos;
    if (usecamera == 0)
    {
        cameraPos = vec3(cameraPosition.x, 70, cameraPosition.z);
        eyeAlt = 70;
    }
    else
    {
        cameraPos = cameraPosition;
        eyeAlt = eyeAltitude;
    }

    vec3 worldDir = normalize(worldPos);
    vec3 totalScattering = vec3(0.0);
    float totalTransmittance = 1.0;
    float b = cloud_depth_half;
    float within_a = smoothstep(float(cloud_height) - b, float(cloud_height), eyeAlt) *
                     (1.0 - smoothstep(cloud_height_max, cloud_height_max + b, eyeAlt));
    float within_b = mix(1.0, 0.0, pow(clamp(distance(eyeAlt, cloud_height_mid), 0, cloud_depth_half) / cloud_depth_half, 8.0));
    float within = within_a;

    bool isBelowVol = eyeAlt < cloud_height_mid;
    bool visibleVol = (worldDir.y > 0.0 && isBelowVol) || (worldDir.y < 0.0 && !isBelowVol);
    float alpha = 0;

    if (visibleVol || within > 0.0)
    {

        float vDotL = dot(normalize(worldPos), normalize(lightPos));

        vec3 bottom = worldDir * ((cloud_height - eyeAlt) / worldDir.y);
        vec3 top = worldDir * ((cloud_height_max - eyeAlt) / worldDir.y);

        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;

        vec3 temp = gbufferModelViewInverse[3].xyz + cameraPos;
        vec3 dpos = temp + worldPos;

        float test = exp2(-max(dpos.y - cloud_height, 0.0));

        if (terrain)
        {
            if (test < 1.0 || isBelowVol)
                end = worldPos;
        }

        float noise = dither * 1.0;
        vec3 rStep = (end - start) / float(PRIMARY_STEP_COUNT);
        // rStep = mix(rStep, worldDir * 256.0 * noise, within);

        vec3 rPos = rStep * noise + start + cameraPos;
        // rPos = mix(rPos, cameraPos, within);

        if (terrain && within > 0.0)
        {
            if (test < 1.0 || isBelowVol)
                rPos = worldPos;
        }

        vec3 scattering = vec3(0.0);
        float transmittance = 1.0;
        float ambientMult = exp(-(-0.3 + 0.24 + 0.8 * rainStrength) * cloudDensity * 75.0);
        float mieDay = max(max(phaseg(vDotL, 0.6), phaseg(vDotL, 0.4 - 1.4 * lightPos.y)), phaseg(vDotL, -0.2));

        float lss = (cloud_depth / SHADOW_STEP_COUNT);
        vec3 skyCol0 = atmosphereUp * ambientMult;
        vec3 sunContribution = mieDay * sunLight;

        float mieDayMulti = phaseg(vDotL, 0.2);
        vec3 sunContributionMulti = mieDayMulti * sunLight;

        vec3 dV_Sun = lightPos * lss;
        float mult_a = length(worldPos) / PI;
        float mult_b = length(worldPos * cloud_depth / worldPos.y / PRIMARY_STEP_COUNT);
        float mult = mix(mult_b, mult_a, within);

        const float expFactor = PI;
        float max_distance = mix(PRIMARY_MAX_DISTANCE, SECONDARY_MAX_DISTANCE, within_b);

        for (int i = 0; i < PRIMARY_STEP_COUNT; ++i)
        {

            vec3 stepPos = calculateStepPosition(i, noise, expFactor, PRIMARY_STEP_COUNT, max_distance, rStep, rPos, cameraPos, worldDir, within);

            float dist = distance(stepPos.xz, cameraPos.xz);
            float atmosFade = exp(-dist * 0.00025);

            float heightFraction = GetHeightFractionForPoint(stepPos.y);
            float gradient = densityHeightGradient(heightFraction, cloudtype);
            if (gradient < 1.0 && stepPos.y < cloud_height && stepPos.y > cloud_height_max && distance(stepPos.xz, cameraPos.xz) > PRIMARY_MAX_DISTANCE && atmosFade < 0.005)
                continue;

            float cloud = getCloudDensity(stepPos, int(4));
            if (cloud <= 0.001)
                continue;

            float muE = cloud * cloudDensity;

            float shadow = calculateSunShadow(sunContribution.g, stepPos.xz, cameraPos.xz, stepPos, dV_Sun, lss);

            scattering += calculateLighting(stepPos.y, muE, shadow, mult, sunContribution, sunContributionMulti, skyCol0, transmittance);

            transmittance *= max(exp(-mult * muE), 0.0);
            alpha += atmosFade;

            if (transmittance < 1e-5)
                break;
        }
        totalScattering = scattering;
        totalTransmittance = transmittance;
    }
    totalScattering = tonemap(totalScattering);

    vec4 result = vec4(clamp(totalScattering, 0.0, 1.0), totalTransmittance);
    result = mix(vec4(0, 0, 0, 1), result, clamp(alpha, 0, 1));
    return result;
}

//////////////////
