/*****************************************************************************
 *                                                                           *
 * flow - Flow is an interactive program for displaying arcascii files,      *
 *        computing flow direction and accumulation on a grid                *
 *                                                                           * 
 * run as: flow                                                              *
 *                                                                           *
 * Author: Jonathan Todd                                                     *
 *                                                                           *
 *****************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include "flow.h"

#ifdef __APPLE__
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif

#define LINE_LENGTH 80
#define BUCKETS 5

//Menu constants
#define EXIT 1
#define FILL 2
#define FLY 3

#define RED 10
#define YELLOW 11
#define MAGENTA 12
#define GREEN 13
#define CYAN 14
#define BLUE 15
#define WHITE 16
#define GRAY 17
#define BLACK 18
#define SAND 19
#define FGREEN 20
#define LGREEN 21
#define DBROWN 22
#define BROWN 23
#define DBLUE 24
#define PURPLE 25

#define GRAYSCALE 4
#define RAINBOW 5
#define CUSTOM 6
#define FLOWACCU 7

GLfloat red[3] = {1.0, 0.0, 0.0};
GLfloat green[3] = {0.0, 1.0, 0.0};
GLfloat blue[3] = {0.0, 0.0, 1.0};
GLfloat black[3] = {0.0, 0.0, 0.0};
GLfloat white[3] = {1.0, 1.0, 1.0};
GLfloat gray[3] = {0.5, 0.5, 0.5};
GLfloat yellow[3] = {1.0, 1.0, 0.0};
GLfloat magenta[3] = {1.0, 0.0, 1.0};
GLfloat cyan[3] = {0.0, 1.0, 1.0};
GLfloat sand[3] = {251.0/255, 224.0/255, 142.0/255};
GLfloat fgreen[3] = {27.0/255, 96.0/255, 17.0/255};
GLfloat lgreen[3] = {46.0/255, 165.0/255, 29.0/255};
GLfloat dbrown[3] = {88.0/255, 67.0/255, 27.0/255};
GLfloat brown[3] = {112.0/255, 112.0/255, 56.0/255};
GLfloat dblue[3] = {51.0/255, 51.0/255, 102.0/255};
GLfloat purple[3] = {128.0/255, 0.0/255, 128.0/255};

GLfloat pos[3] = {0,0,0};
GLfloat theta[3] = {0,0,0};

GLint fillmode = 0;
Grid *grid;             // grid structure for display function
Grid *color_grid;       // grid structure for maping color onto display grid
int resolution = 64;    // 1 is best resolution, 2 is 1/2 the resolution of 1..
int colorStyle = 0;     // Specifies the color to display
float maxZ = 1.0;       // specifies the max value for the z coordinate
int window_width = 500;
int window_height = 500;
int bu = BUCKETS;       // number of buckets
GLfloat *color[BUCKETS+1] = {blue,sand,fgreen,brown,dbrown,white};  // color scheme
float delta = 2.0;      // increment to change theta and pos
float pitch = 0.0;
float yaw = 0.0;
float roll = 0.0;
int window_open = 0;    // glut window open y/n?
float angle=0.0;        // angle for fly through adjustment
int fly = 0;            // fly through mode


// main
//
int main(int argc, char** argv) {

  //---> Setup glut window
  // open a window
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
  glutInitWindowSize(window_width, window_height);
  glutInitWindowPosition(100,100);
  glutCreateWindow(argv[0]);
  glEnable(GL_DEPTH_TEST);

  // camera is at (0,0,0) looking along negative y axis
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(60, 1, 0.01, 10.0); 
  glMatrixMode(GL_MODELVIEW);
  
  // Set background color to black
  glClearColor(0, 0, 0, 0); 
    
  // callback functions
  glutMotionFunc(mouse);
  glutDisplayFunc(display); 
  glutKeyboardFunc(keypress);
  //glutReshapeFunc(reshape);
  glutSpecialFunc(specialkey);
  glutIdleFunc(renderScene);

  // creat menu
  int mainm,colorm,customm, bucketm[bu+1],x;
  
  for(x = 0; x <= bu; x++){
    bucketm[x] = glutCreateMenu(main_menu);
    glutAddMenuEntry("Red",RED+(x*100));
    glutAddMenuEntry("Yellow",YELLOW+(x*100));
    glutAddMenuEntry("Magenta",MAGENTA+(x*100));
    glutAddMenuEntry("Green",GREEN+(x*100));
    glutAddMenuEntry("Cyan",CYAN+(x*100));
    glutAddMenuEntry("Blue",BLUE+(x*100));
    glutAddMenuEntry("White",WHITE+(x*100));
    glutAddMenuEntry("Gray",GRAY+(x*100));
    glutAddMenuEntry("Black",BLACK+(x*100));
    glutAddMenuEntry("Sand",SAND+(x*100));
    glutAddMenuEntry("Forest Green",FGREEN+(x*100));
    glutAddMenuEntry("Light Green",LGREEN+(x*100));
    glutAddMenuEntry("Dark Brown",DBROWN+(x*100));
    glutAddMenuEntry("Brown",BROWN+(x*100));
    glutAddMenuEntry("Dark Blue",DBLUE+(x*100));
    glutAddMenuEntry("Purple",PURPLE+(x*100));
  }
  
  customm = glutCreateMenu(main_menu);
  glutAddMenuEntry("Use custom",CUSTOM);
  for(x = 0; x <= bu; x++){
    char num[8];
    sprintf(num,"Bucket: %d",x);
    glutAddSubMenu(num,bucketm[x]);
  }

  colorm = glutCreateMenu(main_menu);
  glutAddMenuEntry("Gray Scale",GRAYSCALE);
  glutAddMenuEntry("Rainbow",RAINBOW);
  glutAddMenuEntry("Flow-coloring",FLOWACCU);  
  glutAddSubMenu("Custom",customm);
  
  mainm = glutCreateMenu(main_menu);
  glutAddMenuEntry("Fill/Outline", FILL);
  glutAddMenuEntry("Fly Through", FLY);
  glutAddSubMenu("Color",colorm);
  glutAddMenuEntry("Quit", EXIT);
  glutAttachMenu(GLUT_LEFT_BUTTON);

  //---> Start interactive prompt
  interact();

  return 0;
}


// interact - creates interactive prompt with user
//
void interact(){
// Interactive prompt
  while(1){
    char line[LINE_LENGTH],
    str1[LINE_LENGTH]="",
    str2[LINE_LENGTH]="",
    str3[LINE_LENGTH]="";

    printf("flow>");
    if( fgets(line, LINE_LENGTH, stdin) != NULL){
      sscanf(line,"%s %s %s",str1,str2,str3);
      // help
      if(strcmp(str1,"help")==0){
	printf("\n\'flow\' provides and interactive way to compute flow "\
	       "direction, flow accumumlation, display arcascii grids and "\
	       "fly through the terrain.\n\n"\
	       "Commands:\n"\
	       "flow>flowdir <input_file> <output_file>\n"\
	       "flow>flowaccu <input_file> <output_file>\n"\
	       "flow>display <input_file> [input_color_file]\n"\
	       "flow>info <input_file>\n\n");
      }

      // exit
      else if(strcmp(str1,"exit")==0){
	printf("Done.\n");
	exit(0);
      }
      // info
      else if(strcmp(str1,"info")==0){
	// Validate input
	if(strcmp(str2,"")==0)
	  printf("format:info <input_file>\n");
	else
	  printGrid(importGrid(str2));
      }

      // display
      else if(strcmp(str1,"display")==0){
	// Validate input
	if(strcmp(str2,"")==0)
	  printf("format:display <input_file> [input_color_file]\n");
	
	else{
	  // import grid
	  grid = importGrid(str2);
	  
	  // color grid?
	  if(strcmp(str3,"") != 0){
	    color_grid = importGrid(str3);
	  }
	  else
	    color_grid = grid;

	  // Validate that color_grid and grid are same size
	  if(color_grid->ncols != grid->ncols ||
	     color_grid->nrows != grid->nrows)
	    printf("flow: the color grid and input grid files are "\
		   "not the same size.\n");
	  
	  // Grids match in size, now decide how to display
	  else{
	    // Open glut window if not already open otherwise break
	    printf("Type: 't' in the glut window to get interactive prompt back.\n");
	    if(window_open == 0){  
	      window_open = 1;
	      // event handler
	      glutMainLoop();
	    }
	    else{
	      // give control back to mainloop
	      break;
	    }
	  }
	}
      }

      // flowdir
      else if(strcmp(str1,"flowdir")==0){
	// Validate input
	if(strcmp(str2,"")==0 || strcmp(str3,"")==0)
	  printf("format:flowdir <input_file> <output_file>\n");
	else
	  writeGrid(flow_direction(importGrid(str2)),str3);
      }

      // flowaccu
      else if(strcmp(str1,"flowaccu")==0){
	// Validate input
	if(strcmp(str2,"")==0 || strcmp(str3,"")==0)
	  printf("format:flowaccu <input_file> <output_file>\n");
	else{
	  writeGrid(flow_accu(importGrid(str2)),str3);
	}
      }

      // null
      else if(strcmp(str1,"")==0){
	printf("Type \"help\" for more information.\n");
      }

      // unknown
      else
	printf("Invalid command. Type \"help\" for more information.\n");
    }
  }
}


// flow_accu - given a flow direction grid returns a flow accumulation grid
//
Grid *flow_accu(Grid *fd){
  int i,j;
  // Allocate flow accumulation grid
  Grid *fa = (Grid *) malloc(sizeof(Grid));
  fa->name = fd->name;
  fa->ncols = fd->ncols;
  fa->nrows = fd->nrows;
  fa->x = fd->x;
  fa->y = fd->y;
  fa->cellsize = fd->cellsize;
  fa->nodata = fd->nodata;
  if((fa->data = (int**) malloc((fa->nrows)*sizeof(int*)))==NULL){
     printf("flow: insufficient memory for FA grid.");
     exit(1);
  }
  for(i=0;i<fa->nrows;i++){
    if((fa->data[i] = (int*) malloc((fa->ncols)*sizeof(int)))==NULL){
      printf("flow: insufficient memory for FA grid.");
      exit(1);
    }
  }

  // Initialize to 0
  for(i=0;i<fa->nrows;i++){
    for(j=0;j<fa->ncols;j++){
      if(fd->data[i][j] == fd->nodata)
	fa->data[i][j]= fa->nodata;
      else
	fa->data[i][j]=0;
    }
  }
  
  // Compute flow direction
  for(i=0;i<fa->nrows;i++){
    for(j=0;j<fa->ncols;j++){
      // If nodata then disregard
      if(fd->data[i][j] != fd->nodata)
	fa->data[i][j] = compFA(fd,fa,i,j);
    }
  }
  
  return fa;
}


// compFA - recusive function for computing flow accumulation
//
int compFA(Grid *fd, Grid *fa, int i, int j){
  int x,y;
  if(fa->data[i][j] != 0)
    return fa->data[i][j];
  else{
    // Cell starts at 1
    fa->data[i][j]=1;
    //go through the 8 surrounding cells
    for(x=i-1;x<=i+1;x++){
      for(y=j-1;y<=j+1;y++){
	// make sure values are in range
	if(x < fa->nrows && x >= 0 &&       // x not out of range
	   y < fa->ncols && y >=0 &&        // y not out of range
	   !(x==i && y==j)){               // don't look at the center
	  
	  // if neighbor flows into current cell then add it's FA
	  if(fd->data[x][y] == pow(2,(3*(2-(x-i+1)))+2-(y-j+1)))
	    fa->data[i][j] += compFA(fd,fa,x,y);
	}
      }
    }
  }
  return fa->data[i][j];
}


// min - return the min of two ints
//
int min(int x, int y){
  return (x<y)?x:y;
}


// flow_direction - compute single most flow direction for each cell
//
Grid *flow_direction(Grid *g){
  int i,j,x,y;

  // Allocate flow direction grid
  Grid *fd = (Grid *) malloc(sizeof(Grid));
  fd->name = g->name;
  fd->ncols = g->ncols;
  fd->nrows = g->nrows;
  fd->x = g->x;
  fd->y = g->y;
  fd->cellsize = g->cellsize;
  fd->nodata = g->nodata;
  if((fd->data = (int**) malloc((fd->nrows)*sizeof(int*)))==NULL){
     printf("flow: insufficient memory");
     exit(1);
  }
  for(i=0;i<fd->nrows;i++){
    if((fd->data[i] = (int*) malloc((fd->ncols)*sizeof(int)))==NULL){
      printf("flow: insufficient memory");
      exit(1);
    }
  }
  
  // Compute FD for each cell
  for(i=0;i<g->nrows;i++){
    for(j=0;j<g->ncols;j++){
      if(g->data[i][j]!=g->nodata){
	int max = 0;
	fd->data[i][j]=0;
	//go through the 8 surrounding cells
	for(x=i-1;x<=i+1;x++){
	  for(y=j-1;y<=j+1;y++){
	    // make sure values are in range
	    if(x < g->nrows && x >= 0 && // x not out of range
	       y < g->ncols && y >=0 ){  // y not out of range
	      // greatest down slope?
	      if((g->data[i][j]-g->data[x][y]) > max && 
		 g->data[x][y] != g->nodata ){
		max = g->data[i][j]-g->data[x][y];
		fd->data[i][j] = pow(2,(3*(x-i+1))+(y-j+1));
	      }
	    }
	  }
	}
      }
      else
	fd->data[i][j]=fd->nodata;
    } 
  }
  return fd;
}


// draw_grid - given a valid grid, this will draw the elevation grid in 2D
//
void draw_grid(){
  float interval,x,y,x1,y1;
  int i,j,mincols,minrows;
  
  // Decide whether to fill the triangles
  if (fillmode) {
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
  } 
  else {
    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
  }
  
  // set x1 and y1 to the upper left corner of the data
  if(grid->nrows > grid->ncols){
    interval = 2.0/grid->nrows;
    x1 = -((interval*grid->ncols)/2);
    y1 = 1;
  }
  else{
    interval = 2.0/grid->ncols;
    x1 = -1;
    y1 = ((interval*grid->nrows)/2);
  }
  x = x1;
  y = y1;
  
  // loop through the data and draw 2 triangles per cell marked by top left 
  // corner (x,y). Drawn in clockwise rotation
  for(i=0;i<grid->nrows;i+=resolution){
    for(j=0;j<grid->ncols;j+=resolution){
      mincols=min(grid->ncols-1-j,resolution);
      minrows=min(grid->nrows-1-i,resolution);
      // draw upper left triangle
      glBegin(GL_POLYGON);
      color_map(color_grid->data[i][j],color_grid);
      glVertex3f(x,y,z_map(grid->data[i][j],grid));
      color_map(color_grid->data[i][j+mincols],color_grid);
      glVertex3f(x+interval*mincols,y,z_map(grid->data[i][j+mincols],grid));
      color_map(color_grid->data[i+minrows][j],color_grid);
      glVertex3f(x,y-interval*minrows,z_map(grid->data[i+minrows][j],grid));
      glEnd();
      // draw lower right triangle
      glBegin(GL_POLYGON);
      color_map(color_grid->data[i][j+mincols],color_grid);
      glVertex3f(x+interval*mincols,y,z_map(grid->data[i][j+mincols],grid));
      color_map(color_grid->data[i+minrows][j+mincols],color_grid);
      glVertex3f(x+interval*mincols,y-interval*minrows,
		 z_map(grid->data[i+minrows][j+mincols],grid));
      color_map(color_grid->data[i+minrows][j],color_grid);
      glVertex3f(x,y-interval*minrows,z_map(grid->data[i+minrows][j],grid));
      glEnd();
      // move to the next column
      x+=interval*resolution;
    }
    // move back to the first column and the next row
    x = x1;
    y-=interval*resolution;
  }
}


// z_map - map height into a viewable z range (ie -1.0 to -10.0)
//
float z_map(int h, Grid *g){
  if(h == g->nodata)
    return 0;
  else
    return (maxZ/(g->max - g->min))*(h - g->min);
}


// color_map - set the color based on a given height
//
void color_map(int h, Grid *g){
  int i,j;
  float x,k;
  GLfloat fc[3];
  float flow_buckets[]={g->min, g->min+5.0, g->min+30.0, 
			g->min+100.0, g->min+1000.0, g->max};

  // map h from 0 -> 1
  x = (1.0/(g->max - g->min))*(h - g->min); 

  // If nodata color white
  if (h == g->nodata)
    glColor3f(0,0,0);

  // Greyscale
  else if(colorStyle == GRAYSCALE)
    glColor3f(x,x,x);

  // Rainbow
  else if(colorStyle == RAINBOW){
    if(x < (1.0/5))
      glColor3f(1-(5*x),1,0);
    if(x < 2.0/5 && x >= 1.0/5)
      glColor3f(0,1,(x-(1.0/5))*5);
    if(x < 3.0/5 && x >= 2.0/5)
      glColor3f(0,((3.0/5)-x)*5,1);
    if(x < 4.0/5 && x >= 3.0/5)
      glColor3f((x-(3.0/5))*5,0,1);
    if(x <= 5.0/5 && x >= 4.0/5)
      glColor3f(1,0,(1-x)*5);
  }

  // Custom
  else {

    // Regular size buckets
    if(colorStyle != FLOWACCU){
      // check each bucket
      for(i = 0; i < bu; i++){
	if(x <= (float)(i+1)/bu && x >= (float)i/bu){
	  // map bucket from 0 -> 1
	  k = (x-((float)i/bu))/(1.0/bu);
	  // map each color (rgb) from the color scheme
	  for(j = 0; j < 3; j++){
	    // Low to high map (ie 0.2 -> 0.8)
	    if(color[i][j] <= color[i+1][j])
	      fc[j]=((color[i+1][j] - color[i][j]) * k) + color[i][j];
	    // high to low map
	    else
	      fc[j]=((color[i][j] - color[i+1][j]) * (1.0 - k)) + color[i+1][j];
	  }
	}
      }
    }

    // Set bucket size for flow accumulation
    else{
      for(i = 0; i < 5; i++){
	//printf("h: %f buk: %f %f \n",x,flow_buckets[i],flow_buckets[i+1]);
	if(h >= flow_buckets[i] && h <= flow_buckets[i+1]){
	  // map bucket from 0 -> 1
	  k = (h-flow_buckets[i])/(flow_buckets[i] - flow_buckets[i+1]);
	  // map each color (rgb) from the color scheme
	  for(j = 0; j < 3; j++){
	    // Low to high map (ie 0.2 -> 0.8)
	    if(color[i][j] <= color[i+1][j])
	      fc[j]=((color[i+1][j] - color[i][j]) * k) + color[i][j];
	    // high to low map
	    else
	      fc[j]=((color[i][j] - color[i+1][j]) * (1.0 - k)) + color[i+1][j];

	    //printf("x: %f k: %f fc[j]: %f color[i][j]: %f color[i+1][j]:"\
		//"%f\n",x,k,fc[j],color[i][j],color[i+1][j]);
	  }
	}
      }
    }
    
    glColor3f(fc[0],fc[1],fc[2]);
  }
}

// renderScene - for animation
//
void renderScene(void) {
  
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glPushMatrix();
  
  // Rotate around z
  glRotatef(angle,0.0,0.0,1.0);
  
  draw_grid();
  glPopMatrix();
  glutSwapBuffers();
  
  if(fly == 1){
    // increase the angle for the next frame
    angle++;
  }
}

// display - controls what is drawn to the screen
//
void display(void) {
  // clear window
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity(); 
  
  // adjust yaw and pitch
  glRotatef(pitch,1,0,0);
  glRotatef(yaw,0,1,0);
  glRotatef(roll,0,0,1);

  // move terrain
  glTranslatef(pos[0], pos[1], pos[2]);

  // Look at terrain from a better perspective
  gluLookAt(0, -2, 2,
	    0, 0, 0,
	    0, 0, 1);

  // rotate terrain
  glRotatef(theta[0], 1,0,0);
  glRotatef(theta[1], 0,1,0);
  glRotatef(theta[2], 0,0,1); 

  draw_grid();
 
  glFlush();
}


// mouse - assigns action to mouse events
//
void mouse(int x, int y){
  pos[0]=(x-(window_width/2.0))/(window_width/2.0);
  pos[1]=-(y-(window_height/2.0))/(window_height/2.0);
  glutPostRedisplay();
}


// reshape - assigns action to resizing the window
//
void reshape(int w, int h){
  window_width = w;
  window_height = h;
}


// specialkey - handle special keyboard inputs
//
void specialkey(int key, int x, int y) {
  switch(key) {
  case GLUT_KEY_DOWN:
    pitch += 1.0;
    glutPostRedisplay();
    break;
  case GLUT_KEY_UP:
    pitch -= 1.0;
    glutPostRedisplay();
    break;
  case GLUT_KEY_LEFT:
    yaw -= 1.0;
    glutPostRedisplay();
    break;
  case GLUT_KEY_RIGHT:
    yaw += 1.0;
    glutPostRedisplay();
    break;
    }
}


// keypress - assigns actioins to keypresses
//
void keypress(unsigned char key, int x, int y) {
  float i = .05;
  switch(key) {
     /* various object rotations: */
  case 't':
    interact();
    glutPostRedisplay();
    break;
  case 'e':
    maxZ -= 0.1; 
    glutPostRedisplay();
    break;
  case 'E':
    maxZ += 0.1;
    glutPostRedisplay();
    break;
  case 'x':
    theta[0] += delta;
    glutPostRedisplay();
    break;
  case 'y':
    theta[1] += delta;
    glutPostRedisplay();
    break;
  case 'z':
    theta[2] += delta;
    glutPostRedisplay();
    break;
  case 'X':
    theta[0] -= delta; 
    glutPostRedisplay();
    break;
  case 'Y':
    theta[1] -= delta; 
    glutPostRedisplay();
    break;
  case 'Z':
    theta[2] -= delta; 
    glutPostRedisplay();
    break;
  case 'b':
    pos[2] -= delta/8; 
    glutPostRedisplay();
    break;
  case 'f':
    pos[2] += delta/8; 
    glutPostRedisplay();
    break;
  case 'u': 
     pos[1] += delta/8; 
    glutPostRedisplay();
    break;
  case 'd': 
    pos[1] -= delta/8; 
    glutPostRedisplay();
    break;
  case 'r':
    pos[0] += delta/8; 
    glutPostRedisplay();
    break;
  case 'l':
    pos[0] -= delta/8; 
    glutPostRedisplay();
    break;

  case '-':
    if(delta >= 0.5)
      delta -= 0.5;
    printf("Delta: %f\n",delta);
    break;
  case '+':
    if(delta <= 360)
      delta += 0.5; 
    printf("Delta: %f\n",delta);
    break;
  case 'q':
    exit(0);
    break;
  case '>':  // increase resolution
    if(resolution > 1){
      resolution /= 2;
      glutPostRedisplay();
    }
    break;
  case '<':  // decrease resolution
    if(resolution < grid->max){
      resolution *= 2;
      glutPostRedisplay();
    }
    break;
  } 
}


