//rain/snow - localized
//temperature - changes over position, time, elevation, possibly latitude.
//fault slippage
//mudslide/slippage
//erosion - water flow
//ice freezing - breaks rock
//snow accumulation
//evaporation

//wind
//impact craters
//volcanic upwell/lava flow
//rectangular or spherical coords

#include <GL/glut.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <time.h>

#define ROTATETOOL 0
#define TRANSLATETOOL 1
#define ZOOMTOOL 2
#define WIRE 3

//Globals
int viewwidth, viewheight;
int mainwindow;
int wireframe = 0;

int curtool = ROTATETOOL;
int mainmenu;
int lastx, lasty;

double vx = 0, vy = 0;
double vdir = 0.;

	int width = 129;
	int height = 129;
	int numsteps = 100000;
	double timescale = 1.0;

	double snowlevel = 5.;
	double minwater = 0.2;
	double minwaterdraw = 0.3;
	double sedrate = 0.05; //how much dirt is picked up by water
	double rocksed = 0.3; //how much rock can be cut away by water.
	double annualrain = 0.01;
	double chanceofrain = 0.1;
	double evaporate = annualrain * chanceofrain / 3.;
	double evapend = annualrain * chanceofrain * 1.5;

	int timetoevap = 100;
	int t = 0;

	double *rock;
	double *dirt;
	double *water;



double gaussianRand(double m, double var)
{
	double r1 = ((double)(rand() + 1) / (double)(RAND_MAX + 1));
	double r2 = ((double)rand() / (double)RAND_MAX);
	return m + sqrt(var) * sqrt(-2. * log(r1)) * cos(2 * 3.1415926535 * r2);
}


void rain(double *w, int width, int height, double timescale, double precipitation)
{
	for(int y = 0; y < height; y++)
	{
		for(int x = 0; x< width; x++)
		{
			double r = gaussianRand(timescale * precipitation, timescale * precipitation / 3.);
			if(r < 0) r = 0;
			w[y * width + x] += r;
		}
	}
}

void evap(double *w, int width, int height, double timescale, double rate)
{
	for(int y = 0; y < width * height; y++)
	{
		w[y] -= rate * timescale;
		if(w[y] < 0) w[y] = 0;
	}
}


void genFractal(double *r, int width, int height, double dim, int pitch)
{
	//dimensions must be equal and a power of 2 + 1
	assert(width == height);
	if(width < 3) return;
	int nextw = (width - 1) / 2 + 1;
	
	double scale = (double)width * pow(.5, 2 * dim);

	double a = r[0];
	double b = r[width - 1];
	double c = r[(height - 1) * pitch];
	double d = r[(height - 1) * pitch + width - 1];

	//middle
	r[(nextw - 1) + (nextw - 1) * pitch] = (a + b + c + d) / 4. + gaussianRand(0, 1) * scale;
	//top
	if(r[nextw - 1] == 100000.)
	{
		r[nextw - 1] = (a + b) / 2. + gaussianRand(0, 1) * scale;
	}
	//left
	if(r[(nextw - 1) * pitch] == 100000.)
	{
		r[(nextw - 1) * pitch] = (a + c) / 2. + gaussianRand(0, 1) * scale;
	}
	//right
	r[(nextw - 1) * pitch + width - 1] = (b + d) / 2. + gaussianRand(0, 1) * scale;
	//bottom
	r[nextw - 1 + (width - 1) * pitch] = (c + d) / 2. + gaussianRand(0, 1) * scale;

	//recurse
	genFractal(r, nextw, nextw, dim, pitch);
	genFractal(r + nextw - 1, nextw, nextw, dim, pitch);
	genFractal(r + (nextw - 1) * pitch, nextw, nextw, dim, pitch);
	genFractal(r + nextw - 1 + (nextw - 1) * pitch, nextw, nextw, dim, pitch);
}

void initHeightField(double **r, double **d, double **w, int width, int height)
{
	*r = new double[width * height];
	*d = new double[width * height];
	*w = new double[width * height];
	for (int y = 0; y < height; y++)
	{
		for(int x = 0; x < width; x++)
		{
			(*r)[y * width + x] = 100000.;
			(*d)[y * width + x] = 0;
			(*w)[y * width + x] = 0;

		}
	}
	(*r)[0] = 0.0;
	(*r)[width - 1] = 0.0;
	(*r)[(height - 1) * width] = 0.0;
	(*r)[(height - 1) * width + width - 1] = 0.0;

	genFractal(*r, width, height, 1.5, width);
}

