/*
   Program to read molecular information and display the molecule so
   it can be manipulated in simple ways, including rotation in three
   axes, zooming in and out, changing the size of the atoms in the
   molecule, and making the atoms more or less transparent.
   
   One of these is a highly simplified MDL molecular specification
   file.  Reference:
      CTFile formats, MDL Information Systems, Inc.
      San Leandro, CA 94577, 1999, http://www.mdli.com/

   The other is the Protien Data Bank (PDB) file format.  Reference:
      Protein Data Bank Contents Guide: Atomic Coordinate Entry Format 
      Description, Protein Data Bank, Brookhaven National Laboratory,
      Upton, NY 11973-5000 USA, http://www.pdb.bnl.gov
      
   A menu will allow the user to select from a small number of
   visualization options, and some of these, as well as rotations,
   are available from the keyboard.

   Proof of concept for chemistry project for intro graphics course
   
   (c) 1999, Steve Cunningham
*/

#include "glut.h"
#include "molecules.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <ctype.h>

// definitions of parameters for program

#define AMAX 512	// maximum number of atoms in molecule
#define BMAX 512	// maximum number of bonds in molecule
#define ZOOM 0.5	// size of zoom in - zoom out step
#define GRAIN 25	// resolution of glu models
#define ANGLE 2.0	// rotation angle in degrees

// file declaration for molecular description input file
FILE	*fp;
int		argc;
char	**argv;

// data for storing molecular description

int		natoms, nbonds;
float 	sizeMult = 1.0, alphaAdd = 0.0, t = 0.0;

typedef struct atomdata {
					float x, y, z;
					char  name[5];
					int   colindex;
					} atomdata;
atomdata atoms[AMAX];
typedef struct bonddata{
					int first, second, bondtype;
					} bonddata;
bonddata bonds[BMAX];

//	data for interaction
int Tag1=1, Tag2=2, Tag3=3, Tag4=4, Tag5=5, Tag6=6, Tag7=7, Tag8=8,
	Tag9=9, selection;
float	ep; // eye point variable
char *fname, *mainName, *ext;
static char ch; // character from keyboard that controls the rotation
static GLfloat saveState[16] = {1.0,0.0,0.0,0.0,
                                0.0,1.0,0.0,0.0,
                                0.0,0.0,1.0,0.0,
                                0.0,0.0,0.0,1.0},
               viewProj[16]; // hold transforms
typedef GLfloat point3[3];

// function prototypes
char *getFileName(void);
void separateName(char *);
void readMolFile(char *);
void readPdbFile(char *);
void myinit(void);
void molecule(void);
void display(void);
void reshape(int, int);
void keyboard(unsigned char, int, int);
void animate(void);
void options_menu(int);
void giveHelp(void);
void giveAbout(void);
int  lookUpName(int);

char *getFileName(void)
{
	char *myFileName = "";
	
	printf("enter name of molecule file to be displayed > ");
	scanf("%s",myFileName);
	return(myFileName);
}

void separateName(char *fname)
{
	int i=0, j;
	
	mainName = (char *)malloc(20);
	ext = (char *)malloc(4);
	//	loop to find the last period in the file name; save at value j
	//	and separate the file name to extract the extension, which is
	//	used to determine which file reading function is used
	for (i=0; fname[i]!='\0'; i++)
		if (fname[i]=='.')
			j = i;
	for (i=0; i<j; i++)
		mainName[i]=fname[i];
	mainName[j]='\0';
	for (i=0; i<3; i++)
		ext[i]=fname[i+j+1];
	ext[3]='\0';
	return;
}