// main_menu - creates a menu
//
void main_menu(int value)
{
  switch (value){
  case FILL: 
    /* toggle outline/fill */
    fillmode = !fillmode;
    break;
  case FLY: 
    /* toggle fly-through */
    fly = !fly;
    break;
  case GRAYSCALE: 
    colorStyle = GRAYSCALE;
    break;
  case RAINBOW:
    colorStyle = RAINBOW;
    break;
  case CUSTOM:
    colorStyle = CUSTOM;
    color[0] = blue;
    color[1] = sand;
    color[2] = fgreen;
    color[3] = brown;
    color[4] = dbrown;
    color[5] = white;
    break;
  case FLOWACCU:
    colorStyle = FLOWACCU;
    color[0] = white;
    color[1] = yellow;
    color[2] = cyan;
    color[3] = blue;
    color[4] = dblue;
    color[5] = purple;
    break;
  case EXIT: 
    exit(0);
  }
  int x;
  for(x = 0; x <= bu; x++){
    if(value == RED+(x*100))
      color[x] = red;
    if(value == YELLOW+(x*100))
      color[x] = yellow;
    if(value == MAGENTA+(x*100))
      color[x] = magenta;
    if(value == GREEN+(x*100))
      color[x] = green;
    if(value == CYAN+(x*100))
      color[x] = cyan;
    if(value == BLUE+(x*100))
      color[x] = blue;
    if(value == WHITE+(x*100))
      color[x] = white;
    if(value == GRAY+(x*100))
      color[x] = gray;
    if(value == BLACK+(x*100))
      color[x] = black;
    if(value == SAND+(x*100))
      color[x] = sand;
    if(value == FGREEN+(x*100))
      color[x] = fgreen;
    if(value == LGREEN+(x*100))
      color[x] = lgreen;
    if(value == DBROWN+(x*100))
      color[x] = dbrown;
    if(value == BROWN+(x*100))
      color[x] = brown;
    if(value == DBLUE+(x*100))
      color[x] = dblue;
    if(value == PURPLE+(x*100))
      color[x] = purple;
  }
  glutPostRedisplay();
}


