/*
 * A Cascading Reaction-Diffusion simulator for PLY objects. (rdcascade)
 *
 *		for usage info, run program without any parameters
 *
 */

#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include "rply.h"

/* R-D constants and arrays */
GLfloat DA = 0.04199;
GLfloat DI = 0.229;
GLfloat S = 0.003;
GLfloat S2 = 0.02;

GLfloat *BETA;
GLfloat *cA;
GLfloat *cI;
GLfloat *oldCA;
GLfloat *oldCI;
GLfloat *activator;

/* 3d object data arrays */
GLfloat **vertices;
GLfloat **normals;
GLfloat **diffusion;
long **faces;
long **neighbours;

/* just some vars for extracting data from PLY (could be done better) */
long vertex=0, axis=0, face=0, fvertex=0, home=0, neighb=0;
GLfloat maxv = -100;
GLfloat minv = 100;

/* angles and stuff for rotating and zooming the object */
GLfloat xangle = 0.0;
GLfloat yangle = 0.0;
GLfloat zangle = 0.0;
GLfloat zoom = 1.5;
GLfloat atX = 0.0;
GLfloat atY = 0.0;

/* global for number of vertices and faces */
long nvertices,nfaces;
e_ply_type vertex_t, face_length_t, face_value_t;

/* GLUT window identifier and viewport size (so we see whole object)*/
static int win;
GLfloat view;

/*** extra variables for cascading */
int cascaded=0;
GLfloat convergence=0, oldconvergence;
int convcounter=0,convcount=0,convstandby=0;


/* input PLY filename */
char *filename;
char *nbsfile;

/* just some global variables for stuff that should be implemented better */
int plyHasNeighbs=0, plyHasNormals=0, plyHasDiff=0;
int runningRD = 1;


/*
 *
 * Display usage info
 *
 */
void printUsage()
{
	printf("Reaction-Diffusion Simulator (Cascading)\n\nusage:\trdcascade <plyfile> [-p <paramfile>]\n\n");
	printf("\tOptions:\n\t-p\t\t- Use R-D parameters in <paramfile>\n\t\t\t  [If not specified defaults are used]\n\n");
	printf("\tControls:\n\t\tRotate on X axis: 'a' and 'd'\n\t\tRotate on Y axis: 's' and 'w'\n\t\tRotate on Z axis: 'z' and 'c'\n\t\tZoom in and out: '+', '-'\n\t\tMove object: Left, Right, Up and Down arrows\n\t\tReset view: 'r'\n\t\tCascade: Spacebar\n\t\tPause/resume simulation: 'p'\n\t\tExport concentration levels: 'e'\n\t\tQuit: 'q'\n\n");
	printf("\tParameter file format:\n\t\t<DA value>\t\t- Activator diffusion constant (~0.04)\n\t\t<DI value>\t\t- Inhibitor diffusion constant (~0.2)\n\t\t<S value>\t\t- Initial S value (~0.002)\n\t\t<S2 value>\t\t- 2nd run S value (cascading) (~0.01)\n\n");
	printf("If the PLY input file does not contain additional data about vertex normals and neighbour cells, this will be generated before R-D simulation begins. In this case a PLY file with the extra data is stored at <plyfile>_nbs.ply, and this file should be used in future to save on computation time.\n\n");
}


/*
 *
 * The Reaction-Diffusion Jazz
 *
 */