void readMolFile(char *fname)
{
	int i,j;
	char ignore=' ', myLine[80];

//	read data from file for the model
	if ( (fp = fopen(fname, "r")) == NULL )
		{
			fprintf(stderr,"Unable to open file!!!\n" );
			exit(1);
		}
//	begin by skipping the three-line header block
	for (i=0; i<3; i++) {
		fgets(myLine, 80, fp);
	}
//	now read the counts followed by the atom and bond data
	fscanf(fp,"%d %d%[^\n]\n", &natoms, &nbonds, myLine);
	for (i=0; i<natoms; i++)	//	read atom location data and get table index
		{
			fscanf(fp, "%f %f %f", &atoms[i].x, &atoms[i].y, &atoms[i].z );
			fscanf(fp, "%c", &ignore);
			for (j=0; j<4; j++)
				fscanf(fp,"%c", &atoms[i].name[j]);
			atoms[i].name[4]='\0';
			atoms[i].colindex = lookUpName(i);
			fgets(myLine, 80, fp);
		}
	for (i=0; i<nbonds; i++)	//	get bond data
		{
			fscanf(fp,"%d %d %d", &bonds[i].first, &bonds[i].second,
				&bonds[i].bondtype);
			fgets(myLine, 80, fp);
		}
	fclose(fp);
}

void readPdbFile(char *fname)
{
	int i=0,j=0, k, m, n=0, atomflag = 0, conectflag = 0;
	int source, from, destination, to;
	float x, y, z;
	char ignore=' ', *myLine, *token, *atomName;

	for (k=0; k<BMAX; k++)
		bonds[k].bondtype = 0;
	token = (char *)malloc(8);
	myLine = (char *)malloc(80);
	atomName = (char *)malloc(3);
//	read data from file for the model
	if ( (fp = fopen(fname, "r")) == NULL )
		{
			fprintf(stderr,"Unable to open file!!!\n" );
			exit(1);
		}
	while (fgets(token,7,fp) != NULL) {
		if (strcmp(token,"ATOM  ")==0) {
			fgets(myLine,7,fp);
			fgets(atomName,3,fp);
			fgets(myLine,15,fp);
			fscanf(fp,"%f %f %f",&x,&y,&z);
			// now change atomName into standard form
			// extract first letter of atomName
			k = 0;
			while (atomName[k]==' ') k++;
			atoms[i].name[0] = atomName[k];
			k++;
			if (atomName[k] != ' ') {
				if ((atomName[k] >='a') && (atomName[k] <= 'z'))
					atoms[i].name[1] = atomName[k];
				else 
					if ((atomName[k] >='A') && (atomName[k] <= 'Z'))
						atoms[i].name[1] = tolower(atomName[k]);
					else atoms[i].name[1] = ' ';
				}
			atoms[i].name[2] = ' ';
			atoms[i].name[3] = ' ';
			atoms[i].name[4] = '\0';
			// now look up atom in the table and set data in atom list
			atoms[i].colindex = lookUpName(i);
			atoms[i].x = x;
			atoms[i].y = y;
			atoms[i].z = z;
			i++;
			}
/*
	NOTE: double bonds are not indicated in PDB files...
*/
		if (strcmp(token,"CONECT")==0){
			fgets(token,6,fp);
			source = atoi(token);
			k = 0;
			while (!feof(fp) && (k<4)) {
				fgets(token,6,fp);
				destination = atoi(token);
				if (destination != 0) {
					if (destination < source) {
						to = source;
						from = destination;
					}
					else {
						to = destination;
						from = source;
					}			
					for (m=0;m<n;m++) { // is this bond is already read?
						if ((bonds[m].first == from)&&(bonds[m].second == to)) {
							m = n+1;	// set flag to not enter the bond
							break;
						}
					}
					if (m == n) {
						bonds[n].first = from;
						bonds[n].second = to;
						bonds[n].bondtype++;
						n++;
					}
				}
				k++;
			}
			j++;
		}
		fgets(myLine,80,fp);
	}
	natoms = i;
	nbonds = n;
}

