uniform sampler2D 	depthTex,
					colorTex;
uniform sampler3D 	shapeCloudTexture;
uniform sampler3D 	detailCloudTexture;
uniform sampler2D	weatherTexture;

uniform vec3		sundir;
varying vec2 		texcoord;
varying vec2		VPOS;

uniform vec2 		screenSize;

uniform mat4		eyeToWorld;
uniform vec3		campos;

uniform float		wFar;
uniform float 		wNear;

uniform vec3 		earth_center;
uniform float 		earth_radius;
uniform float 		atmo_bottom_radius;
uniform float 		atmo_top_radius;
uniform float 		length_unit;
uniform float		inScatter;

uniform vec3 		suncolor;
uniform float		cDensity;
uniform vec2		cloudsMat;// x=density, y=anisotropy

vec4 Wpos=vec4(0.0,0.0,0.0,1.0);

/*#define GetSolarRadiance GetSolarLuminance
#define GetSkyRadiance GetSkyLuminance
#define GetSkyRadianceToPoint GetSkyLuminanceToPoint
#define GetSunAndSkyIrradiance GetSunAndSkyIlluminance
#define GetSunLuminance GetSunIrradiance
#define GetSkyLuminanceToAtmosphere	GetSkyRadianceToAtmosphere
#define GetIndirectIrradiance GetIndirectIrradiance
*/

vec3 GetSolarRadiance();
vec3 GetSkyRadiance(vec3 camera, vec3 view_ray, float shadow_length,
    vec3 sun_direction, out vec3 transmittance);
vec3 GetSkyRadianceToPoint(vec3 camera, vec3 point, float shadow_length,
    vec3 sun_direction, out vec3 transmittance);
vec3 GetSunAndSkyIrradiance(
    vec3 p, vec3 normal, vec3 sun_direction, out vec3 sky_irradiance);
vec3 GetTrasmittanceTexture(vec2 uv);
vec3 GetScatteringTexture(vec3 uv);
vec3 GetSunIrradiance(vec3 p, vec3 sun_direction);
vec3 GetSkyRadianceToAtmosphere(vec3 camera, vec3 point, float shadow_length, vec3 sun_direction, out vec3 transmittance);
vec3 GetIndirectIrradiance(vec3 p, vec3 sun_direction);

////////////////////////////
const float PI=3.1415926535;
vec3 g_weather;

#define MAX_STEPS 128//128
#define MIN_STEPS 32//32

#define SHADOWS_STEPS	8

#define TOP_BOUND 		atmo_top_radius
#define BOTTOM_BOUND	atmo_bottom_radius

float blend1D(float f1, float f2, float a1, float a2, float alpha)
{
	float depth=0.5;
	float ma=max(a1+(1.0-alpha),a2+alpha)-depth;
	
	float b1=max(a1 + (1.0-alpha) - ma,0.0);
	float b2=max(a2 + alpha - ma,0.0);
	
	return (f1*b1 + f2*b2)/(b1+b2);
	//return (a1+(1.0-alpha)) > (a2+alpha) ? f1 : f2;
}

float bayer2(vec2 a){
    a = floor(a);
    return fract( dot(a, vec2(.5, a.y * .75)) );
}

#define bayer4(a)   (bayer2( .5*(a))*.25+bayer2(a))
#define bayer8(a)   (bayer4( .5*(a))*.25+bayer2(a))
#define bayer16(a)  (bayer8( .5*(a))*.25+bayer2(a))
#define bayer32(a)  (bayer16(.5*(a))*.25+bayer2(a))
#define bayer64(a)  (bayer32(.5*(a))*.25+bayer2(a))
#define bayer128(a) (bayer64(.5*(a))*.25+bayer2(a))


float powder(float od)
{
	return 1.0 - exp2(-od * 2.0);
}

float beerPowder(float od, float cos)
{
	return exp(-od)*(1.0-(cos+1)/2.0*exp(-od*2.0));//(-exp(-(od-2.2)*2)*0.1) + exp(-(od-2.2))*0.9;
}