void updateConc()
/*
 * updateConc(): calculate the new values of activator (cA) and inhibitor (cI) 
 * concentrations for each cell, using the R-D equations.
 */
{
	long i;
	GLfloat iR, aR, iD, aD;
	
	/* run through each triangular face (cell) */
	for (i=1;i<nfaces-1;i++){
		if (!((cascaded) && (cA[i] == 4.0) && (cI[i] == 4.0))) {
			/* reaction */
			iR=S*(16-oldCI[i]*oldCA[i]);
			aR=S*(oldCI[i]*oldCA[i]-oldCA[i]-BETA[i]);
			
			/* diffusion */
			iD=DI*((diffusion[i][0] * oldCI[neighbours[i][0]]) + (diffusion[i][1] * oldCI[neighbours[i][1]]) + (diffusion[i][2] * oldCI[neighbours[i][2]]) -3*oldCI[i]);
			aD=DA*((diffusion[i][0] * oldCA[neighbours[i][0]]) + (diffusion[i][1] * oldCA[neighbours[i][1]]) + (diffusion[i][2] * oldCA[neighbours[i][2]]) -3*oldCA[i]);
			
			/* diffusion without diffusion weightings 
			iD=DI*(oldCI[neighbours[i][0]] + oldCI[neighbours[i][1]] + oldCI[neighbours[i][2]] -3*oldCI[i]);
			aD=DA*(oldCA[neighbours[i][0]] + oldCA[neighbours[i][1]] + oldCA[neighbours[i][2]] -3*oldCA[i]);
			*/
			
			/* all put together */
			cI[i]+= iR + iD;
			cA[i]+= aR + aD;
			
			/* test for convergence */
			if (!cascaded) convergence += (aR+aD)*(aR+aD);
			
			if (cA[i]<0.0)
				cA[i]=0.0;
		}
	}
}

void saveConc()
/*
 * saveConc(): store old chemical concentration levels for use in the next iteration
 * of updateConc.
 */
{
	long i;
	for (i=0;i<nfaces;i++){
		oldCA[i]=cA[i];
		oldCI[i]=cI[i];
		activator[i]=cA[i]/6.0;
	}
}

void doCascade()
/*
 * doCascade(): cascade, set new value for S parameter, and give dark cells activator
 * and inhibitor levels of 4.0.
 */
{
	long i;
	
	printf("Cascading...\n");
	cascaded=1;
	S=S2;
	for (i=0;i<nfaces;i++)
		if (cA[i] < 2.4) cA[i] = cI[i] = 4.0;
}


/*
 *
 * Output the concentration levels to a file. This is for use in rdviewer.
 *
 */
void outputConc(void)
{
	long i;
	char outfilename[20];
	FILE *fp;
	
	printf("Enter filename for pattern export:\n\t");
	scanf("%s",outfilename);
	
	if (fp = fopen(outfilename,"w")) {
		printf("Exporting concentration levels to %s... ",outfilename);
		if ((!plyHasNeighbs) || (!plyHasNormals))
			fprintf(fp, "%s\n", nbsfile);
		else
			fprintf(fp, "%s\n", filename);
		
		for (i=0;i<nfaces;i++)
			fprintf(fp, "%f\n", activator[i]);
		
		fclose(fp);
		printf("Done\n");
	} else
		printf("Error writing to file pattern.cnc.\n");
}


/*
 *
 * GLUT routines, callbacks and whatnot
 *
 */
void update(void)
/*
 * update(): GLUT callback function. called every refresh. We use it to run R-D
 * simulation and test for convergence, then call for a redisplay.
 */
{	 
	if (runningRD) {
		updateConc();
		saveConc();
		
		if (!cascaded) {
			if (convcounter++ >= convcount) {
				//printf("convergence: %f\n", (float)convergence);
				convcounter=0;
				/* test if change is less than previous iterations (need threshold because it drops slightly at beginning */
				if ((oldconvergence - convergence) >= 2) {
					convstandby=1;
					printf("Change is dropping, standby to cascade...\n");
				}
				/* see if we are starting to rise again after dropping a bit, if so trigger cascade */
				if ((convstandby) && (oldconvergence < convergence)) doCascade();
				oldconvergence=convergence;
				convergence=0.0;
			}
		}
	}
	
	glutPostRedisplay();
}

void face3(long face)
/*
 * face3(face): Render the triangular polygon indexed to in the faces array by
 * input parameter face. Colour is set according to activator[] value for the face.
 * The glMaterial and glNormal calls are required for lighting.
 */
{
	GLfloat colour[] = {activator[face], activator[face], activator[face]}; // - colour - {activator[face]/1.4, activator[face]/2.0, 0};
	GLfloat material[] = {activator[face], activator[face], activator[face], 1.0}; // - colour - {activator[face]/1.4, activator[face]/2.0, 0, 1.0};
	glColor3fv(colour);
	glMaterialfv(GL_FRONT, GL_SPECULAR, material);
	glMaterialfv(GL_FRONT, GL_AMBIENT, material);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, material);
	glBegin(GL_POLYGON);
		glNormal3fv(normals[faces[face][0]]);
		glVertex3fv(vertices[faces[face][0]]);
		glNormal3fv(normals[faces[face][1]]);
		glVertex3fv(vertices[faces[face][1]]);
		glNormal3fv(normals[faces[face][2]]);
		glVertex3fv(vertices[faces[face][2]]);
	glEnd();
}

