/***********************************************/
/*           Copyright (c) 2025 Belmu          */
/*             All Rights Reserved             */
/***********************************************/

const float halfSector = bladeSectorAngle * 0.5;

bool isInsideAperture(vec2 coords, float aperture) {
    coords /= aperture;

    float angle = atan(coords.y, coords.x);
          angle = mod(angle + halfSector, bladeSectorAngle) - halfSector;

    return length(coords) <= cos(halfSector) / cos(angle);
}

bool intersectApertureStop(inout vec3 position, vec3 direction, float opticalAxisZ, LensElement element, out vec3 normal) {
    // Ray-Plane intersection
    float t = (opticalAxisZ - position.z) / direction.z;

    position += direction * t;

    if (element.isApertureStop && !isInsideAperture(position.xy, element.aperture))
        return false;

    normal = vec3(0.0, 0.0, -sign(direction.z));

    return t >= 0.0;
}

bool intersectSphericalLensElement(inout vec3 position, vec3 direction, float opticalAxisZ, LensElement element, out vec3 normal) {
    // Ray-Sphere intersection
    vec3 center = vec3(0.0, 0.0, opticalAxisZ + element.radius);
    vec2 sphere = intersectSphere(position - center, direction, element.radius);

    float t = direction.z * element.radius > 0.0 ? sphere.x : sphere.y;

    position += direction * t;

    if (length(position.xy) > element.aperture)
        return false;
    
    normal  = normalize(position - center);
    normal *= -sign(dot(direction, normal));

    return t >= 0.0;
}

float getSurfaceSagitta(vec3 position, LensElement element) {
    float curvature = 1.0 / element.radius;
    float height    = length(position.xy);
    float heightSq  = height * height;

    float sagitta = heightSq * curvature / (1.0 + sqrt(1.0 - heightSq * curvature * curvature));

    for (int i = 0; i < 6; i++) {
        sagitta += element.asphericalCoefficients[i] * pow(height, 2 * (i + 2));
    }

    return sagitta;
}

vec3 getSurfaceSagittaNormal(vec3 position, LensElement element) {
    float epsilon = 1e-6;
    vec3 dx = vec3(epsilon, 0.0, 0.0);
    vec3 dy = vec3(0.0, epsilon, 0.0);

    float dzdx = (getSurfaceSagitta(position + dx, element) - getSurfaceSagitta(position - dx, element)) / (2.0 * epsilon);
    float dzdy = (getSurfaceSagitta(position + dy, element) - getSurfaceSagitta(position - dy, element)) / (2.0 * epsilon);

    return normalize(vec3(-dzdx, -dzdy, 1.0));
}

bool intersectAsphericalLensElement(inout vec3 position, vec3 direction, float opticalAxisZ, float prevOpticalAxisZ, LensElement element, out vec3 normal) {
    vec3 center = vec3(0.0, 0.0, opticalAxisZ + element.radius);
    vec2 sphere = intersectSphere(position - center, direction, element.radius);

    float upperBound = direction.z * element.radius > 0.0 ? sphere.x : sphere.y;
    float lowerBound = ((prevOpticalAxisZ - position.z) / direction.z) * sign(direction.z);

    if (upperBound < 0.0 || lowerBound < 0.0) return false;

    vec3 origin = position;

    for (int i = 0; i < 128; i++) {
        float midPoint = (lowerBound + upperBound) * 0.5;

        position = origin + direction * midPoint;

        float sagitta  = getSurfaceSagitta(position, element);
        float surfaceZ = opticalAxisZ + sagitta;

        if (position.z - surfaceZ > 0.0)
            upperBound = midPoint;
        else
            lowerBound = midPoint;

        if (abs(upperBound - lowerBound) < 1e-7) {
            normal  = getSurfaceSagittaNormal(position, element);
            normal *= -sign(dot(direction, normal));
            return true;
        }
    }

    return false;
}