// exponential integral
float Ei( float z )
{
	return 0.5772156649015328606065 + log( 1e-4 + abs(z) ) + z * (1.0 + z * (0.25 + z * ( (1.0/18.0) + z * ( (1.0/96.0) + z * (1.0/600.0) ) ) ) ); // For x!=0
}

float hgPhase(float x, float g)
{
    float g2 = g*g;
	return 0.25 * ((1.0 - g2) * pow(1.0 + g2 - 2.0*g*x, -1.5));
}

float phase2Lobes(float x)
{
    const float m = 0.6;
    const float gm = cloudsMat.y;//0.8;
    
	float lobe1 = hgPhase(x, 0.8 * gm);
    float lobe2 = hgPhase(x, -0.5 * gm);
    
    return mix(lobe2, lobe1, m);
}

float getHeightSignal(vec3 pos, vec2 weather)
{
	float height=atmo_top_radius-atmo_bottom_radius;
	float h=(pos.y-atmo_bottom_radius)/height;
	h=h*(10.0*(weather.x*weather.x)+2.0)-1.0;
	return -(h*h)+1.0;
	
	/*float dAlt=pos.y-atmo_bottom_radius;
	return dAlt*(dAlt-height)*oneOverHeight*oneOverHeight*4.0;*/
}

float getGradient(vec3 pos)
{
	float height=atmo_top_radius-atmo_bottom_radius;
	return (pos.y-atmo_bottom_radius)/height;
}

float getDensity(vec3 pos)
{
	vec3 npos=normalize(pos);
	float dz = length(pos)/(atmo_top_radius-atmo_bottom_radius);
	vec3 coords=vec3(pos.x/50.0,pos.y/25.0,pos.z/50.0);
	
	vec4 shape = texture3D(shapeCloudTexture,coords); 
	vec4 shape2 = texture3D(shapeCloudTexture,coords*0.5); 
	vec4 detail = texture3D(detailCloudTexture,coords*10.0); 
	g_weather = texture2D(weatherTexture,coords.xz/vec2(5.0));
	float density=g_weather.r-(1.0-cDensity);//clamp(blend1D(0.0, 1.0, 0.8, g_weather.r, cDensity),0.0,1.0);	
	density=max(density,0.0);
	//density=density*(1.0/cDensity);
	//density=sqrt(density);
	//density=clamp((density-0.2)*1.25,0.0,1.0);
	//shape2.a=clamp((shape2.a-0.5)*2.0,0.0,1.0);
	
	density*=getHeightSignal(pos, g_weather.gb);
	density*=clamp((shape2.a*shape.x*shape.y*shape.z),0.0,1.0);
	
	
	density-=((1.0-detail.x)*0.03);//*detail.y*detail.z*0.2);
	density-=((1.0-detail.y)*0.02);
	density-=(detail.z*0.01);
	density=clamp(density,0.0,1.0);	
	density*=getGradient(pos);
	
	//shape2.a=clamp((shape2.a-0.5)*2.0,0.0,1.0);
	
	//float density=getSignal(pos);
	//density*=getSignal(pos);
	/*density*=getGradient(pos);
	density*=1.0-getGradient(pos);
	density*=clamp((shape2.a*shape.x*shape.y*shape.z),0.0,1.0);
	density-=((1.0-detail.x)*0.01);//*detail.y*detail.z*0.2);
	density-=((1.0-detail.y)*0.005);
	density-=(detail.z*0.01);
	*/
//	density=clamp(density,0.0,1.0);	
//	density=blend1D(0.0, 1.0, 0.8, density, cDensity);	
	
	return density;
}

float getSunVisibility(vec3 pos, vec3 sundir)
{
	float rSteps = (TOP_BOUND-BOTTOM_BOUND) / float(SHADOWS_STEPS);
    
    vec3 increment = sundir * rSteps;
    vec3 position = increment * 0.5 + pos;
    
    float transmittance = 0.0;
    
    for (int i = 0; i < SHADOWS_STEPS; i++, position += increment)
		transmittance += getDensity(position)*cloudsMat.x*rSteps;
    
    return clamp(exp2(-transmittance * rSteps),0.0,1.0);
}

