version 7
const float SunFadeDist=100000.0;
const float hFadeDistance= 100.0;
const float hdrExposure = 0.4;
const float IGRID1_SIZE = 1.0/5488.0;
const float IGRID2_SIZE = 1.0/392.0;
const float IGRID3_SIZE = 1.0/28.0;
const float IGRID4_SIZE = 1.0/2.0;
const float M_PI        = 3.141592657;
const vec3  seaColor    = vec3(0.004, 0.0157,0.047);
const float SUN_INTENSITY = 100.0;
const float earthRadius = 6360010.0;
const float Rg = 6360000.0;
const float Rt = 6420000.0;

uniform vec3 Range;  // scr_err*fov/screen_width,amplitude/2, OceanMaxDistance
uniform sampler2D skyIrradianceSampler;
#ifdef mask
in  vec2  texcoord; // texture coordinates
uniform sampler2D maskSampler;
#endif
#ifdef mask1
in  vec2  texcoord1; // texture coordinates
uniform sampler2D maskSampler1;
#endif
uniform sampler2D transmittanceSampler;
uniform sampler2D skySampler;
uniform sampler3D slopeVarianceSampler;
uniform sampler2DArray fftWavesSampler;

uniform vec3  lightvec; // sun direction in world space
uniform vec2  acolor;  //  alpha, sun_reflected_power

in  vec2  u; // coordinates in world space used to compute P(u)
in  vec3  P; // vector [wave point P(u) in world space,campos]
in  float distance2cam;
flat in float sunpower;

out vec4 fragColor;

// ---------------------------------------------------------------------------
// REFLECTED SUN RADIANCE
// ---------------------------------------------------------------------------

/*
_________________________________________________________*/
float Lambda(float cosTheta, float sigmaSq)
{
	float v = cosTheta / sqrt((1.0 - cosTheta * cosTheta) * (2.0 * sigmaSq));
  return exp(-v * v) / (2.0 * v * sqrt(M_PI));
}

/* L, V, N, Tx, Ty in world space
_________________________________________________________*/
float reflectedSunRadiance(vec3 L, vec3 V, vec3 N, vec3 Tx, vec3 Ty, vec2 sigmaSq)
{
    vec3 H = normalize(L + V);
    float zH = max(.1,dot(H, N));  // cos of facet normal zenith angle
    float zetax = dot(H, Tx) / zH;
    float zetay = dot(H, Ty) / zH;

    float zL  = max(dot(L, N),0.01); // cos of source zenith angle
    float zV  = max(dot(V, N),0.01); // cos of receiver zenith angle

    float p = exp(-0.5 * (zetax * zetax / sigmaSq.x + zetay * zetay / sigmaSq.y)) / (2.0 * M_PI * sqrt(sigmaSq.x * sigmaSq.y));

    vec2 csVLTy,csVLTx;
    csVLTy.x = dot(V, Ty);
    csVLTx.x = dot(V, Tx);
    csVLTy.y = dot(L, Ty);
    csVLTx.y = dot(L, Tx);
    vec2 tanVL = atan(csVLTy,csVLTx);

    vec2  cosVL2  = 1.0 / (1.0 + tanVL * tanVL);
    vec2 sigmaVL2 = sigmaSq.y * (1.0 - cosVL2) + sigmaSq.x * cosVL2;

    float fresnel = 0.02 + 0.98 * pow(1.0 - dot(V, H), 5.0);
    return fresnel * p / ((1.0 + Lambda(zL, sigmaVL2.y) + Lambda(zV, sigmaVL2.x)) * zV * pow(zH,4.0) * 4.0);
}

// ---------------------------------------------------------------------------
// REFLECTED SKY RADIANCE
// ---------------------------------------------------------------------------


/* V, N, Tx, Ty in world space
_________________________________________________________*/
vec2 U(vec2 zeta, vec3 V, vec3 N, vec3 Tx, vec3 Ty)
{
    vec3 f = normalize(vec3(-zeta, 1.0)); // tangent space
    vec3 F = f.x * Tx + f.y * Ty + f.z * N; // world space
    vec3 R = 2.0 * dot(F, V) * F - V;
    return R.xy / (1.0 + R.z);
}

/*
_________________________________________________________*/
float meanFresnel(float cosThetaV, float sigmaV)
{
	return pow(1.0 - cosThetaV, 5.0 * exp(-2.69 * sigmaV)) / (1.0 + 22.7 * pow(sigmaV, 1.5));
}

/* V, N in world space
_________________________________________________________*/
float meanFresnel(vec3 V, vec3 N, vec2 sigmaSq)
{
    vec2 v = V.xy; // view direction in wind space
    vec2 t = v * v / (1.0 - V.z * V.z); // cos^2 and sin^2 of view direction
    float sigmaV2 = dot(t, sigmaSq); // slope variance in view direction
    return meanFresnel(dot(V, N), sqrt(sigmaV2));
}

/* V, N, Tx, Ty in world space;
_________________________________________________________*/
vec3 meanSkyRadiance(vec3 V, vec3 N, vec3 Tx, vec3 Ty, vec2 sigmaSq)
{
    const float eps = 0.001;
    vec2 u0 = U(vec2(0.0), V, N, Tx, Ty);
    vec2 dux = 2.0 * (U(vec2(eps, 0.0), V, N, Tx, Ty) - u0) / eps * sqrt(sigmaSq.x);
    vec2 duy = 2.0 * (U(vec2(0.0, eps), V, N, Tx, Ty) - u0) / eps * sqrt(sigmaSq.y);

    return textureGrad(skySampler, u0 * (0.5 / 1.1) + 0.5, dux * (0.5 / 1.1), duy * (0.5 / 1.1)).rgb;
}

