#define IMPROVED_INTEGRATION

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

uniform vec3		moondir;
uniform vec3		moonlightcolor;
uniform vec3 		zenithColor;

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

uniform vec2		wind;
uniform vec2		ditherRot;		

uniform mat4		reprojTM;
uniform mat4		projTM;

uniform vec3 		cloudsColor;
uniform float		waterDensityMul;

uniform float		scatterMul;

uniform float 		cloudsCoordsScale;

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

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;
vec4 g_weather;
float s_cloudAtmosphereInvSize=1.0/(atmo_top_radius-atmo_bottom_radius);
float s_posZ= texture2D(depthTex, texcoord.xy).r;

uniform float minSteps; //#define MIN_STEPS 				32//8//32
uniform float maxSteps; //#define MAX_STEPS 				64//64//128

#define SHADOWS_STEPS			4//8

#define TOP_BOUND 				atmo_top_radius
#define BOTTOM_BOUND			atmo_bottom_radius
#define CLOUDS_RANGE_MIN		4900000.0
#define CLOUDS_RANGE_MAX		5000000.0
#define TEMPORAL_FILTER_ALPHA	0.7

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 hash(float n) { return fract(sin(n) * 1e4); }
float hash(vec2 p) { return fract(1e4 * sin(17.0 * p.x + p.y * 0.1) * (0.1 + abs(sin(p.y * 13.0 + p.x)))); }

float noise(vec2 x) {
    vec2 i = floor(x);
    vec2 f = fract(x);

	// Four corners in 2D of a tile
	float a = hash(i);
    float b = hash(i + vec2(1.0, 0.0));
    float c = hash(i + vec2(0.0, 1.0));
    float d = hash(i + vec2(1.0, 1.0));

    // Simple 2D lerp using smoothstep envelope between the values.
	// return vec3(mix(mix(a, b, smoothstep(0.0, 1.0, f.x)),
	//			mix(c, d, smoothstep(0.0, 1.0, f.x)),
	//			smoothstep(0.0, 1.0, f.y)));

	// Same code, with the clamps in smoothstep and common subexpressions
	// optimized away.
    vec2 u = f * f * (3.0 - 2.0 * f);
	return mix(a, b, u.x) + (c - a) * u.y * (1.0 - u.x) + (d - b) * u.x * u.y;
}


vec2 getScreenPos(vec4 pos)
{
	vec4 ppos;
	vec2 spos;
		
	// compute texture space pos
	ppos=projTM*pos;
	ppos=ppos/ppos.w;
	spos=vec2(1.0)-ppos.xy;
	return spos;
}

vec4 getEyePos(vec2 coords)
{
	vec4 pos;
	pos.z = s_posZ;//texture2D(depthTex, coords.xy).r;
	pos.xy=VPOS*-pos.z;
	pos.w=1.0;
	return pos;
}

vec2 reproject(vec2 coords)
{
	vec4 eye=getEyePos(coords.st);
	vec4 prevEye=reprojTM*eye;
	return getScreenPos(prevEye);
}

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

