#version 120
#extension GL_EXT_gpu_shader4 : enable

const int MAX_MARCHING_STEPS    = 64;
const float CLOUD_START_ALT     = 1500.0;
const float CLOUD_END_ALT       = 3000.0;
const float CUTOFF_DISTANCE     = 100000.0;

const float COVERAGE_SCALE      = 0.001;

const int LIGHT_MARCHING_STEPS  = 6;

const float FORWARD_SCATTERING  = 0.8;
const float BACKWARD_SCATTERING = -0.5;
const float SCATTERING_MIX      = 0.5;

const vec4 NO_CLOUD = vec4(0.0, 0.0, 0.0, 1.0);

const vec4 STRATUS_GRADIENT = vec4(0.02f, 0.05f, 0.09f, 0.11f);
const vec4 STRATOCUMULUS_GRADIENT = vec4(0.02f, 0.2f, 0.48f, 0.625f);
const vec4 CUMULUS_GRADIENT = vec4(0.01f, 0.0625f, 0.78f, 1.0f);


uniform float BASE_SCALE;
uniform float EROSION_SCALE;
uniform float CLOUD_EROSION_STRENGTH;
uniform float test_value;

uniform unsigned int osg_FrameNumber;

uniform vec2 fg_ViewportSize;
uniform vec3 fg_CameraPositionCart;
uniform vec3 fg_CameraPositionGeod;
uniform mat4 fg_PrevViewMatrix;
uniform mat4 fg_PrevViewMatrixInverse;
uniform mat4 fg_PrevProjectionMatrix;
uniform mat4 fg_ViewMatrix;
uniform mat4 fg_ViewMatrixInverse;
uniform mat4 fg_ProjectionMatrixInverse;
uniform vec3 fg_LightDirection;


uniform sampler2D prevframe_tex;
uniform sampler2D depth_tex;
uniform sampler3D noise_base_tex;
uniform sampler3D noise_erosion_tex;

uniform float avisibility;
uniform float cloud_self_shading;
uniform float eye_alt;
uniform float ground_scattering;
uniform float hazeLayerAltitude;
uniform float moonlight;
uniform float overcast;
uniform float scattering;
uniform float terminator;
uniform float terrain_alt;
uniform float visibility;
uniform float air_pollution;
uniform float snowlevel;
uniform float snow_thickness_factor;

// constants needed by the light and fog computations ##########################
const float EarthRadius = 5800000.0;
const float terminator_width = 200000.0;

vec4 cloudPositionWorldSpace;


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


float Noise2D(in vec2 coord, in float wavelength);
float Noise3D(in vec3 coord, in float wavelength);

vec3 rayleigh_out_shift(in vec3 color, in float outscatter);
vec3 get_hazeColor(in float lightArg);
vec3 moonlight_perception(in vec3 light);


////////////////////////////////////////////////////////////////////////////////
// ALS functions

float light_func (in float x, in float a, in float b, in float c, in float d, in float e)
{
    if (x > 30.0) {return e;}
    if (x < -15.0) {return 0.0;}
    return e / pow((1.0 + a * exp(-b * (x-c)) ),(1.0/d));
}

void getALSLight(in vec3 pos, out vec3 ambient, out vec3 diffuse)
{
    vec3 shadedFogColor = vec3(0.55, 0.67, 0.88);
    vec3 moonLightColor = vec3 (0.095, 0.095, 0.15) * moonlight;

    vec3 lightHorizon = fg_LightDirection - normalize(pos) *
        dot(normalize(pos), fg_LightDirection);
    float yprime = -dot(pos, lightHorizon);
    float yprime_alt = yprime - sqrt(2.0 * EarthRadius * eye_alt);
    float lightArg = (terminator-yprime_alt)/100000.0;

    float earthShade = 0.6 * (1.0 - smoothstep(-terminator_width + terminator, terminator_width + terminator, yprime_alt)) + 0.4;

    vec3 light_diffuse;
    vec3 light_ambient;
    float intensity;

    light_diffuse.b = light_func(lightArg, 1.330e-05, 0.264, 2.227, 1.08e-05, 1.0);
    light_diffuse.g = light_func(lightArg, 3.931e-06, 0.264, 3.827, 7.93e-06, 1.0);
    light_diffuse.r = light_func(lightArg, 8.305e-06, 0.161, 3.827, 3.04e-05, 1.0);

    light_ambient.r = light_func(lightArg, 0.236, 0.253, 1.073, 0.572, 0.33);
    light_ambient.g = light_ambient.r * 0.4/0.33;
    light_ambient.b = light_ambient.r * 0.5/0.33;

    if (earthShade < 0.5) {
        intensity = length(light_ambient.rgb);
        light_ambient.rgb = intensity *
            normalize(mix(light_ambient.rgb, shadedFogColor,
                          1.0 - smoothstep(0.1,0.8, earthShade) ));
        light_ambient.rgb = light_ambient.rgb + moonLightColor *
            (1.0 - smoothstep(0.4, 0.5, earthShade));

        intensity = length(light_diffuse.rgb);
        light_diffuse.rgb = intensity *
            normalize(mix(light_diffuse.rgb, shadedFogColor,
                          (1.0 - smoothstep(0.5,0.9, cloud_self_shading))));
    }


    ambient = light_ambient;
    diffuse = light_diffuse;
}