void display(void)
/*
 * display(): GLUT display callback, it sets up the view, and rotation of the scene,
 * then it renders all the polygons, and swaps buffers (we're using double-buffering).
 */
{
	long i;
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	/* Set up the modelview transformation. */
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//gluLookAt(6.0, 4.0, 14.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	gluLookAt(0, 0, zoom * view, atX, atY, 0.0, 0.0, 1.0, 0.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);

	
	for (i=0;i<nfaces;i++)
		face3(i);
	
	glutSwapBuffers();
}

void reshape(int w, int h)
/*
 * reshape(w, h): this is needed for when the window is resized, it keeps the
 * shape and perspective of the viewport correct.
 */
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(40.0, (GLdouble)w/(GLdouble)h, 0.01, 2000.0);
	glViewport(0, 0, (GLsizei)w, (GLsizei)h);
}

void keyb(unsigned char key, int x, int y)
/*
 * keyb(key, x, y): GLUT keyboard callback, we use it to rotate the scene, zoom in
 * and out, perform other user functions, and finally to free all the malloced memory.
 */
{
	long i;
	
	switch(key) {
		case 'a':
			yangle -= 10;
			break;
		case 'd':
			yangle += 10;
			break;
		case 'w':
			xangle -= 10;
			break;
		case 's':
			xangle += 10;
			break;
		case 'c':
			zangle -= 10;
			break;
		case 'z':
			zangle += 10;
			break;
		case '+':
		case '=':
			zoom-=0.2;
			break;
		case '-':
			zoom+=.2;
			break;
		case 'e':
			outputConc();
			break;
		case 'r':
			xangle = yangle = zangle = atX = atY = 0.0;
			zoom =1.5;
			break;
		case 'p':
			if (!runningRD) {
				runningRD = 1;
				printf("R-D simulation running...\n");
			} else {
				runningRD = 0;
				printf("R-D simulation paused. Press 'p' to resume.\n");
			}
			break;
		case ' ':
			if (!cascaded) doCascade();	
			break;
		case 'q':
		case 'Q':	
			printf("Freeing memory\n");
			for(i=0;i<nvertices;i++) {
				free(vertices[i]);
				free(normals[i]);
			}
			free(vertices);
			free(normals);
			for(i=0;i<nfaces;i++) {
				free(faces[i]);
				free(neighbours[i]);
				free(diffusion[i]);
			}
			free(faces);
			free(neighbours);
			free(diffusion);
			
			free(BETA);
			free(cA);
			free(cI);
			free(oldCA);
			free(oldCI);
			free(activator);
			
			printf("All done.\n");
			glutDestroyWindow(win);
			exit(0);
	}
}

void special(int key, int x, int y)
/*
 * special(key, x, y): GLUT special keyboard callback, used for arrow keys to move scene.
 */
{
	switch(key) {
		case GLUT_KEY_RIGHT:
			atX=atX-0.02*view;
			break;
		case GLUT_KEY_LEFT:
			atX=atX+0.02*view;
			break;
		case GLUT_KEY_UP:
			atY=atY-0.02*view;
			break;
		case GLUT_KEY_DOWN:
			atY=atY+0.02*view;
			break;
		default:
			return;
	}
}
/*
 *
 * All the PLY model handling routines
 *
 */
