#type compute
#version 450 core

const float PI = 3.141592;

layout(binding = 0, rgba32f) restrict writeonly uniform imageCube o_CubeMap;

uniform vec3 u_TurbidityAzimuthInclination;

#define PI 3.14159265359

vec3 GetCubeMapTexCoord()
{
    vec2 st = gl_GlobalInvocationID.xy / vec2(imageSize(o_CubeMap));
    vec2 uv = 2.0 * vec2(st.x, 1.0 - st.y) - vec2(1.0);

    vec3 ret;
    if (gl_GlobalInvocationID.z == 0)      ret = vec3(  1.0, uv.y, -uv.x);
    else if (gl_GlobalInvocationID.z == 1) ret = vec3( -1.0, uv.y,  uv.x);
    else if (gl_GlobalInvocationID.z == 2) ret = vec3( uv.x,  1.0, -uv.y);
    else if (gl_GlobalInvocationID.z == 3) ret = vec3( uv.x, -1.0,  uv.y);
    else if (gl_GlobalInvocationID.z == 4) ret = vec3( uv.x, uv.y,   1.0);
    else if (gl_GlobalInvocationID.z == 5) ret = vec3(-uv.x, uv.y,  -1.0);
    return normalize(ret);
}


float saturatedDot( in vec3 a, in vec3 b )
{
    return clamp(dot(a, b), 0.0, 1.0);
}

vec3 YxyToXYZ( in vec3 Yxy )
{
    float Y = Yxy.r;
    float x = Yxy.g;
    float y = Yxy.b;

    float X = x * ( Y / y );
    float Z = ( 1.0 - x - y ) * ( Y / y );

    return vec3(X,Y,Z);
}

vec3 XYZToRGB( in vec3 XYZ )
{
    // CIE/E
    mat3 M = mat3
    (
    2.3706743, -0.9000405, -0.4706338,
    -0.5138850,  1.4253036,  0.0885814,
    0.0052982, -0.0146949,  1.0093968
    );

    return XYZ * M;
}


vec3 YxyToRGB( in vec3 Yxy )
{
    vec3 XYZ = YxyToXYZ( Yxy );
    vec3 RGB = XYZToRGB( XYZ );
    return RGB;
}

void calculatePerezDistribution( in float t, out vec3 A, out vec3 B, out vec3 C, out vec3 D, out vec3 E )
{
    A = vec3(  0.1787 * t - 1.4630, -0.0193 * t - 0.2592, -0.0167 * t - 0.2608 );
    B = vec3( -0.3554 * t + 0.4275, -0.0665 * t + 0.0008, -0.0950 * t + 0.0092 );
    C = vec3( -0.0227 * t + 5.3251, -0.0004 * t + 0.2125, -0.0079 * t + 0.2102 );
    D = vec3(  0.1206 * t - 2.5771, -0.0641 * t - 0.8989, -0.0441 * t - 1.6537 );
    E = vec3( -0.0670 * t + 0.3703, -0.0033 * t + 0.0452, -0.0109 * t + 0.0529 );
}

vec3 calculateZenithLuminanceYxy( in float t, in float thetaS )
{
    float chi  	 	= ( 4.0 / 9.0 - t / 120.0 ) * ( PI - 2.0 * thetaS );
    float Yz   	 	= ( 4.0453 * t - 4.9710 ) * tan( chi ) - 0.2155 * t + 2.4192;

    float theta2 	= thetaS * thetaS;
    float theta3 	= theta2 * thetaS;
    float T 	 	= t;
    float T2 	 	= t * t;

    float xz =
    ( 0.00165 * theta3 - 0.00375 * theta2 + 0.00209 * thetaS + 0.0)     * T2 +
    (-0.02903 * theta3 + 0.06377 * theta2 - 0.03202 * thetaS + 0.00394) * T +
    ( 0.11693 * theta3 - 0.21196 * theta2 + 0.06052 * thetaS + 0.25886);

    float yz =
    ( 0.00275 * theta3 - 0.00610 * theta2 + 0.00317 * thetaS + 0.0)     * T2 +
    (-0.04214 * theta3 + 0.08970 * theta2 - 0.04153 * thetaS + 0.00516) * T +
    ( 0.15346 * theta3 - 0.26756 * theta2 + 0.06670 * thetaS + 0.26688);

    return vec3( Yz, xz, yz );
}

vec3 calculatePerezLuminanceYxy( in float theta, in float gamma, in vec3 A, in vec3 B, in vec3 C, in vec3 D, in vec3 E )
{
    float cosTheta = max(cos(theta), 1e-6);
    return ( 1.0 + A * exp( B / cos( theta ) ) ) * ( 1.0 + C * exp( D * gamma ) + E * cos( gamma ) * cos( gamma ) );
}

vec3 calculateSkyLuminanceRGB( in vec3 s, in vec3 e, in float t )
{
    vec3 A, B, C, D, E;
    calculatePerezDistribution( t, A, B, C, D, E );

    float thetaS = acos(clamp(dot(s, vec3(0,1,0)), 0.0, 1.0));
    float thetaE = acos(clamp(dot(e, vec3(0,1,0)), 0.0, 1.0));
    float gammaE = acos( saturatedDot( s, e )		   );

    vec3 Yz = calculateZenithLuminanceYxy( t, thetaS );

    vec3 fThetaGamma = calculatePerezLuminanceYxy( thetaE, gammaE, A, B, C, D, E );
    vec3 fZeroThetaS = calculatePerezLuminanceYxy( 0.0,    thetaS, A, B, C, D, E );

    vec3 Yp = Yz * ( fThetaGamma / fZeroThetaS );

    return YxyToRGB( Yp );
}

layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;
void main()
{
    vec3 cubeTC = GetCubeMapTexCoord();

    float turbidity     = u_TurbidityAzimuthInclination.x;
    float azimuth       = u_TurbidityAzimuthInclination.y;;
    float inclination   = u_TurbidityAzimuthInclination.z;
    vec3 sunDir     	= normalize( vec3( sin(inclination) * cos(azimuth), cos(inclination), sin(inclination) * sin(azimuth) ) );
    vec3 viewDir  		= cubeTC;

    const float SUN_ANGULAR_RADIUS = 0.03465;
    const float SUN_INTENSITY = 100.0;
    vec3 skyLuminance;

    if (viewDir.y < 0.0) {
        skyLuminance = vec3(0.02);
    } else {
        skyLuminance = calculateSkyLuminanceRGB(sunDir, viewDir, turbidity);
    }

    float cosAngle = dot(viewDir, sunDir);
    float angle = acos(cosAngle);
    if (angle < SUN_ANGULAR_RADIUS) {
        skyLuminance = vec3(SUN_INTENSITY);
    } else {
        float haloWidth = 0.1;
        if (angle < SUN_ANGULAR_RADIUS + haloWidth) {
            float t = (angle - SUN_ANGULAR_RADIUS) / haloWidth;
            float haloFactor = 1.0 - smoothstep(0.0, 1.0, t);
            skyLuminance += vec3(SUN_INTENSITY * 0.1 * haloFactor);
        }
    }

    vec4 color = vec4(skyLuminance * 0.05, 1.0);
    imageStore(o_CubeMap, ivec3(gl_GlobalInvocationID), color);
}