highp float GetNoHSquared(highp float radiusTan, highp float NoL, highp float NoV, highp float VoL){
highp float radiusCos = 1.0 / sqrt(1.0 + radiusTan * radiusTan);  
highp float RoL = 2.0 * NoL * NoV - VoL;
if(RoL >= radiusCos) return 1.0;
highp float rOverLengthT = radiusCos * radiusTan / sqrt(1.0 - RoL * RoL);
highp float NoTr = rOverLengthT * (NoV - RoL * NoL);
highp float VoTr = rOverLengthT * (2.0 * NoV * NoV - 1.0 - RoL * VoL);
highp float triple = sqrt(clamp(1.0 - NoL * NoL - NoV * NoV - VoL * VoL + 2.0 * NoL * NoV * VoL, 0.0, 1.0));
highp float NoBr = rOverLengthT * triple, VoBr = rOverLengthT * (2.0 * triple * NoV);
highp float NoLVTr = NoL * radiusCos + NoV + NoTr, VoLVTr = VoL * radiusCos + 1.0 + VoTr;
highp float p = NoBr * VoLVTr, q = NoLVTr * VoLVTr, s = VoBr * NoLVTr;
highp float xNum = q * (-0.5 * p + 0.25 * VoBr * NoLVTr);
highp float xDenom = p * p + s * ((s - 2.0 * p)) + NoLVTr * ((NoL * radiusCos + NoV) * VoLVTr * VoLVTr + q * (-0.5 * (VoLVTr + VoL * radiusCos) - 0.5));
highp float twoX1 = 2.0 * xNum / (xDenom * xDenom + xNum * xNum);
highp float sinTheta = twoX1 * xDenom;
highp float cosTheta = 1.0 - twoX1 * xNum;
NoTr = cosTheta * NoTr + sinTheta * NoBr;
VoTr = cosTheta * VoTr + sinTheta * VoBr;
highp float newNoL = NoL * radiusCos + NoTr;
highp float newVoL = VoL * radiusCos + VoTr;
highp float NoH = NoV + newNoL;
highp float HoH = 2.0 * newVoL + 2.0;
return max(NoH * NoH / HoH, 0.0);}
highp float GGX(highp vec3 V, highp vec3 L, highp vec3 N, highp float r, highp float f0, float sunSize){

if(r < 0.05) r = 0.05;
highp float a2 = pow(r, 4.0);

highp vec3 H = normalize(L + V);

highp float HoL = clamp(dot(H, L), 0.0, 1.0);

highp float NoL = clamp(dot(N, L), 0.0, 1.0);
highp float NoH = GetNoHSquared(sunSize, NoL, dot(N, V), dot(V, L));

highp float F = (1.0 - f0) * pow(1.0 - HoL, 5.0) + f0;

highp float denom = pow(NoH, 2.0) * (a2 - 1.0) + 1.0;

const highp float pi = 3.14159265358979;

return NoL * a2 / (pi * pow(denom, 2.0)) * F / (pow(HoL, 2.0) * (1.0 - a2) + a2);
}

highp float Fresnel(highp vec3 N, highp vec3 V, vec3 spec){
highp float NdotV = 1.0 - saturate(dot(N, normalize(V)));
return mix((1.0 - spec.r) * pow(NdotV,3.0), 1.0, spec.g);}