// importGrid - brings grid file into an array
//
Grid *importGrid(char *path)
{
  FILE *inputf;
  int i,j,value;

  // Validate input file
  if ((inputf = fopen(path, "r"))== NULL){
     printf("flow: can't open %s\n",path);
     exit(1);
  }

  // Input file info into the grid structure
  Grid *g = (Grid *) malloc(sizeof(Grid));
  g->name = path;
  int scan = fscanf(inputf,"%*s%d%*s%d%*s%f%*s%f%*s%d%*s%d",
	 &g->ncols,&g->nrows,&g->x,
	 &g->y,&g->cellsize,&g->nodata);
  if(scan==0){
    printf("flow: trouble reading input file!\n");
    exit(1);
  }

  // Dynamically allocate array for data
  if((g->data = (int**) malloc((g->nrows)*sizeof(int*)))==NULL){
     printf("flow: insufficient memory");
     exit(1);
  }
  for(i=0;i<g->nrows;i++){
    if((g->data[i] = (int*) malloc((g->ncols)*sizeof(int)))==NULL){
      printf("flow: insufficient memory");
      exit(1);
    }
  }
  
  // Put data into the array and calculate max & min
  g->min = 9999999;
  g->max = 0;
  for(i=0;i<g->nrows;i++){
    for(j=0;j<g->ncols;j++){
      if(fscanf(inputf,"%d",&value)!=EOF){
	if(value < g->min && value != g->nodata)
	  g->min = value;
	if(value > g->max)
	  g->max = value;
	g->data[i][j]=value;
      }
      else{
	printf("flow: data file is corrupt");
	exit(1);
      }
    }
  }

  // Set resolution as a function of gridsize
  resolution = min(g->ncols,g->nrows)/8;

  // Close file
  fclose(inputf);

  // Return grid
  return g;
}

