
######


variable.float.lcR1 =  lightRadiance * exp(-5.22 * lstv) + skycolnight.r
variable.float.lcR0 =  lightRadiance * exp(-7.09 * lstv) + skycolnight.r
variable.float.lcR4 =  lightRadiance * exp(-5.32 * lstv) + skycolnight.r

variable.float.lcG1 =  lightRadiance * exp(-9.715 * lstv) + skycolnight.g
variable.float.lcG0 =  lightRadiance * exp(-13.21 * lstv) + skycolnight.g
variable.float.lcG4 =  lightRadiance * exp(-9.715 * lstv) + skycolnight.g

variable.float.lcB1 =  lightRadiance * exp(-18.59 * lstv) + skycolnight.b
variable.float.lcB0 =  lightRadiance * exp(-25.28 * lstv) + skycolnight.b
variable.float.lcB4 =  lightRadiance * exp(-18.96 * lstv) + skycolnight.b


####
variable.float.R01X = 0.070552 * lcR1
variable.float.R01Y = 0.157153 * lcG1
variable.float.R01Z = 0.341242 * lcB1

variable.float.R00X = 0.784334 * lcR0
variable.float.R00Y = 0.972254 * lcG0
variable.float.R00Z = 0.999842 * lcB0

variable.float.R04X = 0.151315 * lcR4
variable.float.R04Y = 0.318453 * lcG4
variable.float.R04Z = 0.607804 * lcB4



####

variable.float.M01X = 0.025206 * lcR1
variable.float.M01Y = 0.025206 * lcG1
variable.float.M01Z = 0.025200 * lcB1

variable.float.M00X = 0.740012 * lcR0
variable.float.M00Y = 0.740012 * lcG0
variable.float.M00Z = 0.739924 * lcB0

variable.float.M04X = 0.058328 * lcR4
variable.float.M04Y = 0.058328 * lcG4
variable.float.M04Z = 0.058314 * lcB4



#############
variable.float.costh_up_r = 0.12 * (1.0 + lightPos.y * lightPos.y)
variable.float.costh_up_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * lightPos.y) * 1.72 - 1.7 * lightPos.y)

variable.float.costh_left_r = 0.12 * (1.0 + lightPos.x * lightPos.x)
variable.float.costh_left_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * lightPos.x) * 1.72 - 1.7 * lightPos.x)

variable.float.costh_left_50_up_r = 0.12 * (1.0 + (-0.8944 * lightPos.x + 0.4472 * lightPos.y) * (-0.8944 * lightPos.x + 0.4472 * lightPos.y))
variable.float.costh_left_50_up_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * (-0.8944 * lightPos.x + 0.4472 * lightPos.y)) * 1.72 - 1.7 * (-0.8944 * lightPos.x + 0.4472 * lightPos.y))

variable.float.costh_forward_r = 0.12 * (1.0 + lightPos.z * lightPos.z)
variable.float.costh_forward_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * lightPos.z) * 1.72 - 1.7 * lightPos.z)

variable.float.costh_forward_50_up_r = 0.12 * (1.0 + (0.4472 * lightPos.y + 0.8944 * lightPos.z) * (0.4472 * lightPos.y + 0.8944 * lightPos.z))
variable.float.costh_forward_50_up_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * (0.4472 * lightPos.y + 0.8944 * lightPos.z)) * 1.72 - 1.7 * (0.4472 * lightPos.y + 0.8944 * lightPos.z))

variable.float.costh_right_r = 0.12 * (1.0 + lightPos.x * lightPos.x)
variable.float.costh_right_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * lightPos.x) * 1.72 - 1.7 * lightPos.x)

variable.float.costh_right_50_up_r = 0.12 * (1.0 + (0.8944 * lightPos.x + 0.4472 * lightPos.y) * (0.8944 * lightPos.x + 0.4472 * lightPos.y))
variable.float.costh_right_50_up_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * (0.8944 * lightPos.x + 0.4472 * lightPos.y)) * 1.72 - 1.7 * (0.8944 * lightPos.x + 0.4472 * lightPos.y))

variable.float.costh_backward_r = 0.12 * (1.0 + lightPos.z * lightPos.z)
variable.float.costh_backward_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * lightPos.z) * 1.72 - 1.7 * lightPos.z)

variable.float.costh_backward_50_up_r = 0.12 * (1.0 + (-0.8944 * lightPos.z + 0.4472 * lightPos.y) * (-0.8944 * lightPos.z + 0.4472 * lightPos.y))
variable.float.costh_backward_50_up_m = 0.135 / (12.56 * sqrt(1.72 - 1.7 * (-0.8944 * lightPos.z + 0.4472 * lightPos.y)) * 1.72 - 1.7 * (-0.8944 * lightPos.z + 0.4472 * lightPos.y))

###########
uniform.vec3.kRCosthLeft50UpR = vec3(R04X * costh_left_50_up_r, R04Y * costh_left_50_up_r, R04Z * costh_left_50_up_r)
uniform.vec3.kRCosthRight50UpR = vec3(R04X * costh_right_50_up_r, R04Y * costh_right_50_up_r, R04Z * costh_right_50_up_r)
uniform.vec3.kRCosthForward50UpR = vec3(R04X * costh_forward_50_up_r, R04Y * costh_forward_50_up_r, R04Z * costh_forward_50_up_r)
uniform.vec3.kRCosthBackward50UpR = vec3(R04X * costh_backward_50_up_r, R04Y * costh_backward_50_up_r, R04Z * costh_backward_50_up_r)
uniform.vec3.kRCosthUpR = vec3(R01X * costh_up_r, R01Y * costh_up_r, R01Z * costh_up_r)
uniform.vec3.kRCosthForwardR = vec3(R04X * costh_forward_r, R04Y * costh_forward_r, R04Z * costh_forward_r)
uniform.vec3.kRCosthBackwardR = vec3(R04X * costh_backward_r, R04Y * costh_backward_r, R04Z * costh_backward_r)
uniform.vec3.kRCosthRightR = vec3(R00X * costh_right_r, R00Y * costh_right_r, R00Z * costh_right_r)
uniform.vec3.kRCosthLeftR = vec3(R00X * costh_left_r, R00Y * costh_left_r, R00Z * costh_left_r)

uniform.vec3.kMCosthLeft50UpM = vec3(M04X * costh_left_50_up_m, M04Y * costh_left_50_up_m, M04Z * costh_left_50_up_m)
uniform.vec3.kMCosthRight50UpM = vec3(M04X * costh_right_50_up_m, M04Y * costh_right_50_up_m, M04Z * costh_right_50_up_m)
uniform.vec3.kMCosthForward50UpM = vec3(M01X * costh_forward_50_up_m, M01Y * costh_forward_50_up_m, M01Z * costh_forward_50_up_m)
uniform.vec3.kMCosthBackward50UpM = vec3(M01X * costh_backward_50_up_m, M01Y * costh_backward_50_up_m, M01Z * costh_backward_50_up_m)
uniform.vec3.kMCosthUpM = vec3(M01X * costh_up_m, M01Y * costh_up_m, M01Z * costh_up_m)
uniform.vec3.kMCosthForwardM = vec3(M01X * costh_forward_m, M01Y * costh_forward_m, M01Z * costh_forward_m)
uniform.vec3.kMCosthBackwardM = vec3(M01X * costh_backward_m, M01Y * costh_backward_m, M01Z * costh_backward_m)
uniform.vec3.kRCosthRightM = vec3(M00X * costh_right_m, M00Y * costh_right_m, M00Z * costh_right_m)
uniform.vec3.kRCosthLeftM = vec3(M00X * costh_left_m, M00Y * costh_left_m, M00Z * costh_left_m)

variable.vec3.mult4 =  vec3(0.4,0.4,0.4)
variable.vec3.mult15 = vec3(0.15,0.15,0.15)

##########
uniform.vec3.R_UP = (kRCosthUpR * mult4 + kRCosthLeft50UpR * mult15 + kRCosthRight50UpR * mult15 + kRCosthForward50UpR * mult15 + kRCosthBackward50UpR * mult15)
uniform.vec3.M_UP = (kMCosthUpM * mult4 + kMCosthLeft50UpM * mult15 + kMCosthRight50UpM * mult15 + kMCosthForward50UpM * mult15 + kMCosthBackward50UpM * mult15)
uniform.vec3.RM_UP = (R_UP + M_UP) * vec3(sunpow,sunpow,sunpow)



vec4 textureAniso(sampler2D T, vec2 p) {
    vec2 R = vec2(textureSize(T, 0));

    mat2 J = inverse(mat2(dFdx(p), dFdy(p)));       // dFdxy: pixel footprint in texture space
    J = transpose(J) * J;                            // quadratic form
    float d = determinant(J), t = J[0][0] + J[1][1], // find ellipse: eigenvalues, max eigenvector
    D = sqrt(abs(t * t - 4. * d)),                 // abs() fix a bug: in weird view angles 0 can be slightly negative
    V = (t - D) / 2., v = (t + D) / 2.,                     // eigenvalues. ( ATTENTION: not sorted )
    M = 1. / sqrt(V), m = 1. / sqrt(v), l = log2(m * R.y); // = 1./radii^2
    vec2 A = M * normalize(vec2(-J[0][1], J[0][0] - V)); // max eigenvector = main axis
    vec4 O = vec4(0);
    for(float i = -7.5; i < 8.; i++)                       // sample x16 along main axis at LOD min-radius
        O += textureLod(T, p + (i / 16.) * A, l);
    return O / 16.;
}

uint PackVec2ToR16UI(vec2 data) {
    uint x = uint(data.x * 65535.0) & 0xFFFF;  // First component
    uint y = uint(data.y * 65535.0) & 0xFFFF;  // Second component

    // No need for int-to-uint casting here as all values are already uint
    return (x << 0) | (y << 16);  // Packing x, y into 32 bits
}

vec2 UnpackVec2FromR16UI(uint packedData) {
    float x = float(packedData & 0xFFFF) / 65535.0;          // Extract and normalize
    float y = float((packedData >> 16) & 0xFFFF) / 65535.0;  // Extract and normalize

    return vec2(x, y);
}
uvec3 PackVec3ToR16UI(vec3 data) {
    uvec3 packed;
    packed.r = uint(data.r * 65535.0);  // Pack first component
    packed.g = uint(data.g * 65535.0);  // Pack second component
    packed.b = uint(data.b * 65535.0);  // Pack third component

    return packed;
}
vec3 UnpackVec3FromR16UI(uvec3 packedData) {
    vec3 unpacked;
    unpacked.r = float(packedData.r) / 65535.0;  // Unpack and normalize
    unpacked.g = float(packedData.g) / 65535.0;  // Unpack and normalize
    unpacked.b = float(packedData.b) / 65535.0;  // Unpack and normalize

    return unpacked;
}