void myinit(void)
{
	int i;
	float max = 0.0;
	
//	set up overall light data, including specular=ambient=light colors
	GLfloat light_position[]={ 10.0, 10.0, -10.0, 1.0 };
	GLfloat light_color[]={ 1.0, 1.0, 1.0, 1.0 };
	GLfloat ambient_color[]={ 0.2, 0.2, 0.2, 1.0 };
	GLfloat mat_specular[]={ 1.0, 1.0, 1.0, 1.0 };

	glClearColor( 0.7, 0.7, 1.0, 0.0 );
	glShadeModel(GL_SMOOTH);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_specular );
	glLightfv(GL_LIGHT0, GL_POSITION, light_position );
	glLightfv(GL_LIGHT0, GL_AMBIENT, ambient_color );
	glLightfv(GL_LIGHT0, GL_SPECULAR, light_color );
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light_color );

//	attributes
	glEnable(GL_LIGHTING);   // so lighting models are used
	glEnable(GL_LIGHT0);     // we'll use LIGHT0
	glEnable(GL_DEPTH_TEST); // allow z-buffer display
	glEnable(GL_BLEND);      // enable alpha-channel blending        
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	glEnable(GL_LINE_SMOOTH);// allow line smoothing

//	scan range of atom locations to set initial eye point
	for (i=0; i<natoms; i++) {
		if (fabs(atoms[i].x) > max) max = fabs(atoms[i].x);
		if (fabs(atoms[i].y) > max) max = fabs(atoms[i].y);
		if (fabs(atoms[i].z) > max) max = fabs(atoms[i].z);
	}
	ep = 2.0 * max;
}

void molecule(void)
{
//	Create the molecular model in a way that fits a display list.
//	Note that this means that anything that will be changed must be
//	referenced indirectly, because anything defined directly cannot
//	be changed when the display list is executed.  Thus the colors
//	and sizes are defined by references into external (global) arrays

#define ANGTOAU(ang)    ( (ang)*0.529177 )
#define DBGAP 0.05

	int i, j, k, startindex, endindex, bondtype;
	GLUquadric *atomSphere;

    GLfloat myColor[4];
    GLfloat color1[]={1.0, 0.0, 0.0, 0.7};
    GLfloat color3[]={0.0, 1.0, 1.0, 1.0};
    GLfloat mat_shininess[]={ 50.0 };

    glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );
    
//  use the location and link material that was read from the file
//	to display the molecule

//	note that we probably should change the color of the bond for dbl, trpl
    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, bondBlack );
    glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, bondBlack );
	for (i=0; i<nbonds; i++)
		{	// draw the bonds - note that the file is 1-origin indexed
			// while our tables are 0-origin indexed (language is C) --
			startindex = bonds[i].first-1;	//so we decrease index by 1
			endindex   = bonds[i].second-1;
			bondtype   = bonds[i].bondtype;
			if (bondtype == 1) {
				glLineWidth(5.0);
				glBegin(GL_LINE_STRIP);
					glVertex3f(atoms[startindex].x,atoms[startindex].y,
						atoms[startindex].z);
					glVertex3f(atoms[endindex].x,atoms[endindex].y,
						atoms[endindex].z);
				glEnd();
			}
			if (bondtype == 2) {
				glLineWidth(3.0);
				glBegin(GL_LINE_STRIP);
					glVertex3f(atoms[startindex].x-DBGAP,atoms[startindex].y-DBGAP,
						atoms[startindex].z-DBGAP);
					glVertex3f(atoms[endindex].x-DBGAP,atoms[endindex].y-DBGAP,
						atoms[endindex].z-DBGAP);
				glEnd();
				glBegin(GL_LINE_STRIP);
					glVertex3f(atoms[startindex].x+DBGAP,atoms[startindex].y+DBGAP,
						atoms[startindex].z+DBGAP);
					glVertex3f(atoms[endindex].x+DBGAP,atoms[endindex].y+DBGAP,
						atoms[endindex].z+DBGAP);
				glEnd();
			}
		}

	for (i=0; i<natoms; i++)
		{	// draw the atoms
			glPushMatrix();
			atomSphere=gluNewQuadric();
			j = atoms[i].colindex;	// index of color for atom i
			for (k=0; k<4; k++) // copy atomColors[j] and adjust alpha by alphaMult
			{
				myColor[k] = atomColors[j][k];
			}
			if (j==5) myColor[3] += alphaAdd;
			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, myColor );
			glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, myColor );
			glTranslatef(atoms[i].x, atoms[i].y, atoms[i].z);
			if (j==5)
				gluSphere(atomSphere, sizeMult*ANGTOAU(atomSizes[j]), GRAIN, GRAIN);
			else
				gluSphere(atomSphere, ANGTOAU(atomSizes[j]), GRAIN, GRAIN);
			glPopMatrix();
		}
    glutSwapBuffers();
 }

