//
//  liquidSurf
//
//  Created by Bill on 25/04/2007.
//
#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <time.h>
using namespace std;

#define kWindowWidth 800
#define kWindowHeight 600

int nx, ny, numNormals;
clock_t ctim;

GLfloat** mesh;
GLfloat** oldmesh;
GLfloat** velocities;
GLfloat** normals;
GLfloat k1, k2, damp, localMean, meanDist, stime, dt, cdt, nextDrip;
GLfloat sx, sy, dx, dy;
GLfloat ax, ay, az, bx, by, bz, cx, cy, cz;
GLfloat lightAmbient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
GLfloat lightDiffuse[] = { 1.0f, 1.0f, 1.0f, 1.0f };
GLfloat lightSpecular[] = { 0.2f, 0.2f, 0.2f, 1.0f };
GLfloat lightPosition[] = { 20.0f, 10.0f, 5.0f, 1.0f };
GLfloat materialAmbient[] = { 0.0f, 0.0f, 0.0f, 1.0f };
GLfloat materialDiffuse[] = { 0.5f, 0.6f, 1.0f, 1.0f };
GLfloat materialSpecular[] = { 0.1f, 0.1f, 0.1f, 1.0f };

GLvoid InitGL(GLvoid);
GLvoid DrawGLScene(GLvoid);
GLvoid ReSizeGLScene(int Width, int Height);
void initSurface();

int main(int argc, char** argv)
{
  glutInit(&argc, argv);
  glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize (kWindowWidth, kWindowHeight);
  glutInitWindowPosition (100, 100);
  glutCreateWindow (argv[0]);
  
  InitGL();

  stime = 0.0f;
  dt = 0.2f;
  cdt = CLOCKS_PER_SEC/25;

  initSurface();

  glutDisplayFunc(DrawGLScene);
  glutReshapeFunc(ReSizeGLScene);
  
  glutMainLoop();
  
  return 0;
}

void initSurface()
{
  sx = 20;
  sy = 20;
  nx = 120;
  ny = 120;
  dx = sx/(nx-1);
  dy = sy/(ny-1);
  k1 = 5;
  k2 = 20;
  damp = 0.99;
  srand(1);
  nextDrip = 50*dt + 50*dt*float(rand())/RAND_MAX;

  // declare velocity array
  velocities = new GLfloat*[nx];
  for (int i = 0; i < nx; i++)
    velocities[i] = new GLfloat[ny];

  // declare vertices arrays
  mesh = new GLfloat*[nx];
  for (int i = 0; i < nx; i++)
    mesh[i] = new GLfloat[ny];
  oldmesh = new GLfloat*[nx];
  for (int i = 0; i < nx; i++)
    oldmesh[i] = new GLfloat[ny];

  // fill arrays
  // velocities
  for (int j = 0; j < ny; j++)
    for (int i = 0; i < nx; i++)
	  velocities[i][j] = 0;
  // meshes
  for (int j = 0; j < ny; j++)
    for (int i = 0; i < nx; i++)
	  mesh[i][j] = 0;
  for (int j = 0; j < ny; j++)
    for (int i = 0; i < nx; i++)
	  oldmesh[i][j] = 0;

  // declare normals array
  numNormals = (nx-1)*(ny-1)*2;
  normals = new GLfloat*[numNormals];
  for (int i = 0; i < numNormals; i++)
    normals[i] = new GLfloat[3];
}

void calculateVertex(int i, int j)
{
  meanDist = localMean-oldmesh[i][j];
  velocities[i][j] = damp*velocities[i][j] - k1*dt*oldmesh[i][j] + k2*dt*meanDist;
  mesh[i][j] = oldmesh[i][j] + dt*velocities[i][j];
}

