#define VL_CLOUDS_SHADOWS
#define CLOUDS_SHADOWS
#define CLOUDS_SHADOWS_STRENGTH 1.0 //[0.1 0.125 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.9 1.0]
#define VL_SAMPLES2 6 //[4 6 8 10 12 14 16 20 24 30 40 50]
#define Ambient_Mult 1.0 //[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.5 2.0 3.0 4.0 5.0 6.0 10.0]
#define SEA_LEVEL 70 //[0 10 20 30 40 50 60 70 80 90 100 110 120 130 150 170 190]	//The volumetric light uses an altitude-based fog density, this is where fog density is the highest, adjust this value according to your world.
#define ATMOSPHERIC_DENSITY 1.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 4.0 5.0 7.5 10.0 12.5 15.0 20.]
#define fog_mieg1 0.40 //[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.75 0.8 0.85 0.9 0.95 1.0]
#define fog_mieg2 0.10 //[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.75 0.8 0.85 0.9 0.95 1.0]
#define fog_coefficientRayleighR 5.8 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
#define fog_coefficientRayleighG 1.35 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
#define fog_coefficientRayleighB 3.31 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]

#define fog_coefficientMieR 3.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
#define fog_coefficientMieG 3.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
#define fog_coefficientMieB 3.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
 
#define Cloudy_Fog_Density 5.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]
#define Uniform_Fog_Density 1.0 //[0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0]

#define uniformfog_fade 10 // [5 10 20 30 40 50 60 70 80  90 100]
#define cloudyfog_fade 10 // [5 10 20 30 40 50 60 70 80  90 100]

float luma(vec3 color) {
	return dot(color,vec3(0.299, 0.587, 0.114));
}

float phaseRayleigh(float cosTheta) {
	const vec2 mul_add = vec2(0.1, 0.28) /acos(-1.0);
	return cosTheta * mul_add.x + mul_add.y; // optimized version from [Elek09], divided by 4 pi for energy conservation
}

float fog_densities_atmospheric = max((1.0-TimeOfDayFog)*100, 24); // this is seperate from the cloudy and uniform fog.