void display( void )
{
//  display the molecule with the rotations you specified from the keyboard
//	and the parameters you set in the menu

//  The GL_MODELVIEW_MATRIX starts out including the viewing transformation,
//  so we'll first take that off, then construct the modeling transformation,
//  then take *THAT* off, then put them back together.  Sounds a little like
//  Humpty-Dumpty, but...

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    //          eye point   center of view      up
    gluLookAt(0.0, 0.0, ep, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	//	save the viewing projection for later use
	glGetFloatv( GL_MODELVIEW_MATRIX, viewProj );

    //	Put identity onto modelview stack to start viewing transform
    glLoadIdentity();

	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':
				glRotatef( ANGLE, 1.0, 0.0, 0.0); break;
			case 'w':
			case 'p':
				glRotatef(-ANGLE, 1.0, 0.0, 0.0); break;
        	case 'a':
        	case 'k':
            	glRotatef( ANGLE, 0.0, 1.0, 0.0); break;
        	case 's':
        	case 'l':
            	glRotatef(-ANGLE, 0.0, 1.0, 0.0); break;
        	case 'z':
        	case 'm':
            	glRotatef( ANGLE, 0.0, 0.0, 1.0); break;
        	case 'x':
        	case ',':
            	glRotatef(-ANGLE, 0.0, 0.0, 1.0); break;
	    }
    }
    //  NOW apply the rest of the modeling transformation by POST-multiplying
    //  by the saved matrix, and then save that and put the identity back on
    //  the stack.  *whew*
    glMultMatrixf( saveState );
    glGetFloatv( GL_MODELVIEW_MATRIX, saveState );
    glLoadIdentity();
    
    //  Finally rebuild the overall modelview matrix by multiplying by the
    //  viewing transformation, saved from the reshape() function, and then
    //	by the modeling transformation created and saved above.
    glMultMatrixf( viewProj );
    glMultMatrixf( saveState );

	ch = ' ';	// ensure any residual rotation is killed
	// rotation finished, if any; now call display list containing molecule
	//	call list disabled because we need to change the model with time
	molecule();	//	glCallList(1);
	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,4.0*ep);
}

void keyboard(unsigned char key, int x, int y)
{
    ch = ' ';
    switch (key)
    {
        case 'q' :	// rotate around X; q,o = positive, w,p = negative
        case 'o' :
        case 'w' :
        case 'p' :
        case 'a' :	// rotate around Y; a,k = positive, s,l = negative
        case 'k' :
        case 's' :
        case 'l' :
        case 'z' :	// rotate around Z; z,m = positive, x,, = negative
        case 'm' :
        case 'x' :
        case ',' :
        	ch = key; break;	// save the input key to the global ch
        case 'f' :
        	ep -= ZOOM; break;	// move forward 'f' or back 'b'
        case 'b' :
        	ep += ZOOM; break;
    }
    glutPostRedisplay(); // perform display again
}

void animate(void)
{
#define TWOPI 6.28318

	t += 0.1;
	if (t > TWOPI) t -= TWOPI;
	sizeMult = (1.0+0.5*sin(t));
	alphaAdd = 0.2*cos(t);
	glutPostRedisplay();
}

