/*
Real time PBR Volumetric Clouds by robobo1221.
Single scattering
Also includes volumetric light.
http://shadertoy.com/user/robobo1221

First ever somewhat PBR effect I decided to work on.
It uses the same algorithm to calculate the worldPosition as in: https://www.shadertoy.com/view/lstfR7

Feel free to fork and edit it. (Credit me please.)
Hope you enjoy!
*/

/*#define VOLUMETRIC_LIGHT
#define SPHERICAL_PROJECTION

#define cameraMode 2 					//1 is free rotation, 2 is still camera but free sun rotation

#define cloudSpeed 0.2
#define cloudHeight 1600.0
#define cloudThickness 500.0
#define cloudDensity 0.03

#define fogDensity 0.00003

#define volumetricCloudSteps 16			//Higher is a better result with rendering of clouds.
#define volumetricLightSteps 8			//Higher is a better result with rendering of volumetric light.

#define cloudShadowingSteps 8			//Higher is a better result with shading on clouds.
#define volumetricLightShadowSteps 4

#define rayleighCoeff (vec3(0.25, 0.5, 1.25) * 0.5e-5)	//Not really correct
#define mieCoeff vec3(0.9e-6)							//Not really correct

const float sunBrightness = 5.0;

varying vec2 		texcoord;

uniform vec2 		screenSize;
uniform sampler2D 	noiseTex;
uniform float		iTime;

uniform vec2		sunPos;
uniform float 		cDensity;

//////////////////////////////////////////////////////////////////

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))

//////////////////////////////////////////////////////////////////

#define cloudMinHeight cloudHeight
#define cloudMaxHeight (cloudThickness + cloudMinHeight)

#define sunPosition vec3(1.0, 1.0, 0.0)

const float pi=3.14159265359;
//const float pi = acos(-1.0);
const float rPi = 1.0 / pi;
const float hPi = pi * 0.5;
const float tau = pi * 2.0;

mat3 rotationMatrix(vec3 axis, float angle)
{
    axis = normalize(axis);
    float s = sin(angle);
    float c = cos(angle);
    float oc = 1.0 - c;
    
    float xx = axis.x * axis.x;
    float yy = axis.y * axis.y;
    float zz = axis.z * axis.z;
    
    float xy = axis.x * axis.y;
    float xz = axis.x * axis.z;
    float zy = axis.z * axis.y;
    
    return mat3(oc * xx + c, oc * xy - axis.z * s, oc * xz + axis.y * s,
                oc * xy + axis.z * s, oc * yy + c, oc * zy - axis.x * s, 
                oc * xz - axis.y * s, oc * zy + axis.x * s, oc * zz + c);
}

struct positionStruct
{
	vec2 texcoord;
    vec2 mousecoord;
    vec3 worldPosition;
    vec3 worldVector;
    vec3 sunVector;
} pos;

vec3 sphereToCart(vec3 sphere) {
    float cp = cos(sphere.x);
    float sp = sin(sphere.x);
    float ct = cos(sphere.y);
    float st = sin(sphere.y);
    float r = sphere.z;
    return vec3(r * cp * ct, r * st, r * sp * ct);
}

vec3 calculateWorldSpacePosition(vec2 p)
{
	p = p * 2.0 - 1.0;
    
    vec3 worldSpacePosition =  vec3(p.x, p.y, 1.0);
    
    #ifdef SPHERICAL_PROJECTION
		worldSpacePosition = sphereToCart(worldSpacePosition * vec3(pi, hPi, 1.0));
	#endif
    
    return worldSpacePosition;
}

void gatherPositions(inout positionStruct pos, vec2 fragCoord, vec2 mouseCoord, vec2 screenResolution)
{
	pos.texcoord = fragCoord;//*vec2(1.0,0.5)+vec2(0.0,0.5);// / screenResolution;
    pos.mousecoord = mouseCoord;// / screenResolution;
    
    pos.mousecoord = pos.mousecoord.x < 0.001 ? vec2(0.4, 0.64) : pos.mousecoord;
    
    pos.worldPosition = calculateWorldSpacePosition(pos.texcoord);
    pos.sunVector = normalize(calculateWorldSpacePosition(pos.mousecoord));
    
    pos.worldVector = normalize(pos.worldPosition);
}

///////////////////////////////////////////////////////////////////////////////////

#define d0(x) (abs(x) + 1e-8)
#define d02(x) (abs(x) + 1e-3)

const vec3 totalCoeff = rayleighCoeff + mieCoeff;

vec3 scatter(vec3 coeff, float depth){
	return coeff * depth;
}

vec3 absorb(vec3 coeff, float depth){
	return exp2(scatter(coeff, -depth));
}

float calcParticleThickness(float depth){
   	
    depth = depth * 2.0;
    depth = max(depth + 0.01, 0.01);
    depth = 1.0 / depth;
    
	return 100000.0 * depth;   
}

float calcParticleThicknessH(float depth){
   	
    depth = depth * 2.0 + 0.1;
    depth = max(depth + 0.01, 0.01);
    depth = 1.0 / depth;
    
	return 100000.0 * depth;   
}

float calcParticleThicknessConst(const float depth){
    
	return 100000.0 / max(depth * 2.0 - 0.01, 0.01);   
}

float rayleighPhase(float x){
	return 0.375 * (1.0 + x*x);
}

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 miePhaseSky(float x, float depth)
{
 	return hgPhase(x, exp2(-0.0000003 * depth));
}

vec3 calcAtmosphericScatter(positionStruct pos, out vec3 absorbLight){
    const float ln2 = log(2.0);
    
    float lDotW = dot(pos.sunVector, pos.worldVector);
    float lDotU = dot(pos.sunVector, vec3(0.0, 1.0, 0.0));
    float uDotW = dot(vec3(0.0, 1.0, 0.0), pos.worldVector);
    
	float opticalDepth = calcParticleThickness(uDotW);
    float opticalDepthLight = calcParticleThickness(lDotU);
    
    vec3 scatterView = scatter(totalCoeff, opticalDepth);
    vec3 absorbView = absorb(totalCoeff, opticalDepth);
    
    vec3 scatterLight = scatter(totalCoeff, opticalDepthLight);
         absorbLight = absorb(totalCoeff, opticalDepthLight);
    	 
    vec3 absorbSun = abs(absorbLight - absorbView) / d0((scatterLight - scatterView) * ln2);
    
    vec3 mieScatter = scatter(mieCoeff, opticalDepth) * miePhaseSky(lDotW, opticalDepth);
    vec3 rayleighScatter = scatter(rayleighCoeff, opticalDepth) * rayleighPhase(lDotW);
    
    vec3 scatterSun = mieScatter + rayleighScatter;
    
    vec3 sunSpot = smoothstep(0.9999, 0.99993, lDotW) * absorbView * sunBrightness;
    
    return (scatterSun * absorbSun + sunSpot) * sunBrightness;
}

vec3 calcAtmosphericScatterTop(positionStruct pos){
    const float ln2 = log(2.0);
    
    float lDotU = dot(pos.sunVector, vec3(0.0, 1.0, 0.0));
    
	float opticalDepth = calcParticleThicknessConst(1.0);
    float opticalDepthLight = calcParticleThickness(lDotU);
    
    vec3 scatterView = scatter(totalCoeff, opticalDepth);
    vec3 absorbView = absorb(totalCoeff, opticalDepth);
    
    vec3 scatterLight = scatter(totalCoeff, opticalDepthLight);
    vec3 absorbLight = absorb(totalCoeff, opticalDepthLight);
    
    vec3 absorbSun = d02(absorbLight - absorbView) / d02((scatterLight - scatterView) * ln2);
    
    vec3 mieScatter = scatter(mieCoeff, opticalDepth) * 0.25;
    vec3 rayleighScatter = scatter(rayleighCoeff, opticalDepth) * 0.375;
    
    vec3 scatterSun = mieScatter + rayleighScatter;
    
    return (scatterSun * absorbSun) * sunBrightness;
}

float Get3DNoise(vec3 pos) 
{
    float p = floor(pos.z);
    float f = pos.z - p;
    
    const float invNoiseRes = 1.0 / 64.0;
    
    float zStretch = 20.0 * invNoiseRes;
    
    vec2 coord = pos.xy * invNoiseRes + (p * zStretch);
    
	vec2 noise = vec2(texture2D(noiseTex, coord).x-cDensity, texture2D(noiseTex, coord + zStretch).x-cDensity);
    
    return mix(noise.x, noise.y, f);
}

float getClouds(vec3 p)
{
    if (p.y < cloudMinHeight || p.y > cloudMaxHeight)
        return 0.0;
    
    float time = iTime * cloudSpeed;
    vec3 movement = vec3(time, 0.0, time);
    
    vec3 cloudCoord = (p * 0.001) + movement;
    
	float noise = Get3DNoise(cloudCoord) * 0.5;
    	  noise += Get3DNoise(cloudCoord * 2.0 + movement) * 0.25;
    	  noise += Get3DNoise(cloudCoord * 7.0 - movement) * 0.125;
    	  noise += Get3DNoise((cloudCoord + movement) * 16.0) * 0.0625;
    
    const float top = 0.004;
    const float bottom = 0.01;
    
    float horizonHeight = p.y - cloudMinHeight;
    float treshHold = (1.0 - exp2(-bottom * horizonHeight)) * exp2(-top * horizonHeight);
    
    float clouds = smoothstep(0.55, 0.6, noise);
          clouds *= treshHold;
    
    return clouds * cloudDensity;
}
    
float getCloudShadow(vec3 p, positionStruct pos)
{
	const int steps = volumetricLightShadowSteps;
    const float rSteps = cloudThickness / float(steps);
    
    vec3 increment = pos.sunVector * rSteps / pos.sunVector.y;
    vec3 position = pos.sunVector * (cloudMinHeight - p.y) / pos.sunVector.y + p;
    
    float transmittance = 0.0;
    
    for (int i = 0; i < steps; i++, position += increment)
    {
		transmittance += getClouds(position);
    }
    
    return exp2(-transmittance * rSteps);
}

float getSunVisibility(vec3 p, positionStruct pos)
{
	const int steps = cloudShadowingSteps;
    const float rSteps = cloudThickness / float(steps);
    
    vec3 increment = pos.sunVector * rSteps;
    vec3 position = increment * 0.5 + p;
    
    float transmittance = 0.0;
    
    for (int i = 0; i < steps; i++, position += increment)
    {
		transmittance += getClouds(position);
    }
    
    return exp2(-transmittance * rSteps);
}

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

vec3 getVolumetricCloudsScattering(float opticalDepth, float phase, vec3 p, vec3 sunColor, vec3 skyLight, positionStruct pos)
{
    float intergal = exp2(-opticalDepth * log(2.0));
    
	vec3 sunlighting = (sunColor * getSunVisibility(p, pos)) * (opticalDepth * phase) * hPi * sunBrightness;
    vec3 skylighting = skyLight * (opticalDepth * intergal) * 0.25 * rPi;
    
    return (sunlighting + skylighting) * pi;
}

float getHeightFogOD(float height)
{
	const float falloff = 0.001;
    
    return exp2(-height * falloff) * fogDensity;
}

vec3 getVolumetricLightScattering(float opticalDepth, float phase, vec3 p, vec3 sunColor, vec3 skyLight, positionStruct pos)
{
    float intergal = exp2(-opticalDepth * log(2.0));
    
	vec3 sunlighting = sunColor * (opticalDepth * phase) * hPi * sunBrightness;
         sunlighting *= getCloudShadow(p, pos);
    vec3 skylighting = skyLight * (opticalDepth * intergal) * 0.25 * rPi;
    
    return (sunlighting + skylighting) * pi;
}

vec3 calculateVolumetricLight(positionStruct pos, vec3 color, float dither, vec3 sunColor)
{
    #ifndef VOLUMETRIC_LIGHT
    	return color;
    #endif
    
	const int steps = volumetricLightSteps;
    const float iSteps = 1.0 / float(steps);
    
    vec3 increment = pos.worldVector * cloudMinHeight / clamp(pos.worldVector.y, 0.1, 1.0) * iSteps;
    vec3 rayPosition = increment * dither;
    
    float stepLength = length(increment);
    
    vec3 scattering = vec3(0.0);
    vec3 transmittance = vec3(1.0);
    
    float lDotW = dot(pos.sunVector, pos.worldVector);
    float phase = hgPhase(lDotW, 0.8);
    
    vec3 skyLight = calcAtmosphericScatterTop(pos);
    
    for (int i = 0; i < steps; i++, rayPosition += increment)
    {
        float opticalDepth = getHeightFogOD(rayPosition.y) * stepLength;
        
        if (opticalDepth <= 0.0)
            continue;
        
		scattering += getVolumetricLightScattering(opticalDepth, phase, rayPosition, sunColor, skyLight, pos) * transmittance;
        transmittance *= exp2(-opticalDepth);
    }
    
    return color * transmittance + scattering;
}

vec3 calculateVolumetricClouds(positionStruct pos, vec3 color, float dither, vec3 sunColor)
{
	const int steps = volumetricCloudSteps;
    const float iSteps = 1.0 / float(steps);
    
    if (pos.worldVector.y < 0.0)
        return color;
    
    vec3 startPosition = (pos.worldPosition * cloudMinHeight) / pos.worldPosition.y;
    
    vec3 increment = pos.worldVector * cloudThickness / clamp(pos.worldVector.y, 0.1, 1.0) * iSteps;
    vec3 cloudPosition = increment * dither + startPosition;
    
    float stepLength = length(increment);
    
    vec3 scattering = vec3(0.0);
    float transmittance = 1.0;
    
    float lDotW = dot(pos.sunVector, pos.worldVector);
    float phase = phase2Lobes(lDotW);
    
    vec3 skyLight = calcAtmosphericScatterTop(pos);
    
    for (int i = 0; i < steps; i++, cloudPosition += increment)
    {
        float opticalDepth = getClouds(cloudPosition) * stepLength;
        
        if (opticalDepth <= 0.0)
            continue;
        
		scattering += getVolumetricCloudsScattering(opticalDepth, phase, cloudPosition, sunColor, skyLight, pos) * transmittance;
        transmittance *= exp2(-opticalDepth);
    }
    
    return mix(color * transmittance + scattering, color, clamp(length(startPosition) * 0.00001, 0.0, 1.0));
}

vec3 robobo1221Tonemap(vec3 color)
{
    #define rTOperator(x) (x / sqrt(x*x+1.0))

    float l = length(color);

    color = mix(color, color * 0.5, l / (l+1.0));
    color = rTOperator(color);

    return color;
}

void main()
{
	gatherPositions(pos, texcoord, sunPos, screenSize);
    
    float dither = bayer16(texcoord*screenSize);
    
    vec3 lightAbsorb = vec3(0.0);

    vec3 color = vec3(0.0);
    color = calcAtmosphericScatter(pos, lightAbsorb);
	color = calculateVolumetricClouds(pos, color, dither, lightAbsorb);
    //    color = calculateVolumetricLight(pos, color, dither, lightAbsorb);
    gl_FragColor = vec4(color, 1.0);
}

varying vec2 		texcoord;

uniform vec2 		screenSize;
uniform sampler2D 	noiseTex;
uniform float		iTime;

uniform vec2		sunPos;
uniform float 		cDensity;

#define float2 vec2
#define float3 vec3
#define float4 vec4

#define PI 3.1415926535f
#define PI_2 (3.1415926535f * 2.0)

#define EPSILON 1e-5

#define SAMPLES_NUMS 16

float saturate(float x){ return clamp(x, 0.0, 1.0); }

struct ScatteringParams
{
    float sunRadius;
	float sunRadiance;

	float mieG;
	float mieHeight;

	float rayleighHeight;

	float3 waveLambdaMie;
	float3 waveLambdaOzone;
	float3 waveLambdaRayleigh;

	float earthRadius;
	float earthAtmTopRadius;
	float3 earthCenter;
};

float3 ComputeSphereNormal(float2 coord, float phiStart, float phiLength, float thetaStart, float thetaLength)
{
	float3 normal;
	normal.x = -sin(thetaStart + coord.y * thetaLength) * sin(phiStart + coord.x * phiLength);
	normal.y = -cos(thetaStart + coord.y * thetaLength);
	normal.z = -sin(thetaStart + coord.y * thetaLength) * cos(phiStart + coord.x * phiLength);
	return normalize(normal);
}

float2 ComputeRaySphereIntersection(float3 position, float3 dir, float3 center, float radius)
{
	float3 origin = position - center;
	float B = dot(origin, dir);
	float C = dot(origin, origin) - radius * radius;
	float D = B * B - C;

	float2 minimaxIntersections;
	if (D < 0.0)
	{
		minimaxIntersections = float2(-1.0, -1.0);
	}
	else
	{
		D = sqrt(D);
		minimaxIntersections = float2(-B - D, -B + D);
	}

	return minimaxIntersections;
}

float3 ComputeWaveLambdaRayleigh(float3 lambda)
{
	const float n = 1.0003;
	const float N = 2.545E25;
	const float pn = 0.035;
	const float n2 = n * n;
	const float pi3 = PI * PI * PI;
	const float rayleighConst = (8.0 * pi3 * pow(n2 - 1.0,2.0)) / (3.0 * N) * ((6.0 + 3.0 * pn) / (6.0 - 7.0 * pn));
	return rayleighConst / (lambda * lambda * lambda * lambda);
}

float ComputePhaseMie(float theta, float g)
{
	float g2 = g * g;
	return (1.0 - g2) / pow(1.0 + g2 - 2.0 * g * saturate(theta), 1.5) / (4.0 * PI);
}

float ComputePhaseRayleigh(float theta)
{
	float theta2 = theta * theta;
	return (theta2 * 0.75 + 0.75) / (4.0 * PI);
}

float ChapmanApproximation(float X, float h, float cosZenith)
{
	float c = sqrt(X + h);
	float c_exp_h = c * exp(-h);

	if (cosZenith >= 0.0)
	{
		return c_exp_h / (c * cosZenith + 1.0);
	}
	else
	{
		float x0 = sqrt(1.0 - cosZenith * cosZenith) * (X + h);
		float c0 = sqrt(x0);

		return 2.0 * c0 * exp(X - x0) - c_exp_h / (1.0 - c * cosZenith);
	}
}

float GetOpticalDepthSchueler(float h, float H, float earthRadius, float cosZenith)
{
	return H * ChapmanApproximation(earthRadius / H, h / H, cosZenith);
}

float3 GetTransmittance(ScatteringParams setting, float3 L, float3 V)
{
	float ch = GetOpticalDepthSchueler(L.y, setting.rayleighHeight, setting.earthRadius, V.y);
	return exp(-(setting.waveLambdaMie + setting.waveLambdaRayleigh) * ch);
}

float2 ComputeOpticalDepth(ScatteringParams setting, float3 samplePoint, float3 V, float3 L, float neg)
{
	float rl = length(samplePoint);
	float h = rl - setting.earthRadius;
	float3 r = samplePoint / rl;

	float cos_chi_sun = dot(r, L);
	float cos_chi_ray = dot(r, V * neg);

	float opticalDepthSun = GetOpticalDepthSchueler(h, setting.rayleighHeight, setting.earthRadius, cos_chi_sun);
	float opticalDepthCamera = GetOpticalDepthSchueler(h, setting.rayleighHeight, setting.earthRadius, cos_chi_ray) * neg;

	return float2(opticalDepthSun, opticalDepthCamera);
}

void AerialPerspective(ScatteringParams setting, float3 start, float3 end, float3 V, float3 L, bool infinite, out float3 transmittance, out float3 insctrMie, out float3 insctrRayleigh)
{
	float inf_neg = infinite ? 1.0 : -1.0;

	float3 sampleStep = (end - start) / float(SAMPLES_NUMS);
	float3 samplePoint = end - sampleStep;
	float3 sampleLambda = setting.waveLambdaMie + setting.waveLambdaRayleigh + setting.waveLambdaOzone;

	float sampleLength = length(sampleStep);

	float3 scattering = float3(0.0);
	float2 lastOpticalDepth = ComputeOpticalDepth(setting, end, V, L, inf_neg);

	for (int i = 1; i < SAMPLES_NUMS; i++, samplePoint -= sampleStep)
	{
		float2 opticalDepth = ComputeOpticalDepth(setting, samplePoint, V, L, inf_neg);

		float3 segment_s = exp(-sampleLambda * (opticalDepth.x + lastOpticalDepth.x));
		float3 segment_t = exp(-sampleLambda * (opticalDepth.y - lastOpticalDepth.y));
		
		transmittance *= segment_t;
		
		scattering = scattering * segment_t;
		scattering += exp(-(length(samplePoint) - setting.earthRadius) / setting.rayleighHeight) * segment_s;

		lastOpticalDepth = opticalDepth;
	}

	insctrMie = scattering * setting.waveLambdaMie * sampleLength;
	insctrRayleigh = scattering * setting.waveLambdaRayleigh * sampleLength;
}

float ComputeSkyboxChapman(ScatteringParams setting, float3 eye, float3 V, float3 L, out float3 transmittance, out float3 insctrMie, out float3 insctrRayleigh)
{
	bool neg = true;

	float2 outerIntersections = ComputeRaySphereIntersection(eye, V, setting.earthCenter, setting.earthAtmTopRadius);
	if (outerIntersections.y < 0.0) return 0.0;

	float2 innerIntersections = ComputeRaySphereIntersection(eye, V, setting.earthCenter, setting.earthRadius);
	if (innerIntersections.x > 0.0)
	{
		neg = false;
		outerIntersections.y = innerIntersections.x;
	}

	eye -= setting.earthCenter;

	float3 start = eye + V * max(0.0, outerIntersections.x);
	float3 end = eye + V * outerIntersections.y;

	AerialPerspective(setting, start, end, V, L, neg, transmittance, insctrMie, insctrRayleigh);

	bool intersectionTest = innerIntersections.x < 0.0 && innerIntersections.y < 0.0;
	return intersectionTest ? 1.0 : 0.0;
}

float4 ComputeSkyInscattering(ScatteringParams setting, float3 eye, float3 V, float3 L)
{
	float3 insctrMie = float3(0.0);
	float3 insctrRayleigh = float3(0.0);
	float3 insctrOpticalLength = float3(1.0);
	float intersectionTest = ComputeSkyboxChapman(setting, eye, V, L, insctrOpticalLength, insctrMie, insctrRayleigh);

	float phaseTheta = dot(V, L);
	float phaseMie = ComputePhaseMie(phaseTheta, setting.mieG);
	float phaseRayleigh = ComputePhaseRayleigh(phaseTheta);
	float phaseNight = 1.0 - saturate(insctrOpticalLength.x * EPSILON);

	float3 insctrTotalMie = insctrMie * phaseMie;
	float3 insctrTotalRayleigh = insctrRayleigh * phaseRayleigh;

	float3 sky = (insctrTotalMie + insctrTotalRayleigh) * setting.sunRadiance;

	float angle = saturate((1.0 - phaseTheta) * setting.sunRadius);
	float cosAngle = cos(angle * PI * 0.5);
	float edge = ((angle >= 0.9) ? smoothstep(0.9, 1.0, angle) : 0.0);
                         
	float3 limbDarkening = GetTransmittance(setting, -L, V);
	limbDarkening *= pow(float3(cosAngle), float3(0.420, 0.503, 0.652)) * mix(vec3(1.0), float3(1.2,0.9,0.5), edge) * intersectionTest;

	sky += limbDarkening;

	return float4(sky, phaseNight * intersectionTest);
}

float3 TonemapACES(float3 x)
{
	const float A = 2.51f;
	const float B = 0.03f;
	const float C = 2.43f;
	const float D = 0.59f;
	const float E = 0.14f;
	return (x * (A * x + B)) / (x * (C * x + D) + E);
}

float noise(float2 uv)
{
	return fract(dot(sin(uv.xyx * uv.xyy * 1024.0), float3(341896.483, 891618.637, 602649.7031)));
}

uniform vec3 earth_center;

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);
*/

