
#include "/lib/cubemap.glsl"


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(rayleighCoeff, R - planetRadius, scaleHeights.x) * chapmanApproximation(R / scaleHeights.x, cosTheta);
    opticalDepth += upOpticalDepth(mieCoeff, R - planetRadius, scaleHeights.y) * chapmanApproximation(R / 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));
}

vec2 getUvFromValues(vec2 values, const vec2 res) {
	return values * (1.0 - 1.0 / res) + (0.5 / res);
}
vec2 getValuesFromUv(vec2 uv, const vec2 res) {
	return (uv - 0.5 / res) / (1.0 - 1.0 / res);
}

// Faster alternative to acos
// Source: https://seblagarde.wordpress.com/2014/12/01/inverse-trigonometric-functions-gpu-optimization-for-amd-gcn-architecture/#more-3316
// Max relative error: 3.9 * 10^-4
// Max absolute error: 6.1 * 10^-4
// Polynomial degree: 2
float acosApprox(float x) {
	const float C0 = 1.57018;
	const float C1 = -0.201877;
	const float C2 = 0.0464619;

	float res = (C2 * abs(x) + C1) * abs(x) + C0; // p(x)
	res *= sqrt(1.0 - abs(x));

	return (x >= 0) ? res : PI - res; // Undo range reduction
}

const vec2 skyCaptureRes = vec2(2560, 1280);

// Sky capture parameterization from https://sebh.github.io/publications/egsr2020.pdf
// Implementation from Photon shader by SixthSurge
vec2 projectSky(vec3 viewDir) {
	float centerDist = planetRadius;

	vec2 projectedDir = normalize(viewDir.xz);

	float horizonCos    = sqrt(pow2(centerDist) - pow2(planetRadius)) / centerDist;
	float horizonAngle  = acosApprox(horizonCos);
	float azimuthAngle  = PI + atan(projectedDir.x, -projectedDir.y);
	float altitudeAngle = horizonAngle - acosApprox(viewDir.y);

	vec2 coord;
	coord.x = azimuthAngle * rcp(TAU);
	coord.y = 0.5 + 0.5 * sign(altitudeAngle) * sqrt(2.0 * RCPPI * abs(altitudeAngle)); // Section 5.3

	return getUvFromValues(coord, skyCaptureRes);
}

vec3 unprojectSky(vec2 coord) {
	coord = getValuesFromUv(coord, skyCaptureRes);

	// Non-linear mapping of altitude angle. See section 5.3 of the paper
	coord.y = (coord.y < 0.5)
		? -pow2(1.0 - 2.0 * coord.y)
		:  pow2(2.0 * coord.y - 1.0);

	float centerDist = planetRadius;
	float horizonCos = sqrt(pow2(centerDist) - pow2(planetRadius)) / centerDist;

	float horizonAngle  = acosApprox(horizonCos) - 0.5 * PI;
	float azimuthAngle  = coord.x * TAU - PI;
	float altitudeAngle = coord.y * 0.5 * PI - horizonAngle;

	float altitudeCos = cos(altitudeAngle);
	float altitudeSin = sin(altitudeAngle);
	float azimuthCos  = cos(azimuthAngle);
	float azimuthSin  = sin(azimuthAngle);

	return normalize(vec3(altitudeCos * azimuthSin, altitudeSin, -altitudeCos * azimuthCos));
}


vec4 samplePrecomputedSky(vec3 worldDir) {


    // ivec2 sampleCoord = ivec2(projectSky(worldDir)*skyCaptureRes);
    // vec4 col = texelFetch(colortex3, sampleCoord, 0);
    vec4 col = texture(colortex3, projectSky(worldDir));
    return vec4(linearToAP1(col.rgb), col.a);
}

vec3 sunDisk(vec3 sunDir, vec3 viewDir, const float sizeMulti) {
    return vec3(1.0)*float(1.0 - dot(sunDir, viewDir) < 0.00017*sizeMulti*SUN_SIZE) * 1200.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))));
}