#if !defined ATMOS_LOOKUP_LIB
#define ATMOS_LOOKUP_LIB
    bool atmosphereRayIntersectsLowerLimit(float R, float Mu) {
        return Mu < 0.0f && R * R * (Mu * Mu - 1.0f) + atmosphereLowerLimitSquared >= 0.0f;
    }

    float atmosphereDistanceToUpperLimit(float R, float Mu) {
        float discriminant = R * R * (Mu * Mu - 1.0f) + atmosphereRadiusSquared;
        return -R * Mu + sqrt(discriminant);
    }
    float atmosphereDistanceToLowerLimit(float R, float Mu) {
        float discriminant = R * R * (Mu * Mu - 1.0f) + atmosphereLowerLimitSquared;
        return -R * Mu - sqrt(discriminant);
    }

    float addUvMargin(float uv, int resolution) {
        return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution);
    }

    vec2 addUvMargin(vec2 uv, ivec2 resolution) {
        return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution);
    }

    vec3 addUvMargin(vec3 uv, ivec3 resolution) {
        return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution);
    }

    vec4 addUvMargin(vec4 uv, ivec4 resolution) {
        return uv * (1.0 - 1.0 / resolution) + (0.5 / resolution);
    }

    vec2 atmosphereTransmittanceLookupUv(float R, float Mu) {
        const float H = sqrt(atmosphereRadiusSquared - atmosphereLowerLimitSquared);

        float rho = sqrt(R * R - atmosphereLowerLimitSquared);

        float d = atmosphereDistanceToUpperLimit(R, Mu);
        float dMin = atmosphereRadius - R;
        float dMax = rho + H;

        float uvMu = (d - dMin) / (dMax - dMin);
        float uvR = rho / H;

        return vec2(uvMu, uvR);
    }

    vec4 atmosphereScatteringLookupUv(float R, float Mu, float MuS, float V) {
        const float H = sqrt(atmosphereRadiusSquared - atmosphereLowerLimitSquared);

        float rho = sqrt(R * R - atmosphereLowerLimitSquared);
        float uvR = addUvMargin(rho / H, resR);

        float RMu = R * Mu;
        float discriminant = RMu * RMu - R * R + atmosphereLowerLimitSquared;
        float uvMu = 0.5;
        if (atmosphereRayIntersectsLowerLimit(R, Mu)) {
            float d = -RMu - sqrt(discriminant);
            float dMin = R - atmosphereLowerLimit;
            float dMax = rho;
            uvMu -= 0.5 * addUvMargin(dMax == dMin ? 0.0 : (d - dMin) / (dMax - dMin), resMu / 2);
        } else {
            float d = -RMu + sqrt(discriminant + H * H);
            float dMin = atmosphereRadius - R;
            float dMax = rho + H;
            uvMu += 0.5 * (addUvMargin((d - dMin) / (dMax - dMin), resMu / 2) * (1.0 - 3.0/resMu));
        }

        float d = atmosphereDistanceToUpperLimit(atmosphereLowerLimit, MuS);
        float dMin = atmosphereRadius - atmosphereLowerLimit;
        float dMax = H;
        float a = (d - dMin) / (dMax - dMin);
        float A = -2.0 * atmosphere_MuS_min * atmosphereLowerLimit / (dMax - dMin);
        float uvMuS = addUvMargin(max(1.0 - a / A, 0.0) / (1.0 + a), resMuS);

        float uvV = addUvMargin((V + 1.0) / 2.0, resV);

        return vec4(uvV, uvMuS, uvMu, uvR);
    }

    vec2 atmosphere4DLookupUv(vec4 coord) {
        const ivec4 resolution = ivec4(resV, resMuS, resMu, resR);

        vec2 xy = coord.xy * resolution.xy;
        ivec2 zw = ivec2(floor(coord.zw * resolution.zw)) * resolution.xy;
        return xy + zw;
    }
#endif