/*
   Simulation of a semi-permeable membrane

   (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?

typedef GLfloat point3[3];

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

// 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 heavy and light molecules
		for (i=0; i<NHEAVY; i++) {
			heavy[i][0] = (float)(random()%GRAN)/(float)GRAN;
			heavy[i][1] = (float)(random()%GRAN)/(float)GRAN;
			heavy[i][2] = (float)(random()%GRAN)/(float)GRAN;
		}
		for (i=0; i<NLIGHT; i++) {
			light[i][0] = (float)(random()%GRAN)/(float)GRAN;
			light[i][1] = (float)(random()%GRAN)/(float)GRAN;
			light[i][2] = (float)(random()%GRAN)/(float)GRAN;
		}

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

// draw the scene
void show(void)
    {
		int i;

		//	draw the "cage" : -1 to 1 in X, 0 to 1 in each of Y and Z
		glColor4f(1.0,1.0,1.0,1.0);
		glBegin(GL_LINE_STRIP);
			glVertex3f( 1.0, 1.0, 1.0);
			glVertex3f(-1.0, 1.0, 1.0);
			glVertex3f(-1.0, 0.0, 1.0);
			glVertex3f( 1.0, 0.0, 1.0);
			glVertex3f( 1.0, 1.0, 1.0);
		glEnd();
		glBegin(GL_LINE_STRIP);
			glVertex3f( 1.0, 1.0, 0.0);
			glVertex3f(-1.0, 1.0, 0.0);
			glVertex3f(-1.0, 0.0, 0.0);
			glVertex3f( 1.0, 0.0, 0.0);
			glVertex3f( 1.0, 1.0, 0.0);
		glEnd();
		glBegin(GL_LINES);
			glVertex3f( 1.0, 1.0, 1.0);
			glVertex3f( 1.0, 1.0, 0.0);
			glVertex3f(-1.0, 1.0, 1.0);
			glVertex3f(-1.0, 1.0, 0.0);
			glVertex3f( 1.0, 0.0, 1.0);
			glVertex3f( 1.0, 0.0, 0.0);
			glVertex3f(-1.0, 0.0, 1.0);
			glVertex3f(-1.0, 0.0, 0.0);
		glEnd();

		//	draw the membrane and its boundary
		glColor4f(0.1,0.1,0.1,1.0);
		glBegin(GL_POLYGON);
			glVertex3f(0.0,1.0,1.0);
			glVertex3f(0.0,1.0,0.0);
			glVertex3f(0.0,0.0,0.0);
			glVertex3f(0.0,0.0,1.0);
		glEnd();
		glColor4f(1.0,1.0,1.0,1.0);
		glBegin(GL_LINE_STRIP);
			glVertex3f(0.0,1.0,1.0);
			glVertex3f(0.0,1.0,0.0);
			glVertex3f(0.0,0.0,0.0);
			glVertex3f(0.0,0.0,1.0);
			glVertex3f(0.0,1.0,1.0);
		glEnd();
		
		//	draw the points
		glBegin(GL_POINTS);
			glPointSize(4.0);
			glColor4f(1.0,0.0,0.0,1.0);
			for ( i=0; i<=NHEAVY; i++ ) {
				glVertex3fv(heavy[i]);
			}
			glPointSize(2.0);
			glColor4f(0.0,1.0,1.0,1.0);
			for ( i=0; i<=NLIGHT; i++ ) {
				glVertex3fv(light[i]);
			}
		glEnd();
    }
    
// update the positions of all the molecules
void animate(void)
{
	#define HEAVYLEFT  0.8
	#define HEAVYRIGHT 0.9
	#define LIGHTLEFT  0.3
	#define LIGHTRIGHT 0.1
	
	#define LEFT  0
	#define RIGHT 1
	int i,sign,whichside;
	
	for(i=0; i<NHEAVY; i++) {
		if (heavy[i][0] < 0.0) whichside = LEFT;
		else whichside = RIGHT;
		
		sign = random()%2; if (sign == 0) sign = -1;
		heavy[i][0] += sign*(float)(random()%GRAN)/(float)(16*GRAN);
		sign = random()%2; if (sign == 0) sign = -1;
		heavy[i][1] += sign*(float)(random()%GRAN)/(float)(16*GRAN);
		sign = random()%2; if (sign == 0) sign = -1;
		heavy[i][2] += sign*(float)(random()%GRAN)/(float)(16*GRAN);

		if (heavy[i][0] > 1.0) heavy[i][0] = 1.0 - heavy[i][0];

		if ((whichside == RIGHT)&&(heavy[i][0] < 0.0))		// does it cross right to left?
			if ( (float)(random()%GRAN)/(float)(GRAN) >= HEAVYLEFT)
				heavy[i][0] = -heavy[i][0];

		if ((whichside == LEFT)&&(heavy[i][0] > 0.0))		// does it cross left to right?
			if ( (float)(random()%GRAN)/(float)(GRAN) >= HEAVYRIGHT)
				heavy[i][0] = -heavy[i][0];

		if (heavy[i][0] <-1.0) heavy[i][0] =-1.0 - heavy[i][0];
		if (heavy[i][1] > 1.0) heavy[i][1] = 1.0 - heavy[i][1];
		if (heavy[i][1] < 0.0) heavy[i][1] = -heavy[i][1];
		if (heavy[i][2] > 1.0) heavy[i][2] = 1.0 - heavy[i][2];
		if (heavy[i][2] < 0.0) heavy[i][2] = -heavy[i][2];
	}
	for(i=0; i<NLIGHT; i++) {
		if (light[i][0] < 0.0) whichside = LEFT;
		else whichside = RIGHT;

		sign = random()%2; if (sign == 0) sign = -1;
		light[i][0] += sign*(float)(random()%GRAN)/(float)(16*GRAN);
		sign = random()%2; if (sign == 0) sign = -1;
		light[i][1] += sign*(float)(random()%GRAN)/(float)(16*GRAN);
		sign = random()%2; if (sign == 0) sign = -1;
		light[i][2] += sign*(float)(random()%GRAN)/(float)(16*GRAN);

		if (light[i][0] > 1.0) light[i][0] = 1.0 - light[i][0];

		if ((whichside == RIGHT)&&(light[i][0] < 0.0))		// does it cross right to left?
			if ( (float)(random()%GRAN)/(float)(GRAN) >= LIGHTLEFT)
				light[i][0] = -light[i][0];

		if ((whichside == LEFT)&&(light[i][0] > 0.0))		// does it cross left to right?
			if ( (float)(random()%GRAN)/(float)(GRAN) >= LIGHTRIGHT)
				light[i][0] = -light[i][0];

		if (light[i][0] <-1.0) light[i][0] =-2.0 - light[i][0];
		if (light[i][1] > 1.0) light[i][1] = 1.0 - light[i][1];
		if (light[i][1] < 0.0) light[i][1] = -light[i][1];
		if (light[i][2] > 1.0) light[i][2] = 1.0 - light[i][2];
		if (light[i][2] < 0.0) light[i][2] = -light[i][2];
	}
	glutPostRedisplay();
}

void display( void )
{   

	static GLfloat camX, camY, camZ;

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

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

    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)
{
        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;
    	    case 'x' :
				camRadius += 0.5; break;
    	    case 'z' :
    			camRadius -= 0.5;
				if(camRadius < 1.0)
					camRadius = 1.0;
				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 all the molecules
void tally(void)
{
	int hl, hr, ll, lr, i;
	float hratio, lratio;
	
	hl = hr = ll = lr = 0;
	for (i=0;i<NHEAVY;i++)
		if (heavy[i][0]<0.0) hl++;
	else hr++;
	for (i=0;i<NHEAVY;i++)
		if (light[i][0]<0.0) ll++;
	else lr++;
	hratio = (float)hl/(float)hr;
	lratio = (float)ll/(float)lr;
	printf("Heavy: left %d, right %d L/R ratio %f\nLight: left %d, right %d L/R ratio %f\n\n",
			hl,hr,hratio,ll,lr,lratio);
}

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("Semi-permeable Membrane Simulation");
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutKeyboardFunc(keyboard);
        glutIdleFunc(animate);
        glutMouseFunc(mouse);
        glutMotionFunc(motion);
		
        myinit();
        glutMainLoop();

        return 0;
}