/////



//// PENUMBRA SHADOWS


        const int iterations = SHADOW_SAMPLES;
        const vec3 filterDepth = vec3(2.0, 40, 32.0); // x=min radius, y=max radius, z=max depth
        const float filterDepthRange = filterDepth.y - filterDepth.x;
        ProjShadowPos = ProjShadowPos * vec3(0.5, 0.5, 0.0833333333333333) + 0.5;


        shading = 0.0;
        avgDepth = 0;
        pdepth = 0;

        float blockerCount = 0.0;
        float avgBlockerDepth0 = 0.0;




        #ifdef SSS_ALT
        diffuseSun = mix(diffuseSun, max(diffuseSun, 0.5), sssAmount > 0.5);
        float scaleFac = 1.0 / (mix(distortFac, 0.5 + noise * 10, float(sssAmount > 0.0)) * 6000.0);
        #else
        float scaleFac = 1.0 / (distortFac * 6000.0);
        #endif

        ///
        float diffthresh = ((sqrt(1.0 - diffuseSun * diffuseSun) / diffuseSun + 0.7) * scaleFac) * max(shadowResxDist, 0.95);

        ///
        float d01k = (d0 + ((length(fragpos) / far) * 0.15)) * kdistort * distortFac;
        float rdMul0 = d01k / shadowMapResolution;
        float RadiusMax = rdMul0 * filterDepth.y;
        float diffthreshM = diffthresh * filterDepth.y * d01k;
        ///



        int counter = (frameCounter % 40000) * int(iterations);

        for(int i = 0; i < iterations; i++) {
            vec2 offset = tapLocationBlueNoise(i, iterations, 1.0, (R2_samples3(counter + i) + noise2.rg)); // Directly apply RadiusMax

            vec2 shadowCoord = ProjShadowPos.xy + offset * RadiusMax; // Calculate shadow coordinate
            float d = texelFetch(shadowtex0, ivec2(shadowCoord.xy * shadowMapResolution), 0).x; // Sample shadow texture
            float b = clamp(1.0 - (d - ProjShadowPos.z) / diffthreshM, 0.0, 1.0);
            blockerCount += b;
            avgBlockerDepth0 += d * b;

            float blockerCountInv = 1.0 / blockerCount; // Compute inverse once if blockerCount > 0 to avoid division multiple times.
            float isBlocker = step(0.9, blockerCount); // Compute step function once to use in subsequent calculations.

            float avgBlockerDepth = mix(ProjShadowPos.z, avgBlockerDepth0 * blockerCountInv, isBlocker);

            float t0 = clamp((ProjShadowPos.z - avgBlockerDepth) / filterDepth.z, 0.0, 1.0) * 1500.0;
            pdepth = (t0 * filterDepthRange + filterDepth.x) * isBlocker * rdMul0;
            avgDepth += max(ProjShadowPos.z - d, 0.0); // Update average depth¨

         

            float shadowMapFactor = shadowMapResolution * (pdepth / iterations);

            float weight = 1.0 + noise * shadowMapFactor;

            vec3 shadowCoordWeighted = ProjShadowPos + vec3(offset * pdepth, -diffthresh * weight);

            shading += texture(shadowtex1, shadowCoordWeighted); // Accumulate shadow intensity
            avgDepth += max(ProjShadowPos.z - shading, 0.0); // Update average depth¨
            shading = mix(1.0, shading, nightblendShadows);

        }
        shading /= iterations;

        avgDepth /= iterations;
        avgDepth = sqrt(avgDepth * 100 * SSS_MULTI);
        avgDepth = mix(1, avgDepth, nightblendShadows);

        if(sssAmount > 0) {
            float extinction = 1.0 - luma(albedo) * 0.85;
            float angleFactor = dot(normalize(p3), lightPos);

            float extinctionFactor = exp(-avgDepth * 11.0 * extinction);
            float sssBase = extinctionFactor + 3.0 * pow(extinctionFactor, 1.0 / 3.0);

            float phaseFactor = phaseg(angleFactor, 0.85);
            float scattering = clamp((0.7 + 0.3 * PI * phaseFactor) * 1.5 / 4.0 * sssAmount, 0.0, 1.0);
            diffuseSun *= 1.0 - sssAmount;

            SSS = clamp(sssBase * scattering, 0, 1);

        }
    
//// FANCY FILTERING

        shading = 0.0;
        float pdepth = 0.001;
        const int iterations = SHADOW_SAMPLES;
        const vec3 filterDepth = vec3(2.0, 40, 32.0); // x=min radius, y=max radius, z=max depth
        const float filterDepthRange = filterDepth.y - filterDepth.x;



        ProjShadowPos = ProjShadowPos * vec3(0.5, 0.5, 0.0833333333333333) + 0.5;

        diffuseSun = mix(diffuseSun, 0.5, sssAmount);




        float d01 = (d0 + ((length(fragpos) / far) * 0.15));
        float d01k = d01 * kdistort * distortFac;
        float rdMul0 = d01k / shadowMapResolution;


        float scaleFac = 1.0 / (mix(distortFac, 0.5 + noise * 10, float(sssAmount > 0.0)) * 6000.0);
        float diffthresh = ((sqrt(1.0 - diffuseSun * diffuseSun) / diffuseSun + 0.7) * scaleFac) * max(shadowResxDist, 0.95);


        int counter = (frameCounter % 40000) * int(iterations);
        float shadowMapFactor = shadowMapResolution * (pdepth / iterations);
        float weight = 1.0 + noise * shadowMapFactor;


        for(int i = 0; i < iterations; i++) {
            vec2 offset = tapLocationBlueNoise(i, iterations, 1.0, (R2_samples3(counter + i) + noise2.rg)) * pdepth;
            vec3 shadowCoordWeighted = ProjShadowPos + vec3(offset, -diffthresh * weight);

            shading += texture(shadowtex1, shadowCoordWeighted);

        }

        shading /= iterations;
        shading = mix(1, shading, nightblendShadows);


    

//// BASIC FILTERING
        shading = 0.0;
        ProjShadowPos = ProjShadowPos * vec3(0.5, 0.5, 0.0833333333333333) + 0.5;
        diffuseSun = mix(diffuseSun, 0.5, sssAmount);
        float scaleFac = 1.0 / (mix(distortFac, 0.5 + noise * 10, float(sssAmount > 0.0)) * 6000.0);

        float term1 = (sqrt(1.0 - diffuseSun * diffuseSun) / diffuseSun) * scaleFac;

        float diffthresh = term1 * max(shadowResxDist, 0.95);
        mat2 noiseM = mat2(cos(prec), -sin(prec), sin(prec), cos(prec)) * 0.001;
        float weight = 1.0 + noise * shadowMapResolution * 0.00016;
        float dweight = -diffthresh * weight;

        for(int i = 0; i < 6; i++) {
            vec2 offset = noiseM * shadowOffsets[i];

            vec3 shadowCoordWeighted = ProjShadowPos + vec3(offset, dweight);

            shading += texture(shadow, shadowCoordWeighted);

        }
        shading *= 0.1666;



//// NO FILTERING
        shading = 0.0;
        ProjShadowPos = ProjShadowPos * vec3(0.5, 0.5, 0.0833333333333333) + 0.5;
        diffuseSun = mix(diffuseSun, 0.5, sssAmount);
        float scaleFac = 1.0 / (mix(distortFac, 0.5 + noise * 10, float(sssAmount > 0.0)) * 6000.0);

        float term1 = (sqrt(1.0 - diffuseSun * diffuseSun) / diffuseSun) * scaleFac;

        float diffthresh = term1 * max(shadowResxDist, 0.95);
        ProjShadowPos.z -= diffthresh;
        shading = texture(shadow, ProjShadowPos);

        

        /////


vec3 mapUVToRoundedSquare(vec2 uv) {

    // Normalize UV coordinates to [-1, 1]
    vec2 centeredUV = uv * 2.0 - 1.0;
    float r = max(abs(centeredUV.x), abs(centeredUV.y));
    r = r * 2.0;
    float distSquared = dot(centeredUV, centeredUV);

    // Mix between circular and square coordinates
    vec2 squareUV = sign(centeredUV) * sqrt(abs(centeredUV));
    vec2 mixedUV = mix(centeredUV, squareUV, 0.5);

    // Rescale mixedUV back to [0, 1] range
    mixedUV = (mixedUV + 1.0) / 2.0;

    // Standard spherical coordinates calculation with a rotation correction
    float u = mixedUV.x - 0.5;
    float v = 0.5 - mixedUV.y;
    float theta = atan(v, u) - PI / 2.0; // Subtract PI / 2 to rotate

    // Make sure theta wraps around correctly
    if(theta < 0.0) {
        theta += 2.0 * PI;
    }

    float phi0 = PI * length(vec2(u, v));
    float phi1 = r * PI * 0.5;
    float phi = mix(phi0, phi1, 0.5);
    phi = phi;

    // Convert spherical coordinates to Cartesian coordinates
    float x = sin(phi) * cos(theta);
    float y = cos(phi);
    float z = sin(phi) * sin(theta);

    return vec3(x, y, z);
}
vec3 mapUVToSphere2(vec2 uv) {
    // Convert uv coordinates so that (0.5, 0.5) is the north pole
    float u = uv.x - 0.5;
    float v = 0.5 - uv.y;

    float theta = atan(v, u);
    float phi = PI * length(vec2(u, v));

    float x = sin(phi) * cos(theta);
    float y = cos(phi);
    float z = sin(phi) * sin(theta);

    return vec3(x, y, z);
}



vec3 mapUVToSphere3(vec2 uv) {
    // Convert uv coordinates so that (0.5, 0.5) is the north pole
    float u = uv.x - 0.5;
    float v = 0.5 - uv.y;
    float r = max(abs(u), abs(v));
    r = r * 2.0;
    float theta = atan(v, u);
    float phi = r * PI * 0.5;

    float x = sin(phi) * cos(theta);
    float y = cos(phi);
    float z = sin(phi) * sin(theta);

    return vec3(x, y, z);
}
vec2 EncodeNormal3(vec3 normal) {
    float p = sqrt(normal.z * 8.0 + 8.0);
    return vec2(normal.xy / p + 0.5);
}

