/*             
   Project 2:  Gas Chromatography
   Author:     Mike Dibley
   User name:  MDIBLEY
   Course:     CS3600 Computer Graphics I
   Instructor: Dr. Cunningham
   Due date:   October 7, 2002
   File:       mjd_p2.c
	
   Comments:   In my original design, I planned on using data that was as true to the real 
               world as possible. In reality, it was necessary  to tweak numbers somewhat 
			   to produce a nice visual.
               
               The top half of the screen is the representation of 3 gases (analytes)
               traveling through a GC column. The bottom half of the screen shows the 
               resulting chromatogram. The chromatogram shows the number of molecules
               exiting the column per unit time.
*/

#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>

#define COLUMNLENGTH  100    // meters
#define COLUMNRADIUS  265    // micrometers
#define XYSCALE       40     // this number is used to scale the column in xy plane. The
                             // COLUMNRADIUS is devided by this number
#define NUMVERTICES   100    // number of points in a circle from which the column is made
#define PI            3.14159
#define MAXDATAPOINTS 2000   // maximum number of points on the GC chromatogram
#define NUMMOLECULES  2000   // number of molecules in each gas
#define FLOWRATE      0.0025 // a number that produces a good looking result. This value
                             // has no actual tie to reality

typedef GLfloat point2D[2];
typedef GLfloat point3D[3];
typedef point3D gas[NUMMOLECULES];

gas      g_analyte1;
gas      g_analyte2;
gas      g_analyte3;
point2D  g_circle[NUMVERTICES];
point2D  g_data[MAXDATAPOINTS];  // data for the GC chromatogram
int      g_numdatapoints = 0;    // number of points on the GC chromatogram, current

/* function: Populate()  
   comments: randomly places the molecules for a given gas within a unit circle. The
             function places all molecules at z=0 (the left end of the column).
*/
void Populate(gas anal)
{
    int i=0;
    GLfloat randx, randy;
    
    while(i<NUMMOLECULES)
    {
        randx = 2 * ((GLfloat)rand()/RAND_MAX) -1;
        randy = 2 * ((GLfloat)rand()/RAND_MAX) -1;
        
 // if not true, then (x,y) lies outside of the unit circle, so do not keep
        if( 1.0 > (randx * randx + randy * randy) ) 
        {
            anal[i][0] = randx;
            anal[i][1] = randy;
            anal[i][2] = 0;
            i++;
        }
    } 
} 

/* function: MyInit()  
   comments: creates a unit circle from a number of vertices specified above as:
             NUMVERTICES. g_circle is then used by Cylinder() to create a cylinder.
             MyInit() also calls Populate for the 3 analytes (gases).
*/
void MyInit()
{
    int i;
    GLfloat angle;
    glEnable(GL_DEPTH_TEST);
    glClearColor(0.0, 0.0, 0.0, 1.0);
    
// create an array of vertices (a circle) for function cylinder()    
    for(i=0; i<NUMVERTICES; i++)
    {
        angle = 2 * PI * i / NUMVERTICES; 
        g_circle[i][0] = cos(angle);
        g_circle[i][1] = sin(angle);
    }
    
    Populate(g_analyte1);
    Populate(g_analyte2);
    Populate(g_analyte3);  
}

/* function: Cylinder()  
   comments: creates a cylinder of length=1 and radius=1. Note that the ends are NOT
             closed.
*/
void Cylinder()
{
    int i;
    glBegin (GL_QUAD_STRIP);
    
    for(i=0; i<NUMVERTICES; i++)
    {
        glVertex3f(g_circle[i][0], g_circle[i][1], 0);
        glVertex3f(g_circle[i][0], g_circle[i][1], 1);
    }
    glVertex3f(g_circle[0][0], g_circle[0][1], 0);
    glVertex3f(g_circle[0][0], g_circle[0][1], 1);
    glEnd();
}

/* function: Update()  
   comments: updates the position of all molecules in a gas. Note that only the 
             z-coordinate is changed. For a given single call to Update() not all
             molecules will move. The number of molecules that will move is determined
             by the partition coefficient, K, where:
             
             K = (molecules absorbed to stationary phase) / (molecules in mobile phase)

NOTE: I HAD THE ABOVE EQUATION REVERSED IN THE ORIGINAL SUBMISSION. IT IS CORRECT AS IT IS NOW.
             
             The function returns an integer equal to the number of molecules whose
             z-position has exceed 1.0 during the current function call. This is used
             to count the molecules exiting the column.
*/
int Update(gas anal, float part_coef)
{
    int i, detected = 0;
    GLfloat move;
    for(i=0; i<NUMMOLECULES; i++)
    {
        move = (GLfloat)rand()/RAND_MAX;
        
        if( move > (part_coef / (1 + part_coef)) )
        {
            anal[i][2] += FLOWRATE;
            if (anal[i][2] >= 1.0 && anal[i][2] < (1.0 + FLOWRATE) )
                detected++;
        }   
    } 
    return detected;
}