vec3 getSunScattering(vec3 sunLight, float opticalLen, vec3 pos, float phase, vec3 vdir)
{
	//float integral = exp2(-opticalLen * log(2.0));
    float beersPowder = beerPowder(opticalLen,dot(sundir,vdir));//powder(opticalLen * log(2.0));
	vec3 sunlighting = sunLight*beersPowder*getSunVisibility(pos, sundir)*phase;//*opticalLen;//*opticalLen*phase*getSunVisibility(pos, sundir);//(sunColor */*vec3(1.0,1.0,1.0) */ getSunVisibility(pos, sundir) * beersPowder);
	
	return sunlighting;//*integral*PI;
}

float DistanceToCloudBoundary(vec3 r0, vec3 rd, vec3 s0, float sr) 
{
	float a= dot(rd,rd);
	float b= dot(rd,r0)*2.0;
	float c= dot(r0,r0)-sr*sr;
	
	float discriminant=b*b-4.0*a*c;
	if(discriminant<0.0)
		return -1;
	
	if(discriminant == 0) 
		return - 0.5 * b / a; 
		
	float q = (b > 0) ? -0.5 * (b + sqrt(discriminant)) : -0.5 * (b - sqrt(discriminant)); 		
	float x0 = q / a; 
    float x1 = c / q; 
	return max(x0,x1);
	//return max(0.0, max(-b-sqrt(discriminant)/(2.0*a),-b + sqrt(discriminant)/(2.0*a)));
}

float DistanceToCloudBoundary2(vec3 r0, vec3 rd, vec3 s0, float sr, out float x0, out float x1) 
{
	float a= dot(rd,rd);
	float b= dot(rd,r0)*2.0;
	float c= dot(r0,r0)-sr*sr;
	
	float discriminant=b*b-4.0*a*c;
	if(discriminant<0.0)
		return -1;
	
	if(discriminant == 0) 
		return - 0.5 * b / a; 
		
	float q = (b > 0) ? -0.5 * (b + sqrt(discriminant)) : -0.5 * (b - sqrt(discriminant)); 		
	x0 = q / a; 
    x1 = c / q; 
	if(x0<x1)
	{
		float tmp=x0;
		x0=x1;
		x1=tmp;
	}
	return max(x0,x1);
	//return max(0.0, max(-b-sqrt(discriminant)/(2.0*a),-b + sqrt(discriminant)/(2.0*a)));
}

float raySphereIntersect(vec3 r0, vec3 rd, vec3 s0, float sr) 
{
    // - r0: ray origin
    // - rd: normalized ray direction
    // - s0: sphere center
    // - sr: sphere radius
    // - Returns distance from r0 to first intersecion with sphere,
    //   or -1.0 if no intersection.
    float a = dot(rd, rd);
    vec3 s0_r0 = r0 - s0;
    float b = 2.0 * dot(rd, s0_r0);
    float c = dot(s0_r0, s0_r0) - (sr * sr);
	
    if (b*b - 4.0*a*c < 0.0) 
		return -1.0;
    
	return (-b - sqrt((b*b) - 4.0*a*c))/(2.0*a);
}

vec3 getCloudColor(vec3 color, vec3 rayOrigin, vec3 vdir)
{
	float dist=0.0;
	float density=0.0;
	vec3 start=rayOrigin;
	vec3 end=rayOrigin;
	
	float ds;
	float de;
	if(length(rayOrigin)>TOP_BOUND)
	{
		dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), TOP_BOUND, de, ds);
		if(dist<0.0)
			return color;	
		
		start=rayOrigin + vdir*ds;
