
#include "/lib/cubemap.glsl"
#include "/lib/atmos/constants.glsl"

const float RAYLEIGH_SCALE = 8.5e3;
const float MIE_SCALE      = 1.2e3;
const float OZONE_HEIGHT   = 35e3;

vec2 scaleHeights = vec2(RAYLEIGH_SCALE, MIE_SCALE);

vec3 approx_r_ho(vec3 rd) {
    // Approximate densities 
    // x - rayleigh
    // y - mie
    // z - ozone
    float y = rd.y * atmosphere_planetRadius;
    vec4 atmosRadius = vec4(RAYLEIGH_SCALE, MIE_SCALE, OZONE_HEIGHT, OZONE_HEIGHT + RAYLEIGH_SCALE) + atmosphere_planetRadius;
    vec4 z = sqrt(atmosRadius*atmosRadius + y*y - atmosphere_planetRadius*atmosphere_planetRadius) - y;

    return vec3(z.xy, (z.w - z.z) * 18.0);
}

vec3 calcSunTransmitted(vec3 dirToSun, vec3 L_0) {
     vec3 sun_r_ho = approx_r_ho(dirToSun);
     vec3 airdensity_InScatter = atmosphere_coefficientRayleigh*sun_r_ho.x + 
                                 atmosphere_coefficientMie*sun_r_ho.y + 
                                 atmosphere_coefficientOzone*sun_r_ho.z; // The probability of a scatter event 
     vec3 sunExtinction = expf(-airdensity_InScatter);
     return L_0 * sunExtinction;
}

float chapmanApproximation(float ratio, float cosTheta) {
    float sqrtRatio = sqrt(ratio);
    float additive = 1.0 / (1.0 + sqrtRatio);

    if (cosTheta >= 0.0) {
        return sqrtRatio / (1.0 + sqrtRatio * cosTheta) + additive;
    } else {
        float sinTheta = sqrt(clamp01(1.0 - cosTheta * cosTheta));
        return 2.0 * sqrt(sinTheta * ratio) * expf(ratio * (1.0 - sinTheta)) - sqrtRatio / (1.0 - sqrtRatio * cosTheta) + additive;
    }
}

vec3 upOpticalDepth(vec3 attenuationCoefficient, float height, float scaleHeight) {
    return attenuationCoefficient * scaleHeight * expf(-height / scaleHeight);
}

vec3 atmosphereTransmittanceApproximation(vec3 rayPosition, vec3 rayDirection) {
    float R = length(rayPosition);
    float cosTheta = dot(rayPosition, rayDirection) / R;
    vec3 opticalDepth = vec3(0.0);
    opticalDepth += upOpticalDepth(atmosphere_coefficientRayleigh, R - atmosphere_planetRadius, atmosphere_scaleHeights.x) * chapmanApproximation(R / atmosphere_scaleHeights.x, cosTheta);
    opticalDepth += upOpticalDepth(atmosphere_coefficientMie, R - atmosphere_planetRadius, atmosphere_scaleHeights.y) * chapmanApproximation(R / atmosphere_scaleHeights.y, cosTheta);
    return clamp01(expf(-opticalDepth));
}

const float[3] cubemapKernel = float[3](0.19858494730562265, 0.6028301053887546, 0.19858494730562265);

const ivec4[4] cubemapXOffsets = ivec4[4](
//class:   0  1  2  3
	ivec4( 0, 1, 1, 0),
	ivec4( 2, 3, 3, 2),
	ivec4( 0, 1, 1, 0),
	ivec4( 2, 3, 3, 2)
);
const ivec4[4] cubemapYoffsets = ivec4[4](
//class:   0  1  2  3
	ivec4( 0, 1, 0, 1),
	ivec4( 0, 1, 0, 1),
	ivec4( 2, 3, 2, 3),
	ivec4( 2, 3, 2, 3)
);

// vec4 upscaleCubemap(sampler2D tex, ivec2 coord, ivec2 fullResCoord) {
//     ivec2 coordInTile = fullResCoord % 2;
// 	int set = coordInTile.x + 2 * coordInTile.y;

//     vec4 sum = vec4(0.0);
//     float weightSum = 1e-5;
//     int index = getFaceIndex(coord);
//     for(int i = 0; i < 4; i++) {
//         int x = cubemapXOffsets[i][set]-2;
// 		int y = cubemapYoffsets[i][set]-2;

// 		ivec2 offsetCoord = coord + ivec2(x,y);

//         float weight = float(getFaceIndex(offsetCoord) == index && 
//                 clamp(offsetCoord, ivec2(0), ivec2(CUBEMAP_RES*6, CUBEMAP_RES)-1) == offsetCoord);
//         sum += texelFetch(tex, offsetCoord, 0).rgba * weight;
//         weightSum += weight;

//     }
//     return sum /= weightSum;

// }

