#version 130
#define IMPROVED_INTEGRATION

precision highp float;

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= -10000000.0;//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;///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(vec2 pos)
{
	vec2 coords=vec2(pos.x/50.0,pos.y/50.0);
	
	g_weather = texture2D(weatherTexture,(coords.xy/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=clamp(density,0.0,1.0);

	return density;
}

float DistanceToCloudBoundary(vec3 r0, vec3 rd, 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);
}

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;
	float BOUND=(TOP_BOUND+BOTTOM_BOUND)*0.5;
	
	dist=DistanceToCloudBoundary(rayOrigin, vdir, BOUND, de, ds);
	if(dist<0.0)
		return color;	
	
	start=rayOrigin + vdir*de;
	if(start.y<rayOrigin.y)
		return color;
	
	density=getDensity(start.xz);
	float opticalLen=density*cloudsMat.x*0.025;
	transmittance=exp(-opticalLen);//1.0-density;
	
	vec3 sunLight =  GetSunIrradiance(end, sundir)*suncolor;
	float lDotW = dot(sundir, vdir);
	float phase = phase2Lobes(lDotW);
	float beersPowder = beerPowder(opticalLen, g_weather.b*100.0*waterDensityMul);
	
	float moon_lDotW = dot(moondir, vdir);
	float moon_phase = phase2Lobes(moon_lDotW);
	
	vec3 skyTransmittance=vec3(0.0,0.0,0.0);
	vec3 skylight=GetSkyRadianceToAtmosphere(rayOrigin, end, 0.0, sundir, skyTransmittance);
			
	// DAY
	vec3 scattering=sunLight*phase*opticalLen*transmittance*beersPowder*4.0;//*sunvis;
	scattering += sunLight/(2.0*PI)*clamp(dot(sundir,vec3(0.0,1.0,0.0)),0.0,1.0)*scatterMul;
	scattering += skylight*(2.0*PI)*cloudsColor;
	
	// NIGHT
	scattering += moonlightcolor*moon_phase*opticalLen*transmittance*beersPowder*4.0;//*sunvis;
	scattering += moonlightcolor/(2.0*PI)*0.25;
	scattering += zenithColor*(2.0*PI)*cloudsColor;
	
	vec3 cloudColor=color*transmittance+scattering;
	
	vec3 CurrentT;
	scattering=GetSkyRadianceToPoint(rayOrigin, start, 0.0, sundir, CurrentT)*inScatter;
	cloudColor = mix(cloudColor*CurrentT + scattering,cloudColor,transmittance);
	
	return cloudColor; 
}

void main()
{
	vec4 pos;	
	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);
	vec4 PreScattering;
	float transmittance;
	finalColor=getCloudColor(finalColor.rgb, (campos.xyz*length_unit)-earth_center.xyz, camdir, transmittance, PreScattering);
	
	gl_FragColor.xyz=finalColor.rgb;//*(1.0 - step(1.0,transmittance));
	gl_FragColor.w=1.0 - transmittance;
}