void calculateSurface()
{
  GLfloat** temp;
  temp = oldmesh;
  oldmesh = mesh;
  mesh = temp;

  // see if it's time for a new drip
  if (stime > nextDrip)
  {
    int dripx = rand() % (nx-2) + 1;
	int dripy = rand() % (ny-2) + 1;
	GLfloat newV = -(rand() % 250 + 50);
	newV = newV/50;
    velocities[dripx][dripy] = newV;
	velocities[dripx+1][dripy] = newV;
	velocities[dripx][dripy+1] = newV;
	velocities[dripx+1][dripy+1] = newV;
	nextDrip = stime + 50*dt + 50*dt*float(rand())/RAND_MAX;
  }

  // calculate bulk vertices
  for (int j = 1; j < ny-1; j++)
    for (int i = 1; i < nx-1; i++)
	{
	  localMean = (oldmesh[i-1][j] + oldmesh[i+1][j] + oldmesh[i][j-1] + oldmesh[i][j+1]
	               + 0.25f*oldmesh[i-1][j+1] + 0.25f*oldmesh[i+1][j+1]
				   + 0.25f*oldmesh[i-1][j-1] + 0.25f*oldmesh[i+1][j-1])/5;
      calculateVertex(i, j);
    }
  // calculate top / bottom edges
  for (int i = 1; i < nx-1; i++)
  {
    localMean = (oldmesh[i-1][0] + oldmesh[i+1][0] + 2*oldmesh[i][1])/4;
    calculateVertex(i, 0);
	localMean = (oldmesh[i-1][ny-1] + oldmesh[i+1][ny-1] + 2*oldmesh[i][ny-2])/4;
	calculateVertex(i, ny-1);
  }
  // calculate left / right edges
  for (int j = 1; j < ny-1; j++)
  {
    localMean = (oldmesh[0][j-1] + oldmesh[0][j+1] + 2*oldmesh[1][j])/4;
	calculateVertex(0, j);
	localMean = (oldmesh[nx-1][j-1] + oldmesh[nx-1][j+1] + 2*oldmesh[nx-2][j])/4;
	calculateVertex(nx-1, j);
  }
  // calculate corners
  localMean = (oldmesh[1][0] + oldmesh[0][1])/2;
  calculateVertex(0, 0);
  localMean = (oldmesh[1][ny-1] + oldmesh[0][ny-2])/2;
  calculateVertex(0, ny-1);
  localMean = (oldmesh[nx-2][0] + oldmesh[nx-1][1])/2;
  calculateVertex(nx-1, 0);
  localMean = (oldmesh[nx-2][ny-1] + oldmesh[nx-1][ny-2])/2;
  calculateVertex(nx-1, ny-1);

  // calculate new normals
  int ni = 0;
  for (int j = 0; j < (ny-1); j++)
    for (int i = 0; i < (nx-1); i++)
	{
	  ax = dx*1.0f;
	  ay = mesh[i+1][j] - mesh[i][j];
	  az = 0.0f;
	  bx = dx*1.0f;
	  by = mesh[i+1][j+1] - mesh[i][j];
	  bz = dy*1.0f;
	  cx = 0.0f;
	  cy = mesh[i][j+1] - mesh[i][j];
	  cz = dy*1.0f;
      normals[ni][0] = -ay*bz;// + az*by;
	  normals[ni][1] = /*-az*bx +*/ ax*bz;
	  normals[ni][2] = -ax*by + ay*bx;
	  normals[ni+1][0] = -by*cz + bz*cy;
	  normals[ni+1][1] = /*-bz*cx +*/ bx*cz;
	  normals[ni+1][2] = -bx*cy;// + by*cx;
	  ni += 2;
    }
}

GLvoid InitGL(GLvoid)
{
  glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
  glClearDepth(1.0);
  glDepthFunc(GL_LESS);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_NORMALIZE);
  glShadeModel(GL_SMOOTH);

  glLightfv(GL_LIGHT1, GL_AMBIENT, lightAmbient);
  glLightfv(GL_LIGHT1, GL_DIFFUSE, lightDiffuse);
  glLightfv(GL_LIGHT1, GL_SPECULAR, lightSpecular);
  glLightfv(GL_LIGHT1, GL_POSITION, lightPosition);
  glEnable(GL_LIGHT1);
  glEnable(GL_LIGHTING);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();

  gluPerspective(45.0f,(GLfloat)kWindowWidth/(GLfloat)kWindowHeight,0.1f,100.0f);
  glMatrixMode(GL_MODELVIEW);
}

GLvoid DrawGLScene(GLvoid)
{
  ctim = clock();

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();

  calculateSurface();

  glTranslatef(-sx/2, 7.5f, -27.5f-sy/2);
  glRotatef(30.0f, 1.0f, 0.0f, 0.0f);
  glBegin(GL_TRIANGLES);
  glColor3f(1.0f, 1.0f, 1.0f);
  glMaterialfv(GL_FRONT, GL_AMBIENT, materialAmbient);
  glMaterialfv(GL_FRONT, GL_DIFFUSE, materialDiffuse);
  glMaterialfv(GL_FRONT, GL_SPECULAR, materialSpecular);
  int ni = 0;
  for (int j = 0; j < (ny-1); j++)
    for (int i = 0; i < (nx-1); i++)
	{
	  glNormal3f(normals[ni][0], normals[ni][1], normals[ni][2]);
	  glVertex3f(dx*i, mesh[i][j], dy*j);
	  glVertex3f(dx*(i+1), mesh[i+1][j], dy*j);
	  glVertex3f(dx*(i+1), mesh[i+1][j+1], dy*(j+1));
	  glNormal3f(normals[ni+1][0], normals[ni+1][1], normals[ni+1][2]);
	  glVertex3f(dx*i, mesh[i][j], dy*j);
	  glVertex3f(dx*(i+1), mesh[i+1][j+1], dy*(j+1));
	  glVertex3f(dx*i, mesh[i][j+1], dy*(j+1));
	  ni += 2;
	}
  for (int i = 0; i < (nx-1); i++)
  {
    glNormal3f(0.0f, 0.0f, 1.0f);
	glVertex3f(dx*i, -2.0f, dy*(ny-1));
	glVertex3f(dx*i, mesh[i][ny-1], dy*(ny-1));
	glVertex3f(dx*(i+1), mesh[i+1][ny-1], dy*(ny-1));
	glNormal3f(0.0f, 0.0f, 1.0f);
	glVertex3f(dx*i, -2.0f, dy*(ny-1));
	glVertex3f(dx*(i+1), mesh[i+1][ny-1], dy*(ny-1));
	glVertex3f(dx*(i+1), -2.0f, dy*(ny-1));
  }
  glEnd();

  glFlush();
  glutSwapBuffers();
  glutPostRedisplay();

  stime += dt;
  while (GLfloat(clock()-ctim) < cdt);
}

GLvoid ReSizeGLScene(int Width, int Height)
{
  glViewport (0, 0, (GLsizei) Width, (GLsizei) Height);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();

  gluPerspective(45.0, (GLfloat) Width / (GLfloat) Height, 0.1, 100.0);
   
  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
}