//*****************************************************************************************
// Disclaimer, Authorship and License
/*
	This code has been written to serve as a learning tool for the course MAE 574, 
	Virtual Reality Applications and Research, Spring 2009.
	
	Large parts of the code has been adapted from the examples given in the OpenGL 
	Programming Guide by Woo et. al and also from Computer Graphics Using OpenGL by FS Hill
	Some parts of the code has also been adapted from various resources available through 
	out the internet.I thank all those whose code and resources I have used to write these 
	examples. You may use this code for any non-commerical purpose. Feel free to include 
	and use parts of this code for your own application, but remember that this code is not
	entirely bug free, use it at your own risk. If you plan on using the code for your any 
	academic purpose please drop me an email. I would like to attach a link to your course 
	from my home page.

	The Author of this code is Govindarajan Srimathveeravalli, Dept. of Mech and Aero Eng. ,
	University at Buffalo.
*/
//******************************************************************************************

// This code gives a quick demonstration of two things, collision detection (through ray-plane
// intersection and subsequent point in polygon testing) and giving simple, but dynamically acc
// urate physics response. All objects move under two forces, gravity and a air friction - disspative
// force. Collisions are almost elastic (direction vectors are not completely preserved). 
// Try the following
// 1. Obviously, the collision detection is not fool proof. After a number of iterations, it fails, and the balls
//    slip through the plane. How to fix this?
// 2. The spheres do not interact with each other. Add a sphere-sphere intersection algorithm to handle this.
//    How will you modify the struct-functions for ball to handle collision forces coming in from other spheres?

#define EPSILON 0.001
#define TWOPI	6.28
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <glut.h>
#include <vector>
#include <time.h>

#include "geom.h"

using namespace std;
#define G 9.81

// for the plane on which the ball is going to bounce
cVector			planeNormal;
float			A, B, C, D;
cVector			vertex1, vertex2, vertex3, vertex4;

// White directional light
GLfloat white_light[] = { 1., 1., 1., 1. };
GLfloat light_posn[] = { 1., 1., 1., 0. }; 

// Blue material for the object
GLfloat mat_ambient[] = { 0., 0., 1., 1.0 };
GLfloat mat_diffuse[] = { 0.0, 0., 1., 1.0 };
GLfloat mat_specular[] = { 0.0, 0.0, 0.8, 1.0 };

// -------------------- code for the particle physics
float delT = 0.1;

struct balls{
		cVector posn;
		cVector vel;
		float mass;
};

balls b1;
balls allballs[100];

// Gravitation Force = -m*g*U
cVector evalGravitationalForce( void )
{
	cVector gravity;
	gravity = cVector( 0., -1., 0.)*G; // all objects have unit mass
	return gravity;
}

cVector evalExtnlFrictionForce( cVector currentState )
{
	// like dissipation force, in the opp direction of current velocity
	cVector dissipForce;	
	dissipForce = -currentState*0.1; // constant reduction of 10% for every timestep
	return dissipForce;
}

cVector evalForces( cVector currentState )
{
	cVector accumalator;
	accumalator += evalGravitationalForce();
	accumalator += evalExtnlFrictionForce( currentState);
	return accumalator;
}

void rungeKutta(balls &b)
{
	cVector k1 = new cVector( 0., 0., 0. );
	k1 = evalForces(b.vel)*delT;

	cVector k2 = new cVector( 0., 0., 0. );
	k2 = evalForces(b.vel+k1/2)*delT; 

	cVector k3 = new cVector( 0., 0., 0. );
	k3 = evalForces(b.vel+k2/2)*delT; 

	cVector k4 = new cVector( 0., 0., 0. );
	k4 = evalForces(b.vel+k3)*delT; 

	b.vel += (k1*1/6 + k2*1/3 + k3*1/3 + k4*1/6)*1/b.mass;
	b.posn += b.vel*delT;
}

