/*
   Simulation of gas laws in a bounded volume

   (c) February 2001, Steve Cunningham & Ken Brown
*/
#define OSX

#ifdef MAC
#include "glut.h"
#endif
#ifdef OSX
#include <GLUT/glut.h>
#endif
#ifdef UNIX
#include <GL/glut.h>
#endif
#ifdef PC
#include <GL/glut.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define TRUE 1
#define FALSE 0

// variables to keep track of the current rotation on each axis
GLfloat camPitch = 0.0, camYaw = 0.0, camRadius = 3.0;
int mx,my; // used to calculate mouse displacement
int mdown; // is a button down?
int buttonDown; // if so, which button?

//static char ch;
float bound = 1.0, distance = .06;
int bounce;
typedef GLfloat point3[3];

//	definitions and data for the molecules in the space
#define NMOLS 1000	// numbers of heavy and light molecules
#define GRAN  1000	// granularity of positions and motions in the initial space
point3 mols[NMOLS];

//	trail for one molecule and number of points in that trail
#define TRAILSIZE 50
int npoints = 0;
point3 trail[TRAILSIZE];

// function prototypes
void myinit(void);
void show(void);
void display(void);
void reshape(int,int);
void keyboard(unsigned char,int,int);
void animate(void);
void tally();

// initialize OpenGL and all necessary global variables
void myinit(void)
{
		int i;

//		initialize random number generator
		srand( clock() );
		
//		initialize the locations of the molecules
		for (i=0; i<NMOLS; i++) {
			mols[i][0] = -1.0 + 2.0*(float)(random()%GRAN)/(float)GRAN;
			mols[i][1] = -1.0 + 2.0*(float)(random()%GRAN)/(float)GRAN;
			mols[i][2] = -1.0 + 2.0*(float)(random()%GRAN)/(float)GRAN;
		}

        glClearColor( 0.0, 0.0, 1.0, 0.0 );
        glEnable(GL_DEPTH_TEST);
}

// draw all entities in the scene
void show(void)
{
	int i, j, k;

	//	draw the "cage" : initially -1 to 1 in each of X, Y, and Z but bound can vary
	glColor4f(1.0,1.0,1.0,1.0);
	glBegin(GL_LINE_STRIP);
		glVertex3f( bound, bound, bound);
		glVertex3f(-bound, bound, bound);
		glVertex3f(-bound,-bound, bound);
		glVertex3f( bound,-bound, bound);
		glVertex3f( bound, bound, bound);
	glEnd();
	glBegin(GL_LINE_STRIP);
		glVertex3f( bound, bound,-bound);
		glVertex3f(-bound, bound,-bound);
		glVertex3f(-bound,-bound,-bound);
		glVertex3f( bound,-bound,-bound);
		glVertex3f( bound, bound,-bound);
	glEnd();
	glBegin(GL_LINES);
		glVertex3f( bound, bound, bound);
		glVertex3f( bound, bound,-bound);
		glVertex3f(-bound, bound, bound);
		glVertex3f(-bound, bound,-bound);
		glVertex3f( bound,-bound, bound);
		glVertex3f( bound,-bound,-bound);
		glVertex3f(-bound,-bound, bound);
		glVertex3f(-bound,-bound,-bound);
	glEnd();
		
	//	draw the points
	glPointSize(4.0);
	glColor4f(1.0,0.0,0.0,1.0);
	glBegin(GL_POINTS);
		for ( i=0; i<=NMOLS; i++ ) {
			glVertex3fv(mols[i]);
			if (i == 0) {
				if (npoints < TRAILSIZE) npoints++;
				for (j=TRAILSIZE-1; j>0; j--)
					for (k=0; k<3; k++)
						trail[j][k]=trail[j-1][k];
				for (k=0; k<3; k++){
					trail[0][k] = mols[i][k];
				}
			}
		}
	glEnd();
	glColor4f(1.0,1.0,0.0,1.0);
	glBegin(GL_LINE_STRIP);
		for (j=0; j<npoints; j++) {
			glVertex3fv(trail[j]);
		}
	glEnd();
}

// update the position of all the molecules
void animate(void)
{
	int i,j,sign;
	
	bounce = 0;
	for(i=0; i<NMOLS; i++) {
		for (j=0; j<3; j++) {
			sign = -1+2*(random()%2);
			mols[i][j] += (float)sign*distance*((float)(random()%GRAN)/(float)(GRAN));
			if (mols[i][j] > bound)
				{mols[i][j] = 2.0*bound - mols[i][j]; bounce++;}
			if (mols[i][j] <-bound)
				{mols[i][j] = 2.0*(-bound) - mols[i][j]; bounce++;}
		}
	}
	glutPostRedisplay();
}