static int vertex_cb(p_ply_argument argument)
/*
 * vertex_cb(argument): RPly vertex callback. This is called everytime RPly encounters
 * a vertex as it parses the PLY file. We store the vertex coordinates and (rather crudely)
 * find the max and min coordinates, for use in setting the viewport correctly.
 */
{
	long eol;
	p_ply_property property;
	
	ply_get_argument_user_data(argument, NULL, &eol);
	ply_get_argument_property(argument, &property, NULL, NULL);
	ply_get_property_info(property, NULL, &vertex_t, NULL, NULL);
	
	vertices[vertex][axis] = (GLfloat)ply_get_argument_value(argument);

	if (vertices[vertex][axis] > maxv) maxv = vertices[vertex][axis];
	if (vertices[vertex][axis] < minv) minv = vertices[vertex][axis];
	
	axis++;
	if (eol) {
		vertex++;
		axis=0;
	}
	return 1;
}

static int normal_cb(p_ply_argument argument)
/*
 * normal_cb(argument): RPly normal callback. Store the normal for each vertex.
 */
{
	long eol;
	
	ply_get_argument_user_data(argument, NULL, &eol);
	normals[vertex-1][axis] = (GLfloat)ply_get_argument_value(argument);
	
	axis++;
	if (eol) {
		//vertex++;
		axis=0;
	}
	return 1;
}

static int face_cb(p_ply_argument argument)
/*
 * face_cb(argument): RPly face callback. Store vertex indices for each face.
 */
{
	long length, value_index;
	p_ply_property property;
	
	ply_get_argument_property(argument, &property, &length, &value_index);
	ply_get_property_info(property, NULL, NULL, &face_length_t, &face_value_t);
  
  switch (value_index) {
		case 0:
		case 1: 
			faces[face][fvertex++] = (long)ply_get_argument_value(argument);
			break;
		case 2:
			faces[face][fvertex++] = (long)ply_get_argument_value(argument);
			face++;
			fvertex=0;
			break;
		default: 
			break;
	}
	return 1;
}

static int neighbours_cb(p_ply_argument argument)
/*
 * neighbours_cb(argument): RPly neighbour callback. Retrieves the neighbour indices.
 * home variable indexs current cell, neighb indexs its neighbours. Bit of a cludge to
 * cover the case where neighbours ar found but not diffusion weightings.
 */
{
	long value_index;
	
	ply_get_argument_property(argument, NULL, NULL, &value_index);
  
  switch (value_index) {
		case 0:
		case 1: 
			neighbours[home][neighb++] = (long)ply_get_argument_value(argument);
			break;
		case 2:
			neighbours[home][neighb++] = (long)ply_get_argument_value(argument);
			if (plyHasDiff == 0) home++;
			neighb=0;
			break;
		default:
			break;
	}
	return 1;
}

static int diffusion_cb(p_ply_argument argument)
/*
 * diffusion_cb(argument): RPly diffusion weightings callback.
 */
{
	long value_index;
	
	ply_get_argument_property(argument, NULL, NULL, &value_index);
  
  switch (value_index) {
		case 0:
		case 1: 
			diffusion[home][neighb++] = (GLfloat)ply_get_argument_value(argument);
			break;
		case 2:
			diffusion[home][neighb++] = (GLfloat)ply_get_argument_value(argument);
			home++;
			neighb=0;
			break;
		default:
			break;
	}
	return 1;
}