double slide(double *r, double *d, double *w, int x1, int x2)
{
	double minslope = 1;
	double slidepct = .03;
	double snowslide = 1;
	double dirtslide = .9;
	double rockslide = .02;


	double h1 = r[x1] + d[x1];
	double hd = (r[x1] + d[x1]) - (r[x2] + d[x2]);
	if(h1 > snowlevel)
	{
		//snow slides too...
		hd += w[x1] - w[x2];
	}
	if (hd <= gaussianRand(minslope, minslope / 10.)) return d[x2] + r[x2] + w[x2];;

	//snow first
	double slideamt = hd * slidepct;
	if(h1 > snowlevel)
	{
		if(w[x1] / snowslide >= slideamt)
		{
			w[x1] -= slideamt * snowslide;
			w[x2] += slideamt * snowslide;
			return d[x2] + r[x2] + w[x2];;
		}
		slideamt -= w[x1] / snowslide;
		w[x2] += w[x1];
		w[x1] = 0;
	}

	//dirt
	if(d[x1] / dirtslide >= slideamt)
	{
		d[x1] -= slideamt * dirtslide;
		d[x2] += slideamt * dirtslide;
		return d[x2] + r[x2] + w[x2];;
	}
	slideamt -= d[x1] / dirtslide;
	d[x2] += d[x1];
	d[x1] = 0;

	//rock->dirt
	d[x2] += slideamt * rockslide;
	r[x1] -= slideamt * rockslide;
	return d[x2] + r[x2] + w[x2];
}

void erode(double *r, double *d, double *w, int width, int height, double timescale)
{
/*
	double *nr = new double[width * height];
	double *nd = new double[width * height];
	double *nw = new double[width * height];
*/
	//could do either multiple or single flow directions.
	//rock always erodes at a low value
	//dry dirt erodes quickly, medium dirt erodes slowly, and wet dirt erodes quickly.
	//for w < 2 er = (w - 1)^2 for w > 2 er = log(w + (e - 2))   
	//(both + base rock erosion rate) / 5, maxed to 0.5
	//this is for all lower spots, 
	//how to avoid problem of multiple sources/sinks - check the footprint paper
	//each point will only feed one other - treat as an inflow problem....
	//poll neighbors and see if I am its lowest neighbor, then collect all inflows at once.
	//do in two passes - first create the table of lowest neighbors, then pull in.

	//slides... should do this in order of steepness
	for(int y = 1; y < height - 1; y++)
	{
		for(int x = 1; x < width - 1; x++)
		{
			int lowest = x - 1 + (y - 1) * width;
			double lh = slide(r, d, w, x + y * width, x - 1 + (y - 1) * width);
			double th = slide(r, d, w, x + y * width, x + (y - 1) * width);
			if(th < lh)
			{
				lowest = x + (y - 1) * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x + 1 + (y - 1) * width);
			if(th < lh)
			{
				lowest = x + 1 + (y - 1) * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x - 1 + y * width);
			if(th < lh)
			{
				lowest = x - 1 + y * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x + 1 + y * width);
			if(th < lh)
			{
				lowest = x + 1 + y * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x - 1 + (y + 1) * width);
			if(th < lh)
			{
				lowest = x - 1 + (y + 1) * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x + (y + 1) * width);
			if(th < lh)
			{
				lowest = x + (y + 1) * width;
				lh = th;
			}
			th = slide(r, d, w, x + y * width, x + 1 + (y + 1) * width);
			if(th < lh)
			{
				lowest = x + 1 + (y + 1) * width;
				lh = th;
			}

			int mypt = x + y * width;
			if((r[mypt] + d[mypt] <= snowlevel) && (w[mypt] >= minwater))
			{

				int ht = r[mypt] + d[mypt] + w[mypt];
				
				if(ht > lh)
				{


					double wmove;
					if((ht - lh) / 2. >= (w[mypt] - minwater))
					{
						wmove = w[mypt] - minwater;
					}
					else
					{
						wmove = (ht - lh) / 2.;
					}
					w[lowest] += wmove;
					w[mypt] -= wmove;

					//how much dirt to move
					double dmove = wmove * sedrate;
					//dirt
					if(d[mypt] >= dmove)
					{
						d[mypt] -= dmove;
						d[lowest] += dmove;
					}
					else
					{
						dmove -= d[mypt];
						d[lowest] += d[mypt];
						d[mypt] = 0;

						//rock->dirt
						d[lowest] += dmove * rocksed;
						r[mypt] -= dmove * rocksed;
					}
				}
			}
		}
	}
}

