
vec3 calculate_oren_nayar_model(in vec3 n, in vec3 l, in vec3 e, in float roughness, in float specularPower)
{
	float EdotN = dot(e, n);
	float LdotN = dot(l, n);
	vec3 hv = normalize(l + e);

	float gamma = max( dot(e - n*dot(e, n), l - n*dot(l, n)), 0.0);
	float rough_sq = roughness * roughness;
	float A = 1.0 - 0.5 * (rough_sq / (rough_sq + 0.57));
	float B = 0.45 * (rough_sq / (rough_sq + 0.09));

	float acosEdotN = acos(EdotN);
	float acosLdotN = acos(LdotN);
	float alpha = max(acosEdotN, acosLdotN);
	float beta  = min(acosEdotN, acosLdotN);
	float C = sin(alpha) * tan(beta);
	
	float diffTerm = max(LdotN, 0.0)*(A + B*C);

	float specTerm = pow(max(dot(n, hv), 0.0), specularPower)*diffTerm;

	return vec3(diffTerm, specTerm, 1.0);
}


vec3 calculate_oren_nayar_model_abs(in vec3 n, in vec3 l, in vec3 e, in float roughness, in float specularPower)
{
	float EdotN = dot(e, n);
	float LdotN = abs(dot(l, n));
	vec3 hv = normalize(l + e);

	float gamma = max( dot(e - n*dot(e, n), l - n*dot(l, n)), 0.0);
	float rough_sq = roughness * roughness;
	float A = 1.0 - 0.5 * (rough_sq / (rough_sq + 0.57));
	float B = 0.45 * (rough_sq / (rough_sq + 0.09));

	float acosEdotN = acos(EdotN);
	float acosLdotN = acos(LdotN);
	float alpha = max(acosEdotN, acosLdotN);
	float beta  = min(acosEdotN, acosLdotN);
	float C = sin(alpha) * tan(beta);
	
	float diffTerm = max(LdotN, 0.0)*(A + B*C);

	float specTerm = pow(max(dot(n, hv), 0.0), specularPower)*diffTerm;

	return vec3(diffTerm, specTerm, 1.0);
}
