// most of the following code are wrapped from PBRT: physically-based rendering. 
// I copied them here just to make the evaluator simpler. 

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;


inline float Clamp(float val, float low, float high) 
{
	if (val < low) return low;
	else if (val > high) return high;
	else return val;
}

inline float FrDiel(float cosi, float cost, const float &etai, const float &etat) 
{
		float Rparl = ((etat * cosi) - (etai * cost)) /
						((etat * cosi) + (etai * cost));
		float Rperp = ((etai * cosi) - (etat * cost)) /
						((etai * cosi) + (etat * cost));
		return (Rparl*Rparl + Rperp*Rperp) / 2.f;
}


class FresnelDielectric{
public:
	// FresnelDielectric Public Methods
	FresnelDielectric(float ei, float et) {
		eta_i = ei;
		eta_t = et;
	}

	float Evaluate(float cosi) const
	{
		// Compute Fresnel reflectance for dielectric
		cosi = Clamp(cosi, -1.f, 1.f);
		// Compute indices of refraction for dielectric
		bool entering = cosi > 0.;
		float ei = eta_i, et = eta_t;
		if (!entering) std::swap(ei, et);		
		float eta = ei / et;
		float sint2 = eta * eta * std::max(0.f, 1.f-cosi*cosi);
		if (sint2 > 1.) {
			// Handle total internal reflection
			return 1.;
		}
		else {
			float cost = sqrtf(std::max(0.f, 1.f-sint2));
			return FrDiel(fabsf(cosi), cost, ei, et);
		}
	}
public:
	// FresnelDielectric Private Data
	float eta_i, eta_t;
};


// Geometry Declarations
class Vector {
public:
	// Vector Public Methods
	Vector(float _x=0, float _y=0, float _z=0)
		: x(_x), y(_y), z(_z) {
	}
	Vector operator+(const Vector &v) const {
		return Vector(x + v.x, y + v.y, z + v.z);
	}
	
	Vector& operator+=(const Vector &v) {
		x += v.x; y += v.y; z += v.z;
		return *this;
	}
	Vector operator-(const Vector &v) const {
		return Vector(x - v.x, y - v.y, z - v.z);
	}
	
	Vector& operator-=(const Vector &v) {
		x -= v.x; y -= v.y; z -= v.z;
		return *this;
	}
	bool operator==(const Vector &v) const {
		return x == v.x && y == v.y && z == v.z;
	}
	Vector operator*(float f) const {
		return Vector(f*x, f*y, f*z);
	}
	
	Vector &operator*=(float f) {
		x *= f; y *= f; z *= f;
		return *this;
	}
	Vector operator/(float f) const {
		assert(f!=0);
		float inv = 1.f / f;
		return Vector(x * inv, y * inv, z * inv);
	}
	
	Vector &operator/=(float f) {
		assert(f!=0);
		float inv = 1.f / f;
		x *= inv; y *= inv; z *= inv;
		return *this;
	}
	Vector operator-() const {
		return Vector(-x, -y, -z);
	}
	float operator[](int i) const {
		assert(i >= 0 && i <= 2);
		return (&x)[i];
	}
	
	float &operator[](int i) {
		assert(i >= 0 && i <= 2);
		return (&x)[i];
	}
	// Vector Public Data
	float x, y, z;
};
inline ostream &operator<<(ostream &os, const Vector &v) {
	os << v.x << ", " << v.y << ", " << v.z;
	return os;
}
inline Vector operator*(float f, const Vector &v) {
	return v*f;
}
inline float Dot(const Vector &v1, const Vector &v2) {
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

inline float CosTheta(const Vector &w) { return w.z; }
inline float SinTheta(const Vector &w) {
	return sqrtf(std::max(0.f, 1.f - w.z*w.z));
}
inline float SinTheta2(const Vector &w) {
	return 1.f - CosTheta(w)*CosTheta(w);
}

inline float PhaseHG(const Vector &w, const Vector &wp, float g) {
	float costheta = Dot(w, wp);
	return 1.f / (4.f * M_PI) * (1.f - g*g) /
		powf(1.f + g*g - 2.f * g * costheta, 1.5f);
}