vec3 DecodeNormal3(vec2 enc) {
    vec2 fenc = enc * 4.0 - 2.0;
    float f = dot(fenc, fenc);
    float g = sqrt(1.0 - f / 4.0);
    vec3 normal;
    normal.xy = fenc * g;
    normal.z = 1.0 - f / 2.0;
    return normal;
}

vec4 PackFloatRGBA(float v) {
    vec4 enc = vec4(1.0, 255.0, 65025.0, 16581375.0) * v;
    enc = fract(enc);
    enc -= enc.yzww * vec4(1.0 / 255.0, 1.0 / 255.0, 1.0 / 255.0, 0.0);
    return enc;
}

float UnpackFloatRGBA(vec4 rgba) {
    return dot(rgba, vec4(1.0, 1.0 / 255.0, 1.0 / 65025.0, 1.0 / 16581375.0));
}

// float PackTwo16BitTo32Bit(float a, float b)
// {
// 	a = clamp(a, 0.0, 1.0);
// 	b = clamp(b, 0.0, 1.0);

// 	a *= 65536.0;
// 	b *= 65536.0;

// 	int ai = int(a);
// 	int bi = int(b);

// 	int data = ai << 16;
// 	data += bi & 0x0000FFFF;

// 	float dataf = float(data) / 0xFFFFFFFF;

// 	return dataf;
// }

// vec2 UnpackTwo16BitFrom32Bit(float value)
// {
// 	int data = int(value * 0xFFFFFFFF);

// 	int ai = data >> 16;
// 	int bi = data & 0x0000FFFF;

// 	float a = float(ai) / 65536.0;
// 	float b = float(bi) / 65536.0;

// 	return vec2(a, b);
// }

#if 0

float PackTwo16BitTo32Bit(float a, float b) {
    float data;

    a = clamp(a, 0.0, 255.0 / 256.0);
    b = clamp(b, 0.0, 255.0 / 256.0);

    a *= 65535.0;
    b *= 65535.0;

    a = floor(a);
    b = floor(b);

    data = a * exp2(16.0);
    data += b;

    data /= exp2(32.0) - 1;

    return data;

	// float data;

	// a = clamp(a, 0.01, 2047.0 / 2048.0);
	// b = clamp(b, 0.01, 2047.0 / 2048.0);

	// a *= 2047.0;
	// b *= 2047.0;

	// a = floor(a);
	// b = floor(b);

	// data = a * exp2(11.0);
	// data += b;

	// data /= exp2(22.0) - 1;
	// data += 1.0;

	// return data;
}

vec2 UnpackTwo16BitFrom32Bit(float value) {
    vec2 data;

    value *= exp2(32.0) - 1;

    data.x = floor(value / exp2(16.0));
    data.y = mod(value, exp2(16.0));

    data.x /= 65535.0;
    data.y /= 65535.0;

    return data;

	// vec2 data;

	// value -= 1.0;
	// value *= exp2(22.0) - 1;

	// data.x = floor(value / exp2(11.0));
	// data.y = mod(value, exp2(11.0));

	// data.x /= 2047.0;
	// data.y /= 2047.0;

	// return data;
}

#else

float PackTwo16BitTo32Bit(float a, float b) {
    vec2 v = vec2(a, b);
	// v = clamp(v, vec2(0.0), vec2(1.0));
    return dot(floor(v * 8191.9999), vec2(1. / 8192., 1.));
}
vec2 UnpackTwo16BitFrom32Bit(float v) {
    return vec2(fract(v) * (8192. / 8191.), floor(v) / 8191.);
}

#endif

vec3 rgb2hsv(vec3 c) {
    vec4 K = vec4(0.0, -1.0 / 3.0, 2.0 / 3.0, -1.0);
    vec4 p = mix(vec4(c.bg, K.wz), vec4(c.gb, K.xy), step(c.b, c.g));
    vec4 q = mix(vec4(p.xyw, c.r), vec4(c.r, p.yzx), step(p.x, c.r));

    float d = q.x - min(q.w, q.y);
    float e = 1.0e-10;
    return vec3(abs(q.z + (q.w - q.y) / (6.0 * d + e)), d / (q.x + e), q.x);
}