void draw(double *r, double *d, double *w, int width, int height)
{

	float *nr;
	nr = new float[viewwidth * viewheight * 3];

	memset(nr, 0, viewwidth * viewheight * 3 * sizeof(float));

	delete []nr;


	double *heights = new double[width * height];
	for(int i = 0; i < width * height; i++)
	{
		heights[i] = r[i] + d[i] + w[i];
	}

	
	double *normalx = new double[width * height];
	double *normaly = new double[width * height];

	for(int y = 1; y < height - 1; y++)
	{
		for(int x = 1; x < width - 1; x++)
		{
			double h = heights[x + y * width];
			double hd1 = heights[x - 1 + y * width] - h;
			double hd2 = heights[x + (y - 1) * width] - h;
			double hd3 = heights[x + 1 + y * width] - h;
			double hd4 = heights[x + (y + 1) * width] - h;
			normalx[x + y * width] = (hd1 - hd3) / 2.;
			normaly[x + y * width] = (-hd4 + hd2) / 2.;
		}
	}

	glMatrixMode(GL_PROJECTION);
	glPushMatrix();

	double ht = 0;
	if((vx > -.5) && (vx < .5) && (vy < .5) && (vy > - .5))
	{
		int vpx = (vx + .5) * width;
		int vpy = (vy + .5) * height;
		ht = heights[vpx + vpy * width];
	}

	glRotated(-90, 1, 0, 0);
	glRotated(vdir / 3.1415926536 * 180., 0, 0, 1);
	glTranslated(-vx, -vy, -(ht + 3) / width);
	
//	nr = new float[width * height * 3];

	for(y = 0; y < height - 1; y++)
	{
		glBegin(GL_QUAD_STRIP);
		for(int x = 0; x < width; x++)
		{
			int spot = x + (y + 1) * width;

			double theight = heights[spot];

			glNormal3d(normalx[spot], normaly[spot], 1);
			float color[3];
			if(w[spot] > minwaterdraw)
			{
				if(r[spot] + d[spot] > snowlevel)
				{
					//white snow
					color[0] = 1.;
					color[1] = 1.;
					color[2] = 1.;
				}
				else
				{
					//blue water
					color[0] = 0.;
					color[1] = 0.;
					color[2] = 1.;
				}
			}
			else if(d[spot] > 0)
			{
				if(w[spot] > 0)
				{
					//green plants
					color[0] = 0.;
					color[1] = 1.;
					color[2] = 0.;
				}
				else
				{
					//brown dirt
					color[0] = .5;
					color[1] = .2;
					color[2] = .2;
				}
			}
			else
			{
				//grey rock
				color[0] = .5;
				color[1] = .4;
				color[2] = .6;
			}

/*			nr[spot * 3] = color[0];
			nr[spot * 3 + 1] = color[1];
			nr[spot * 3 + 2] = color[2];
*/
			theight /= (double)width;

			glMaterialfv(GL_FRONT, GL_AMBIENT, color);
			glVertex3d((double)x / (double)width - .5, (double)(y + 1) / (double) height - .5, theight);

			spot = x + y * width;

			theight = heights[spot] / (double)width;
			if(w[spot] > minwaterdraw)
			{
				if(r[spot] + d[spot] > snowlevel)
				{
					//white snow
					color[0] = 1.;
					color[1] = 1.;
					color[2] = 1.;
				}
				else
				{
					//blue water
					color[0] = 0.;
					color[1] = 0.;
					color[2] = 1.;
				}
			}
			else if(d[spot] > 0)
			{
				if(w[spot] > 0)
				{
					//green plants
					color[0] = 0.;
					color[1] = 1.;
					color[2] = 0.;
				}
				else
				{
					//brown dirt
					color[0] = .5;
					color[1] = .2;
					color[2] = .2;
				}
			}
			else
			{
				//grey rock
				color[0] = .5;
				color[1] = .4;
				color[2] = .6;
			}
			glNormal3d(normalx[spot], normaly[spot], 1);
			glVertex3d((double)x / (double)width - .5, (double)y / (double) height - .5, theight);


		}
		glEnd();
	}
	glPopMatrix();
	delete []normalx;
	delete []normaly;
	delete []heights;
//	glDrawPixels(width, height, GL_RGB, GL_FLOAT, nr);
//	delete []nr;
		
}