////////////////////////////////////////////////////////////////////////////////
// General utility functions

float map(float s, float a1, float a2, float b1, float b2)
{
    return b1+(s-a1)*(b2-b1)/(a2-a1);
}

vec2 cartToGeoc(vec3 cart)
{
    vec2 geoc = vec2(0.0);
    geoc.x = atan(cart.y, cart.x);
    float nxy = sqrt(cart.x*cart.x + cart.y*cart.y);
    geoc.y = atan(cart.z, nxy);
    return geoc;
}

vec3 getFragmentWorldPos(vec2 uv)
{
    vec3 posCS = vec3(uv, texture2D(depth_tex, uv).r) * 2.0 - 1.0;
    vec4 pos = fg_ViewMatrixInverse * fg_ProjectionMatrixInverse *
        vec4(posCS, 1.0);
    return pos.xyz / pos.w;
}

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

// Henyey-Greenstein phase function is used instead of the much more complicated
// Mie phase function to approximate the angular distribution of scattered light
float HenyeyGreenstein(float costheta, float g)
{
    float gg = g * g;
    return (1.0 - gg) * pow(1.0 + gg - 2.0 * g * costheta, -1.5) * 0.25;
}

vec3 getAmbientLight(float altitude)
{
    // TODO: Actually return the cloud ambient light for this altitude
	return mix(
		vec3(0.5f, 0.67f, 0.82f),
		vec3(1.0f, 1.0f, 1.0f),
		altitude);
}

float densityHeightGradient(float heightFrac, float cloudType)
{
    float stratus = 1.0 - clamp(cloudType * 2.0, 0.0, 1.0);
	float stratocumulus = 1.0 - abs(cloudType - 0.5) * 2.0;
	float cumulus = clamp(cloudType - 0.5, 0.0, 1.0) * 2.0;
	vec4 cloudGradient =
        STRATUS_GRADIENT * stratus +
        STRATOCUMULUS_GRADIENT * stratocumulus +
        CUMULUS_GRADIENT * cumulus;
	return smoothstep(cloudGradient.x, cloudGradient.y, heightFrac) -
        smoothstep(cloudGradient.z, cloudGradient.w, heightFrac);
}

// Get the cloud density in a given sky position and altitude
float sampleDensity(vec3 pos, float alt)
{
    // Get the low frequency noise that defines the base shape of the cloud
    float baseCloud = texture3D(noise_base_tex, pos / BASE_SCALE).r;

    // Calculate the cloud coverage value (this could be a CPU generated
    // weather texture instead)
    vec2 pos2d = cartToGeoc(pos);
    float coverage = Noise2D(pos2d, COVERAGE_SCALE) * 0.625 +
        Noise2D(pos2d, COVERAGE_SCALE / 2.0) * 0.250 +
        Noise2D(pos2d, COVERAGE_SCALE / 4.0) * 0.125;
    coverage = smoothstep(0.5, 1.0, coverage);
    //coverage *= weather_coverage;

    float baseCloudWithCoverage = map(baseCloud, 1.0 - coverage, 1.0, 0.0, 1.0);
    baseCloudWithCoverage *= coverage;

    baseCloudWithCoverage *= densityHeightGradient(alt, 1.0);

    // Apply some erosion to add details
    float highFreqNoise = texture3D(noise_erosion_tex, pos / EROSION_SCALE).r;
    float highFreqNoiseModifier = mix(1.0 - highFreqNoise, // whispy
                                      highFreqNoise,
                                      alt);
    float density = map(baseCloudWithCoverage,
                        highFreqNoiseModifier * CLOUD_EROSION_STRENGTH,
                        1.0, 0.0, 1.0);

    return clamp(density, 0.0, 1.0);
}

