/*
	Estimation of a volume by Monte Carlo techniques

	Source file to be used with
	Cunningham, Computer Graphics: Programming in OpenGL for Visual Communication, Prentice-Hall, 2007

	Intended for class use only
*/

#include <GLUT/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

typedef GLfloat point3[3];

//	definitions and data for the points in the space
#define NPTS 2000	// numbers of heavy and light molecules -- larger for screen viewing
#define GRAN  1000	// granularity of positions and motions in the initial space
point3 pts[NPTS];

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

//	the spheres
#define NSPHERES 50
typedef struct theSphere {
	float x, y, z, r;
	} theSphere;
theSphere aSphere[NSPHERES];
int incount=0, outcount=0, counting=0;
float proportion;

#define ANGLE 2.0

static char ch;
float       ep = 3.0, bound = 1.0, distance = .06;
int			bounce;
GLfloat		xAngle = 0., yAngle = 0., zAngle = 0.;

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

void myinit(void)
{
		int i;
		float maxR;

//		initialize random number generator
		srand( clock() );
		
//		initialize the locations of the molecules
		for (i=0; i<NPTS; i++) {
			pts[i][0] = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
			pts[i][1] = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
			pts[i][2] = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
		}
//		initialize the locations of the spheres
		for ( i=0; i<NSPHERES; i++ ) {
			aSphere[i].x = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
			aSphere[i].y = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
			aSphere[i].z = -1.0 + 2.0*(float)(rand()%GRAN)/(float)GRAN;
			maxR = 1.0;
			if (fabs(aSphere[i].x-1)<maxR) maxR = fabs(aSphere[i].x-1);
			if (fabs(aSphere[i].x+1)<maxR) maxR = fabs(aSphere[i].x+1);
			if (fabs(aSphere[i].y-1)<maxR) maxR = fabs(aSphere[i].y-1);
			if (fabs(aSphere[i].y+1)<maxR) maxR = fabs(aSphere[i].y+1);
			if (fabs(aSphere[i].z-1)<maxR) maxR = fabs(aSphere[i].z-1);
			if (fabs(aSphere[i].z+1)<maxR) maxR = fabs(aSphere[i].z+1);
			aSphere[i].r = maxR*(float)(rand()%GRAN)/(float)GRAN;
		}

        glClearColor( 0.4, 0.4, 0.4, 0.0 );
        glDisable(GL_DEPTH_TEST);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
        glEnable(GL_BLEND);
}

void show(void)
    {
		int i;
		GLUquadric *S;

		//	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 spheres in the space
		glColor4f(1.0, 1.0, 0.0, 0.1);
		for ( i=0; i<NSPHERES; i++ ) {
			S = gluNewQuadric();
			glPushMatrix();
			glTranslatef(aSphere[i].x,aSphere[i].y,aSphere[i].z);
			gluSphere(S,aSphere[i].r,20,20);
			glPopMatrix();
		}
		
		//	draw the points
		glPointSize(2.0);
		glBegin(GL_POINTS);
			for ( i=0; i<=NPTS; i++ ) {
				//	red if inside, green if not inside
				if (inside(i)) {
					glColor4f(1.0,0.0,0.0,1.0);
					glPointSize(8.0);
					if (counting == 0) incount++;
				}
				else {
					glColor4f(0.0,1.0,0.0,1.0);
					glPointSize(2.0);
					if (counting == 0) outcount++;
				}
				glVertex3fv(pts[i]);
			}
		glEnd();
		proportion = 8.0*(float)incount/(float)(incount+outcount);
		if (counting == 0) {
			counting = 1;
			printf("Counts: : %d in, %d out, volume = %f\n",incount,outcount,proportion);
		}
    }

int inside(int j)
{
	int i;
	
	for (i=0; i<NSPHERES; i++) {
		if ((pts[j][0]-aSphere[i].x)*(pts[j][0]-aSphere[i].x)+
			(pts[j][1]-aSphere[i].y)*(pts[j][1]-aSphere[i].y)+
			(pts[j][2]-aSphere[i].z)*(pts[j][2]-aSphere[i].z) < aSphere[i].r*aSphere[i].r)
				return 1;
	}
	return 0;
}

void display( void )
{    
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//           eye point      center of view      up
	gluLookAt( 0.0, 1.2*ep, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);

	glRotatef(xAngle, 1.0, 0.0, 0.0);
	glRotatef(yAngle, 0.0, 1.0, 0.0);
	glRotatef(zAngle, 0.0, 0.0, 1.0);
	if (ch != ' ') {	// keyboard indicated a rotation...
	//	rotation control uses two keypads embedded in standard keyboard
	//    q  w                     o  p  -- rotate around x-axis + -
	//     a  s                   k  l   -- rotate around y-axis + -
	//      z  x                 m  ,    -- rotate around z-axis + -
	switch(ch) {
			case 'q':
			case 'o':
				xAngle += ANGLE;
				glRotatef( ANGLE, 1.0, 0.0, 0.0); break;
			case 'w':
			case 'p':
				xAngle -= ANGLE;
				glRotatef(-ANGLE, 1.0, 0.0, 0.0); break;
        	case 'a':
        	case 'k':
				yAngle += ANGLE;
            	glRotatef( ANGLE, 0.0, 1.0, 0.0); break;
        	case 's':
        	case 'l':
				yAngle -= ANGLE;
            	glRotatef(-ANGLE, 0.0, 1.0, 0.0); break;
        	case 'z':
        	case 'm':
				zAngle += ANGLE;
            	glRotatef( ANGLE, 0.0, 0.0, 1.0); break;
        	case 'x':
        	case ',':
				zAngle -= ANGLE;
            	glRotatef(-ANGLE, 0.0, 0.0, 1.0); break;
	    }
    }
	ch = ' ';	// ensure any residual rotation is killed

    show();

    glutSwapBuffers();
 }

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

void keyboard(unsigned char key, int x, int y)
{
	ch = ' ';
	switch (key)
	{
		case 'w' :
		case 'q' :
		case 's' :
		case 'a' :
		case 'x' :
		case 'z' :
			ch = key; break;
		case 'f' :
			ep = ep - 0.2; break;
		case 'b' :
			ep = ep + 0.2; break;
	}
	glutPostRedisplay(); /* perform display again */
}

int main(int argc, char** argv)
{
/* Standard GLUT initialization */
	glutInit(&argc,argv);
	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(70,70);
	glutCreateWindow("Monte Carlo volume estimation");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);
		
	myinit();
	glutMainLoop();
}