// copy from https://www.shadertoy.com/view/4l2GzW
float r(float n)
{
 	return fract(cos(n*89.42)*343.42);
}
vec2 r(vec2 n)
{
 	return vec2(r(n.x*23.62-300.0+n.y*34.35),r(n.x*45.13+256.0+n.y*38.89)); 
}
float worley(vec2 n,float s)
{
    float dis = 2.0;
    for(int x = -1;x<=1;x++)
    {
        for(int y = -1;y<=1;y++)
        {
            vec2 p = floor(n/s)+vec2(x,y);
            float d = length(r(p)+vec2(x,y)-fract(n/s));
            if (dis>d)
            {
             	dis = d;   
            }
        }
    }
    return 1.0-dis;
	
}

// copy from https://www.shadertoy.com/view/4sc3z2

#define MOD3 vec3(.1031,.11369,.13787)

vec3 hash33(vec3 p3)
{
	p3 = fract(p3 * MOD3);
    p3 += dot(p3, p3.yxz+19.19);
    return -1.0 + 2.0 * fract(vec3((p3.x + p3.y)*p3.z, (p3.x+p3.z)*p3.y, (p3.y+p3.z)*p3.x));
}

float perlin_noise(vec3 p)
{
    vec3 pi = floor(p);
    vec3 pf = p - pi;
    
    vec3 w = pf * pf * (3.0 - 2.0 * pf);
    
    return 	mix(
        		mix(
                	mix(dot(pf - vec3(0, 0, 0), hash33(pi + vec3(0, 0, 0))), 
                        dot(pf - vec3(1, 0, 0), hash33(pi + vec3(1, 0, 0))),
                       	w.x),
                	mix(dot(pf - vec3(0, 0, 1), hash33(pi + vec3(0, 0, 1))), 
                        dot(pf - vec3(1, 0, 1), hash33(pi + vec3(1, 0, 1))),
                       	w.x),
                	w.z),
        		mix(
                    mix(dot(pf - vec3(0, 1, 0), hash33(pi + vec3(0, 1, 0))), 
                        dot(pf - vec3(1, 1, 0), hash33(pi + vec3(1, 1, 0))),
                       	w.x),
                   	mix(dot(pf - vec3(0, 1, 1), hash33(pi + vec3(0, 1, 1))), 
                        dot(pf - vec3(1, 1, 1), hash33(pi + vec3(1, 1, 1))),
                       	w.x),
                	w.z),
    			w.y);
}