int outputPLY(char *filename)
/*
 * outputPLY(filename): calls all the RPly output functions in correct order to write
 * a PLY file with all additional data (neighbours, normals, diffusion weightings).
 */
{
	int i,j;
	p_ply ply;

	ply = ply_create(filename, PLY_ASCII, NULL);
	
	/* add vertex element header info */
	if (!ply_add_element(ply, "vertex", nvertices)) printf("cant add vertex element\n");
	ply_add_scalar_property(ply, "x", vertex_t);
	ply_add_scalar_property(ply, "y", vertex_t);
	ply_add_scalar_property(ply, "z", vertex_t);
	ply_add_scalar_property(ply, "xnorm", vertex_t);
	ply_add_scalar_property(ply, "ynorm", vertex_t);
	ply_add_scalar_property(ply, "znorm", vertex_t);

	/* add faces element info */
	if (!ply_add_element(ply, "face", nfaces)) printf("cant add face element\n");
	ply_add_list_property(ply, "vertex_indices", face_length_t, face_value_t);
	
	/* add neighbours element info */
	if (!ply_add_element(ply, "neighbours", nfaces)) printf("cant add neighbours element\n");
	ply_add_list_property(ply, "vertex_indices", face_length_t, face_value_t);
	ply_add_list_property(ply, "diffusion_weighting", face_length_t, vertex_t);
	
	/* add comment saying created by our prog and added neighbours */
	ply_add_comment(ply, "Extra neighbour and vertex normal data generated for Reaction-Diffusion simulations");
	
	/* write the header */
	if (!ply_write_header(ply)) printf("cant write header\n");
	
	/* write the properties one by one */
	for (i=0;i<nvertices;i++) {
		ply_write(ply, vertices[i][0]);
		ply_write(ply, vertices[i][1]);
		ply_write(ply, vertices[i][2]);
		ply_write(ply, normals[i][0]);
		ply_write(ply, normals[i][1]);
		ply_write(ply, normals[i][2]);
	}
	
	for (i=0;i<nfaces;i++) {
		ply_write(ply, 3);
		ply_write(ply, faces[i][0]);
		ply_write(ply, faces[i][1]);
		ply_write(ply, faces[i][2]);
	}
	
	for (i=0;i<nfaces;i++) {
		ply_write(ply, 3);
		ply_write(ply, neighbours[i][0]);
		ply_write(ply, neighbours[i][1]);
		ply_write(ply, neighbours[i][2]);
		ply_write(ply, 3);
		ply_write(ply, diffusion[i][0]);
		ply_write(ply, diffusion[i][1]);
		ply_write(ply, diffusion[i][2]);
	}
	
	ply_close(ply);
	printf("PLY file with neighbour and normal data written to %s\n",filename);
}

/*
 *
 * Neighbour finding and normal vector algorithms. Damn slow for big models.
 *
 */
void findNeighbours(void)
{
	long i,j,k,p,v,v1,v2,v3,v1st,v2nd,matches;
	GLfloat diff, norm;
	
	printf("Finding neighbours (may take some time)...\n");
	for(i=0;i<nfaces;i++) {
		v1 = faces[i][0];
		v2 = faces[i][1];
		v3 = faces[i][2];
		p = 0;
		for(j=0;j<nfaces;j++) {
			matches = 0;
			for(k=0;k<3;k++) {
				v = faces[j][k];
				if ((v == v1) || (v == v2) || (v == v3)) {
					if (matches == 0) v1st = v;
					if (matches == 1) v2nd = v;
					matches++;
				}
			}
			if (matches == 2) {
				neighbours[i][p] = j;
				diff = sqrt(((vertices[v1st][0] - vertices[v2nd][0]) * (vertices[v1st][0] - vertices[v2nd][0])) + ((vertices[v1st][1] - vertices[v2nd][1]) * (vertices[v1st][1] - vertices[v2nd][1])) + ((vertices[v1st][2] - vertices[v2nd][2]) * (vertices[v1st][2] - vertices[v2nd][2])));
				diffusion[i][p++] = diff;
			}
		}
		norm = 3/(diffusion[i][0] + diffusion[i][1] + diffusion[i][2]);
		diffusion[i][0] = norm * diffusion[i][0];
		diffusion[i][1] = norm * diffusion[i][1];
		diffusion[i][2] = norm * diffusion[i][2];
	}
	printf("Done.\n");
}


GLfloat faceNorm(long face, int axis)
{
	/* [v1 - v2] x [v2 - v3] */
	
	int i;
	GLfloat lhs[3];
	GLfloat rhs[3];
	
	for(i=0;i<3;i++) {
		lhs[i] = vertices[faces[face][0]][i] - vertices[faces[face][1]][i];
		rhs[i] = vertices[faces[face][1]][i] - vertices[faces[face][2]][i];
	}
	
	switch (axis) {
		case 0:
			return (lhs[1] * rhs[2] - lhs[2] * rhs[1]);
		case 1:
			return (lhs[2] * rhs[0] - lhs[0] * rhs[2]);
		case 2:
			return (lhs[0] * rhs[1] - lhs[1] * rhs[0]);
	}
}