float sampleDensityAlongLightRay(vec3 p, float alt, vec3 lightDir)
{
    float stepDelta = 10.0;
    float t = stepDelta;
    float shadow = 1.0;
    for(int i = 0; i < LIGHT_MARCHING_STEPS; ++i) {
        vec3 samplePoint = p + lightDir * t;
        float density = sampleDensity(samplePoint, alt);
        shadow *= exp(-density * stepDelta);
        stepDelta *= 1.3;
        t += stepDelta;
    }
    return shadow;
}

vec3 sampleLighting(vec3 p, float alt, vec3 lightDir, float phase)
{
    vec3 ambient, diffuse;
    getALSLight(p, ambient, diffuse);
    ambient = getAmbientLight(alt);
    float shadow = sampleDensityAlongLightRay(p, alt, lightDir);

    // float loddedDensity = texture3DLod(noise_base_tex, p / BASE_SCALE, 3.0).r;
    // float depthProbability = 0.05 +
    //     pow(loddedDensity, map(alt, 0.3, 0.85, 0.5, 2.0));

    return ambient + diffuse * phase * shadow;
}

bool rayAtmosphereIntersection(in vec3 ro, in vec3 rd, in float alt,
                               out float t1, out float t2)
{
    float radius = length(ro) - fg_CameraPositionGeod.z + alt;
    float a = dot(rd, rd) * 2.0;
    float b = dot(rd, ro) * 2.0;
    float c = dot(ro, ro) - radius * radius;
    float discriminant = b * b - 2.0 * a * c;
    if (discriminant < 0.0)
        return false;

    t1 = max(0.0, (-b - sqrt(discriminant)) / a);
    t2 = max(0.0, (-b + sqrt(discriminant)) / a);
    return true;
}

float radicalInverse_VdC(unsigned int bits)
{
	bits = (bits << 16u) | (bits >> 16u);
	bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
	bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
	bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
	bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
	return float(bits) * 2.3283064365386963e-10; // / 0x100000000
}

vec4 rayMarch(vec3 ro, vec3 rd, float depth, vec3 fragWorldPos)
{
    vec3 lightDir = normalize(fg_LightDirection);
    float costheta = dot(rd, lightDir);

    float phase = mix(HenyeyGreenstein(costheta, FORWARD_SCATTERING),
                      HenyeyGreenstein(costheta, BACKWARD_SCATTERING),
                      SCATTERING_MIX);
    phase = clamp(phase, 0.0, 1.0);


    float stepDelta = depth / float(MAX_MARCHING_STEPS);
    float t = stepDelta * radicalInverse_VdC(osg_FrameNumber);

    // The RGB part contains the scattered light color and the alpha value contains
    // the transmittance, both along the ray
    vec4 result = vec4(0.0, 0.0, 0.0, 1.0);

    for (int i = 0; i < MAX_MARCHING_STEPS; ++i) {
        vec3 samplePoint = ro + rd * t;

        // Skip samples behind other objects
        float fragToCamera = distance(fragWorldPos, fg_CameraPositionCart);
        float sampleToCamera = distance(samplePoint, fg_CameraPositionCart);
        if (fragToCamera < sampleToCamera)
            continue;

        float earthRadius = length(fg_CameraPositionCart) - fg_CameraPositionGeod.z;
        float altitude = length(samplePoint) - earthRadius;
        float normalizedAlt = clamp((altitude - CLOUD_START_ALT) /
                                    (CLOUD_END_ALT - CLOUD_START_ALT), 0.0, 1.0);

        float density = sampleDensity(samplePoint, normalizedAlt);

        // Only evaluate lighting for samples inside clouds
        if (density > 0.0) {
            float transmittance = exp(-density * stepDelta);

            vec3 S = sampleLighting(samplePoint, normalizedAlt, lightDir, phase)
                * density;
            vec3 Sint = (S - S * transmittance) / density;

            result.rgb += result.a * Sint;
            result.a *= transmittance;
        }

        // Early exit if we reached maximum density
        if (result.a < 0.01)
            break;

        t += stepDelta;
    }

    return result;
}