/* function: Display()  
   comments: displays the column. Note that GL_CLIP_PLANE1 is enabled before the modeling
             transformations and GL_CLIP_PLANE2 is enabled after. If I had to do again,
			 I think I would enable both after the transformation. This keeps all 
			 calculations about their location in unit space.
*/
void Display(void)
{
    int i;

// this clip plane allows us to see inside of the column
    GLdouble clip_eq1[] = { 0.0, -1.0, 0.0, 0.75*(COLUMNRADIUS/XYSCALE) };

// this clip plane hides all of the molecules that have exited the column
    GLdouble clip_eq2[] = { 0.0, 0.0, -1.0, 1.0 };
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);  
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
            
    gluLookAt(0.0, 100.0, 0.0,    // eye point
              0.0, 0.0, 0.0,      // center of view
              1.0, 0.0, 0.0);     // up	
    glPushMatrix();
    
    glEnable(GL_CLIP_PLANE1); // this clip plane used to open the column
    glClipPlane(GL_CLIP_PLANE1, clip_eq1);
    
    glColor3f(0.8, 0.8, 0.8);
    glTranslatef(15.0, 0.0, -COLUMNLENGTH/2);
    glScalef(COLUMNRADIUS/XYSCALE, COLUMNRADIUS/XYSCALE, COLUMNLENGTH);
    Cylinder();
    glDisable(GL_CLIP_PLANE1);
    
    glEnable(GL_CLIP_PLANE2); // this clip plane hides molecules that have exited
    glClipPlane(GL_CLIP_PLANE2, clip_eq2);
    
    glColor3f(1.0, 0.0, 0.0);
    glBegin(GL_POINTS);
    for (i=0; i<NUMMOLECULES; i++)
        glVertex3fv(g_analyte1[i]);
    glEnd();
    
    glColor3f(0.0, 1.0, 0.0);
    glBegin(GL_POINTS);
    for (i=0; i<NUMMOLECULES; i++)
        glVertex3fv(g_analyte2[i]);
    glEnd();
    
    glColor3f(0.0, 0.0, 0.1);
    glBegin(GL_POINTS);
    for (i=0; i<NUMMOLECULES; i++)
        glVertex3fv(g_analyte3[i]);
    glEnd();
    
    glDisable(GL_CLIP_PLANE2);
    
    glPopMatrix();
    glPushMatrix();
    
// now, add the chromatogram
    glColor3f(1.0, 1.0, 0.0);
    glLineWidth(2.0);
    glTranslatef(-40.0, 0.0, -COLUMNLENGTH/2);
    glScalef(1.0, 1.0, 0.05);
    
// the column was centered about the z-axis, so the chromatogram needs to be in the xz-plane
// (x is number of molecules seen by detector, z is time)
    glBegin(GL_LINE_STRIP);  
    for(i=0; i<g_numdatapoints; i++) 
        glVertex3f(g_data[i][1], 0.0, g_data[i][0]);
    glEnd();
    
    glDisable(GL_CLIP_PLANE2);

    glutSwapBuffers();
    glPopMatrix();
}

/* function: Reshape()  
   comments: standard reshape functionality
*/
void Reshape(int w,int h)
{
    glViewport(0,0,(GLsizei)w,(GLsizei)h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(45.0, (float)w/(float)h, 1.0, 300.0);
}

/* function: Animate()  
   comments: updates the position of all molecules in each analyte. For each call a 
             new data point is added to g_data, where x represents time (number of
             data points) and y represents number of molecules detected leaving 
             column.
*/
void Animate(void)
{       
    int detected = 0;
    
// note that the second number in the following function calls are partition coefficients.
// In the real world these numbers would be much larger, but larger values in this program
// produced too much spreading of the peaks.
    detected += Update(g_analyte1, 0.5);  
    detected += Update(g_analyte2, 1.0);
    detected += Update(g_analyte3, 2.5);
    
    g_data[g_numdatapoints][0] = g_numdatapoints;
    g_data[g_numdatapoints][1] = detected;
    g_numdatapoints++;
    
    if (g_numdatapoints == MAXDATAPOINTS)
        g_numdatapoints = 0;
	
// finished updating molecules; now do the display
    glutPostRedisplay();
}

int main(int argc, char** argv)
{
// Standard GLUT initialization
    glutInit(&argc,argv);
    glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH); 
    glutInitWindowSize(600,400);
    glutInitWindowPosition(50,50);
    glutCreateWindow("Gas Chromatography");
    glutDisplayFunc(Display);
    glutReshapeFunc(Reshape);
    glutIdleFunc(Animate);
		
    MyInit();
    glutMainLoop(); // enter event loop 
    return 0;
}