vec3 R2(const float n)
{
    return fract(float(n) * vec3(0.245122333753, 0.430159709002, 0.145122333753)) * 2 - 1;
}
vec4 GetNoise(sampler3D noiseSampler, vec3 position)
{
    return texture(noiseSampler, fract(position));
}
vec3 calculateAbsorbance(vec3 scatterCoeff, float attenuationDepth)
{
    return clamp(exp(-scatterCoeff * attenuationDepth), 0.0, 1.0);
}

float bellFunction(float x_y)
{
    const float peak = 70.0;         // Peak of the bell curve
    const float start = 50.0;        // Approximate start of the bell curve
    const float end = 128.0;         // Approximate end of the bell curve
    const float width = end - start; // Total range
    const float sigma = width / 5.0;
    return exp(-pow((x_y - peak) / sigma, 2.0));
}

float cloud_density2(in vec3 x, float fogAmount)
{
#ifdef OVERWORLD
    // Calculate uniform coverage
    float unifCov = bellFunction(x.y);
    unifCov = clamp(unifCov * fogAmount, 0, 1);
    if (unifCov <= 0.00)
        return 0.0;
#else

    float unifCov = 1;

#endif

    // Scale coordinates
    x *= vec3(12.0, 24.0, 12.0);

    // Initialize noise parameters
    float result = 0.0;
    float w = 1.0;
    float w_sum = 0.0;
    vec3 x_oct = x / 8000.0;

    result += dot(GetNoise(colortex15, x_oct * 1.0 + 1.8 * (frameTimeCounter * 0.001)).xy, vec2(0.5, 2.0)) * 1.05;
    result *= dot(GetNoise(colortex15, x_oct * 4.0 + 1.8 * (frameTimeCounter * 0.002)).xy, vec2(0.5, 2.0)) * 1.01;
    result *= unifCov;
    // Apply bias and adjust result
    float bias_noise = max(0.0, 0.5 - 0.1 * exp(-(x.y - 400.0) / 8192.0));
    float bias = 0.7 * (1 - pow(bias_noise, 1.5));
    result = max(0.0, (result - bias) / (1.0 - bias));

    // Final adjustment and return
    return clamp(result * unifCov, 0.0, 1.0);
}

void AtmosphereLoop(
    float rayLength,
    vec3 start,
    vec3 dV,
    float dY,
    float estEyeDepth,
    vec3 refractedSunVec,
    vec3 ambient,
    vec3 sunColor,
    inout float sh0,
    inout vec3 absorbance_w,
    inout vec3 vL_w,
    inout float density,
    vec3 dVWorld,
    float fogAmount)
{
    float dither = noise_standard(gl_FragCoord.xy);

    float SdotV = dot(sunPosWorld, normalize(dVWorld));
    vec3 rL_Water = water_coefficientRayleigh;
    vec3 mieScatterFactor_Water = dirtScatterCoef;
    vec3 ScatterCoeff_water = (totEpsilon + mieScatterFactor_Water);
    float samples = VL_SAMPLES;
    float expFactorInvSamples = log(11.0) / samples;
    float expFactorBase = pow(11.0, dither / samples);
    float sh = 1.0;
    float mie_water = phaseg(SdotV, 0.4);
    float Rayleigh = phaseRayleigh(SdotV);
    sh0 = 0.0;

    vec3 temp = gbufferModelViewInverse[3].xyz + cameraPosition;
    float sunCaustics = 0.25; // Compensates for the caustics we removed

    for (int i = 0; i < samples; i++)
    {
        float exponent = float(i) / samples;
        float d = (expFactorBase * pow(11.0, exponent) - 1.0 / 11.0) / 10.0;
        float dd = expFactorBase * pow(11.0, exponent) * expFactorInvSamples / 10.0;
        vec3 progress = start + dV * d;

        // Apply distortion to the XY components
        vec3 pos = vec3(progress.xy * calcDistort(progress.xy), progress.z);

#ifdef VL_NOISE_DENSITY
        vec3 dpos = temp + d * dVWorld;
        float densitySample = cloud_density2(dpos, fogAmount);
#endif

#ifdef OVERWORLD
        // Check if position is within bounds
        if (abs(pos.x) < 1.0 - 0.5 / 2048.0 && abs(pos.y) < 1.0 - 0.5 / 2048.0)
        {
            pos = pos * vec3(0.5, 0.5, 0.5 / 6.0) + 0.5;
            sh = clamp(shadow2D(shadow, pos).x, rainStrength, 1.0);
        }
#endif

        float attenuationDepth = dd * rayLength;

        // Calculate light contributions
        vec3 sun = sunColor * sh * sunCaustics;
        vec3 ambientMul_water = exp(-max(estEyeDepth - dY * d, 0.0) * rL_Water);
        vec3 sunMul_water = exp(-max(estEyeDepth - dY * d, 0.0) / abs(refractedSunVec.y) * totEpsilon);
        vec3 mie_water_mult = mie_water * sunMul_water;
        vec3 rayleighContribution_Water = Rayleigh * rL_Water;
        vec3 mieScatteringContribution_Water = mieScatterFactor_Water * mie_water_mult;
        vec3 totalScatter_Water = rayleighContribution_Water + mieScatteringContribution_Water;
        vec3 sunContribution_Water = sun * totalScatter_Water;
        vec3 ambientContribution_Water = ambient * ambientMul_water * (rL_Water + mieScatterFactor_Water);
        vec3 light_water = sunContribution_Water + ambientContribution_Water;

        // Accumulate lighting and absorbance
        vL_w += (light_water - light_water * exp(-ScatterCoeff_water * attenuationDepth)) /
                (ScatterCoeff_water + 1e-8) * absorbance_w;
        absorbance_w *= calculateAbsorbance(ScatterCoeff_water, attenuationDepth);

#ifdef VL_NOISE_DENSITY
        density += densitySample;
#endif
        sh0 += sh;
    }

    // Average shadow and density over samples
    sh0 /= samples;

#ifdef VL_NOISE_DENSITY
    density /= samples;
#else
    density = fogAmount;
#endif
}