float noise_sum(vec3 p)
{
    float f = 0.0;
    p = p * 4.0;
   // f += 2.0000 * perlin_noise(p); p = 2.0 * p;
    f += 1.0000 * perlin_noise(p); p = 2.0 * p;
    f += 0.5000 * perlin_noise(p); p = 2.0 * p;
	f += 0.2500 * perlin_noise(p); p = 2.0 * p;
	f += 0.1250 * perlin_noise(p); p = 2.0 * p;
	f += 0.0625 * perlin_noise(p); p = 2.0 * p;
    
    return f*0.5+0.5;
}

float worley_sum(vec2 n,float s)
{
    float f = 0.0;
    f += 1.0000 * worley(n,s);
    f += 0.5000 * worley(n*2.0,s);
	f += 0.2500 * worley(n*4.0,s);
	f += 0.1250 * worley(n*8.0,s);
	f += 0.0625 * worley(n*16.0,s);
    
    return f*0.5;
}

/*
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec3 pos = vec3(fragCoord.xy*0.001 * vec2(iResolution.x/iResolution.y, 1.0), iTime * 0.05);
    float f1= worley_sum(fragCoord.xy*0.5, 32.0);
    float f2= noise_sum(pos);
	fragColor = vec4(vec3(f1*0.5*f2),1.0);
}
*/