void findNormals(void)
{
	long i,j;
	int count=0;
	GLfloat normal[3];
	GLfloat magnitude;
	
	printf("Calculating normals (may take some time)...\n");
	
	for(i=0;i<nvertices;i++) {
		count=0;
		normal[0] = normal[1] = normal[2] = 0.0;
		for(j=0;j<nfaces;j++) {
			if ((faces[j][0] == i) || (faces[j][1] == i) || (faces[j][2] == i)) {
				count++;
				normal[0] += faceNorm(j, 0);
				normal[1] += faceNorm(j, 1);
				normal[2] += faceNorm(j, 2);
			}
		}
		magnitude = sqrt((normal[0] * normal[0]) + (normal[1] * normal[1]) + (normal[2] * normal[2]));
		if (magnitude == 0) {
			normals[i][0] = normal[0];
			normals[i][1] = normal[1];
			normals[i][2] = normal[2];
		} else { 
			normals[i][0] = normal[0]/magnitude;
			normals[i][1] = normal[1]/magnitude;
			normals[i][2] = normal[2]/magnitude;
		}
	}
	printf("Done.\n");
}


/*
 *
 * Lots of initialisation stuff. Loading the object, etc.
 *
 */
void init(void)
{
	long i,j;
	p_ply ply;
	GLfloat mat_specular[] = { 0.8, 0.8, 0.8, 1.0 };
	GLfloat mat_shininess[] = { 10.0 };
	GLfloat light_ambient[] = {0.1, 0.1, 0.1, 0.5};
  GLfloat light_position[] = { 1.0, 1.0, 1.0, 0.0 };
	 
	/* some GL initialisation */
	glClearColor(0.0, 0.0, 0.5, 1.0);
	glEnable(GL_DEPTH_TEST);
	
	/* load PLY file */
	ply = ply_open(filename, NULL);
	if ((!ply) || (!ply_read_header(ply))) {
		printf("Error opening PLY file or reading header.\n");
   	exit;
	}
  
	/* set callbacks for PLY reading and get number of vertices */
  nvertices = ply_set_read_cb(ply, "vertex", "x", vertex_cb, NULL, 0);
  ply_set_read_cb(ply, "vertex", "y", vertex_cb, NULL, 0);
  ply_set_read_cb(ply, "vertex", "z", vertex_cb, NULL, 1);
	
	/* test if PLY file contains vertex normals, if so set callbacks */
	if (ply_set_read_cb(ply, "vertex", "xnorm", normal_cb, NULL, 0)) {
		printf("This PLY file defines normals to vertices. Good work.\n");
		plyHasNormals = 1;
	}
	ply_set_read_cb(ply, "vertex", "ynorm", normal_cb, NULL, 0);
	ply_set_read_cb(ply, "vertex", "znorm", normal_cb, NULL, 1);
	
	/* set face callback and get number of faces */
  nfaces = ply_set_read_cb(ply, "face", "vertex_indices", face_cb, NULL, 0);
	/* try again with element label "vertex_index" if it fails */
  if (nfaces == 0) nfaces = ply_set_read_cb(ply, "face", "vertex_index", face_cb, NULL, 0);
	
	/* test for neighbour and diffusion weighting data and register callbacks */
	if (ply_set_read_cb(ply, "neighbours", "vertex_indices", neighbours_cb, NULL, 0)) {
		printf("The PLY file contains neighbour data. Hoorah!\n");
		plyHasNeighbs = 1;
	}
	if (ply_set_read_cb(ply, "neighbours", "diffusion_weighting", diffusion_cb, NULL, 0))
		plyHasDiff = 1;
	printf("Vertices: %ld\nFaces: %ld\n", nvertices, nfaces);
   
  /* mallocate memory for arrays */
  /* 3d structure data */
  vertices = malloc(nvertices * sizeof(GLfloat *));
	normals = malloc(nvertices * sizeof(GLfloat *));
  for(i=0;i<nvertices;i++) {
  	vertices[i] = malloc(3 * sizeof(GLfloat));
		normals[i] = malloc(3 * sizeof(GLfloat));
	}
  faces = malloc(nfaces * sizeof(long *));
  neighbours = malloc(nfaces * sizeof(long *));
	diffusion = malloc(nfaces * sizeof(GLfloat *));
  for(i=0;i<nfaces;i++) {
		faces[i] = malloc(3 * sizeof(long));
  	neighbours[i] = malloc(3 * sizeof(long));
		diffusion[i] = malloc(3 * sizeof(GLfloat));
  }
   
  /* r-d arrays */
  BETA = malloc(nfaces * sizeof(GLfloat));
  cA = malloc(nfaces * sizeof(GLfloat));
  cI = malloc(nfaces * sizeof(GLfloat));
  oldCA = malloc(nfaces * sizeof(GLfloat));
  oldCI = malloc(nfaces * sizeof(GLfloat));
  activator = malloc(nfaces * sizeof(GLfloat));
  
  printf("Done mallocating memory.\n");
  
	/* a little bit of error checking */
  if (!ply_read(ply)) printf("Error extracting vertices.\n");
  ply_close(ply);
  
  /* making sure we get everything in view (very basic may need to be improved) */
  view = maxv - minv;
  //printf("max: %.2f %.2f\nView: %.2f\n",maxv,minv,view);
	
	/* setup for lighting */
	glShadeModel (GL_SMOOTH);
	glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
	light_position[0] = light_position[1] = light_position[2] = -view;
  glLightfv(GL_LIGHT0, GL_POSITION, light_position);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
	
	/* if we dont have neighbours or normals find them then save the ply file */
	if ((!plyHasNeighbs) || (!plyHasNormals)) {
		printf("PLY file does not contain neighbours and/or normals data.\n");
		if (!plyHasNeighbs) findNeighbours();
		if (!plyHasNormals) findNormals();
		
		/* set filename for outputting with nbs data (this is stupidly complicated) */
		nbsfile = (char *)calloc(strlen(filename)+5, sizeof(char));
		strncpy(nbsfile,filename,strlen(filename)-4);
		strcat(nbsfile, "_nbs.ply");
		
		/* output a PLY file */
		outputPLY(nbsfile);
	}
	
	/* initialise r-d arrays and save as old concetrations */
	convcount=ceil(nfaces/1000);
	//printf("convcount: %d\n",convcount);
	for(i=0;i<nfaces;i++) {
		cA[i] = cI[i] = 4.0;
		BETA[i] = (0.1*rand()/(RAND_MAX+1.0))+11.95;
	}
	saveConc();
	
	printf("Initialisation done.\n");
	printf("Controls:\n\tRotate on X axis: 'a' and 'd'\n\tRotate on Y axis: 's' and 'w'\n\tRotate on Z axis: 'z' and 'c'\n\tZoom in and out: '+', '-'\n\tMove object: Left, Right, Up and Down arrows\n\tReset view: 'r'\n\tCascade: Spacebar\n\tPause/resume simulation: 'p'\n\tExport concentration levels: 'e'\n\tQuit: 'q'\n");
}


