/*
   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) February 2001, Steve Cunningham & Ken Brown
*/

#include <GL/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 FILE_MAX 128 // maximum length of a file name for this program
// for determineFileType
#define FILE_TYPE_UNKNOWN 0
#define FILE_PDB 1
#define FILE_MOL 2

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

// data for storing molecular description

int		natoms, nbonds;
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
enum inputTags { Tag1=1, Tag2, Tag3, Tag4, Tag5, Tag6, Tag7, Tag8 };
float	ep; // eye point variable
char fname[FILE_MAX],
	 mainName[FILE_MAX],
	 ext[4];

#define ROT_AMT 2.0
GLfloat rotx = 0.0, roty = 0.0, rotz = 0.0;

typedef GLfloat point3[3];

// texture data
float D1, D2;
float texParms[4];
static GLuint texName;
float ramp[256][3];

// function prototypes
int determineFileType(char *fname);
void readMolFile(char *);
void readPdbFile(char *);
void makeRamp(void);
void hsv2rgb(float, float, float, float *, float *, float *);
void myinit(void);
void molecule(void);
void display(void);
void reshape(int, int);
void keyboard(unsigned char, int, int);
void options_menu(int);
void giveHelp(void);
int  lookUpName(int);

int determineFileType(char *fname) {
	int i;
	// first search for the period indicating the extension
	// if we can't find it we'll return unknown file type
	for(i = 0; i < FILE_MAX-3; i++) {
		if(fname[i] == '.') {
			i++;

			// check for our two supported types
			if(fname[i] == 'p' && fname[i+1] == 'd' && fname[i+2] == 'b')
				return FILE_PDB;
			if(fname[i] == 'm' && fname[i+1] == 'o' && fname[i+2] == 'l')
				return FILE_MOL;

			// it's some other extension
			return FILE_TYPE_UNKNOWN;
		}
	}

	// if we didn't find an extension, it's not a supported type
	return FILE_TYPE_UNKNOWN;
}

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 makeRamp(void)
{
	int i;
	float h, s, v, r, g, b;

	// Make color ramp for 1D texture: starts at 0, ends at 240, 256 steps
	for (i=0; i<256; i++) {
		h = (float)i*240.0/255.0;
		s = 1.0; v = 1.0;
		hsv2rgb( h, s, v, &r, &g, &b );
		ramp[i][0] = r;ramp[i][1] = g;ramp[i][2] = b;
		}
}

void hsv2rgb( float h, float s, float v, float *r, float *g, float *b)
{
//	conversion from Foley et.al., fig. 13.4
//	don't worry about all possible cases because we know s = v = 1.0
//	and 0 <= h <=240

	float f, p, q, t;
	int   k;

	h = h/60.0;
	k = (int)h;
	f = h - (float)k;
	p = v * (1.0 - s);
	q = v * (1.0 - (s * f));
	t = v * (1.0 - (s * (1.0 - f)));
	switch (k) {
		case 0: *r = v; *g = t; *b = p; break;
		case 1: *r = q; *g = v; *b = p; break;
		case 2: *r = p; *g = v; *b = t; break;
		case 3: *r = p; *g = q; *b = v; break;
		case 4: *r = t; *g = p; *b = v; break;
		case 5: *r = v; *g = p; *b = q; break;
	}
}

// initialize OpenGL and all necessary global variables
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.2, 0.2, 0.2, 0.0 );
	glShadeModel(GL_SMOOTH); // use gourand shading
	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
	
//	set up texture information
	makeRamp();
	glPixelStorei( GL_UNPACK_ALIGNMENT, 1 );
	glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL );
	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP );
 	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR );
	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR );
	glTexImage1D( GL_TEXTURE_1D, 0, 3, 256, 0, GL_RGB, GL_FLOAT, ramp );
	glEnable( GL_TEXTURE_GEN_S );
	glEnable( GL_TEXTURE_1D );

//	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;
	
//	define display list for the molecule
	glNewList(1, GL_COMPILE);
		molecule();
	glEndList();
}

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, startindex, endindex, bondtype;
	GLUquadric *atomSphere;

    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
    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
			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, atomColors[j] );
			glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, atomColors[j] );
			glTranslatef(atoms[i].x, atoms[i].y, atoms[i].z);
			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

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    
//	Now build the texture information before the transformations are made
    glEnable(GL_TEXTURE_1D);
    glEnable( GL_TEXTURE_GEN_S );
    glTexGeni( GL_S, GL_TEXTURE_GEN_MODE, GL_EYE_LINEAR );
    D1 = ep - 6.3; D2 = ep + 6.3;
    texParms[0] = texParms[1] = 0.0;
    texParms[2] = -1.0/(D2-D1);
    texParms[3] = -D1/(D2-D1);
    glTexGenfv( GL_S, GL_EYE_PLANE, texParms);
    glBindTexture(GL_TEXTURE_1D, texName);

    //          eye point   center of view      up
    gluLookAt(0.0, 0.0, ep, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);

	glRotatef(rotx, 1.0,0.0,0.0);
	glRotatef(roty, 0.0,1.0,0.0);
	glRotatef(rotz, 0.0,0.0,1.0);

	// rotation finished, if any; now call display list containing molecule
	glCallList(1);
	glutSwapBuffers();
}

void reshape(int w,int h)
{
    glViewport(0,0,(GLsizei)w,(GLsizei)h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0,(GLfloat)w/(GLfloat)h,1.0,4.0*ep);
}

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

void options_menu(int input)
{
	int i;

	switch(input) {
	case 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();
		}	
		break;
	case 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();
		}
		break;
	case 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();
		}
		break;
	case 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();
		}
		break;
	case Tag5: ep *= .97; break;	// move in
	case Tag6: ep *= 1.03; break;	// move out
	case Tag8: giveHelp(); break;
	}
	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("*\tKeyboard controls are:\n");
	printf("\t\t - keys q/w/o/p rotate around X-axis,\n");
	printf("\t\t - keys a/s/k/l rotate around Y-axis,\n");
	printf("\t\t - keys z/x/m/, rotate around Z-axis,\n");
	printf("\t\t - keys f/b move the clipping plane.\n\n\n");
}

int main(int argc, char** argv)
{
//	Standard GLUT initialization
    glutInit(&argc,argv);
//	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);

//	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
	glutAttachMenu(GLUT_RIGHT_BUTTON); // give menu a name

//	input molecule data from file
	printf("enter name of molecule file to be displayed > ");
	scanf("%s",fname);
	printf("%s\n",fname);
	switch(determineFileType(fname)) {
	case FILE_PDB:
		readPdbFile(fname);
		break;
	case FILE_MOL:
		readMolFile(fname);
		break;
	case FILE_TYPE_UNKNOWN:
		fprintf(stderr,"File type unknown.\n" );
		exit(0);
	}
	
    myinit();
    glutMainLoop();

	return 0;
}