float P22(in vec2 slope, in vec2 alpha) {
    if (any(equal(alpha, vec2(0.0)))) { return 0.0; }
	float tmp = 1.0 + square(slope.x)/square(alpha.x) + square(slope.y)/square(alpha.y);
	float value = 1.0 / (pi * alpha.x * alpha.y) / (tmp * tmp);
	return value;
}

float smith_Lambda(in vec3 wi, in vec2 alpha) {
    if(wi.z > 0.9990) {
        return 0.0;
    }
    if(wi.z < -0.9990) {
        return -1.0;
    }

    float theta_I = acos(wi.z);
    float a = 1.0 / tan(theta_I) / alpha_I(wi, alpha);

    float value = 0.5 * (-1.0 + sign(a) * sqrt(1.0 + 1.0/square(a)));
    
    return value;
}

float projectedArea(in vec3 wi, in vec2 alpha) {
    if(wi.z > 0.9990) {
        return 1.0;
    }
    if(wi.z < -0.9990) {
        return 0.0;
    }

    float theta_I = acos(wi.z);
    float sinTheta_I = sin(theta_I);

    float alphaI = alpha_I(wi, alpha);

    float value = 0.5 * (wi.z + sqrt(square(wi.z) + square(sinTheta_I) * square(alphaI)));

    return value;
}

vec2 sampleP22_11(in float theta_I, in vec2 U) {
    vec2 slope;

    if(theta_I < 0.0001) {
        float r = sqrt(U.x/(1.0-U.x));
        float phi = tau * U.x;
        slope = vec2(r * cos(phi), r * sin(phi));

        return slope;
    }

    float sinTheta_I = sin(theta_I);
    float cosTheta_I = cos(theta_I);
    float tanTheta_I = tan(theta_I);

    float slope_I = cosTheta_I/sinTheta_I;

    float projectedArea = 0.5 * (cosTheta_I + 1.0);
    if(projectedArea < 0.0001 || projectedArea != projectedArea) {
        return vec2(0.0);
    }

    float c = 1.0 / projectedArea;

    float A = 2.0 * U.x / cosTheta_I / c - 1.0;
    float B = tanTheta_I;
    float tmp = 1.0 / (square(A) - 1.0);

    float D = sqrt(max(0.0, square(B) * square(tmp) - (square(A) - square(B)) * tmp));
    vec2 slopeX = vec2(B*tmp - D, B*tmp + D);
    slope.x = (A < 0.0 || slopeX.y > 1.0 / tanTheta_I) ? slopeX.x : slopeX.y;

    float U2;
    float S;

    if(U.y > 0.5) {
        S = 1.0;
        U2 = 2.0 * (U.y - 0.5);
    } else {
        S = -1.0;
        U2 = 2.0  * (0.5 - U.y);
    }
    float z = (U2 * (U2 * (U2 * 0.27385 - 0.73369)+0.46341)) / (U2 * (U2 * (U2 * 0.093073+0.309420) - 1.000000) + 0.597999);
    slope.y = S * z * sqrt(1.0 + square(slope.x));

    return slope;
}