// ------------------------------------  Plane interesection stuff ----//
// We use parametric equations that give a response to plane - line intersections and use that to perform bounc
// ing of the balls
void makePlane(void)
{
	A = vertex1.y*(vertex2.z - vertex3.z) + 
		vertex2.y*(vertex3.z - vertex1.z) + 
		vertex3.y*(vertex1.z - vertex2.z);

	B = vertex1.z*(vertex2.x - vertex3.x) + 
		vertex2.z*(vertex3.x - vertex1.x) + 
		vertex3.z*(vertex1.x - vertex2.x);

	C = vertex1.x*(vertex2.y - vertex3.y) + 
		vertex2.x*(vertex3.y - vertex1.y) + 
		vertex3.x*(vertex1.y - vertex2.y);

	D = -vertex1.x*(vertex2.y*vertex3.z - vertex3.y*vertex2.z) -
		 vertex2.x*(vertex3.y*vertex1.z - vertex1.y*vertex3.z) - 
		 vertex3.x*(vertex1.y*vertex2.z - vertex2.y*vertex1.z);
}

int intersection( cVector origin, const cVector& distance )
{
	// the line is in parametric form, t is parameter at which intersection
	// takes place
	double t;

	t = -( (A*origin.x) + (B*origin.y) + (C*origin.z) + D ) / ( (A*distance.x) + (B*distance.y) + (C*distance.z) );
	
	// if the t value is between 0 and 1 we know that the point intersects
	// the plane made by that particular triangle
	if(t < 0. || t > 1.0)
		return 0;
	else
		return 1;
}

// an example of a generic point in polygon test, can also be used to evaluate if the balls collide with
// the plane (note: the test uses only 3 vertices from 4, hence, it works only for the triangle making
// half the drawn quadrilateral
int PointInPoly( cVector test )
{
	double anglesum = 0. , costheta = 0.;
	double m1 = 0., m2 = 0.;

	// the 
	cVector T1 = cVector(	(vertex1.x - test.x), 
							(vertex1.y - test.y), 
							(vertex1.z - test.z) ) ;

	cVector T2 = cVector(	(vertex2.x - test.x),
							(vertex2.y - test.y), 
							(vertex2.z - test.z) ) ;

	cVector T3 = cVector(	(vertex3.x - test.x), 
							(vertex3.y - test.y), 
							(vertex3.z - test.z) ) ;

	// 1-2 pair
	m1 = sqrt(T1.x*T1.x + T1.y*T1.y + T1.z*T1.z);
	m2 = sqrt(T2.x*T2.x + T2.y*T2.y + T2.z*T2.z);
	
	if (m1*m2 <= EPSILON || m1*m2 == 0)
	{
		std::cout << "We are on a node, consider this inside 1 " << std::endl;
		return 0; 
	}
	else
	{
		if( m1 != 0 && m2 != 0 )
			costheta = (T1.x*T2.x + T1.y*T2.y + T1.z*T2.z) / (m1*m2);
	}
	
	if( abs(costheta) != 1)
		anglesum += acos(costheta);

	m1 = 0., m2 = 0.;
	costheta = 0.;

	// 2-3 pair

	m1 = sqrt(T2.x*T2.x + T2.y*T2.y + T2.z*T2.z);
	m2 = sqrt(T3.x*T3.x + T3.y*T3.y + T3.z*T3.z);
	
	if (m1*m2 <= EPSILON || m1*m2 == 0)
	{
		std::cout << "We are on a node, consider this inside 1 " << std::endl;
		return 0; 
	}
	else
	{
		if( m1 != 0 && m2 != 0 )
			costheta = (T3.x*T2.x + T3.y*T2.y + T3.z*T2.z) / (m1*m2);
	}
	
	if(abs(costheta) != 1)
		anglesum += acos(costheta);
	
	m1 = 0., m2 = 0.;
	costheta = 0.;

	// 3-1 pair
	m1 = sqrt(T1.x*T1.x + T1.y*T1.y + T1.z*T1.z);
	m2 = sqrt(T3.x*T3.x + T3.y*T3.y + T3.z*T3.z);
	
	if (m1*m2 <= EPSILON || m1*m2 == 0)
	{
		std::cout << "We are on a node, consider this inside 1 " << std::endl;
		return 0; 
	}
	else
	{
		if( m1 != 0 && m2 != 0 )
			costheta = (T1.x*T3.x + T1.y*T3.y + T1.z*T3.z) / (m1*m2);
	}
		
	
	if( abs(costheta) != 1)
		anglesum += acos(costheta);

	m1 = 0., m2 = 0.;
	costheta = 0.;

	if( abs(anglesum - TWOPI) >= 0.1 )
		return 1;
	else
		return 0;
}