varying vec2 		texcoord;
uniform vec2 		screenSize;
//uniform sampler2D 	noiseTex;
uniform float		iTime;

void main()
{
	vec2 uv = gl_FragCoord.xy;//fragCoord.xy / iResolution.xy;
    
   /* vec2 mouse = sunPos;//iMouse.xy / iResolution.xy;
    if (mouse.x == 0.0 && mouse.y == 0.0)
        mouse = float2(0.5, 0.7);
    
    float3 V = ComputeSphereNormal(uv, 0.0, PI_2, 0.0, PI);
    float3 L = ComputeSphereNormal(float2(mouse.x, mouse.y), 0.0, PI_2, 0.0, PI);
    
	vec3 transmittance;
	vec3 scattering;
	float4 sky;

	sky.rgb=GetTrasmittanceTexture(uv);
	gl_FragColor = float4(sky.rgb, sky.a);*/
	vec3 pos = vec3(uv*0.0025 * vec2(screenSize.x/screenSize.y, 1.0), iTime*0.1);
    float f1= worley_sum(uv.xy*0.5, 32.0);
    float f2= noise_sum(pos);
	gl_FragColor = vec4(vec3(f1*f2),1.0);//vec4(vec3(f1*0.5+0.5/*+f2*/), 1.0);
}

