#include <GL/glut.h>
#include <iostream.h>
#include <fstream.h>
#include <math.h>

const double pi = 3.1415926535;

double look_az=0, look_el=0;
double look_zoom=100;
double mirror_az=pi/2, mirror_el=pi/4;
int res = 100;	// Num triangles around equator of polar function.
int intensityScale = 40;

double param = 1, cone = pi/6;

double coeff[] = {0, 1, 0, 0};

double TestBRDF(double view_az, double view_el)
{	// Global: mirror_az, mirror_el.

	double viewX = -cos(view_el)*sin(view_az);
	double viewY = sin(view_el);
	double viewZ = cos(view_el)*cos(view_az);
	double sourceX = -cos(mirror_el)*sin(mirror_az);
	double sourceY = sin(mirror_el);
	double sourceZ = cos(mirror_el)*cos(mirror_az);
	double diffX = viewX-sourceX;
	double diffY = viewY-sourceY;
	double diffZ = viewZ-sourceZ;
	double dist = sqrt(diffX*diffX+diffY*diffY+diffZ*diffZ);
	return 1*(0.1 + exp(-dist));
}

inline
double Power(double x, int e)
{
	double p = 1;
	for (int i = 0; i < e; ++i)
		p *= x;
	return p;
}

class Vect
{
public:
	double v[3];
	Vect()
		{v[0]=v[1]=v[2]=0;}
	Vect(const Vect& V)
		{ v[0]=V.v[0]; v[1]=V.v[1]; v[2]=V.v[2]; }
	void operator += (const Vect& V)
		{ v[0]+=V.v[0]; v[1]+=V.v[1]; v[2]+=V.v[2]; }
	void operator -= (const Vect& V)
		{ v[0]-=V.v[0]; v[1]-=V.v[1]; v[2]-=V.v[2]; }
	void operator *= (const double c)
		{ v[0]*=c; v[1]*=c; v[2]*=c; }
	double Norm()
		{ return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); }
	double Dot(const Vect& V)
		{ return v[0]*V.v[0]+v[1]*V.v[1]+v[2]*V.v[2]; }
	Vect(double az, double el)
	{
		v[0] = -cos(el)*sin(az);
		v[1] = sin(el);
		v[2] = cos(el)*cos(az);
	}
	double Az()
	{
		double n = Norm();
		if (v[0]>0)
			return 2*pi-acos((v[2]/n)/cos(asin(v[1]/n)));
		else
			return acos((v[2]/n)/cos(asin(v[1]/n)));
		// +++++ This might be inaccurate for some values?
	}
	double El()
	{
		return asin(v[1]/Norm());
	}
};


double Phong(double view_az, double view_el)
{	// Global: mirror_az, mirror_el.

	const static double Iamb = 0.1;
	const static double Idiff = 0.5;
	const static double Ispec = 0.5;
	const static int specexp = 8;

	Vect view (view_az,view_el);
	Vect mirror (mirror_az,mirror_el);
	Vect source (mirror_az,pi-mirror_el);
	Vect norm (0,pi/2);
	double diffuse = norm.Dot(source);
	double cosalpha = view.Dot(mirror);
	if (cosalpha < 0)
		cosalpha = 0;

	return Iamb + Idiff * diffuse + Ispec*Power(cosalpha,specexp);
}

double Permute1(double view_az, double view_el, double (*brdf)(double,double))
{	// Global: mirror_az, mirror_el

	Vect view (view_az, view_el);
	Vect mirror (mirror_az, mirror_el);

	Vect diff = view;
	diff -= mirror;
	double R = diff.Norm();

	if (R < cone)
	{
		double newR = coeff[0] + coeff[1]*R + coeff[2]*R*R + coeff[3]*R*R*R;
		//newR = cone*pow(R/cone,param);

		Vect view2 = diff;
		view2 *= newR/R;
		view2 += mirror;

		return (*brdf)(view2.Az(), view2.El());
	}
	else
		return (*brdf)(view_az, view_el);
}

double Identity(double view_az, double view_el, double (*brdf)(double,double))
{
	return (*brdf)(view_az, view_el);
}

inline
void PolarVertex(double az, double el, double rad)
{
	rad *= intensityScale;
	double x = -cos(el)*sin(az);
	double y = sin(el);
	double z = cos(el)*cos(az);
	glNormal3f(x, y, z);
	x *= rad;
	y *= rad;
	z *= rad;
	glVertex3f(x, y, z);
}

void DrawBRDF(double (&permut)(double,double,double (*)(double,double)), double (&brdf)(double,double))
{
	for (int i = 0; i < res/4; ++i)
	{
		double el = double(i)/res * 2*pi;
		double elNext = el + double(1)/res * 2*pi;

		glBegin(GL_TRIANGLE_STRIP);

		for (int j = 0; j < res; ++j)
		{
			double az = double(j)/res * 2*pi;

			PolarVertex(az, el, permut(az,el, brdf));
			PolarVertex(az, elNext, permut(az,elNext, brdf));
		}

		PolarVertex(0,el, permut(0,el, brdf));
		PolarVertex(0,elNext, permut(0,elNext, brdf));
		glEnd();
	}
}

GLfloat	bright_light[] =	{1.0, 1.0, 1.0, 1.0},
		mid_light[] =		{0.5, 0.5, 0.5, 1.0},
		dim_light[] =		{0.1, 0.1, 0.1, 1.0},
		no_light[] =		{0.0, 0.0, 0.0, 1.0};

GLfloat	top_light_direction[] =	{0.0, 1.0, 0.0, 0.0},
		camera_light_position[] =	{0.0, 0.0, 0.0, 1.0};