float cloudVol(in vec3 pos){
	// various fade away heights
	float cloudyfade =  exp2(-max(pos.y-SEA_LEVEL,0.0)/ mix(cloudyfog_fade, Biome_Fog_Properties[6], Biome_Fog_Properties[0]));
	float uniformfade = exp2(-max(pos.y-SEA_LEVEL,0.0)/ mix(cloudyfog_fade, Biome_Fog_Properties[7], Biome_Fog_Properties[0]));

	// cloudy fog shapes
	vec3 samplePos = pos*vec3(1.0,1./48.,1.0); // +vec3(0,frameTimeCounter/80,0);
	float main_cloudy_fog    = densityAtPos(samplePos*50);
	float cloudy_fog_erosion = 1.0 - densityAtPos(samplePos*500);
	
	// create the shape of the uniform fog and combine the cloudy fog noises
	float cloudyfog = max(mix(main_cloudy_fog*5, cloudy_fog_erosion*5, main_cloudy_fog-0.7) * cloudyfade, 0.0);
	float uniformfog = uniformfade;

	// density modifiers + specific biome fog densities
	cloudyfog  *= mix((1.0-TimeOfDayFog)*Cloudy_Fog_Density  + Weather_properties.g, Biome_Fog_Properties[1], Biome_Fog_Properties[0]);
	uniformfog *= mix((1.0-TimeOfDayFog)*Uniform_Fog_Density + Weather_properties.b, Biome_Fog_Properties[2], Biome_Fog_Properties[0]);

	// combine them all
	float Fog_Final = cloudyfog+uniformfog ;
	
	// rain fog
	Fog_Final +=  mix(0, uniformfade, rainStrength);

	return clamp(Fog_Final,0.0,255.0);
}
mat2x3 getVolumetricRays(float dither,vec3 fragpos) {

	//project pixel position into projected shadowmap space
	vec3 wpos = mat3(gbufferModelViewInverse) * fragpos + gbufferModelViewInverse[3].xyz;
	vec3 fragposition = mat3(shadowModelView) * wpos + shadowModelView[3].xyz;
	fragposition = diagonal3(shadowProjection) * fragposition + shadowProjection[3].xyz;



	//project view origin into projected shadowmap space
	vec3 start = toShadowSpaceProjected(vec3(0.));


	//rayvector into projected shadow map space
	//we can use a projected vector because its orthographic projection
	//however we still have to send it to curved shadow map space every step
	vec3 dV = (fragposition-start);
	vec3 dVWorld = (wpos-gbufferModelViewInverse[3].xyz);

	float maxLength = min(length(dVWorld),far)/length(dVWorld);
	dV *= maxLength;
	dVWorld *= maxLength;

	//apply dither
	vec3 progress = start.xyz;
	vec3 progressW = gbufferModelViewInverse[3].xyz+cameraPosition;
		vec3 vL = vec3(0.);

		float SdotV = dot(sunVec,normalize(fragpos))*lightCol.a;
		float dL = length(dVWorld);
		//Mie phase + somewhat simulates multiple scattering (Horizon zero down cloud approx)
		float mie = 2*mix(phaseg(SdotV,fog_mieg1),	phaseg(SdotV,fog_mieg2),0.5);


		float rayL = phaseRayleigh(SdotV);

		vec3 ambientCoefs = dVWorld/dot(abs(dVWorld),vec3(1.));

		vec3 ambientLight = ambientUp*clamp(ambientCoefs.y,0.,1.);
		ambientLight += ambientDown*clamp(-ambientCoefs.y,0.,1.);
		ambientLight += ambientRight*clamp(ambientCoefs.x,0.,1.);
		ambientLight += ambientLeft*clamp(-ambientCoefs.x,0.,1.);
		ambientLight += ambientB*clamp(ambientCoefs.z,0.,1.);
		ambientLight += ambientF*clamp(-ambientCoefs.z,0.,1.);

		vec3 skyCol0 = ambientLight*max(eyeBrightnessSmooth.y,5)/vec3(240.)*Ambient_Mult*2.0*8./150./3.;
		
		// Makes fog more white idk how to simulate it correctly
		vec3 sunColor = lightCol.rgb * 8./150./3. ;


		// allow biome fogs to do their thing
		vec3 biomefog = ( luma(lightCol.rgb) * Biome_FogColor ) * 8./150./3.;

		sunColor = mix(sunColor, 	biomefog, 	   Biome_Fog_Properties[0]);
		skyCol0 =  mix(skyCol0.rgb, biomefog*0.025, Biome_Fog_Properties[0]);
		
		float biome_mie = Biome_Fog_Properties[3] * mix(phaseg(SdotV,0.1), phaseg(SdotV,Biome_Fog_Properties[5]), 0.5);
		mie += mix(0, biome_mie,  Biome_Fog_Properties[0]);

		vec3 rC = vec3(fog_coefficientRayleighR*1e-6, fog_coefficientRayleighG*1e-5, fog_coefficientRayleighB*1e-5);
		vec3 mC = vec3(fog_coefficientMieR*1e-6, fog_coefficientMieG*1e-6, fog_coefficientMieB*1e-6);

		float mu = 1.0;
		float muS = 1.0*mu;
		vec3 absorbance = vec3(1.0);
		float expFactor = 11.0;
		vec3 WsunVec = mat3(gbufferModelViewInverse) * sunVec * lightCol.a;

		// vec3 CprogressW = gbufferModelViewInverse[3].xyz+cameraPosition + dVWorld;
		float cloudShadow = 1.0;

		float mie_cloud = cloudShadow*( mix(phaseg(SdotV,0.5),	2*phaseg(SdotV,0.7),0.5)) * rainStrength;

		for (int i=0;i<VL_SAMPLES2;i++) {
			float d = (pow(expFactor, float(i+dither)/float(VL_SAMPLES2))/expFactor - 1.0/expFactor)/(1-1.0/expFactor);
			float dd = pow(expFactor, float(i+dither)/float(VL_SAMPLES2)) * log(expFactor) / float(VL_SAMPLES2)/(expFactor-1.0);
			progress = start.xyz + d*dV;
			progressW = gbufferModelViewInverse[3].xyz+cameraPosition + d*dVWorld;
			//project into biased shadowmap space
			float distortFactor = calcDistort(progress.xy);
			vec3 pos = vec3(progress.xy*distortFactor, progress.z);
			float densityVol = cloudVol(progressW);
			float sh = 1.0;
			if (abs(pos.x) < 1.0-0.5/2048. && abs(pos.y) < 1.0-0.5/2048){
				pos = pos*vec3(0.5,0.5,0.5/6.0)+0.5;
				sh =  shadow2D( shadow, pos).x;
			}
			#ifdef VL_Clouds_Shadows
			#ifdef CLOUDS_SHADOWS
				float max_height = clamp(400.0 - progressW.y, 0.0,1.0); // so it doesnt go beyond the height of the clouds
				vec3 campos = (progressW)-319;
				// get cloud position
				vec3 cloudPos = campos*Cloud_Size + WsunVec/abs(WsunVec.y) * (2250 - campos.y*Cloud_Size);
				// get the cloud density and apply it
				cloudShadow = getCloudDensity(cloudPos, 1);
				cloudShadow = exp(-cloudShadow*sqrt(cloudDensity)*25);
				cloudShadow *= max_height; 
			#endif
			#endif

			//Water droplets(fog)
			float density = densityVol*ATMOSPHERIC_DENSITY*mu*300.;

			//Just air
			vec2 airCoef = exp2(-max(progressW.y-SEA_LEVEL,0.0)/vec2(8.0e3, 1.2e3)*vec2(6.,7.0)) * fog_densities_atmospheric;
			
			//Pbr for air, yolo mix between mie and rayleigh for water droplets
			vec3 rL = rC*airCoef.x;
			vec3 m = (airCoef.y+density)*mC;

			vec3 vL0 = ( sunColor*sh * cloudShadow)*(rayL*rL+m*mie) + skyCol0*(rL+m);
			vL += (vL0 - vL0 * exp(-(rL+m)*dd*dL)) / ((rL+m)+0.00000001)*absorbance;
			absorbance *= clamp(exp(-(rL+m)*dd*dL),0.0,1.0);
		}
	return mat2x3(vL,absorbance);



}