/*void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
	vec2 uv = fragCoord.xy / iResolution.xy;
    
    vec2 mouse = iMouse.xy / iResolution.xy;
    if (mouse.x == 0.0 && mouse.y == 0.0)
        mouse = float2(0.5, 0.7);
    
    float3 V = ComputeSphereNormal(uv, 0.0, PI_2, 0.0, PI);
    float3 L = ComputeSphereNormal(float2(mouse.x, mouse.y), 0.0, PI_2, 0.0, PI);
    
	ScatteringParams setting;
	setting.sunRadius = 500.0;
	setting.sunRadiance = 20.0;
	setting.mieG = 0.96;
	setting.mieHeight = 1200.0;
	setting.rayleighHeight = 8000.0;
	setting.earthRadius = 6360000.0;
	setting.earthAtmTopRadius = 6420000.0;
	setting.earthCenter = float3(0, -setting.earthRadius, 0);
	setting.waveLambdaMie = float3(2e-7);
    
    // wavelength with 680nm, 550nm, 450nm
    setting.waveLambdaRayleigh = ComputeWaveLambdaRayleigh(float3(680e-9, 550e-9, 450e-9));
    
    // see https://www.shadertoy.com/view/MllBR2
	setting.waveLambdaOzone = float3(1.36820899679147, 3.31405330400124, 0.13601728252538) * 0.6e-6 * 2.504;
	
    float3 eye = float3(0,1000.0,0);
   	float4 sky = ComputeSkyInscattering(setting, eye, V, L);
    sky.rgb = TonemapACES(sky.rgb * 2.0);
    sky.rgb = pow(sky.rgb, float3(1.0 / 2.2)); // gamma
  //  sky.rgb += noise(uv*iTime) / 255.0; // dither
   
	fragColor = float4(sky.rgb, 1.0);
}*/