GLfloat normal_mat[] = {10, 10, 10, 1},
		mirror_mat[] = {10, 0, 0, 1},
		surface_old_mat[] = {1, 1, 1, .5},
		surface_new_mat[] = {1, 0, 0, .5},
		dark_mat[] = {.1,.1,.1,.1};			// Unused


void Display()
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();

	glTranslatef(0,0,-look_zoom);
	glRotatef(look_el, 1, 0, 0);
	glRotatef(look_az, 0, 1, 0);

	glMaterialfv(GL_BACK, GL_AMBIENT, dark_mat);
	glMaterialfv(GL_BACK, GL_DIFFUSE, dark_mat);
	glMaterialfv(GL_BACK, GL_SPECULAR, dark_mat);

	glMaterialfv(GL_FRONT, GL_AMBIENT, surface_old_mat);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, surface_old_mat);
	glMaterialfv(GL_FRONT, GL_SPECULAR, surface_old_mat);
	glMaterialf(GL_FRONT, GL_SHININESS, 5);
	DrawBRDF(Identity, Phong);
	glMaterialfv(GL_FRONT, GL_AMBIENT, surface_new_mat);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, surface_new_mat);
	glMaterialfv(GL_FRONT, GL_SPECULAR, surface_new_mat);
	glMaterialf(GL_FRONT, GL_SHININESS, 5);
	DrawBRDF(Permute1, Phong);
//	glutSolidTeapot(50);

	// GREEN: Normal.
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, normal_mat);
	glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, normal_mat);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, normal_mat);
	double h = Identity(0,pi/2,Phong);
	double h2 = Permute1(0,pi/2,Phong);
	if (h2 < h)
		h = h2;
	glBegin(GL_LINES);
		PolarVertex(0,pi/2, h);
		PolarVertex(0,pi/2, 10);
	glEnd();
	/*
	// BLUE: Incoming light angle.
	glBegin(GL_LINES);
		PolarVertex(0,0,0);
		PolarVertex(mirror_az,pi-mirror_el,10);
	glEnd();
	*/
	// RED: Reflected light (mirror) angle.
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, mirror_mat);
	glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, mirror_mat);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mirror_mat);
	h = Identity(mirror_az,mirror_el,Phong);
	h2 = Permute1(mirror_az,mirror_el,Phong);
	if (h2 < h)
		h = h2;
	glBegin(GL_LINES);
		PolarVertex(mirror_az,mirror_el,h);
		PolarVertex(mirror_az,mirror_el,10);
	glEnd();

	glutSwapBuffers();
}

void KeyInputHandler(unsigned char Key, int x, int y)
{
	switch(Key)
	{
		case 27:		// Escape quits.
			exit(0);
			break;
		case '2':
			look_el -= 5;
			break;
		case '4':
			look_az -= 10;
			break;
		case '6':
			look_az += 10;
			break;
		case '8':
			look_el += 5;
			break;
		case '7':
			look_zoom *= 8/7.0;
			break;
		case '9':
			look_zoom /= 8/7.0;
			break;
		case 'f':
			mirror_el -= 2.5/180.0*pi;
			break;
		case 'd':
			mirror_az -= 5/180.0*pi;
			break;
		case 'g':
			mirror_az += 5/180.0*pi;
			break;
		case 'r':
			mirror_el += 2.5/180.0*pi;
			break;
		case 'h':
			param *= 8/7.0;
			break;
		case 'k':
			param /= 8/7.0;
			break;
		case 'u':
			cone *= 8/7.0;
			break;
		case 'j':
			cone /= 8/7.0;
			break;
		case ']':
			res += 4;
			break;
		case '[':
			res -= 4;
			break;
		case '1':
			surface_old_mat[3] += 0.1;
			surface_new_mat[3] -= 0.1;
			break;
		case '3':
			surface_old_mat[3] -= 0.1;
			surface_new_mat[3] += 0.1;
			break;
		default:
			;
	}

	coeff[0] = 0;
	coeff[1] = param;
	coeff[2] = -2*(param-1)/cone;
	coeff[3] = (param-1)/cone/cone;
//	coeff[0] = 0;
//	coeff[1] = param;
//	coeff[2] = -(param-1)/cone;
//	coeff[3] = 0;

	Display();
	cout << "slope parameter = " << param << ", \t" << "region extent = " << cone/pi*180 << endl;
}

void main(int argc, char* argv[])
{
	int width=500, height=500;

	glutInitDisplayMode(GLUT_RGB | GLUT_ALPHA | GLUT_DOUBLE | GLUT_DEPTH);
	glutInitWindowPosition(50, 50);
	glutInitWindowSize(width, height);

	glutCreateWindow("BRDF");

	glDisable(GL_DEPTH_TEST);
	glEnable(GL_CULL_FACE);
	glCullFace(GL_BACK);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	glutDisplayFunc(Display);
	glutKeyboardFunc(KeyInputHandler);

	glClearColor(0.0, 0.0, 0.0, 0.0);

	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(90, 1, .1, 400);
	glMatrixMode(GL_MODELVIEW);

	glLoadIdentity();

	glEnable(GL_LIGHTING);
	// Light source at the eyepoint.
	glEnable(GL_LIGHT0);
	glLightfv(GL_LIGHT0, GL_POSITION, camera_light_position);
	glLightfv(GL_LIGHT0, GL_AMBIENT, no_light);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, dim_light);
	glLightfv(GL_LIGHT0, GL_SPECULAR, no_light);
	// Directional light source from the top..
	glEnable(GL_LIGHT1);
	glLightfv(GL_LIGHT1, GL_POSITION, top_light_direction);
	glLightfv(GL_LIGHT1, GL_AMBIENT, dim_light);
	glLightfv(GL_LIGHT1, GL_DIFFUSE, bright_light);
	glLightfv(GL_LIGHT1, GL_SPECULAR, mid_light);

	glutMainLoop();
}