// writeGrid - write grid structure to file path
//
void writeGrid(Grid *g,char *path){
  FILE *outputf;
  int i,j;

  // Validate output file
  if ((outputf = fopen(path, "w"))== NULL){
     printf("mult: can't write to %s\n",path);
     exit(1);
  }
  
  // Write grid info
  fprintf(outputf,"ncols\t\t%u\n",g->ncols);
  fprintf(outputf,"nrows\t\t%u\n",g->nrows);
  fprintf(outputf,"xllcorner\t%f\n",g->x);
  fprintf(outputf,"yllcorner\t%f\n",g->y);
  fprintf(outputf,"cellsize\t%u\n",g->cellsize);
  fprintf(outputf,"NODATA_value\t%d\n",g->nodata);

  // Write data
  for(i=0;i<g->nrows;i++){
    for(j=0;j<g->ncols;j++){
      fprintf(outputf, "%d ",g->data[i][j]);
    }
    fprintf(outputf, "\n");
  }

  // Close file
  fclose(outputf);
}


// printGrid - Print the grid info 
//
void printGrid(Grid *g){
  printf("\nGrid Info:\n\n");
  printf("Name: %s\n",g->name);
  printf("Number of Columns: %u\n",g->ncols);
  printf("Number of Rows: %u\n",g->nrows);
  printf("Lat Lon X Corner: %f\n",g->x);
  printf("Lat Lon Y Corner: %f\n",g->y);
  printf("Cell size: %u\n",g->cellsize);
  printf("No data value: %d\n",g->nodata);
  printf("Min value: %d\n",g->min);
  printf("Max value: %d\n",g->max);
}