const float multiScatterPhase = 0.5;

float getZenithDensity(float density, float x) {
    return density / pow(max(x, 0.35e-2), 0.75);
}

vec3 getSkyAbsorption(vec3 x, float y){
    return exp2(-x * y) * 2.0;;
}

float getRayleigMultiplier(vec3 worldPos, vec3 lightPos){
    return 1.0 + (1.0 - clamp(distance(worldPos, lightPos), 0.0, 1.0)) * PI;
}

vec3 jodieReinhardTonemap(vec3 c){
    float l = dot(c, vec3(0.2126, 0.7152, 0.0722));
    vec3 tc = c / (c + 1.0);

    return mix(c / (l + 1.0), tc, tc);
}

vec3 getAtmosphericScattering(vec3 worldPos, vec3 viewPos, vec3 lightPos){
    //Variables
    vec3 nViewPos = normalize(viewPos);

    float VoSRaw = dot(nViewPos, sunVec);
    float VoURaw = dot(nViewPos, upVec);
    float VoUClamped = clamp(VoURaw, 0.0, 1.0);
    float VoSClamped = clamp(VoSRaw, 0.0, 1.0);
          VoSClamped = pow(VoSClamped, 1.25);
    float sunVisibilitySqrt = sqrt(sunVisibility);

    //Prepare scattering properties
    float skyDensity = exp(-0.75 * pow(VoUClamped, 2.0 - timeBrightness));
    float scatteringWidth = pow(1.0 - VoUClamped, 2.0 - VoSClamped);

    float sunScatteringMIxer = timeBrightness * (1.0 - wetness * 0.5) * min(exp(VoSClamped), 1.0);
          sunScatteringMIxer = clamp(sunScatteringMIxer * scatteringWidth * 0.33, 0.0, 1.0 - timeBrightness);

    float rayleighScatteringMixer = sunVisibilitySqrt * mix(0.15 + VoSClamped * min(pow(VoUClamped, 0.75) * 2.5, 1.0), 0.5 + VoSClamped * min(pow(VoUClamped, 0.75), 1.0) * 0.25, sunVisibility);
          rayleighScatteringMixer *= (1.0 - wetness * 0.75) * (1.0 - timeBrightness * timeBrightness * 0.5);

    //Realistic sky scattering
    float density = 1.25 - timeBrightness * 0.5;
    float zenith = getZenithDensity(density, worldPos.y);
    float sunPointDistMult = clamp(length(max(lightPos.y + multiScatterPhase, 0.0)), 0.0, 1.0);
    
    const vec3 newSkyCol = vec3(0.39, 0.57, 1.0) * 2.5;
    vec3 absorption = getSkyAbsorption(newSkyCol, zenith);
    vec3 sunAbsorption = getSkyAbsorption(newSkyCol, getZenithDensity(density, lightPos.y + multiScatterPhase));
    vec3 sky = newSkyCol * zenith;
    
    vec3 totalSky = mix(sky * absorption, sky / sqrt(sky * sky + 2.0), sunPointDistMult);
         totalSky *= sunAbsorption * 0.5 + 0.5 * length(sunAbsorption);

    //Day, Night and Rain Sky
    vec3 daySky = mix(normalize(skyColor + 0.0001) * 1.5, skyColSqrt, 0.5 - timeBrightness * 0.25);
         daySky = mix(daySky, lightColSqrt, pow2(1.0 - VoUClamped) * 0.2 * (1.0 - wetness));
         daySky = mix(daySky * skyDensity, totalSky, rayleighScatteringMixer);
         daySky = jodieReinhardTonemap(daySky * PI);
         daySky = pow(daySky, vec3(2.2));
         daySky = mix(daySky, lightCol, sunScatteringMIxer);
    vec3 atmosphere = mix(lightNight * skyDensity * 0.75, daySky, sunVisibilitySqrt);
         atmosphere = mix(atmosphere, weatherCol * clamp(sunVisibility, 0.35, 1.0), wetness * 0.5);

    //Underground Sky
    atmosphere = mix(caveMinLightCol, atmosphere, caveFactor);

    return atmosphere;
}