/*
_________________________________________________________*/
vec3 hdr(vec3 L)
{
    L=max(L,0.0);
    vec3 f =step(1.413 / hdrExposure,L);
    vec3 L1=pow(L * 0.38317 * hdrExposure, vec3(1.0 / 2.2));
    vec3 L2=1.0 - exp(-L * hdrExposure);
    L = mix(L1,L2,f);
#ifdef sRGBdecode
    f =step(0.04045,L);
    L1=L*(1.0/12.92);
    L2=pow( L*(1.0/1.055) + 0.055/1.055,vec3(2.4));
    L =mix(L1,L2,f);
#endif
    return L;
}

/*
_________________________________________________________*/
vec2 getIrradianceUV(float r, float muS)
{
    float uR = (r - Rg) / (Rt - Rg);
    float uMuS = (muS + 0.2) / (1.0 + 0.2);
    return vec2(uMuS, uR);
}

/*
_________________________________________________________*/
vec3 irradiance(sampler2D sampler, float r, float muS)
{
    vec2 uv = getIrradianceUV(r, muS);
    return texture(sampler, uv).rgb;
}
/*
_________________________________________________________*/
vec2 getTransmittanceUV(float r, float mu)
{
    float uR, uMu;
    uR = sqrt((r - Rg) / (Rt - Rg));
    uMu = atan((mu + 0.15) / (1.0 + 0.15) * tan(1.5)) / 1.5;
    return vec2(uMu, uR);
}

// transmittance(=transparency) of atmosphere for infinite ray (r,mu)
// (mu=cos(view zenith angle)), intersections with ground ignored
vec3 transmittance(float r, float mu)
{
    vec2 uv = getTransmittanceUV(r, mu);
    return texture(transmittanceSampler, uv).rgb;
}

// transmittance(=transparency) of atmosphere for infinite ray (r,mu)
// (mu=cos(view zenith angle)), or zero if ray intersects ground
vec3 transmittanceWithShadow(float r, float mu)
{
    return mu < -sqrt(1.0 - (Rg / r) * (Rg / r)) ? vec3(0.0) : transmittance(r, mu);
}


// incident sun light at given position (radiance)
// r=length(x)
// muS=dot(x,s) / r
vec3 sunRadiance(float muS)
{
 return transmittanceWithShadow(earthRadius, muS) * SUN_INTENSITY;
}

// incident sky light at given position, integrated over the hemisphere (irradiance)
// r=length(x)
// muS=dot(x,s) / r
vec3 skyIrradiance(float muS)
{
 return irradiance(skyIrradianceSampler, earthRadius, muS) * SUN_INTENSITY;
}


/*
_________________________________________________________*/
void sunRadianceAndSkyIrradiance(vec3 worldS, out vec3 sunL, out vec3 skyE)
{
    sunL = sunRadiance(worldS.z);
    skyE = skyIrradiance(worldS.z);
}



/*
_________________________________________________________*/
void main()
{
 float ka=clamp(1.0-(distance2cam-Range.z)/Range.z,0.0,1.0);
#ifdef mask
 float _ka = texture(maskSampler,texcoord).r;
#ifdef mask1
  _ka = max(_ka,texture(maskSampler1,texcoord1).r);
#endif
  ka*=_ka;
#endif
 if (ka<0.01) discard;
 vec3 V = normalize( P );
 vec2 slopes = texture(fftWavesSampler, vec3(u * IGRID1_SIZE, 1.0)).xy+
               texture(fftWavesSampler, vec3(u * IGRID2_SIZE, 1.0)).zw+
               texture(fftWavesSampler, vec3(u * IGRID3_SIZE, 2.0)).xy+
               texture(fftWavesSampler, vec3(u * IGRID4_SIZE, 2.0)).zw;
 vec3 N = normalize(vec3(-slopes.x, -slopes.y, 1.0));

 float Jxx = dFdx(u.x);
 float Jxy = dFdy(u.x);
 float Jyx = dFdx(u.y);
 float Jyy = dFdy(u.y);
 float A = Jxx * Jxx + Jyx * Jyx;
 float B = Jxx * Jxy + Jyx * Jyy;
 float C = Jxy * Jxy + Jyy * Jyy;
 float ua = pow(A / 10.0, 0.25);
 float ub = 0.5 + 0.5 * B / sqrt(A * C);
 float uc = pow(C / 10.0, 0.25);
 vec2 sigmaSq = texture(slopeVarianceSampler, vec3(ua, ub, uc)).xy;

 sigmaSq = max(sigmaSq, 2e-5);

 vec3 Ty = normalize(vec3(0.0, N.z, -N.y));
 vec3 Tx = cross(Ty, N);
 float fresnel = 0.02 + 0.98 * meanFresnel(V, N, sigmaSq);
 vec3 Lsun,Esky,Lsky,Lsea;
 sunRadianceAndSkyIrradiance(lightvec, Lsun, Esky);

 float sun_fade=pow(clamp(1.0-distance2cam/SunFadeDist,0.0,1.0),4.0); // fading sun blicks 'cos distance
 float lens    = length(Lsun);
 float sunlen  = sunpower*acolor.y*sun_fade * reflectedSunRadiance(lightvec, V, N, Tx, Ty, sigmaSq) * min(20.,lens);
 Lsun = Lsun * (sunlen/max(1e-16,lens));
 Lsky = fresnel * meanSkyRadiance(V, N, Tx, Ty, sigmaSq);
 Lsea = (1.0 - fresnel) * (seaColor * Esky / M_PI);

 fragColor.rgb = hdr(Lsun+Lsea+Lsky);
 fragColor.a   = ka*clamp(acolor.x+distance2cam/hFadeDistance,0.0,1.0);

}