void options_menu(int input)
{
	int i;

	if (input == Tag1) {
		if (atomSizes[0] > 0.1) {
			for (i=0; i<ATSIZE; i++)
				atomSizes[i] *= 0.9;
			//		define display list for the molecule
			glNewList(1, GL_COMPILE);
				molecule();
			glEndList();
			}	
		}
	if (input == Tag2) {
		if (atomSizes[0] < 2.0) {	// 2.0 is an artificial max
			for (i=0; i<ATSIZE; i++)
				atomSizes[i] *= 1.1;
			//		define display list for the molecule
			glNewList(1, GL_COMPILE);
				molecule();
			glEndList();
			}
		}
	if (input == Tag3) {
		if (atomColors[0][3] > 0.1) {
			for (i=0; i<ATSIZE; i++)
				atomColors[i][3] -= 0.1;
			//		define display list for the molecule
			glNewList(1, GL_COMPILE);
				molecule();
			glEndList();
			}
		}
	if (input == Tag4) {
		if (atomColors[0][3] <0.95) {
			for (i=0; i<ATSIZE; i++)
				atomColors[i][3] += 0.1;
			//		define display list for the molecule
			glNewList(1, GL_COMPILE);
				molecule();
			glEndList();
			}
		}
	if (input == Tag5) {	// move in
		ep *= .97;
		}
	if (input == Tag6) {	// move out
		ep *= 1.03;
		}
	if (input == Tag8) giveHelp();
	if (input == Tag9) giveAbout();
	ch = ' ';	// nullify any standing keyboard character
	glutPostRedisplay();
}

int lookUpName(int i)
{
//	return the index in the atomic tables for which the name in the
//	table matches the name in atoms[i]

	int j;
	
	for (j=0; j<ATSIZE; j++) {
		if (atomNames[j][0]==atoms[i].name[0] &&
		    atomNames[j][1]==atoms[i].name[1] )
				return(j);
		}
	return(-1);
}

void giveHelp(void)
{
	printf("help is given here\n");
}

void giveAbout(void)
{
	printf("about this program...\n");
}

//void main(int argc, char** argv)
void main(void)
{
//	Standard GLUT initialization
    glutInit(&argc,argv);
//	Hardcoded name isn't optimal but the test system doesn't do  stdio  right
    fname = "helvetane.pdb";
//	initialize for double buffering, RGB color, and depth tests
    glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
//	define the window parameters and set up callbacks
    glutInitWindowSize(500,500);
    glutInitWindowPosition(70,70);
    glutCreateWindow(mainName);
    glutKeyboardFunc(keyboard);
    glutDisplayFunc(display);
    glutReshapeFunc(reshape);
    glutIdleFunc(animate);

//	create the user menu to choose the visualization
	glutCreateMenu(options_menu);	// create menu and define entries
	glutAddMenuEntry("Atoms are smaller", Tag1);			// 1
	glutAddMenuEntry("Atoms are larger", Tag2);				// 2
	glutAddMenuEntry("Atoms are more transparent", Tag3);	// 3
	glutAddMenuEntry("Atoms are less transparent", Tag4);	// 4
	glutAddMenuEntry("Move in toward molecule", Tag5);		// 5
	glutAddMenuEntry("Move out away from molecule", Tag6);	// 6
	glutAddMenuEntry("---------------", Tag7);				// null choice
	glutAddMenuEntry("Help", Tag8);							// 6
	glutAddMenuEntry("About this program", Tag9);			// 7
	glutAttachMenuName(GLUT_RIGHT_BUTTON, "Options"); // give menu a name

//	input molecule data from file
//	fname = getFileName();
	separateName(fname);
	if (strcmp(ext,"mol")==0)
		readMolFile(fname);
	if (strcmp(ext,"pdb")==0)
		readPdbFile(fname);
	
    myinit();
    glutMainLoop();
}