void display( void )
{
	static GLfloat camX, camY, camZ;

	// convert from spherical coords to cartesian
	camX = cos(camYaw) * cos(camPitch) * camRadius;
	camY = sin(camYaw) * cos(camPitch) * camRadius;
	camZ = sin(camPitch) * camRadius;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//           eye point      center of view      up
	gluLookAt( camX, camY, camZ, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);

    show();

    glutSwapBuffers();
}

void reshape(int w,int h)
{
        glViewport(0,0,(GLsizei)w,(GLsizei)h);
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
        gluPerspective(60.0,(float)w/(float)h,0.5,30.0);
}

void keyboard(unsigned char key, int x, int y)
{
        int i;
        float epsilon = .01;

        switch (key)
        {
			case 'w' :
				camYaw += .08;
				if(camYaw > 6.28)
					camYaw = camYaw-6.28;
				break;
    	    case 'q' :
				camYaw -= .08;
				if(camYaw < 0.0)
					camYaw = 6.28+camYaw;
				break;
    	    case 's' :
				camPitch += .08;
				if(camPitch > 1.57)
					camPitch = 1.57;
				break;
    	    case 'a' :
				camPitch -= .08;
				if(camPitch < -1.57)
					camPitch = -1.57;
				break;
			// change the size of the cube
			case '-' :
				if(bound > 0.1)
					bound = bound - 0.05; 
				for (i=0; i<NMOLS; i++) {
					if (mols[i][0] > bound) mols[i][0] = bound-epsilon;
					if (mols[i][0] <-bound) mols[i][0] =-bound+epsilon;
					if (mols[i][1] > bound) mols[i][1] = bound-epsilon;
					if (mols[i][1] <-bound) mols[i][1] =-bound+epsilon;
					if (mols[i][2] > bound) mols[i][2] = bound-epsilon;
					if (mols[i][2] <-bound) mols[i][2] =-bound+epsilon;
				}
				break;
			case '+' :
			case '=' :
				bound = bound + 0.05; 
				for (i=0; i<NMOLS; i++) {
					if (mols[i][0] > bound) mols[i][0] = bound-epsilon;
					if (mols[i][0] <-bound) mols[i][0] =-bound+epsilon;
					if (mols[i][1] > bound) mols[i][1] = bound-epsilon;
					if (mols[i][1] <-bound) mols[i][1] =-bound+epsilon;
					if (mols[i][2] > bound) mols[i][2] = bound-epsilon;
					if (mols[i][2] <-bound) mols[i][2] =-bound+epsilon;
				}
				break;
			case 't' :
				tally(); break;
        }
}

void mouse(int button, int state, int x, int y) {
	// save the mouse position for the camera rotation
	mx = x;
	my = y;
	if(state == GLUT_DOWN)
		mdown = TRUE;
	else
		mdown = FALSE;

	buttonDown = button; // remember which button is down
}

void motion(int x, int y) {
	if(buttonDown == GLUT_LEFT_BUTTON) {
		// adjust the camera's yaw
		camYaw += 0.02*(GLfloat)(mx-x);

		// wrap the yaw values at 0 and 2*pi
		if(camYaw < 0.0)
				camYaw = 6.28+camYaw;
		else if(camYaw > 6.28)
				camYaw = camYaw-6.28;

		// adjust the camera's pitch
		camPitch += 0.02*(GLfloat)(y-my);

		// make sure we don't pitch too far
		if(camPitch > 1.57)
			camPitch = 1.57;
		else if(camPitch < -1.57)
			camPitch = -1.57;
	}
	else {
		// zoom in and out if right mouse button is down
		camRadius += 0.5*(GLfloat)(y-my);
		if(camRadius < 1.0)
			camRadius = 1.0;
	}

	mx = x;
	my = y;
}

// output statistics on the pressure in the cube
void tally(void)
{
	float pressure, volume;

	pressure = (float)bounce/(bound*bound);	// hits per unit area
	volume   = bound*bound*bound;		// dimension cubed
	printf("pressure %f, volume %f, product %f\n",pressure,volume,pressure*volume);
}

int main(int argc, char** argv)
{
/* Standard GLUT initialization */
        glutInit(&argc,argv);
        glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
        glutInitWindowSize(500,500);
        glutInitWindowPosition(70,70);
        glutCreateWindow("gas laws simulation");
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutKeyboardFunc(keyboard);
        glutMouseFunc(mouse);
        glutMotionFunc(motion);
        glutIdleFunc(animate);
		
        myinit();
        glutMainLoop();

        return 0;
}