
#define cloudDensity CLOUD_DENSITY

#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
const float maxDist2 = SHADOW_MAX_DISTANCE * SHADOW_MAX_DISTANCE;

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

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

float cloudtype = mix(CLOUD_TYPE, 1.0, rainStrength);
vec2 windGenerator(float t)
{
    vec2 p = vec2(sin(2.2 * t) - cos(1.4 * t),
                  cos(1.3 * t) - sin(1.9 * t));

    return clamp((p * vec2(0.001, 0.0015) + 1.0) * 0.5, 0.0, 1.0);
}

float calculateCloudDensity(vec3 pointInSpace, vec3 position, float mipLevel, float currentTime, vec2 windDirection)
{
    float cloud_coverage = mix(CLOUD_COVER, 1.0, rainStrength) * 0.5;

    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 = GetNoise(colortex15, currentPosition.xyz * 0.00011).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(colortex15, 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);
}

float getCloudDensity(in vec3 pos, in int LoD, vec2 wind)
{

    return calculateCloudDensity(pos / 3.0, pos, LoD, cloudTime, wind);
}

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, vec2 wind)
{
    float shadow = 0.0;

    // Use squared distance to avoid the sqrt in distance()
    vec2 diff = stepPosXZ - cameraPositionXZ;
    if (sunContributionG > 1e-5 && dot(diff, diff) < maxDist2)
    {
        float factor = cloudDensity * lss;
        vec3 samplePos = stepPos;
        // Increment samplePos iteratively to avoid multiplication every loop iteration
        for (int j = 1; j < SHADOW_STEP_COUNT; j++)
        {
            samplePos += dV_Sun;
            if (samplePos.y < cloud_height_max)
            {
                shadow += getCloudDensity(samplePos, 1, wind) * factor;
            }
        }
    }
    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, bool usefade)
{
    // ---------------------------------
    // 1) Camera & altitude setup
    // ---------------------------------
    vec3 cameraPos = vec3(cameraPosition.x, 70.0, cameraPosition.z);
    float eyeAlt = 70.0;
    vec2 wind = windGenerator(cloudTime * 0.01);

    if (usecamera != 0.0)
    {
        cameraPos = cameraPosition;
        eyeAlt = eyeAltitude;
    }

    // ---------------------------------
    // 2) Direction, accumulators
    // ---------------------------------
    // Hoist the normalization of lightPos to avoid re-normalizing inside loops
    vec3 lightDir = normalize(lightPos);
    vec3 worldDir = normalize(worldPos);

    // We’ll accumulate scattering/transmittance directly here (no local copies).
    vec3 totalScattering = vec3(0.0);
    float totalTransmittance = 1.0;
    float alpha = 0.0;

    // ---------------------------------
    // 3) Cloud “within” checks
    // ---------------------------------
    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));

    // Used to modulate max_distance below
    float within_b = mix(
        1.0,
        0.0,
        pow(
            clamp(distance(eyeAlt, cloud_height_mid), 0.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));

    // ---------------------------------
    // 4) Primary intersection check
    // ---------------------------------
    // Only proceed if we’re within or looking into the volume
    if (visibleVol || (within > 0.0))
    {
        // Precompute the dot for the phase function
        float vDotL = dot(worldDir, lightDir);

        // Intersect with cloud layer top & bottom
        vec3 bottom = worldDir * ((cloud_height - eyeAlt) / worldDir.y);
        vec3 top = worldDir * ((cloud_height_max - eyeAlt) / worldDir.y);

        // If looking “out” of the volume from this side, zero them out
        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;

        // Combined camera position for some checks
        vec3 temp = gbufferModelViewInverse[3].xyz + cameraPos;
        vec3 dpos = temp + worldPos;

        // This test is used to clamp 'end' if terrain is present, etc.
        float test = exp2(-max(dpos.y - cloud_height, 0.0));

        // Step setup
        vec3 rStep = (end - start) / float(PRIMARY_STEP_COUNT);
        // Start position, offset by noise
        vec3 rPos = rStep * dither + start + cameraPos;

        // If terrain is present, possibly clamp the end
        if (terrain && (test < 1.0 || isBelowVol))
        {
            end = worldPos;
        }

        // Ambient & Mie scattering terms
        float ambientMult = exp(-(0.74 + 0.8 * rainStrength) * cloudDensity * 75.0);
        float mieDay = max(
            phaseg(vDotL, 0.6),
            max(phaseg(vDotL, 0.4 - 1.4 * lightPos.y),
                phaseg(vDotL, -0.2)));
        vec3 skyCol0 = atmosphereUp * ambientMult;
        vec3 sunContribution = mieDay * sunLight;

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

        // Shadow stepping
        float lss = cloud_depth / float(SHADOW_STEP_COUNT);
        vec3 dV_Sun = lightDir * lss;

        // “mult” factor for transmittance
        float mult_a = length(worldPos) / PI;
        float mult_b = length(worldPos * cloud_depth / worldPos.y / float(PRIMARY_STEP_COUNT));
        float mult = mix(mult_b, mult_a, within);

        // Distance-based step fade
        float max_distance = mix(PRIMARY_MAX_DISTANCE, SECONDARY_MAX_DISTANCE, within_b);

        // ---------------------------------
        // 5) Ray-marching through the volume
        // ---------------------------------
        // Initialize our accumulators to zero for safety
        totalScattering = vec3(0.0);
        totalTransmittance = 1.0;

        for (int i = 0; i < PRIMARY_STEP_COUNT; ++i)
        {
            // Position for this sample
            vec3 stepPos = calculateStepPosition(
                i,
                dither,
                PI, // expFactor
                PRIMARY_STEP_COUNT,
                max_distance,
                rStep,
                rPos,
                cameraPos,
                worldDir,
                within);

            // Distance from camera in XZ
            float dist = distance(stepPos.xz, cameraPos.xz);
            float atmosFade = exp(-dist * 0.00025);

            // Optional height-based gradient check
            float heightFraction = GetHeightFractionForPoint(stepPos.y);
            float gradient = densityHeightGradient(heightFraction, cloudtype);

            // Skip tiny-contribution samples
            if (gradient < 1.0 && stepPos.y < cloud_height && stepPos.y > cloud_height_max && dist > PRIMARY_MAX_DISTANCE && atmosFade < 0.005)
            {
                continue;
            }

            // Sample cloud density (e.g. 3D noise or texture)
            float cloud = getCloudDensity(stepPos, 4, wind);
            if (cloud <= 0.001)
                continue;

            float muE = cloud * cloudDensity;

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

            // Lighting contribution
            vec3 lightAccum = calculateLighting(
                stepPos.y,
                muE,
                shadow,
                mult,
                sunContribution,
                sunContributionMulti,
                skyCol0,
                totalTransmittance);

            totalScattering += lightAccum;
            totalTransmittance *= max(exp(-mult * muE), 0.0);

            alpha += atmosFade;

            // Early out if it becomes opaque
            if (totalTransmittance < 1e-5)
                break;
        }
    }

    // ---------------------------------
    // 6) Tonemap and final output
    // ---------------------------------
    totalScattering = tonemap(totalScattering);

    // Combine with alpha fade
    vec4 result = vec4(clamp(totalScattering, 0.0, 1.0), totalTransmittance);
    // Mix with black if alpha is small
    result = mix(vec4(0.0, 0.0, 0.0, 1.0), result, clamp(alpha, 0.0, 1.0));
    return result;
}

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