int main(int argc, char** argv)
{
	FILE *fp;
	
	/* print out usage information */
	if ((argc == 1) || (argv[1][0] == '-')) {
		printUsage();
		return 0;
	}
	
	/* get PLY filename from command-line */
	filename = argv[1];
	printf("%s\n",filename);
	
	/* test if a parameter file is given, otherwise use defaults */
	if (argv[2] != NULL) {
		switch (argv[2][1]) {
			case 'p':
				if (fp = fopen(argv[3],"r")) {
					fscanf(fp, "%f\n%f\n%f\n%f", &DA, &DI, &S, &S2);
					printf("Taking parameters from %s:\n", argv[3]);
					fclose(fp);
					break;
				} else
					printf("Parameter file %s not found. ",argv[3]);
			default:
				printf("Using default parameters:\n");
		}
	} else {
		printf("Using default parameters:\n");
	}
	printf("\tDA: %f\n\tDI: %f\n\tS: %f\n\tS2: %f\n",DA,DI,S,S2);
	
	/* intialise GLUT and register callbacks, then enter main loop */
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGB);
	glutInitWindowSize(250,250);
	glutInitWindowPosition(300,200);
	win = glutCreateWindow("Cascading R-D Objects");
	init();
	glutDisplayFunc(display);
	glutSpecialFunc(special);
	glutReshapeFunc(reshape);
	glutIdleFunc(update);
	glutKeyboardFunc(keyb);
	reshape(250,250);
	glutMainLoop();
	
	return 0;
}