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

float getSunPoint(vec3 p, vec3 lightdirection){
	return smoothstep(0.01, 0.0096, distance(p, lightdirection));
}

float getRayleigMultiplier(vec3 p, vec3 lightdirection){
	return 1.0 + pow(1.0 - clamp(distance(p, lightdirection), 0.0, 1.0), 2.0) * pi * 0.5;
}

float getMie(vec3 p, vec3 lightdirection){
	float disk = clamp(1.0 - pow(distance(p, lightdirection), 0.1), 0.0, 1.0);
	disk += clamp(1.0 - pow(distance(p, lightdirection), 6.9), 0.0, 1.0)*0.1;
	return disk*disk*(3.0 - 2.0 * disk) * 2.0 * pi;
}

vec3 getAtmosphericScattering(vec3 direction, vec3 lightdirection){
		
	float zenith = zenithDensity(direction.y);
	float sunPointDistMult =  clamp(length(max(lightdirection.y + multiScatterPhase, 0.0)), 0.0, 1.0);
	
	float rayleighMult = getRayleigMultiplier(direction, lightdirection);
	
	vec3 absorption = getSkyAbsorption(skyColor, zenith);
    vec3 sunAbsorption = getSkyAbsorption(skyColor, zenithDensity(lightdirection.y + multiScatterPhase));
	vec3 sky = skyColor * zenith * rayleighMult;
	vec3 sun = getSunPoint(direction, lightdirection) * absorption;
	vec3 mie = getMie(direction, lightdirection) * sunAbsorption*vec3(0.5,0.6,0.7);
	
	vec3 totalSky = mix(sky * absorption, sky / sqrt(sky * sky + 2.0), sunPointDistMult);
         totalSky += sun + mie;
	     totalSky *= sunAbsorption * 0.5 + 0.5 * length(sunAbsorption);
	
	return totalSky;
}