float beerPowder(float od, float rain)
{
	return exp(-od*rain)*(1.0-exp(-od*2.0));//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)
{
    float m = 0.5;//0.6;
    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 gradientRemap(float H, float min, float max, float newMin, float newMax)
{
	return newMin + (((H-min)/(max-min))*(newMax-newMin));
}

float getHeightSignal(float H, vec2 weather)
{
	/*float h=(H-atmo_bottom_radius)*s_cloudAtmosphereInvSize;
	float type=weather.x;
	float stratus = gradientRemap(h,0.0,0.07,0.0,1.0)*gradientRemap(h,0.2,0.3,1.0,0.0);
	float stratocumulus = gradientRemap(h,0.07,0.25,0.0,1.0)*gradientRemap(h,0.45,0.65,1.0,0.0);
	float cumulus = gradientRemap(h,0.0,0.1,0.0,1.0)*gradientRemap(h,0.75,0.99,1.0,0.0);
	
	float mixStratus=clamp(mix(stratus, stratocumulus, type*2.0),0.0,1.0);
	return mix(mixStratus, cumulus, clamp(type-0.5,0.0,1.0)*2.0);*/
	float h=(H-atmo_bottom_radius)*s_cloudAtmosphereInvSize;///s_cloudAtmosphereSize;
	h=h*(10.0*(1.0-weather.x)+2.0)-1.0;
	h=h*h;
	return -(h*h)+1.0;
}

float getGradientByType(float H, vec2 weather)
{
	float h=(H-atmo_bottom_radius)*s_cloudAtmosphereInvSize;	
	float type=weather.x;
	float stratus = gradientRemap(h,0.0,0.07,0.0,1.0);
	float stratocumulus = gradientRemap(h,0.07,0.25,0.0,1.0);
	float cumulus = gradientRemap(h,0.0,0.1,0.0,1.0);
	
	float mixStratus=clamp(mix(stratus, stratocumulus, type*2.0),0.0,1.0);
	return mix(mixStratus, cumulus, clamp(type-0.5,0.0,1.0)*2.0);
}

float getGradient(float h)
{
	return (h-atmo_bottom_radius)*s_cloudAtmosphereInvSize;///s_cloudAtmosphereSize;
}

float getDirectionalGradient(vec3 pos, vec3 ray)
{
	float h=max(dot(pos,ray),0.0);
	return clamp((h-atmo_bottom_radius)*s_cloudAtmosphereInvSize,0.0,1.0);///s_cloudAtmosphereSize;
}

float getDensity(vec3 pos)
{
	vec3 npos=normalize(pos);
	float dz = length(pos)*s_cloudAtmosphereInvSize;///s_cloudAtmosphereSize;
	vec3 coords=vec3(pos.x/50.0,pos.y/25.0,pos.z/50.0);
	
	g_weather = texture2D(weatherTexture,(coords.xz/vec2(5.0)*vec2(cloudsCoordsScale))+wind.xy);
	float density=g_weather.x;

	density=max(density-(1.0-cDensity),0.0)*(1.0/(cDensity+0.0001));
	density*=getHeightSignal(pos.y, g_weather.gb);
		
	vec4 shape = texture3D(shapeCloudTexture,coords+vec3(wind.x*2.0,0.0,wind.y*2.0)); 

	density*=clamp(((shape.a*shape.x*shape.y*shape.z)),0.0,1.0);
	vec4 detail = texture3D(detailCloudTexture,coords*10.0+vec3(wind.x*6.0,0.0,wind.y*6.0)); 
	density-=((1.0-detail.x)*0.01);
	density-=((1.0-detail.y)*0.005);
	density-=(detail.z*0.005);
	
	density=clamp(density,0.0,1.0);

	return density;
}

float getSunVisibility(vec3 pos, vec3 dir)
{
	float rSteps = (TOP_BOUND-BOTTOM_BOUND) / float(SHADOWS_STEPS);
    
    vec3 increment = dir * 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 rain)
{
	float beersPowder = beerPowder(opticalLen, rain);
	vec3 sunlighting = sunLight*beersPowder*getSunVisibility(pos, sundir)*phase;
	
	return sunlighting;//*integral*PI;
}

vec3 getMoonScattering(vec3 moonLight, float opticalLen, vec3 pos, float phase, vec3 vdir, float rain)
{
	float beersPowder = beerPowder(opticalLen, rain);
	vec3 moonlighting = moonLight*beersPowder*getSunVisibility(pos, moondir)*phase;
	
	return moonlighting;//*integral*PI;
}


vec3 getMultiScattering(int N, vec3 sunLight, float opticalLen, vec3 pos, float lDotW, vec3 vdir, float rain)
{
	//float integral = exp2(-opticalLen * log(2.0));
	vec3 sunlighting=vec3(0.0,0.0,0.0);
	float 	a=0.5,
			b=0.5,
			c=0.5;
			
	for(int i=0; i<N; i++)
	{
		float phase=phase2Lobes(lDotW*pow(c,i));
		float beersPowder = beerPowder(opticalLen, rain)*pow(a,i);//powder(opticalLen * log(2.0));
		float Lext=getSunVisibility(pos, sundir)*pow(b,i);

		sunlighting += sunLight*beersPowder*Lext*phase;
	}
	
	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, out float transmittance, out vec4 preScattering)
{
	preScattering.rgb = color;
	preScattering.a = 1.0;
	transmittance=1.0;
	
	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;
		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;
	}	
	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;
		}
		
	}
	
	vec3 occluderVec=Wpos.xyz-rayOrigin;
	vec3 endVec=end-rayOrigin;
	vec3 startVec=start-rayOrigin;
	
	if(start.y < earth_radius)
		return color;
		