vec3 hsv2rgb(vec3 c) {
    vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
    vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
    return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

vec3 adjustVibrancy(vec3 color, float vibrancy) {
    vec3 hsv = rgb2hsv(color);
    hsv.y *= vibrancy; // Adjust the saturation
    return hsv2rgb(hsv);
}

vec2 encodeUnitVector(vec3 vector) {
	// Scale down to octahedron, project onto XY plane
    vector.xy /= abs(vector.x) + abs(vector.y) + abs(vector.z);
	// Reflect -Z hemisphere folds over the diagonals
    return vector.z <= 0.0 ? (1.0 - abs(vector.yx)) * vec2(vector.x >= 0.0 ? 1.0 : -1.0, vector.y >= 0.0 ? 1.0 : -1.0) : vector.xy;
}

vec3 decodeUnitVector(vec2 encoded) {
	// Exctract Z component
    vec3 vector = vec3(encoded, 1.0 - abs(encoded.x) - abs(encoded.y));
	// Reflect -Z hemisphere folds over the diagonals
    float t = max(-vector.z, 0.0);
    vector.xy += vec2(vector.x >= 0.0 ? -t : t, vector.y >= 0.0 ? -t : t);
	// Normalize and return
    return normalize(vector);
}
vec2 EncodeNormal(vec3 n) {
    float f = sqrt(n.z * 8.0 + 8.0);
    return n.xy / f + 0.5;
}

vec3 DecodeNormal(vec2 enc) {
    vec2 fenc = enc * 4.0 - 2.0;
    float f = dot(fenc, fenc);
    float g = sqrt(1.0 - f / 4.0);
    vec3 n;
    n.xy = fenc * g;
    n.z = 1.0 - f / 2.0;
    return n;

}//Packing functions
float PackTwo4BitTo8Bit(float a, float b) {
    float data;

    a = clamp(a, 0.0, 15.0 / 16.0);
    b = clamp(b, 0.0, 15.0 / 16.0);

    a *= 15.0;
    b *= 15.0;

    a = floor(a);
    b = floor(b);

    data = a * exp2(4.0);
    data += b;

    data /= exp2(8.0) - 1;

    return data;
}

vec2 UnpackTwo4BitFrom8Bit(float value) {
    vec2 data;

    value *= exp2(8.0) - 1;

    data.x = floor(value / exp2(4.0));
    data.y = mod(value, exp2(4.0));

    data.x /= 15.0;
    data.y /= 15.0;

    return data;
}

float PackTwo8BitTo16Bit(float a, float b) {
    float data;

    a = clamp(a, 0.0, 255.0 / 256.0);
    b = clamp(b, 0.0, 255.0 / 256.0);

    a *= 255.0;
    b *= 255.0;

    a = floor(a);
    b = floor(b);

    data = a * exp2(8.0);
    data += b;

    data /= exp2(16.0) - 1;

    return data;

	// vec2 d = vec2(a, b);
	// d = clamp(d, vec2(0.0), vec2(255.0 / 256.0));
	// d *= 256.0;

	// float data = dot(d, vec2(1.0 / exp2(8.0), 1.0 / exp2(16.0)));

	// return data;
}

float PackTwo5BitTo10Bit(float a, float b) {
    float data;

    // Clamp the inp values to fit in 5 bits
    a = clamp(a, 0.0, 31.0 / 32.0);
    b = clamp(b, 0.0, 31.0 / 32.0);

    // Scale up to the full 5-bit range (0 to 31)
    a *= 31.0;
    b *= 31.0;

    // Truncate to integer values
    a = floor(a);
    b = floor(b);

    // Pack a into the high 5 bits and b into the low 5 bits
    data = a * exp2(5.0) + b;

    // Normalize to range [0, 1] by dividing by the max 10-bit value (1023)
    data /= exp2(10.0) - 1;

    return data;
}

vec2 UnpackTwo5BitFrom10Bit(float value) {
    vec2 data;

    // Scale back from [0, 1] to the full 10-bit range
    value *= exp2(10.0) - 1;

    // Extract the high and low 5-bit values
    data.x = floor(value / exp2(5.0));
    data.y = mod(value, exp2(5.0));

    // Normalize the components to range [0, 1]
    data.x /= 31.0;
    data.y /= 31.0;

    return data;
}

vec2 UnpackTwo8BitFrom16Bit(float value) {
    vec2 data;

    value *= exp2(16.0) - 1;

    data.x = floor(value / exp2(8.0));
    data.y = mod(value, exp2(8.0));

    data.x /= 255.0;
    data.y /= 255.0;

    return data;
}
/*
    Normals encoding and decoding based on Spectrum by Zombye
*/
vec2 encodeNormal(in vec3 normal) {
    normal.xy /= abs(normal.x) + abs(normal.y) + abs(normal.z);
    return (normal.z <= 0.0 ? (1.0 - abs(normal.yx)) * vec2(normal.x >= 0.0 ? 1.0 : -1.0, normal.y >= 0.0 ? 1.0 : -1.0) : normal.xy) * 0.5 + 0.5;
}
vec3 decodeNormal(in vec2 encodedNormal) {
    encodedNormal = encodedNormal * 2.0 - 1.0;
    vec3 normal = vec3(encodedNormal, 1.0 - abs(encodedNormal.x) - abs(encodedNormal.y));
    float t = max(-normal.z, 0.0);
    normal.xy += vec2(normal.x >= 0.0 ? -t : t, normal.y >= 0.0 ? -t : t);
    return normalize(normal);
}// encode normal in two channels (xy),torch(z) and sky lightmap (w)
vec4 encode(vec3 unenc, vec2 lightmaps) {
    unenc.xy = unenc.xy / dot(abs(unenc), vec3(1.0)) + 0.00390625;
    unenc.xy = unenc.z <= 0.0 ? (1.0 - abs(unenc.yx)) * sign(unenc.xy) : unenc.xy;
    // unenc.xy = mix(unenc.xy, (1.0 - abs(unenc.yx)) * sign(unenc.xy), min(unenc.z, 0));
    vec2 encn = unenc.xy * 0.5 + 0.5;
    return vec4((encn), vec2(lightmaps.x, lightmaps.y));
}
vec3 decode(vec2 encn) {
    vec3 unenc = vec3(0.0);
    encn = encn * 2.0 - 1.0;
    unenc.xy = abs(encn);
    unenc.z = 1.0 - unenc.x - unenc.y;
    // unenc.xy = unenc.z <= 0.0 ? (1.0 - unenc.yx) * sign(encn) : encn;
    unenc.xy = mix(encn, (1.0 - unenc.yx) * sign(encn), min(0.0, unenc.z));
    return normalize(unenc.xyz);
}
vec3 decode21(vec2 enc) {
    // Unpack from [0, 1] to [-1, 1]
    vec2 encn = enc * 2.0 - 1.0;

    // Compute the z-component based on the normalization constraint
    float z = 1.0 - length(encn);
    vec3 unenc = vec3(encn, z);

    // Normalize the vector to account for any minor precision issues
    return normalize(unenc);
}

vec3 decode22(vec2 enc) {
    vec2 fenc = enc * 4 - 2;
    float f = dot(fenc, fenc);
    float g = sqrt(1 - f / 4.0);
    vec3 n;
    n.xy = fenc * g;
    n.z = 1 - f / 2;
    return n;
}
vec2 decodeVec2(float a) {
    int bf = int(a * 65535.);
    return vec2(bf - (256 * (bf / 256)), bf >> 8) / 255.;
}

float encodeVec2v2(vec2 a) {
    ivec2 bf = ivec2(a * 255.);
    return float(bf.x | (bf.y << 8)) / 65535.;
}


vec3 clamp10(vec3 value) {
    // Define approximate safe ranges for R11F_G11F_B10F to avoid NaN.
    // These ranges are conservative to avoid edge cases near the limits.
    float minVal = -1024.0; // Minimum representable value
    float maxVal = 1024.0;  // Maximum representable value before risking NaN

    // Clamp values to the safe range
    value = clamp(value, vec3(minVal), vec3(maxVal));

    return value;
}

vec4 clamp10(vec4 value) {
    // Define approximate safe ranges for R11F_G11F_B10F to avoid NaN.
    // These ranges are conservative to avoid edge cases near the limits.
    float minVal = -1024.0; // Minimum representable value
    float maxVal = 1024.0;  // Maximum representable value before risking NaN

    // Clamp values to the safe range
    value = clamp(value, vec4(minVal), vec4(maxVal));

    return value;
}

vec3 peirceQuincuncial(vec2 uv) {
    // Normalize UV coordinates to range [-1, 1]
    vec2 normalizedUV = 2.0 * uv - 1.0;

    // Calculate angles
    float phi = (TWO_PI * normalizedUV.x / 2.0);
    float theta = PI * normalizedUV.y / 2.0;

    // Perform Peirce quincuncial projection
    float x = cos(theta) * sin(phi);
    float y = cos(theta) * cos(phi);
    float z = sin(theta);

    // Return the result as a vec3 in world space
    return vec3(x, y, z);
}
vec2 inversePeirceQuincuncial(vec3 xyz) {
    // Adjust coordinates for the 90-degree tilt around the y-axis
    float x = -xyz.z; // x becomes -z
    float y = xyz.y;  // y remains y
    float z = xyz.x;  // z becomes x

    // Calculate theta and phi from the adjusted coordinates
    float theta = asin(x); // Using adjusted x for theta calculation
    float phi = atan(z, y); // Using adjusted z and y for phi calculation, ensuring correct quadrant handling

    // Convert theta and phi back to normalized UV coordinates, considering the tilt adjustments
    float u = phi / TWO_PI + 0.5;
    float v = (theta / PI) + 0.5;

    // Normalize UV to range [0, 1]
    return vec2(u, v);
}
// Function to map a point in 3D space onto a 2D Lambert Azimuthal Equal Area projection
vec2 InverseLambertAzimuthalEqualArea(vec3 xyz) {
    // Adjust for y-up coordinate system by swapping y (up) and z (depth)
    float x = xyz.x;
    float y = xyz.z; // Adjust for the y-up system where z represents forward/upward.
    float z = xyz.y; // Adjust for the y-up system, treating y as depth.

    // Normalize inp coordinates to ensure they are on the unit sphere
    vec3 normXYZ = normalize(vec3(x, y, z));

    // Calculate spherical coordinates
    float lambda = atan(normXYZ.x, normXYZ.z); // Longitude calculation
    float phi = asin(normXYZ.y); // Latitude calculation, using adjusted y

    // Conversion to Lambert Azimuthal Equal-Area projection
    float k = sqrt(0.5 / (1.0 + cos(phi) * cos(lambda)));
    vec2 mappedUV = k * vec2(cos(phi) * sin(lambda), sin(phi));

    // No need to adjust mappedUV.y for hemisphere; previous attempts might have over-corrected.

    // Transform UV coordinates from [-1, 1] to [0, 1] for texture mapping
    vec2 uv = (mappedUV + vec2(1.0)) / 2.0;

    return uv;
}

vec3 simulateAzimuthalProjectionEffect(vec2 uv) {
    // Normalize UV coordinates to range from -1 to 1
    vec2 centeredUV = uv * 2.0 - 1.0;
    vec2 centeredUVSquared = centeredUV * centeredUV;

    // Calculate radial distance from the center of the UV space
    float radius = pow(length(centeredUV), 1.0);
    if(radius > 1.0)
        radius = 1.0; // Clamp radius to the unit circle

    // Calculate azimuthal angle from UV coordinates
    float angle = atan(centeredUV.y, centeredUV.x);

    // Adjust this factor to control the spread effect
    float simulatedPhi = radius * PI / 1;

    // Use the calculated azimuthal angle directly
    float phi = angle;

    // Convert polar coordinates back to Cartesian for the view vector
    vec3 viewVector = vec3(sin(simulatedPhi) * cos(phi), sin(simulatedPhi) * sin(phi), cos(simulatedPhi));

    // Rotate around the X-axis by π/2 to simulate viewing from the equator
    float cosLatOffset = cos(1.5708); // cos(π/2) = 0
    float sinLatOffset = sin(-1.5708); // sin(π/2) = 1
    vec3 rotatedViewVector = vec3(viewVector.x, viewVector.y * cosLatOffset - viewVector.z * sinLatOffset, viewVector.y * sinLatOffset + viewVector.z * cosLatOffset);

    return rotatedViewVector;
}

vec2 invertAzimuthalProjectionEffect(vec3 viewVector) {
    // Reverse the rotation around the X-axis
    float cosLatOffset = cos(1.5708); // Reverse rotation, so use negative angle
    float sinLatOffset = sin(-1.5708);
    vec3 originalViewVector = vec3(viewVector.x, viewVector.y * cosLatOffset + viewVector.z * sinLatOffset, -viewVector.y * sinLatOffset + viewVector.z * cosLatOffset);

    // Convert from Cartesian back to polar coordinates
    float simulatedPhi = acos(originalViewVector.z);
    float phi = atan(originalViewVector.y, originalViewVector.x);

    // Reverse the spread effect adjustment to find the original radius
    // Assuming PI was used in the forward projection
    float radius = simulatedPhi / PI * 1;

    // Clamp the radius back if needed (here it's for consistency, as in practice it should not exceed 1)
    if(radius > 1.0)
        radius = 1.0;
    float X = 1.0;

    radius = pow(radius, 1.0 / X);
    // Reverse the calculation of radial distance and azimuthal angle to get original UV coordinates
    vec2 centeredUV = vec2(cos(phi) * radius, sin(phi) * radius);

    // Normalize UV coordinates back to the 0 to 1 range
    vec2 uv = (centeredUV + 1.0) / 2.0;
    return uv;
}

vec3 mapUVToSphere(vec2 uv) {

    return spherical_surface(uv);

    float lon = (uv.x - 0.5) * TWO_PI;

    float lat = uv.y * PI;

    float x = sin(lat) * cos(lon);
    float y = cos(lat);
    float z = sin(lat) * sin(lon);

    return normalize(vec3(x, y, z));
}

vec3 peirceQuincuncial(vec2 uv) {
    // Normalize UV coordinates to range [-1, 1]
    vec2 normalizedUV = 2.0 * uv - 1.0;

    // Calculate angles
    float phi = (TWO_PI * normalizedUV.x / 2.0);
    float theta = PI * normalizedUV.y / 2.0;

    // Perform Peirce quincuncial projection
    float x = cos(theta) * sin(phi);
    float y = cos(theta) * cos(phi);
    float z = sin(theta);

    // Return the result as a vec3 in world space
    return vec3(x, y, z);
}
vec2 inversePeirceQuincuncial(vec3 xyz) {
    // Adjust coordinates for the 90-degree tilt around the y-axis
    float x = -xyz.z; // x becomes -z
    float y = xyz.y;  // y remains y
    float z = xyz.x;  // z becomes x

    // Calculate theta and phi from the adjusted coordinates
    float theta = asin(x); // Using adjusted x for theta calculation
    float phi = atan(z, y); // Using adjusted z and y for phi calculation, ensuring correct quadrant handling

    // Convert theta and phi back to normalized UV coordinates, considering the tilt adjustments
    float u = phi / TWO_PI + 0.5;
    float v = (theta / PI) + 0.5;

    // Normalize UV to range [0, 1]
    return vec2(u, v);
}
// Function to map a point in 3D space onto a 2D Lambert Azimuthal Equal Area projection
vec2 InverseLambertAzimuthalEqualArea(vec3 xyz) {
    // Adjust for y-up coordinate system by swapping y (up) and z (depth)
    float x = xyz.x;
    float y = xyz.z; // Adjust for the y-up system where z represents forward/upward.
    float z = xyz.y; // Adjust for the y-up system, treating y as depth.

    // Normalize inp coordinates to ensure they are on the unit sphere
    vec3 normXYZ = normalize(vec3(x, y, z));

    // Calculate spherical coordinates
    float lambda = atan(normXYZ.x, normXYZ.z); // Longitude calculation
    float phi = asin(normXYZ.y); // Latitude calculation, using adjusted y

    // Conversion to Lambert Azimuthal Equal-Area projection
    float k = sqrt(0.5 / (1.0 + cos(phi) * cos(lambda)));
    vec2 mappedUV = k * vec2(cos(phi) * sin(lambda), sin(phi));

    // No need to adjust mappedUV.y for hemisphere; previous attempts might have over-corrected.

    // Transform UV coordinates from [-1, 1] to [0, 1] for texture mapping
    vec2 uv = (mappedUV + vec2(1.0)) / 2.0;

    return uv;
}
vec3 simulateAzimuthalProjectionEffect(vec2 uv) {
    // Normalize UV coordinates to range from -1 to 1
    vec2 centeredUV = uv * 2.0 - 1.0;
    vec2 centeredUVSquared = centeredUV * centeredUV;

    // Calculate radial distance from the center of the UV space
    float radius = pow(length(centeredUV), 1.0);
    if(radius > 1.0)
        radius = 1.0; // Clamp radius to the unit circle

    // Calculate azimuthal angle from UV coordinates
    float angle = atan(centeredUV.y, centeredUV.x);

    // Adjust this factor to control the spread effect
    float simulatedPhi = radius * PI / 1;

    // Use the calculated azimuthal angle directly
    float phi = angle;

    // Convert polar coordinates back to Cartesian for the view vector
    vec3 viewVector = vec3(sin(simulatedPhi) * cos(phi), sin(simulatedPhi) * sin(phi), cos(simulatedPhi));

    // Rotate around the X-axis by π/2 to simulate viewing from the equator
    float cosLatOffset = cos(1.5708); // cos(π/2) = 0
    float sinLatOffset = sin(-1.5708); // sin(π/2) = 1
    vec3 rotatedViewVector = vec3(viewVector.x, viewVector.y * cosLatOffset - viewVector.z * sinLatOffset, viewVector.y * sinLatOffset + viewVector.z * cosLatOffset);

    return rotatedViewVector;
}

vec2 invertAzimuthalProjectionEffect(vec3 viewVector) {
    // Reverse the rotation around the X-axis
    float cosLatOffset = cos(1.5708); // Reverse rotation, so use negative angle
    float sinLatOffset = sin(-1.5708);
    vec3 originalViewVector = vec3(viewVector.x, viewVector.y * cosLatOffset + viewVector.z * sinLatOffset, -viewVector.y * sinLatOffset + viewVector.z * cosLatOffset);

    // Convert from Cartesian back to polar coordinates
    float simulatedPhi = acos(originalViewVector.z);
    float phi = atan(originalViewVector.y, originalViewVector.x);

    // Reverse the spread effect adjustment to find the original radius
    // Assuming PI was used in the forward projection
    float radius = simulatedPhi / PI * 1;

    // Clamp the radius back if needed (here it's for consistency, as in practice it should not exceed 1)
    if(radius > 1.0)
        radius = 1.0;
    float X = 1.0;

    radius = pow(radius, 1.0 / X);
    // Reverse the calculation of radial distance and azimuthal angle to get original UV coordinates
    vec2 centeredUV = vec2(cos(phi) * radius, sin(phi) * radius);

    // Normalize UV coordinates back to the 0 to 1 range
    vec2 uv = (centeredUV + 1.0) / 2.0;
    return uv;
}

vec2 remapY(vec2 inp) {

    float newY;
    float newX;
    if(inp.y < pivot) {
        newY = lowerStart + (inp.y / pivot) * (lowerEnd - lowerStart);
    } else {
        newY = upperStart + ((inp.y - pivot) / (1.0 - pivot)) * (upperEnd - upperStart);
    }

    if(inp.x < pivot) {
        newX = lowerStart + (inp.x / pivot) * (lowerEnd - lowerStart);
    } else {
        newX = upperStart + ((inp.x - pivot) / (1.0 - pivot)) * (upperEnd - upperStart);
    }
    //return inp;
    //return vec2(newX, inp.y);
    return vec2(inp.x, newY);
}
vec2 inverseRemapY(vec2 inp) {

    float newY;
    float newX;
    if(inp.y < lowerEnd) {
        newY = (inp.y - lowerStart) / (lowerEnd - lowerStart) * pivot;
    } else {
        newY = pivot + (inp.y - upperStart) / (upperEnd - upperStart) * (1.0 - pivot);
    }
    if(inp.x < lowerEnd) {
        newX = (inp.x - lowerStart) / (lowerEnd - lowerStart) * pivot;
    } else {
        newX = pivot + (inp.x - upperStart) / (upperEnd - upperStart) * (1.0 - pivot);
    }
    //return inp;

    //return vec2(newX, inp.y);
    return vec2(inp.x, newY);
}
vec2 sphereToCarteImproved(vec3 dir) {

    return spherical_coordinates(normalize(dir));

    // Normalize the direction vector
    dir = normalize(dir);

    // Calculate longitude and latitude
    float lon = atan(dir.z, dir.x);
    float lat = asin(dir.y);

    float u = 0.5 + lon / TWO_PI;

    float v = 0.5 - lat / PI;

    vec2 uv = (vec2(u, 1.0 - v));

    return uv;
}
/*
//////////////////

// Main function
// Transmittance stores atmospheric transmittance in xyz and planet intersection flag in w.
// Fog parameter can be used in cases where you want to fade the world out
// (like in a game with limited far clip plane).

vec3 GetAtmosphere(vec3 rayDir) {

// Config
    const float DENSITY = 1.0; // Atmosphere density. 1 is Earth-like.
    const float AERIAL_SCALE = 1.0; // Higher value = more aerial perspective. A value of 1 is tuned to match reference implementation.
    const float EXPOSURE = 4.0; // Tuned to match reference.

// Math
    const float INFINITY = 3.402823466e38;

// Atmosphere
    const float ATMOSPHERE_HEIGHT = 100000.0;
    const float PLANET_RADIUS = 6371000.0;
    const vec3 PLANET_CENTER = vec3(0, -PLANET_RADIUS, 0);
    const vec3 C_RAYLEIGH = vec3(5.802, 13.558, 33.100) * 1e-6;
    const vec3 C_MIE = vec3(3.996, 3.996, 3.996) * 1e-6;
    const vec3 C_OZONE = vec3(0.650, 1.881, 0.085) * 1e-6;
    const float RAYLEIGH_MAX_LUM = 2.0;
    const float MIE_MAX_LUM = 0.5;

// Magic numbers
    const float M_TRANSMITTANCE = 0.75;
    const float M_LIGHT_TRANSMITTANCE = 1e6;
    const float M_MIN_LIGHT_ELEVATION = -0.4;
    const float M_DENSITY_HEIGHT_MOD = 1e-12;
    const float M_DENSITY_CAM_MOD = 10.0;

	// Planet and atmosphere intersection to get optical depth
    vec3 oc = 0.0 - PLANET_CENTER;
    float b0 = dot(oc, rayDir);

    // Assuming we're on the surface, we're only interested in when the ray exits the atmosphere.
    // Since we're starting from the surface, c0 simplifies significantly.
    float atmosphereRadius = PLANET_RADIUS + ATMOSPHERE_HEIGHT;
    float c0 = dot(oc, oc) - atmosphereRadius * atmosphereRadius;

    // c0 is zero if you're exactly on the surface, leading to simplification.
    float discriminant = b0 * b0 - c0; // This simplifies to b0 * b0 if starting on the surface.

    float sqrtDiscriminant = sqrt(discriminant);
    vec2 t2;
    t2.x = -b0 - sqrtDiscriminant; // This is towards the planet's center and can be ignored if you're on the surface looking outward.
    t2.y = -b0 + sqrtDiscriminant; // This is the exit point through the atmosphere.

    //float altitude = (length(rayStart - PLANET_CENTER) - PLANET_RADIUS);
    float normAltitude = 0.0 / ATMOSPHERE_HEIGHT;

    float opticalDepth = t2.y;

    // Optical depth modulators
    opticalDepth = min(INFINITY, opticalDepth);
    opticalDepth = min(opticalDepth * AERIAL_SCALE, t2.y);

    // Altitude-based density modulators
    float altitudeFactor = 1.0 / (2.0 + (t2.y * t2.y) * M_DENSITY_HEIGHT_MOD);
    float h = pow(1.0 - altitudeFactor, 1.0 + normAltitude * M_DENSITY_CAM_MOD);

    // Use h directly to calculate densities, avoiding redundant square calculation
    float densityR = h * h * DENSITY; // Density for Rayleigh scattering
    float densityM = h * h * h * h * h * DENSITY; // Density for Mie scattering, equivalent to (h^2)^2 * h

    // Apply light transmittance (makes sky red as sun approaches horizon)
    vec3 lightColor = vec3(1 - nightmix) * exp(-(C_RAYLEIGH + C_MIE + C_OZONE) * pow(smoothstep(1.0, M_MIN_LIGHT_ELEVATION, sunPosWorld.y), 8.0) * DENSITY * h * M_LIGHT_TRANSMITTANCE);

    // Approximate marched Rayleigh + Mie scattering with some exp magic.

    vec3 R = (1.0 - exp(-opticalDepth * densityR * C_RAYLEIGH / RAYLEIGH_MAX_LUM)) * RAYLEIGH_MAX_LUM;
    vec3 M = (1.0 - exp(-opticalDepth * densityM * C_MIE / MIE_MAX_LUM)) * MIE_MAX_LUM;

    float costh = dot(rayDir, sunPosWorld);
    float phaseR = (1.0 + costh * costh) * 0.06; // Rayleigh scattering phase

    float g = 0.85; // Asymmetry factor for the Mie scattering phase
    float g2 = g * g;
    float phaseM = (1.0 - g2) / (4.0 * PI * pow(1.0 + g2 - 2.0 * g * costh, 1.5)); // Mie scattering phase

    vec3 A = phaseR * lightColor;
    vec3 B = phaseM * lightColor;
    vec3 C = (C_RAYLEIGH * densityR + C_MIE * densityM + C_OZONE * densityR) * pow(1.0 - normAltitude, 4.0) * M_TRANSMITTANCE;

        // Combined scattering
    vec3 scattering = R * A + M * B;

        // View extinction, matched to reference
    vec3 transmittance = exp(-opticalDepth * C);

    //return vec3(rayDir * 0.1);
    //return vec3(4.0 * PI * pow(1.0 + g2 - 2.0 * g * costh, 1.5));
    return scattering * EXPOSURE;

}

vec3 GetAtmosphereSimplified(vec3 rayDir) {
    vec3 sunPos = -sunPosWorld;

    float b0 = dot(vec3(0, 6371000.0, 0), rayDir);
    float costh = dot(rayDir, sunPos);

    float t2 = -b0 + sqrt(pow(b0, 2.0) + 1284200000000.0);

    float h = 1.0 - (1.0 / (2.0 + pow(t2, 2.0) * 0.000000000001));

    float phaseR = (1.0 + pow(costh, 2.0)) * 0.06;

    float phaseM = 0.00977 * pow(1.0 - 0.987 * costh, -1.5);
    float dynamicFactor = pow(smoothstep(1.0, -0.4, sunPos.y), 8.0) * h;
    float powh2 = pow(h, 2.0);
    float powh5 = pow(h, 5.0);

    vec3 opticalDepth = vec3(t2, t2, t2);
    vec3 lightColor = vec3(nightmixInv, nightmixInv, nightmixInv) * vec3(exp(-10448.0 * dynamicFactor), exp(-19435.0 * dynamicFactor), exp(-37181.0 * dynamicFactor));
    const vec3 vec3one = vec3(1.0, 1.0, 1.0);
    const vec3 vec3two = vec3(1.0, 1.0, 1.0);
    const vec3 vec3half = vec3(1.0, 1.0, 1.0);
    vec3 R = (vec3one - exp(-opticalDepth * vec3(0.000002901 * powh2, 0.000006779 * powh2, 0.00001655 * powh2))) * vec3two;
    vec3 M = (vec3one - exp(-opticalDepth * vec3(0.000007992 * powh5, 0.000007992 * powh5, 0.000007992 * powh5))) * vec3half;

    vec3 A = vec3(phaseR, phaseR, phaseR) * lightColor;
    vec3 B = vec3(phaseM, phaseM, phaseM) * lightColor;

    vec3 scattering = R * A + M * B;
    scattering = scattering * vec3(4.0, 4.0, 4.0);
    return scattering;
}

// Constants that can be pre-defined outside GetAtmosphere2
const float ATMOSPHERE_RADIUS_SQUARED = 1284200000000.0; // Earth's radius squared plus atmosphere height squared
const float LIGHT_SCATTERING_COEFF = 0.0000000000009;
const vec3 RAYLEIGH_SCATTERING_COEFF = vec3(0.000002901, 0.000006779, 0.00001655);
const vec3 MIE_SCATTERING_COEFF = vec3(0.000007992, 0.000007992, 0.00000799);
const vec3 LIGHT_ATTENUATION = vec3(10.448, 19.435, 37.188);
const float EXPOSURE = 4.0;

vec3 GetAtmosphere2(vec3 rayDir) {
    vec3 sunPos = mix(sunPosWorld, -sunPosWorld, nightmix);

    float b0 = 6371000.0 * rayDir.y; // Simplified dot product to direct multiplication
    float costh = dot(rayDir, sunPos);

    float t2 = -b0 + sqrt(b0 * b0 + ATMOSPHERE_RADIUS_SQUARED);
    float h = 1.0 - 1.0 / (2.0 + t2 * t2 * LIGHT_SCATTERING_COEFF);
    float densityR = h * h;
    float densityM = densityR * densityR * h;

    float phaseR = (1.0 + costh * costh) * 0.12;
    float phaseM = (0.2775 / (12.5664 * sqrt(1.7225 - 1.7 * costh) * (1.7225 - 1.7 * costh))) * 0.5;
    float smoothStepVal = smoothstep(1.0, -0.4, sunPos.y);
    vec3 night = mix(vec3(0.002, 0.005, 0.006) * 3.14, vec3(0.0), 1 - nightmix);
    vec3 lightColor = nightmixInv * exp(-LIGHT_ATTENUATION * pow(smoothStepVal, 8.0) * h) + night;
    vec3 R = (1.0 - exp(-t2 * densityR * RAYLEIGH_SCATTERING_COEFF)) * phaseR * lightColor;
    vec3 M = (1.0 - exp(-t2 * densityM * MIE_SCATTERING_COEFF)) * phaseM * lightColor;

    vec3 scattering = (R + M) * EXPOSURE;
    return scattering;
}





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

vec3 calcAtmosphericScatter3(vec3 pos) {

    float plp = clamp(length(pos - sunPosWorld), -1, 1);

    float maxPosY = max(pos.y + 0.04, 0.0001);
    float optDepthDiv = (maxPosY * maxPosY * (3.0 - 2.0 * maxPosY));

    vec3 opticalDepth = atmosphericScaleV3 / vec3(optDepthDiv, optDepthDiv, optDepthDiv);

    vec3 scatterView = (totalCoeff * opticalDepth);

    vec3 absorbView = exp2((totalCoeff * -opticalDepth));

    vec3 absorbDiff = abs(absorbLight - absorbView);
    vec3 scatterDiff = abs((scatterLight - scatterView) * 0.693147);

    vec3 absorbSun = absorbDiff / (scatterDiff + 1e-8);

    vec3 scatterSun = scatterView * (0.375 * (2.0 - 2.0 * plp + plp * plp));
    vec3 final = (scatterSun * absorbSun) * PI;
    float luma = dot(final, luminanceCoeff);
    //final = mix(final, skyColor.rgb * dot(final, luminanceCoeff), 0.75);
    final = mix(final, luma * skyColorNight, nightmix);

    final = mix(final, dot(final, luminanceCoeff) * rainColorMix, rainStrength);
    //final = final + (final - dot(final, vec3(0.299, 0.587, 0.114))) * 0.3;

    //final *= 0.5;

    //return vec3(dot(pos, sunPosWorld));
    return final * 2.0;
}
*/



vec3 sunFunc() {
    const vec3 LIGHT_ATTENUATION = vec3(10.448, 19.435, 37.188);

    float b0 = 6371000.0 * lightPos.y;

    float t2 = -b0 + sqrt(b0 * b0 + 1284200000000.0);

    float h = 1.0 - (1.0 / (2.0 + t2 * t2 * 0.0000000000009));
    float densityR = h * h;
    float densityM = densityR * densityR * h;

    float smoothStepVal = smoothstep(1.0, -0.4, lightPos.y);
    vec3 lightColor = lightRadiance * exp(-LIGHT_ATTENUATION * pow(smoothStepVal, 8.0) * h) + skycolnight;

    vec3 R = (1.0 - exp(-t2 * densityR * vec3(0.000002901, 0.000006779, 0.00001655))) * 0.24 * lightColor;
    vec3 M = (1.0 - exp(-t2 * densityM * vec3(0.000007992, 0.000007992, 0.00000799))) * 3.2715 * lightColor;

    vec3 scattering = (R + M) * 4.0;

    return vec3(smoothStepVal);
}

// Constants for R calculations
uniform vec3 kRCosthLeft50UpR;
uniform vec3 kRCosthRight50UpR;
uniform vec3 kRCosthForward50UpR;
uniform vec3 kRCosthBackward50UpR;
uniform vec3 kRCosthUpR;
uniform vec3 kRCosthForwardR;
uniform vec3 kRCosthBackwardR;
uniform vec3 kRCosthRightR;
uniform vec3 kRCosthLeftR;

// Constants for M calculations
uniform vec3 kMCosthLeft50UpM;
uniform vec3 kMCosthRight50UpM;
uniform vec3 kMCosthForward50UpM;
uniform vec3 kMCosthBackward50UpM;
uniform vec3 kMCosthUpM;
uniform vec3 kMCosthForwardM;
uniform vec3 kMCosthBackwardM;
uniform vec3 kRCosthRightM;
uniform vec3 kRCosthLeftM;
uniform vec3 R_UP;
uniform vec3 M_UP;
uniform vec3 RM_UP;

vec3 fogAtmosphereUp() {

    vec3 scattering = RM_UP;

    scattering = mix(vec3(scattering.b) * skyColorNight * 10, scattering, nightblend5 * 0.9);

    scattering = mix(scattering, vec3(scattering.b) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphereLeft() {

    vec3 R = kRCosthLeft50UpR * 0.4 + kRCosthUpR * 0.15 + kRCosthLeftR * 0.15 + kRCosthForward50UpR * 0.15 + kRCosthBackward50UpR * 0.15;

    vec3 M = kMCosthLeft50UpM * 0.4 + kMCosthUpM * 0.15 + kRCosthLeftM * 0.15 + kMCosthForward50UpM * 0.15 + kMCosthBackward50UpM * 0.15;

    vec3 scattering = (R + M) * sunpow;

    scattering = mix(vec3(scattering.b) * skyColorNight * 10, scattering, nightblend5 * 0.9);

    scattering = mix(scattering, vec3(scattering.b) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphereRight() {

    vec3 R = kRCosthRight50UpR * 0.4 + kRCosthUpR * 0.15 + kRCosthRightR * 0.15 + kRCosthForward50UpR * 0.15 + kRCosthBackward50UpR * 0.15;

    vec3 M = kMCosthRight50UpM * 0.4 + kMCosthUpM * 0.15 + kRCosthRightM * 0.15 + kMCosthForward50UpM * 0.15 + kMCosthBackward50UpM * 0.15;

    vec3 scattering = (R + M) * sunpow;

    scattering = mix(vec3(scattering.b) * skyColorNight * 10, scattering, nightblend5 * 0.9);

    scattering = mix(scattering, vec3(scattering.b) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphereFront() {

    vec3 R = kRCosthForward50UpR * 0.4 + kRCosthUpR * 0.15 + kRCosthForwardR * 0.15 + kRCosthLeft50UpR * 0.15 + kRCosthRight50UpR * 0.15;

    vec3 M = kMCosthForward50UpM * 0.4 + kMCosthUpM * 0.15 + kMCosthForwardM * 0.15 + kMCosthLeft50UpM * 0.15 + kMCosthRight50UpM * 0.15;

    vec3 scattering = (R + M) * sunpow;

    scattering = mix(vec3(scattering.b) * skyColorNight * 10, scattering, nightblend5 * 0.9);

    scattering = mix(scattering, vec3(scattering.b) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphereBack() {

    vec3 R = kRCosthBackward50UpR * 0.4 + kRCosthUpR * 0.15 + kRCosthBackwardR * 0.15 + kRCosthLeft50UpR * 0.15 + kRCosthRight50UpR * 0.15;

    vec3 M = kMCosthBackward50UpM * 0.4 + kMCosthUpM * 0.15 + kMCosthBackwardM * 0.15 + kMCosthLeft50UpM * 0.15 + kMCosthRight50UpM * 0.15;

    vec3 scattering = (R + M) * sunpow;

    scattering = mix(vec3(scattering.b) * skyColorNight * 10, scattering, nightblend5 * 0.9);

    scattering = mix(scattering, vec3(scattering.b) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphereDemo(vec3 rayDir) {

    rayDir = normalize(rayDir);

    float b = 6371000.0 * rayDir.y;

    float t = -b + sqrt(b * b + 1284200000000.0);

    float h = 1.0 - 1.0 / (2.0 + t * t * 0.0000000000009);
    float densityR = pow(h, 2);
    float densityM = pow(h, 5);

    float costh = dot(rayDir, lightPos);

    float phaseR = 0.12 * (1.0 + pow(costh, 2));
    float term = 1.72 - 1.7 * costh;

    float phaseM = 0.135 / (12.5 * sqrt(term) * term);
    float lposTemp = lightRadiance * pow(smoothstep(1.0, -0.4, lightPos.y), 8.0) * h;
    vec3 lightColor = exp(-vec3(10.448, 19.435, 37.188) * lposTemp) + skycolnight;

    vec3 R = (1.0 - exp(-t * densityR * vec3(0.000002901, 0.000006779, 0.00001655))) * phaseR * lightColor;
    vec3 M = (1.0 - exp(-t * densityM * vec3(0.000007992, 0.000007992, 0.00000799))) * phaseM * lightColor;
    vec3 scattering = (R + M) * sunpow;

    return scattering;
}
vec3 fogAtmosphereDemo2(vec3 rayDir) {

    rayDir = normalize(rayDir);

    float b = 6371000.0 * rayDir.y;

    float t = -b + sqrt(b * b + 1284200000000.0);

    float h = 1.0 - 1.0 / (2.0 + t * t * 0.0000000000009);
    float h2 = h * h;
    float h5 = h * h * h * h * h;

    float costh = dot(rayDir, lightPos);

    float phaseR = 0.12 * (1.0 + costh * costh);
    float a = 1.72 - 1.7 * costh;
    float phaseM = 0.0108 * 1 / max(a * a * 1.3 + a * 0.17, 0);
    //float phaseM = 0.0108 * 1 / (pow(a, 1.5));
//0.5 - 0.5*cos(pi*x)
    float lposTemp = lightRadiance * pow(smoothstep(1.0, -0.4, lightPos.y), 8.0) * h;
    vec3 lightColor = exp(-vec3(10.448, 19.435, 37.188) * lposTemp) + skycolnight;

    vec3 R = (1.0 - exp(-t * h2 * vec3(0.000002901, 0.000006779, 0.00001655))) * phaseR * lightColor;
    vec3 M = (1.0 - exp(-t * h5 * vec3(0.000007992, 0.000007992, 0.00000799))) * phaseM * lightColor;
    vec3 scattering = (R + M) * sunpow;

    return scattering;
}

vec3 fogAtmosphere299(vec3 rayDir, sampler2D none) {

    // Normalize the ray direction
    rayDir = normalize(rayDir);

    // Calculate t2
    float b0 = 6371000.0 * rayDir.y;
    float t2 = -b0 + sqrt(b0 * b0 + 1284200000000.0);

    // Calculate costh (cosine of the angle between ray direction and light position)
    float costheta = dot(rayDir, lightPos);

    // Calculate h and densities
    float h = 1.0 - 1.0 / (2.0 + t2 * t2 * 9e-13);
    float densityR = h * h;
    float densityM = h * h * h * h * h;

    // Phase functions
    float phaseR = 0.12 * (1.0 + costheta * costheta);
    float phaseM = (0.0052 + -0.0015 * costheta + -0.00249 * costheta * costheta) / (1.0 + -1.95 * costheta + 0.9502 * costheta * costheta);

    // Light color computation
    //float temp = pow(smoothstep(1.0, -0.4, lightPos.y), 8.0);
    float temp = 1.010 / (1.01 + exp(12.0 * (lightPos.y + 0.14)));

    vec3 lightColor = lightRadiance * exp(-vec3(10.448, 19.435, 37.188) * temp * h) + skycolnight;

    // Rayleigh and Mie scattering
    vec3 R = (1.0 - exp(-t2 * densityR * vec3(2.901e-6, 6.779e-6, 1.655e-5))) * phaseR;
    vec3 M = (1.0 - exp(-t2 * densityM * vec3(7.992e-6))) * phaseM;
    vec3 scattering = (R + M) * lightColor * sunpow;

    // Apply night blend
    scattering = mix(scattering, dot(scattering + 0.01, vec3(0.299, 0.587, 0.114)) * skyColorNight * 10.0, (1.0 - nightblend5) * 0.9);

    // Apply rain color mix
    scattering = mix(scattering, dot(scattering, vec3(0.299, 0.587, 0.114)) * rainColorMix, rainStrength);

    return scattering;
}
vec3 fogAtmosphere999(vec3 rayDir) {

    float b = 6371000.0 * rayDir.y;
    float t = -b + sqrt(b * b + 1284200000000.0);
    float h = 1.0 - 1.0 / (2.0 + t * t * 9e-13);

    float costheta = dot(rayDir, lightPos);
    float costheta2 = costheta * costheta;

    float phaseR = 0.12 * (1.0 + costheta2);
    float phaseM = (0.00704 + (-0.00638 * costheta)) / (1.0 + (-1.9995 * costheta) + costheta2);

    float temp = (1.010 / (1.01 + exp(12.0 * (lightPos.y + 0.14)))) * h;

    vec3 lightColor = lightRadiance * exp(-vec3(10.448, 19.435, 37.188) * temp);

    vec3 R = (1.0 - exp(-t * h * h * vec3(2.901e-6, 6.779e-6, 1.655e-5))) * phaseR;
    vec3 M = (1.0 - exp(-t * h * h * h * h * h * vec3(7.992e-6))) * phaseM;
    vec3 scattering = (R + M) * (lightColor + skycolnight) * sunpow;

    scattering = mix(scattering, scattering.b * skyColorNight * 10.0, (1.0 - nightblend5) * 0.9);

    scattering = mix(scattering, scattering.b * rainColorMix, rainStrength);

    return scattering;
}

vec3 fogAtmosphereTest(vec3 rayDir) {

    float t = -6371000.0 * rayDir.y + sqrt(40589641000000.0 * rayDir.y * rayDir.y + 1284200000000.0);

    float h = 1.0 - 1.0 / (2.0 + t * t * 0.0000000000009);

    float c = dot(rayDir, lightPos);

    float m = (1.0 - exp(-t * pow(h, 5) * 7.992e-6)) * 0.0108 / (pow(1.71 - 1.7 * c, 1.5));

    float temp = (1 / (1 + exp(12 * (lightPos.y + 0.14)))) * h;

    float Sr = ((1.0 - exp(-t * h * h * 2.901e-6)) * 0.12 * (1.0 + c * c) + m) * (lightRadiance * exp(-10.448 * temp) + skycolnight.r) * sunpow;
    float Sg = ((1.0 - exp(-t * h * h * 6.779e-6)) * 0.12 * (1.0 + c * c) + m) * (lightRadiance * exp(-19.435 * temp) + skycolnight.g) * sunpow;
    float Sb = ((1.0 - exp(-t * h * h * 1.655e-5)) * 0.12 * (1.0 + c * c) + m) * (lightRadiance * exp(-37.188 * temp) + skycolnight.b) * sunpow;

    vec3 scattering = vec3(Sr, Sg, Sb);
    scattering = mix(scattering, scattering.b * skyColorNight * 10.0, (1.0 - nightblend5) * 0.9);

    scattering = mix(scattering, scattering.b * rainColorMix, rainStrength);

    return scattering;
}
void calculateLighting(
       out vec3 ambientUp,
       out vec3 ambientLeft,
       out vec3 ambientRight,
       out vec3 ambientB,
       out vec3 ambientF,
       out vec3 ambientDown,
       out float ambientA,
       out vec4 sunColor
) {
       float sunVis = clamp(sunElevation, 0.0, 0.05) / 0.05 * clamp(sunElevation, 0.0, 0.05) / 0.05;
       sunColor = vec4(sunLight, sunVis);

       ambientUp = fogAtmosphereUp();
       ambientDown = fogAtmosphereUp();

       ambientLeft = fogAtmosphereLeft();
       ambientRight = fogAtmosphereRight();

       ambientF = fogAtmosphereFront();
       ambientB = fogAtmosphereBack();

       ambientA = luma((ambientUp + ambientLeft + ambientRight + ambientB + ambientF + ambientDown) / 6);

}


vec2 OffsetDist(float x)
{
    float n = fract(x * 8.0) * pi;
    return vec2(cos(n), sin(n)) * x;
}
vec4 ld4(vec4 dist)
{
    return (2.0 * near) / (far + near - dist * (far - near));
}
vec4 textureGatherOffsetsApprox(sampler2D tex, vec2 uv, ivec2 size, ivec2 offsets[4])
{
    vec4 result = vec4(0.0);
    for (int i = 0; i < 4; i++)
    {
        // Convert UV to texel space
        ivec2 texelPos = ivec2(uv * vec2(size)) + offsets[i];
        // Fetch the texel value using texelFetch for nearest sampling
        result[i] = texelFetch(tex, texelPos, 0).r;
    }
    return result;
}

float ssao3(float dither, float depth)
{
    float ao = 0.0;

    depth = ld(depth);

    vec2 scale = dither * gl_FragCoord.xy * 0.5;
    float mult = 0.5 * (far - near);

    const vec4 steps = 0.4 + vec4(0.0, 0.2, 0.4, 0.6);

#if defined MC_GL_VENDOR_NVIDIA && defined MC_OS_WINDOWS
    ivec2 texoffsets0[4];
    texoffsets0[0] = ivec2(vec2(0.321630, 0.237497) * scale);
    texoffsets0[1] = ivec2(vec2(0.170197, 0.575310) * scale);
    texoffsets0[2] = ivec2(vec2(0.226930, 0.767481) * scale);
    texoffsets0[3] = ivec2(vec2(1.000000, 0.000000) * scale);

    ivec2 texoffsets1[4];
    texoffsets1[0] = -texoffsets0[0];
    texoffsets1[1] = -texoffsets0[1];
    texoffsets1[2] = -texoffsets0[2];
    texoffsets1[3] = -texoffsets0[3];

    vec4 depthgather0 = ld4(textureGatherOffsets(depthtex0, gl_FragCoord.xy / viewResolution, texoffsets0, 0));
    vec4 depthgather1 = ld4(textureGatherOffsets(depthtex0, gl_FragCoord.xy / viewResolution, texoffsets1, 0));

#else

    ivec2 texelcoord = ivec2(gl_FragCoord.xy);

    ivec2 texelCoords0 = ivec2(texelcoord + vec2(0.321630, 0.237497) * scale);
    ivec2 texelCoords1 = ivec2(texelcoord + vec2(0.170197, 0.575310) * scale);
    ivec2 texelCoords2 = ivec2(texelcoord + vec2(0.226930, 0.767481) * scale);
    ivec2 texelCoords3 = ivec2(texelcoord + vec2(1.000000, 0.000000) * scale);

    ivec2 texelCoordsx = ivec2(texelcoord - vec2(0.321630, 0.237497) * scale);
    ivec2 texelCoordsy = ivec2(texelcoord - vec2(0.170197, 0.575310) * scale);
    ivec2 texelCoordsz = ivec2(texelcoord - vec2(0.226930, 0.767481) * scale);
    ivec2 texelCoordsw = ivec2(texelcoord - vec2(1.000000, 0.000000) * scale);

    vec4 depthgather0 = ld4(vec4(texelFetch(depthtex0, texelCoords0, 0).r, texelFetch(depthtex0, texelCoords1, 0).r,
                                 texelFetch(depthtex0, texelCoords2, 0).r, texelFetch(depthtex0, texelCoords3, 0).r));
    vec4 depthgather1 = ld4(vec4(texelFetch(depthtex0, texelCoordsx, 0).r, texelFetch(depthtex0, texelCoordsy, 0).r,
                                 texelFetch(depthtex0, texelCoordsz, 0).r, texelFetch(depthtex0, texelCoordsw, 0).r));

#endif

    vec4 sample0 = (depth - (depthgather0)) * mult / steps;
    vec4 sample1 = (depth - (depthgather1)) * mult / steps;

    vec4 angle;
    angle += max(0.5 - sample0, 0.0);
    angle += max(0.5 - sample1, 0.0);

    vec4 dist;
    dist += max(0.25 * sample0 - 1.0, 0.0);
    dist += max(0.25 * sample1 - 1.0, 0.0);

    vec4 ao2 = clamp(angle + dist, 0.0, 1.0);
    ao = dot(ao2, vec4(1.0)) * 0.25;

    return ao;
}


#include "/lib/rsm.glsl"

vec3 get_GI(vec3 normal, float noise, vec3 p3)
{

    vec3 worldposition = p3;
    worldposition = mat3(shadowModelView) * p3 + shadowModelView[3].xyz;
    vec3 wpos2 = worldposition.xyz;

    worldposition = diagonal3(shadowProjection) * worldposition + shadowProjection[3].xyz;

    float distortFac = calcDistort(worldposition.xy);
    worldposition.xy *= distortFac;

    worldposition = worldposition * vec3(0.5, 0.5, 0.5 / 6.0) + 0.5;

    vec3 t1 = mat3(shadowModelView) * normal;

    vec3 projNormal = t1.rgb;

    vec3 gi = vec3(0.0);

    noise = noise * 2.0 * PI;

    const int nbSamples = 8;

    float sampleRadius = 0.15;
    for (int i = 0; i < nbSamples; i++)
    {

        vec2 offset = tapLocation(i, noise, nbSamples, 90) * 0.1;
        float weight = length(offset);

        // vec2 TC = ((worldposition.xy + offset * sampleRadius) * 2.0 - 1.0) / 1.8 * 0.5 + 0.5;
        vec2 TC = ((worldposition.xy + offset * sampleRadius));
        vec2 TCScaled = TC * shadowMapResolution;

        vec4 shadowSample =
            vec4(texelFetch2D(shadowcolor0, ivec2(TCScaled), 0).rgb, texelFetch2D(shadowtex0, ivec2(TCScaled), 0).x);
        float sampleDepth = shadowSample.a + 1 / 2000.0;
        vec3 shadowNormal = texelFetch2D(shadowcolor1, ivec2(TCScaled), 0).xyz * 2.0 - 1.0;
        vec3 sAlbedo = pow(shadowSample.rgb, vec3(2.2));

        vec3 S1 = vec3((TC * 2.0 - 1.0) * 1.8 * 0.5 + 0.5, sampleDepth);
        vec4 shadowPos = shadowProjectionInverse * ((vec4(S1, 1.0) * 2.0 - 1.0) * vec4(1.0, 1.0, 7.0, 1.0));
        vec3 V = shadowPos.xyz - wpos2;

        float VdotV = dot(V, V);
        float NdotS = clamp(dot(projNormal, V * inversesqrt(VdotV)), 0.0, 1.0) * 0.8 + 0.2;
        float inHemisphere = sqrt(clamp(-dot(shadowNormal, V * inversesqrt(VdotV)), 0.0, 1.0)) * 0.8 + 0.25;
        float Scoef = 1.0 / (1.0 + VdotV);

        gi += inHemisphere * NdotS * Scoef * pow(weight, 1.5) * sAlbedo;
        //  gi = vec3(Scoef);
    }

    // Scale by (Total_area/Area_covered) to compensate for missing missing sample
    // acov = n*dS
    // atot = r*r*pi*dS
    // dS = 1/smres*1/smres*shadowDist*shadowDist
    // shadowDist = approx 576
    // dS = approx 576/512^2
    // r = sampleRadius * shadowDist
    // r*r*pi/n
    const float dS = pow(576 / 512.0, 2.0);
    // return gi;
    return gi * (sampleRadius * sampleRadius * 512.0 * 512.0 / PI / nbSamples) * dS;
}


vec2 calculateSunTexCoord(vec3 spherePoint, vec3 normalizedSunPos, float sunScale)
{
    // Project the spherePoint onto the plane orthogonal to the sun's position to get the relative position
    vec3 rightVec = normalize(cross(normalizedSunPos, vec3(0, 1, 0))); // Assuming Y is up
    vec3 upVec = cross(rightVec, normalizedSunPos);
    // Calculate horizontal and vertical "visibility" factors
    float horizontalFactor = dot(spherePoint, rightVec);
    float verticalFactor = dot(spherePoint, upVec);

    bool isVisible = abs(horizontalFactor) < 0.17 && abs(verticalFactor) < 0.17;

    float frontBackFactor = dot(normalizedSunPos, spherePoint);
    vec2 relativePos = vec2(0, 0);

    if (isVisible && frontBackFactor > 0.972)
    {
        // Calculate the relative position of the spherePoint in the sun's direction plane

        relativePos.x = dot(spherePoint, rightVec) * sunScale + 0.5; // Scale and bias to move to texture space
        relativePos.y = dot(spherePoint, upVec) * sunScale + 0.5;

        // Clamp the coordinates to [0, 1] to avoid wrapping
        relativePos = clamp(relativePos, 0.0, 1.0);
    }

    return relativePos;
}

vec2 calculateMoonTexCoord(vec3 spherePoint, vec3 normalizedMoonPos, float moonScale, int moonPhase)
{
    const int phasesPerRow = 4;                                         // Number of moon phases per row in the texture
    const int totalRows = 2;                                            // Total number of rows
    vec3 rightVec = normalize(cross(normalizedMoonPos, vec3(0, 1, 0))); // Horizontal direction
    vec3 upVec = cross(rightVec, normalizedMoonPos);                    // Vertical direction

    float horizontalFactor = dot(spherePoint, rightVec);
    float verticalFactor = dot(spherePoint, upVec);

    bool isVisible = abs(horizontalFactor) < 0.17 && abs(verticalFactor) < 0.17;
    float frontBackFactor = dot(normalizedMoonPos, spherePoint);
    vec2 relativePos = vec2(0, 0);

    if (isVisible && frontBackFactor > 0.972)
    {
        float phaseWidth = 1.0 / float(phasesPerRow); // Width of each phase in texture coordinates
        float phaseHeight = 1.0 / float(totalRows);   // Height of each phase in texture coordinates

        // Calculate the texture coordinates as if the entire texture is just one phase
        vec2 baseTexCoord;
        baseTexCoord.x = dot(spherePoint, rightVec) * moonScale + 0.5;
        baseTexCoord.y = dot(spherePoint, upVec) * moonScale + 0.5;

        // Determine the row and column of the current phase
        int row = moonPhase / phasesPerRow;
        int column = moonPhase % phasesPerRow;

        // Now adjust for the phase width, phase height, and the phase offset
        relativePos.x = (baseTexCoord.x / float(phasesPerRow)) + (float(column) * phaseWidth);
        relativePos.y = (baseTexCoord.y / float(totalRows)) + (float(row) * phaseHeight);

        // Clamp the coordinates to the boundaries of the current phase section to avoid wrapping
        relativePos.x = clamp(relativePos.x, float(column) * phaseWidth, (float(column) + 1.0) * phaseWidth);
        relativePos.y = clamp(relativePos.y, float(row) * phaseHeight, (float(row) + 1.0) * phaseHeight);
    }
    return relativePos;
}



///

/*
vec2 SinCos(float x)
{
    return vec2(sin(x), cos(x));
}

vec3 GenerateUnitVector(vec2 xy)
{
    xy.x *= radians(360.0);
    xy.y = xy.y * 2.0 - 1.0;
    return vec3(SinCos(xy.x) * sqrt(1.0 - xy.y * xy.y), xy.y);
}
    const ivec2 samples = ivec2(64, 32);

    ambientRight = vec3(0.0);
    ambientUp = vec3(0.0);
    ambientF = vec3(0.0);

    ambientLeft = vec3(0.0);
    ambientDown = vec3(0.0);
    ambientB = vec3(0.0);

    for (int x = 0; x < samples.x; ++x)
    {
        for (int y = 0; y < samples.y; ++y)
        {
            vec3 dir = GenerateUnitVector((vec2(x, y) + 0.5) / samples);
            if (dir.y < 0.0)
            {
                dir.y = -dir.y;
            }

            vec3 skySample = skyFromTex2(dir, colortex6).rgb;

            ambientRight += skySample * clamp(dir.x, 0.0, 1.0);
            ambientUp += skySample * clamp(dir.y, 0.0, 1.0);
            ambientF += skySample * clamp(dir.z, 0.0, 1.0);

            ambientLeft += skySample * clamp(-dir.x, 0.0, 1.0);
            ambientB += skySample * clamp(-dir.z, 0.0, 1.0);
        }
    }

    const float sampleWeight = 2.0 / (samples.x * samples.y);
    ambientRight *= sampleWeight;
    ambientUp *= sampleWeight;
    ambientF *= sampleWeight;

    ambientLeft *= sampleWeight;
    ambientB *= sampleWeight;

    // super simple fake skylight bounce
    const float fakeBounceAlbedo = 0.2;
    ambientRight += ambientUp * fakeBounceAlbedo * 0.5;
    ambientF += ambientUp * fakeBounceAlbedo * 0.5;

    ambientLeft += ambientUp * fakeBounceAlbedo * 0.5;
    ambientDown += ambientUp * fakeBounceAlbedo;
    ambientB += ambientUp * fakeBounceAlbedo * 0.5;

*/