//		end=rayOrigin + vdir*de;

		dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), BOTTOM_BOUND, de, ds);
		if(dist<0.0)
			return color;	
			
		end=rayOrigin + vdir*ds;
			
		/*dist=DistanceToCloudBoundary(rayOrigin, vdir, vec3(0.0,0.0,0.0), 200.0);
		if(dist<0.0)
			return color;	
		start=rayOrigin + vdir*dist;
		
		dist=DistanceToCloudBoundary(rayOrigin+vdir*(dist+2.0), vdir, vec3(0.0,0.0,0.0), 200.0);
		if(dist<0.0)
			return color;
		end=rayOrigin + vdir*dist;*/
	}	
	else
	{
		if(length(rayOrigin)<BOTTOM_BOUND)
		{
			dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), BOTTOM_BOUND, de, ds);
			if(dist<0.0)
				return color;	
		
			start=rayOrigin + vdir*de;
			
			dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), TOP_BOUND, de, ds);
			if(dist<0.0)
				return color;	
			end=rayOrigin + vdir*de;
		}
		else
		{
			dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), BOTTOM_BOUND, de, ds);
			end=rayOrigin + vdir*ds;
			if(dist<0.0)
			{
				dist=DistanceToCloudBoundary2(rayOrigin, vdir, vec3(0.0,0.0,0.0), TOP_BOUND, de, ds);
				end=rayOrigin + vdir*de;
			}
			if(dist<0.0)
				return color;
			start=rayOrigin;
		}
		
	}
	/*else
	{
		if(length(rayOrigin)<atmo_bottom_radius)
		{
			dist=DistanceToCloudBoundary(rayOrigin, vdir, vec3(0.0,0.0,0.0), atmo_bottom_radius);
			if(dist<0.0)
				return color;
			start=rayOrigin + vdir*dist;
			
			dist=DistanceToCloudBoundary(rayOrigin, vdir, vec3(0.0,0.0,0.0), atmo_top_radius);
			if(dist<0.0)
				return color;
			end=rayOrigin + vdir*dist;
		}
		else
		{
			dist=DistanceToCloudBoundary(rayOrigin, vdir, vec3(0.0,0.0,0.0), atmo_top_radius);
			if(dist<0.0)
			{	
				dist=DistanceToCloudBoundary(rayOrigin, vdir, vec3(0.0,0.0,0.0), atmo_bottom_radius);
				if(dist<0.0)
					return color;
			}
			end=rayOrigin + vdir*dist;
		}
	}		*/
	
	vec3 occluderVec=Wpos.xyz-rayOrigin;
	vec3 endVec=end-rayOrigin;
	vec3 startVec=start-rayOrigin;
	
	if(length(occluderVec)<length(startVec))
		return color;
	
	if(length(occluderVec)<length(endVec))
		endVec=occluderVec;
		
	float L=length(end-start);
	vec3 cloudColor=vec3(0.0,0.0,0.0);
	float transmittance=1.0;
	//float partialTransmittance=1.0;
	vec3 CurrentT=vec3(1.0,1.0,1.0);
	vec3 scattering=vec3(0.0,0.0,0.0);
	float meanDist=0.0;
	int countPos=0;
	float heightB=(start.y-earth_radius);
		
	if(L>0.0)
	{
		// get sky radiance to the top atmosphere bound
		vec3 skyTransmittance=vec3(0.0,0.0,0.0);
		vec3 skylight=GetSkyRadianceToAtmosphere(rayOrigin, end, 0.0, sundir, skyTransmittance);
		
	/*	vec3 ambientLight;
		vec3 sunLight =  GetSunAndSkyIrradiance( end, vec3(0.0,1.0,0.0), sundir, ambientLight);
		sunLight*=suncolor;
		ambientLight*=max(exp(-heightB*0.8),0.0);
		*/
		
		//vec3 ambientLight=GetIndirectIrradiance(end,sundir)*max(exp(-heightB*0.8),0.0);
		vec3 sunLight =  GetSunIrradiance(end, sundir)*suncolor;
		
		float lDotW = dot(sundir, vdir);
		float phase = phase2Lobes(lDotW);
	
		int NUM_STEPS=int(min(mix(float(MIN_STEPS),float(MAX_STEPS),1.0-exp(-L*0.005)),float(MAX_STEPS)));
		float stepLen=L/NUM_STEPS;
		vec3 cpos;
		float opticalLen=0.0;
		vec3 partialScattering;
		float t=bayer16(texcoord.st*screenSize)*stepLen;
		//float densityAcc=0.0;
		for(int i=0; i<NUM_STEPS; i++)
		{
			cpos=start+vdir*t;
			density=getDensity(cpos);
			opticalLen=getDensity(cpos)*cloudsMat.x*g_weather.z;
			const float clampedExt=max(0.0000001,opticalLen);
			opticalLen*=stepLen;
			
			if(opticalLen>0.0)
			{
				partialScattering = getSunScattering(sunLight, opticalLen, cpos, phase, vdir); 	// sun
				partialScattering += skylight;													// sky
				//partialScattering += ambientLight;											// ground
				
			//	partialScattering = (partialScattering - partialScattering*transmittance) / clampedExt;
				scattering+=partialScattering*opticalLen*transmittance; 	//vec3(getGradient(cpos)*exp(-opticalLen));
				//partialTransmittance=exp(-opticalLen);
				transmittance*=exp(-opticalLen);					
				/*if(transmittance<=0.001)
				{
					transmittance=0.0;
					break;
				}*/
				
				//if(transmittance<0.9999)
			//	meanDist+=partialTransmittance*t;
				meanDist+=t;
				countPos+=1;
			}
			if(density>=1.0)
				break;
			//densityAcc+=density;
			t+=stepLen;
		}
	}

	//cloudColor=mix(color,cloudColor+vec3(0.2,0.2,0.2),clamp(density,0.0,1.0));
	/*if(transmittance>0.8)
	{
		transmittance=1.0;
		scattering=vec3(0.0,0.0,0.0);
	}*/
	cloudColor = scattering + color*transmittance;//mix(color,max(cloudColor,vec3(1.0))*transmittance/*+scattering*/,clamp(density,0.0,1.0));//;//+scattering;//*transmittance+scattering;//mix(color,vec3(0.5,0.5,0.5)*transmittance,density);
	
	// compute areal perspective to an average depth into the cloud
	if(transmittance<1.0)
	{
		meanDist/=countPos;
		scattering=GetSkyRadianceToPoint(rayOrigin, start+vdir*meanDist, 0.0, sundir, CurrentT)*inScatter;
		cloudColor = (scattering + cloudColor*CurrentT);
		//cloudColor = mix(color, cloudColor, 1.0-max(transmittance,1.0));
	}
	return cloudColor; 
}

void main()
{
	vec4 pos;	
	pos.z = texture2D(depthTex,texcoord.st).r; 
	vec4 color=texture2D(colorTex,texcoord.st);

	// compute position
	pos.xy=VPOS*-pos.z;
	pos.w=1.0;
	vec3 finalColor=color.xyz;
	Wpos=eyeToWorld*pos;
	/*if(-pos.z>500000.0)
	{*/
	vec3 camdir=normalize(Wpos.xyz-campos.xyz);
	Wpos.xyz=(Wpos.xyz*length_unit)-earth_center.xyz;
	
	//GetSkyRadiance((campos.xyz*length_unit)-earth_center.xyz, camdir, 0.0, sundir, absorb);
	finalColor=getCloudColor(color.rgb, (campos.xyz*length_unit)-earth_center.xyz, camdir);
	
	/*}*/
	
/*	if(Wpos.x>-200.0 && Wpos.x<200.0 && Wpos.y>-200.0 && Wpos.y<200.0 && Wpos.z>-200.0 && Wpos.z<200.0)
		finalColor.xyz += texture3D(shapeCloudTexture,(Wpos.xyz/400.0)-vec3(0.5)).a; */
	gl_FragColor.xyz=finalColor;
	gl_FragColor.w=1.0;
}