/*	if(length(start)>50000.0)
		return color;	
	*/
	if(length(occluderVec)<length(startVec))
	{
		float z0=-texture2D(depthTex, texcoord.xy+vec2(2.0/screenSize.x,0.0)).r;
		if(z0<wFar)
		{
			float z1=-texture2D(depthTex, texcoord.xy-vec2(2.0/screenSize.x,0.0)).r;
			if(z1<wFar)
			{
				float z2=-texture2D(depthTex, texcoord.xy+vec2(0.0,2.0/screenSize.y)).r;
				if(z2<wFar)
				{
					float z3=-texture2D(depthTex, texcoord.xy-vec2(0.0,2.0/screenSize.y)).r;
					if(z3<wFar)
						return color;
				}
			}
		}
	}
	
	if(length(occluderVec)<length(endVec))
		endVec=occluderVec;
		
	float L=length(end-start);
	vec3 cloudColor=color;
	vec3 CurrentT;
	vec3 scattering=vec3(0.0,0.0,0.0);
	float meanDist=0.0;
	int countPos=0;
	//float heightB=(start.y-earth_radius);
	int NUM_STEPS=0;
	//float totalDensity=0.0;

	if(L>0.0)
	{
		float deltaSteps=maxSteps/64.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);
		
		// desaturate skylight
		/*vec3 lum=dot(skylight.rgb,vec3(0.3,0.59,0.11));
		skylight.rgb=mix(skylight.rgb,lum,0.8);
		*/
		vec3 skylightDelta = skylight*deltaSteps;
		
		vec3 zenithColorDelta=zenithColor*deltaSteps;
		
		vec3 sunLight =  GetSunIrradiance(end, sundir)*suncolor;
		vec3 sunLightDelta=sunLight*deltaSteps;
		float sunlightLen=length(sunLightDelta);
		
		float lDotW = dot(sundir, vdir);
		float phase = phase2Lobes(lDotW);
		
		vec3 moonLightDelta=moonlightcolor*deltaSteps;
		float moonlightLen=length(moonLightDelta);
		float moon_lDotW = dot(moondir, vdir);
		float moon_phase = phase2Lobes(moon_lDotW);
	
		NUM_STEPS=int(min(mix(maxSteps,minSteps,exp(-L*0.025)),maxSteps));
		float stepLen=L/NUM_STEPS;
		vec3 cpos;
		float opticalLen=0.0;
		vec3 partialScattering;
		float t=bayer32((texcoord.st*screenSize*0.5)+vec2(ditherRot.x*ditherRot.y,ditherRot.x*(1.0-ditherRot.y)))*stepLen;
		
		for(int i=0; i<NUM_STEPS; i++)
		{
			cpos=start+vdir*t;
			density=getDensity(cpos);
			opticalLen=density*cloudsMat.x*stepLen;
			
			if(opticalLen>0.0)
			{
				partialScattering=vec3(0.0,0.0,0.0);
				if(sunlightLen>0.0)
				{
					partialScattering += getSunScattering(sunLightDelta, opticalLen, cpos, phase, vdir, g_weather.b*100.0*waterDensityMul); 				// sun
					partialScattering += pow(getGradient(cpos.y),0.8)*sunLightDelta/(2.0*PI)*clamp(dot(sundir,vec3(0.0,1.0,0.0)),0.0,1.0)*scatterMul;		//		pow(getDirectionalGradient(cpos,sundir),0.8)*sunLight/(2.0*PI)*0.5;				
				}
				if(moonlightLen>0.0)
				{
					partialScattering += getMoonScattering(moonLightDelta, opticalLen, cpos, moon_phase, vdir, g_weather.b*100.0*waterDensityMul); 	// sun
					partialScattering += pow(getGradient(cpos.y),0.8)*moonLightDelta/(2.0*PI)*0.25;
				}						
				partialScattering += pow(getGradient(cpos.y),0.8)*skylightDelta*(2.0*PI)*cloudsColor;									// sky
				partialScattering += pow(getGradient(cpos.y),0.8)*zenithColorDelta*(2.0*PI)*cloudsColor;
				
				#ifdef IMPROVED_INTEGRATION
					float clampedExt=max(0.0000001,density*cloudsMat.x);
					float partialTrans=exp(-opticalLen);				
					partialScattering = (partialScattering - partialScattering*partialTrans) / clampedExt;
					scattering+=partialScattering*opticalLen*transmittance; 				//vec3(getGradient(cpos)*exp(-opticalLen));
					transmittance*=partialTrans;	
				#else
					transmittance*=exp(-opticalLen);				
					scattering+=partialScattering*opticalLen*transmittance; 				//vec3(getGradient(cpos)*exp(-opticalLen));
				#endif
				
				meanDist+=t*transmittance;
				countPos+=1;
			}
			
			if(transmittance<0.01)
				break;
			t+=stepLen;
		}
		
		// last step = cirrus clouds
		cpos=start+vdir*t*100.0;
		float baseCirrus=texture2D(cirrusTex,vec2(cpos.x,cpos.z)*0.0001+wind*2.0).g;
		density=baseCirrus*texture2D(cirrusTex,vec2(cpos.x,cpos.z)*0.002-wind).b;
		density+=(1.0-baseCirrus)*texture2D(cirrusTex,vec2(cpos.x,cpos.z)*0.001+wind*4.0).r;
		
		density=max(density-(1.0-cDensity*0.8),0.0);
		density=density*(1.0-clamp(length(start.xz-cpos.xz)/5000.0,0.0,1.0));
		opticalLen=density*cloudsMat.x*0.0025;
		if(opticalLen>0.0)
		{
			partialScattering = sunLight*phase/(2.0*PI);
			partialScattering += skylight;//(2.0*PI);											// sky
			transmittance*=exp(-opticalLen);				
			
			scattering+=partialScattering*opticalLen*transmittance; 						//vec3(getGradient(cpos)*exp(-opticalLen));
			meanDist+=t*transmittance;
			countPos+=1;
		}
	}
	
	//cloudColor=mix(color,cloudColor+vec3(0.2,0.2,0.2),clamp(density,0.0,1.0));
	preScattering=vec4(scattering,transmittance);
	if(transmittance<1.0)
	{
		vec4 prevST=preScattering;
		vec2 uv=reproject(texcoord);
		if(uv.x>=0 && uv.y>=0 && uv.x<1.0 && uv.y<1.0)
		{
			prevST=texture2D(prevClouds,uv.st);
			/*if(prevST.a<1.0)
			{*/
				scattering=mix(scattering,prevST.rgb,TEMPORAL_FILTER_ALPHA);
				transmittance=mix(transmittance,prevST.a,TEMPORAL_FILTER_ALPHA);
			//}
		}
		preScattering.xyz=scattering;
		preScattering.w=transmittance;
		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
		meanDist/=countPos;
		scattering=GetSkyRadianceToPoint(rayOrigin, start+vdir*meanDist, 0.0, sundir, CurrentT)*inScatter;
		cloudColor = mix(cloudColor*CurrentT + scattering,cloudColor,transmittance);
	//	cloudColor = mix(cloudColor, color, clamp( (-s_posZ-CLOUDS_RANGE_MIN)/(CLOUDS_RANGE_MAX-CLOUDS_RANGE_MIN),0.0,1.0)/*clamp(length(start) * 0.0001, 0.0, 1.0)*/);
	}
	//gl_FragData[1]=PreScattering;
	return cloudColor; 
}

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

	// compute position
	pos=getEyePos(texcoord);
	vec3 finalColor=color.rgb;
	Wpos=eyeToWorld*pos;
	
	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);
	vec4 PreScattering;
	float transmittance;
	finalColor=getCloudColor(finalColor.rgb, (campos.xyz*length_unit)-earth_center.xyz, camdir, transmittance, PreScattering);
	
	gl_FragData[0].xyz=finalColor.rgb;//*(1.0 - step(1.0,transmittance));
	gl_FragData[0].w=1.0 - step(1.0,transmittance);
	
	gl_FragData[1]=PreScattering;
}