vec4 upscaleCubemap(sampler2D t, vec2 x, vec2 textureSize) {
    x *= vec2(textureSize);

    vec2 p = floor(x);
    vec2 f = fract(x);

    // vec4 a = texture2D(t, (p                 ) / vec2(textureSize));
    // vec4 b = texture2D(t, (p + vec2(1.0, 0.0)) / vec2(textureSize));
    // vec4 c = texture2D(t, (p + vec2(0.0, 1.0)) / vec2(textureSize));
    // vec4 d = texture2D(t, (p + vec2(1.0, 1.0)) / vec2(textureSize));
    vec4 a = texelFetch(t, ivec2(p                 ), 0);
    vec4 b = texelFetch(t, ivec2(p + vec2(1.0, 0.0)), 0);
    vec4 c = texelFetch(t, ivec2(p + vec2(0.0, 1.0)), 0);
    vec4 d = texelFetch(t, ivec2(p + vec2(1.0, 1.0)), 0);

    int index = getFaceIndex(ivec2(p));
    float bValid = float(getFaceIndex(ivec2(p + vec2(1.0, 0.0))) == index && 
                clamp(ivec2(p + vec2(1.0, 0.0)), ivec2(0), ivec2(CUBEMAP_RES*6, CUBEMAP_RES)-1) == ivec2(p + vec2(1.0, 0.0)));
    float cValid = float(getFaceIndex(ivec2(p + vec2(0.0, 1.0))) == index && 
                clamp(ivec2(p + vec2(0.0, 1.0)), ivec2(0), ivec2(CUBEMAP_RES*6, CUBEMAP_RES)-1) == ivec2(p + vec2(0.0, 1.0)));
    float dValid = float(getFaceIndex(ivec2(p + vec2(1.0, 1.0))) == index && 
                clamp(ivec2(p + vec2(1.0, 1.0)), ivec2(0), ivec2(CUBEMAP_RES*6, CUBEMAP_RES)-1) == ivec2(p + vec2(1.0, 1.0)));
    
    // b = bValid ? b : a;
    // c = cValid ? c : a;
    // d = dValid ? d : a;
    return mix(mix(a, b, f.x*bValid), mix(c, d, f.x*dValid), f.y*clamp01(cValid + dValid));
}
vec4 samplePrecomputedSkyFiltered(vec3 worldDir, ivec2 fullResCoord, vec3 noise) {

	// float phi = acos(worldDir.y);
    // float theta = atan(worldDir.z, worldDir.x);
    // if(theta < 0.0) {
    //     theta = 2.0*PI + theta;
    // }

	// vec2 sampleUV = vec2(theta / (2.0*PI), phi / PI);
    worldDir = normalize(worldDir + noise*0.0015);
    int index;
    vec2 sampleUV = convert_xyz_to_cube_uv(worldDir*vec3(1,-1,1), index);
        sampleUV.x = (sampleUV.x + float(index))/6.0;

    ivec2 sampleCoord = ivec2(convert_xyz_to_cube_uv(worldDir*vec3(1,-1,1), index)*CUBEMAP_RES);
        sampleCoord.x += index*CUBEMAP_RES;
    

    // vec4 col = upscaleCubemap(colortex3, sampleUV, vec2(CUBEMAP_RES*6.0, CUBEMAP_RES));

    // vec4 col = upscaleCubemap(colortex3, sampleUV, vec2(CUBEMAP_RES*6.0, CUBEMAP_RES));

    vec4 col = texelFetch(colortex3, sampleCoord, 0);

    vec4 sum = vec4(0.0);

    return vec4(linearToAP1(col.rgb), col.a);
    vec2 fullres = floor(vec2(CUBEMAP_RES*6.0, CUBEMAP_RES));
    // return linearToAP1(upscaleCubemap(colortex3, sampleUV, vec4(rcp(fullres), fullres), 0.75).rgb);
    // return linearToAP1(textureLod(colortex3, sampleUV, 0.0).rgb);
}

vec4 samplePrecomputedSky(vec3 worldDir) {

	// float phi = acos(worldDir.y);
    // float theta = atan(worldDir.z, worldDir.x);
    // if(theta < 0.0) {
    //     theta = 2.0*PI + theta;
    // }

	// vec2 sampleUV = vec2(theta / (2.0*PI), phi / PI);
    int index;
    // vec2 sampleUV = convert_xyz_to_cube_uv(worldDir*vec3(1,-1,1), index);
    //     sampleUV.x = (sampleUV.x + float(index))/6.0;

    ivec2 sampleCoord = ivec2(convert_xyz_to_cube_uv(worldDir*vec3(1,-1,1), index)*CUBEMAP_RES);
        sampleCoord.x += index*CUBEMAP_RES;
    vec4 col = texelFetch(colortex3, sampleCoord, 0);
    return vec4(linearToAP1(col.rgb), col.a);
    // return textureLod(colortex3, sampleUV, 0.0).rgb;
}

vec3 sunDisk(vec3 sunDir, vec3 viewDir, const float sizeMulti) {
    return vec3(1.0)*float(1.0 - dot(sunDir, viewDir) < 0.00017*sizeMulti*SUN_SIZE) * 2000.0;
}
vec3 moonDisk(vec3 sunDir, vec3 viewDir, const float sizeMulti) {
    return vec3(1.0)*float(1.0 - dot(sunDir, viewDir) < 0.00017*sizeMulti*SUN_SIZE) * 2000.0;
}

vec3 calcNetherSky(vec3 rdW) {
    return 0.2*mix(netherSkyCol*0.4, netherSkyCol*1.3, abs(dot(rdW, vec3(0,1,0))));
}