void step()
{
	if(rand() < (double)RAND_MAX * chanceofrain) 
		rain(water, width, height, 1, annualrain);
	erode(rock, dirt, water, width, height, 1);
	evap(water, width, height, 1, evaporate);
	glutPostRedisplay();
	t++;
	if(t >= timetoevap) evaporate = evapend;
}


// display: called by glut to redraw
void display(void)
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	draw(rock, dirt, water, width, height);

	glutSwapBuffers();
}

// mousemove: called by glut when the mouse is moved with no button pressed
//  x: new x pos
//  y: new y pos
void mousemove(int _x, int _y)
{
}

// mousedrag: called by glut when the mouse is moved with a button pressed
//  x: new x pos
//  y: new y pos
void mousedrag(int _x, int _y)
{
	double dx = _x - lastx;
	double dy = _y - lasty;

	vx -= dy / 1000. * sin(vdir);
	vy -= dy / 1000. * cos(vdir);

	vdir += dx / 100.;


	glutPostRedisplay();
	lastx = _x;
	lasty = _y;
}

//mouse: called by glut for mouse clicks
void mouse(int button, int state, int _x, int _y)
{
	//we only care about drags
	lastx = _x;
	lasty = _y;
	//transform coords...
	float x = (float)(_x * 2 - viewwidth) / (float)viewwidth;
	float y = -(float)(_y * 2 - viewheight) / (float)viewheight;
	glutPostRedisplay();
}

// reshape: called by GLUT when the viewport changes.
//  width: new width
//  height: new height
void reshape(int width, int height)
{
	viewwidth = width;
	viewheight = height;
	glViewport(0, 0, width, height);
}

void menu(int value)
{
	if(value == WIRE)
	{
		wireframe = !wireframe;
		glutPostRedisplay();
	}
	else
	{
		curtool = value;
	/*	switch(curtool)
		{
			case ROTATETOOL:
			{
				glutSetCursor(GLUT_CURSOR_CYCLE);
				break;
			}
			case TRANSLATETOOL:
			{
				glutSetCursor(GLUT_CURSOR_FULL_CROSSHAIR);
				break;
			}
			case ZOOMTOOL:
			{
				glutSetCursor(GLUT_CURSOR_UP_DOWN);
				break;
			}
		}*/
	}
}

void initmenus()
{
	mainmenu = glutCreateMenu(menu);
	glutAddMenuEntry("Rotate", ROTATETOOL);
	glutAddMenuEntry("Translate", TRANSLATETOOL);
	glutAddMenuEntry("Zoom", ZOOMTOOL);
	glutAddMenuEntry("Wireframe", WIRE);
	glutAttachMenu(GLUT_RIGHT_BUTTON);
}



int main(int argc, char **argv)
{
	glutInit(&argc, argv);
	
	
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
	mainwindow = glutCreateWindow("Fun with 3D!");
	glutDisplayFunc(display);
	glutMouseFunc(mouse);
	glutMotionFunc(mousedrag);
	glutPassiveMotionFunc(mousemove);
	glutReshapeFunc(reshape);
	glutSetCursor(GLUT_CURSOR_CROSSHAIR);
	glutIdleFunc(step);
	
	glDepthFunc(GL_LESS);
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_LIGHTING);
	//glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
	float ambient[4];
	ambient[0] = 1.0;
	ambient[1] = 1.0;
	ambient[2] = 1.0;
	ambient[3] = 1.0;
	//glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient);
	glShadeModel(GL_SMOOTH);
	glEnable(GL_NORMALIZE);
	glEnable(GL_CULL_FACE);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	float pos[3];
	pos[0] = 0;
	pos[1] = 2;
	pos[2] = 3;
	glLightfv(GL_LIGHT0, GL_POSITION, pos);

	glMatrixMode(GL_PROJECTION);
	glFrustum(-.01, .01, -.01, .01, .01, .5);
	glEnable(GL_NORMALIZE);
	initmenus();

	srand(time(NULL));
	//either load or initialize height field
	initHeightField(&rock, &dirt, &water, width, height);


	glutMainLoop();
	return 0;
}