void idle( int dummy )
{
	
	for( unsigned int i=0; i<100; i++)
	{
		if(intersection( allballs[i].posn, allballs[i].vel/10.))
		{
			allballs[i].vel = new cVector( allballs[i].vel.x/2., -allballs[i].vel.y, allballs[i].vel.z/2.);
		}
		else // simulate the system
		{
			rungeKutta( allballs[i] );
		}
	}

	glutTimerFunc( 50, idle, 0);
	glutPostRedisplay();
}


void reshape( int w, int h )
{
	glViewport( 0, 0, (GLsizei)w, (GLsizei)h );
	glMatrixMode( GL_PROJECTION );
	glLoadIdentity();
	glFrustum( -.4, .4, -.3, .3, 1., 2000. );
	glMatrixMode( GL_MODELVIEW );
	glutPostRedisplay();
}

void draw( )
{
	glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
	glLoadIdentity();
	//The first viewing command
	gluLookAt( 120., 120., 120. , 20., 0., 20., 0., 1., 0. );

	glBegin( GL_QUADS );
		glNormal3f( planeNormal.x, planeNormal.y, planeNormal.z );
		glVertex3f( vertex1.x, vertex1.y, vertex1.z);
		glVertex3f( vertex2.x, vertex2.y, vertex2.z);
		glVertex3f( vertex3.x, vertex3.y, vertex3.z);
		glVertex3f( vertex4.x, vertex4.y, vertex4.z);
	glEnd();

	glPushAttrib( GL_LIGHTING );
	glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
    glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
	for( unsigned int i=0; i<100; i++)
	{
		glPushMatrix();
			glTranslatef( allballs[i].posn.x, allballs[i].posn.y, allballs[i].posn.z );
			glutSolidSphere( 1.0, 10, 10 );
		glPopMatrix();
	}
	glPopAttrib();

	glutSwapBuffers();
}

void init( )
{
	glClearColor( .46, .53, .6, 0. );
	srand (time(NULL));

	// We keep our light within the display func so that we can influence its properties
	// Create the first white light
	glLightfv( GL_LIGHT0, GL_POSITION, light_posn );
	glLightfv( GL_LIGHT0, GL_AMBIENT, white_light );
	glLightfv( GL_LIGHT0, GL_DIFFUSE, white_light );
	glLightfv( GL_LIGHT0, GL_SPECULAR, white_light );	

	glEnable( GL_LIGHTING );
	glEnable( GL_LIGHT0 );

	// Enable depth testing
	glEnable( GL_DEPTH_TEST );
	glShadeModel( GL_FLAT );

	// set up a large sloping plane using glQuad
	vertex1 = cVector( -150., 0., 150. );
	vertex2 = cVector(  150., 0., 150. );
	vertex3 = cVector(  150., 10., -150. );
	vertex4 = cVector( -150., 10., -150. );

	// calculate plane parameters for intersection testing
	makePlane();
	planeNormal = cVector( 0., 1., .1 );

	for( unsigned int i=0; i<100; i++)
	{
		allballs[i].posn = cVector( float( rand() % 50), float( rand() % 50), float( rand() % 50) );
		allballs[i].vel = cVector( float( rand() % 10), -1., float( rand() % 10) );
		allballs[i].mass = 1.;
	}
}



int main(int argc, char* argv[])
{
	glutInit( &argc, argv );
	glutInitDisplayMode( GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
	glutInitWindowPosition( 50, 50 );
	glutInitWindowSize( 640, 480 );
	glutCreateWindow( "Bouncer" );
	init( );
	glutDisplayFunc( draw );
	glutReshapeFunc( reshape );
	glutTimerFunc( 50, idle, 0);
	glutMainLoop();

	return 0;
}