vec4 getCloudColor(vec3 rd)
{
    vec3 ro = fg_CameraPositionCart;

    vec3 entry = vec3(0.0);
    float depth = 0.0;
    float t1 = 0.0, t2 = 0.0;
    if (fg_CameraPositionGeod.z < CLOUD_START_ALT) {
        // We are below the clouds
        // Discard fragments below the horizon
        if (dot(normalize(ro), rd) < 0.0)
            return NO_CLOUD;

        rayAtmosphereIntersection(ro, rd, CLOUD_START_ALT, t1, t2);
        entry = ro + rd * t2;
        rayAtmosphereIntersection(ro, rd, CLOUD_END_ALT, t1, t2);
        vec3 exit = ro + rd * t2;
        depth = distance(entry, exit);
    } else if (fg_CameraPositionGeod.z > CLOUD_END_ALT) {
        // We are over the clouds
        if (!rayAtmosphereIntersection(ro, rd, CLOUD_END_ALT, t1, t2))
            return NO_CLOUD;
        entry = ro + rd * t1;
        if (rayAtmosphereIntersection(ro, rd, CLOUD_START_ALT, t1, t2)) {
            vec3 exit = ro + rd * t1;
            depth = distance(entry, exit);
        } else {
            // We didn't intersect with the inner bound, use a fixed cloud depth
            depth = CUTOFF_DISTANCE;
        }
    } else {
        // We are inside the clouds
        entry = ro;
        vec3 exit;
        if (!rayAtmosphereIntersection(ro, rd, CLOUD_END_ALT, t1, t2)) {
            rayAtmosphereIntersection(ro, rd, CLOUD_START_ALT, t1, t2);
            exit = ro + rd * t2;
        } else {
            exit = ro + rd * t1;
        }
        depth = max(distance(entry, exit), CUTOFF_DISTANCE);
    }

    cloudPositionWorldSpace = vec4(entry, 1.0);
    vec3 fragWorldPos = getFragmentWorldPos(gl_TexCoord[0].st);
    return rayMarch(entry, rd, depth, fragWorldPos);
}

const unsigned int bayerMatrix4[16] = unsigned int[](
    0, 8, 2, 10,
    12, 4, 14, 6,
    3, 11, 1, 9,
    15, 7, 13, 5
    );

void main()
{
    vec2 uv = gl_TexCoord[0].st;
    vec2 rayNDS = uv * 2.0 - 1.0;
    vec4 rayCS = vec4(rayNDS, -1.0, 1.0);
    vec4 rayVS = fg_ProjectionMatrixInverse * rayCS;
    rayVS = vec4(rayVS.xy, -1.0, 0.0);

    vec4 rayWS = fg_ViewMatrixInverse * rayVS;
    vec3 rayWSNorm = normalize(rayWS).xyz;

    //vec3 backgroundColor = texture2D(color_tex, uv).xyz;

    vec4 fragColor;

    vec4 cloudColor = getCloudColor(rayWSNorm);

    vec4 rayPrime = fg_PrevProjectionMatrix * fg_PrevViewMatrix *
        cloudPositionWorldSpace;
    rayPrime /= rayPrime.w;
    vec2 prev_uv = rayPrime.xy * 0.5 + 0.5;
    vec4 prevColor = texture2D(prevframe_tex, prev_uv);

    vec3 fragViewSpacePos = (fg_ViewMatrix * vec4(getFragmentWorldPos(gl_TexCoord[0].st), 1.0)).xyz;

    bool isOut = any(greaterThan(abs(prev_uv - 0.5), vec2(0.5)));
    const float alpha = 0.1;
    if (isOut || fragViewSpacePos.z > -500.0) {
        fragColor = cloudColor;
    } else {
        fragColor = cloudColor * alpha + prevColor * (1.0 - alpha);
    }

    // vec2 pixel = gl_TexCoord[0].st * fg_ViewportSize;
    // unsigned int bayerIndex = unsigned int(mod(pixel.x*pixel.y, 16));
    // if (unsigned int(mod(float(osg_FrameNumber), 16)) != bayerMatrix4[bayerIndex]) {
    //     // It's NOT the time to render clouds
    //     vec4 rayPrime = fg_PrevProjectionMatrix * fg_PrevViewMatrix * rayWS;
    //     vec2 prev_uv = rayPrime.xy * 0.5 + 0.5;
    //     vec4 prevColor = texture2D(prevframe_tex, prev_uv);
    //     //vec2 haveInfo = step(vec2(0.0), prev_uv) - step(vec2(1.0), prev_uv);
    //     //fragColor = mix(vec4(0.0, 0.0, 0.0, 1.0), prevColor, haveInfo.x*haveInfo.y);
    //     fragColor = prevColor;
    // } else {
    //     vec4 cloudColor = getCloudColor(rayWSNorm);
    //     fragColor = cloudColor;
    // }

    // TEMPORAL
    //fragColor = mix(backgroundColor, prevColor, haveInfo.x*haveInfo.y);
    // CURRENT
    //fragColor = cloudColor.rgb + backgroundColor.rgb * cloudColor.a;

    gl_FragColor = fragColor;
}