// PAOLO'S CHANGES {
double Fresnel( double phi, double n )
{
	// calculates fresnel coefficient for a surface with an 
	// index of refraction of n, and an angle of phi between
	// m (normal at hit pt) and s (dir from hit pt to light).
	// the equation is given by:

	// F = 1/2 * (g-c)^2 / (g+c)^2 * { 1 + [ (c(g+c)-1) / (c(g-c)+1) ]^2 }
	// where c = cos(phi), g = sqrt( n^2 + c^2 - 1 )

	double F; // the fresnel coefficient
	double c; // cos(phi)
	double g; // sqrt( n^2 + c^2 - 1 )
	double temp;

	c = cos( phi );
	g = sqrt( (n*n) + (c*c) - 1 );
	if(g<0) g=0;

	F = 0.5;
	F *= ((g-c)*(g-c))/((g+c)*(g+c));
	temp = ((c * (g+c)) - 1) / ((c * (g-c)) + 1);
	temp *= temp;
	temp += 1;
	F *= temp;
	return F;
}

double Specular( Vector3& m, Vector3& s, Vector3& v, Material& surface )
{
	// this function computes ps for the Cook-Torrence lighting model
	// m = surface normal
	// s = unit vector pointing at light source 
	// v = unit vector pointing at camera
	// h = vector 'half-way' between s and v
	// ps = (Fresnel(phi, surface refraction)) / PI * (D * G / ((m.s)(m.v)))
	// phi = angle between m and s
	// sigma = angle between h and m
	// D = (1 / (4 * roughness^2 * cos(sigma)^4)) * e^-(tan(sigma)/roughness)^2
	// G = min(1, Gm, Gs)
	// Gm = 2(m.h)(m.s)/(h.s)
	// Gs = 2(m.h)(m.v)/(h.s)

	Vector3 h;	// half-way vector
	double ps;	// this is what we are computing
	double phi;	// angle between m and s;
	double sigma;	// angle  between h and m
	double D;	// Beckmann distribution
	double G;	// geometric term
	double Gm, Gs; // candidate geometric terms
	double temp;	// temporary variable

	m.normalize(); s.normalize(); v.normalize();

	// compute h
	// NOTE: we assume that s and v are normalized!
	h.set( s.x + v.x, s.y + v.y, s.z + v.z );
	h.normalize();

	// compute phi
	// NOTE: we assume that m and s are normalized!
	phi = acos( m.dot(s) );

	// compute sigma
	// NOTE: we assume that m is normalized!
	sigma = acos( h.dot(m) );

	// compute D
	temp = tan(sigma) / surface.roughness;
	temp *= temp;
	D = 1.0 / cos(sigma);
	D *= D; // (1/cos^2)
	D *= D; // (1/cos^4)
	D /= (4.0 * surface.roughness * surface.roughness);
	D *= exp( -temp );

	// compute G
	G = 1.0;
	Gm = 2.0 * m.dot(h) * m.dot(s) / h.dot(s);
	Gs = 2.0 * m.dot(h) * m.dot(v) / h.dot(s);
	if(Gm < G) G = Gm;
	if(Gs < G) G = Gs;

	// compute the specular reflection
	temp = (D * G / ( m.dot(s) * m.dot(v) ) );
	ps = (Fresnel( phi, surface.indexOfRefraction )) / PI * temp;

	return ps;
}

// } PAOLO'S CHANGES