float Square(in float x) { return x * x; }

float Theta(in vec3 w)
{
    return acos(w.z / length(w));
}

float CosTheta(in vec3 w)
{
    return cos(Theta(w));
}

float EricHeitz2018GGXG1Lambda(in vec3 V, in float alpha_x, in float alpha_y)
{
    float Vx2 = Square(V.x);
    float Vy2 = Square(V.y);
    float Vz2 = Square(V.z);
    float ax2 = Square(alpha_x);
    float ay2 = Square(alpha_y);
    return (-1.0 + sqrt(1.0 + (Vx2 * ax2 + Vy2 * ay2) / Vz2)) / 2.0;
}

float EricHeitz2018GGXG1(vec3 V, float alpha_x, float alpha_y)
{
    return 1.0 / (1.0 + EricHeitz2018GGXG1Lambda(V, alpha_x, alpha_y));
}

// wm: microfacet normal in frame
float EricHeitz2018GGXD(in vec3 N, in float alpha_x, in float alpha_y)
{
    float Nx2 = Square(N.x);
    float Ny2 = Square(N.y);
    float Nz2 = Square(N.z);
    float ax2 = Square(alpha_x);
    float ay2 = Square(alpha_y);
    return 1.0 / (PI * alpha_x * alpha_y * Square(Nx2 / ax2 + Ny2 / ay2 + Nz2));
}

float EricHeitz2018GGXG2(in vec3 V, in vec3 L, in float alpha_x, in float alpha_y)
{
    return EricHeitz2018GGXG1(V, alpha_x, alpha_y) * EricHeitz2018GGXG1(L, alpha_x, alpha_y);
}

vec3 EricHeitz2018GGX(in vec3 V, in vec3 L, in float roughness, in float anisotropic, in float ior)
{
    float alpha = roughness * roughness;
    float aspect = sqrt(1.0 - 0.9 * anisotropic);
    float alpha_x = alpha * aspect;
    float alpha_y = alpha / aspect;

    vec3 H = normalize(L + V);
    float NoV = CosTheta(V);
    float NoL = CosTheta(L);
    if (NoV < 0.0 || NoL < 0.0) return vec3(0);

    float VoH = dot(V, H);
    float NoH = CosTheta(H);
    vec3 F0 = vec3(abs((1.0 - ior) / (1.0 + ior)));
    F0 = F0 * F0;

    float D = EricHeitz2018GGXD(H, alpha_x, alpha_y);
    vec3 F = SchlickFresnel(max(NoV, 0.0), F0);
    float G = EricHeitz2018GGXG2(V, L, alpha_x, alpha_y);
    return (F * D * G) / (4.0f * VoH * NoH);
}