#include <stdlib.h>
#include "BaseGL.h"
#include "Water.h"

Water::Water()
{
	nx = 0;
	nz = 0;
	dx = 0;
	dz = 0;
	diffmat[0] = 0;
	diffmat[1] = 0;
	diffmat[2] = 0.8f;
	diffmat[3] = 1;
	c = 15.0f;
	w = 0.01f;
	outer = NULL;
	osize = 0;
	r=0;
	rix=0;
	riz=0;
}

Water::~Water()
{
}

void Water::Init(int sx, int sz, float ddx, float ddz, float r)
{
	nx = sx;
	nz = sz;
	dx = ddx;
	dz = ddz;
	wv = new sv3[nx*nz];
	wn = new sv3[nx*nz];
	wvv = new float[nx*nz];

	for (int i=0; i<nx; i++) {
		for (int j=0; j<nz; j++){
			wv[i*nx + j].set(i*dx,-1,j*dz);
			wn[i*nx + j].set(0,1,0);
			wvv[i*nx + j] = 0.0f;
		}
	}
	//EmbossSphere Initialization
	this->r = r;
	rix = (int)(r/dx+.5f);				// +.5 to round correctly
	riz = (int)(r/dz+.5f);
	osize = (2*rix+1)*(2*riz+1);	//upper bound for # if nodes in sphere-square
	outer = new int[osize];				//index to outer points int EmbossSphere.here: r const
}

void Water::Update(float timestep)
{
	int i,j;
	int r,l,u,d; // (i,j) is center,(l,j) left,(r,j) right, (i,u) up, (i,d) lower vertex
	
	//perform a finite difference step & update vertex normals
	for (i=0; i<nx; i++) {
		r = (i+1)%nx;
		l = (i == 0) ? (nx-1) : (i-1);
		for (j=0; j<nz; j++){
			u = (j+1)%nz;
			d = (j == 0) ? (nz-1) : (j-1);
			//fd-step for all velocities
			wvv[i*nx+j] = wvv[i*nx+j] + timestep*c/(dx*dz)*
  							  (-wvv[i*nx+j]*w + wv[l*nx+j].v[1]+wv[r*nx+j].v[1]+
													 wv[i*nx+u].v[1]+wv[i*nx+d].v[1]-
													 4*wv[i*nx+j].v[1]);
		}
	}

	for (i=0; i<nx; i++) {
		r = (i+1)%nx;
		l = (i == 0) ? (nx-1) : (i-1);
		for (j=0; j<nz; j++){
			u = (j+1)%nz;
			d = (j == 0) ? (nz-1) : (j-1);
			//fd-step to update positions 
			wv[i*nx+j].v[1] = wv[i*nx+j].v[1] + timestep * wvv[i*nx+j];
		}
	}

	for (i=0; i<nx; i++) {
		r = (i+1)%nx;
		l = (i == 0) ? (nx-1) : (i-1);
		for (j=0; j<nz; j++){
			u = (j+1)%nz;
			d = (j == 0) ? (nz-1) : (j-1);

			wn[i*nx+j].v[0] = (wv[r*nx+j].v[1] - wv[i*nx+j].v[1]) / dx;
			wn[i*nx+j].v[1] = 1;
			wn[i*nx+j].v[2] = (wv[i*nx+u].v[1] - wv[i*nx+j].v[1]) / dz;
		}
	}
}

void Water::Render()
{
	int i,j,k=0;
	glEnable(GL_LIGHTING);
	glBegin(GL_TRIANGLES);
	glColor3f(0,0,1);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, diffmat);
	for (i=0; i<nx-1; i++) {
		for (j=0; j<nz-1; j++){
			k = i*nx + j;
			glNormal3fv(wn[k].v);
			glVertex3fv(wv[k].v);
			glNormal3fv(wn[k+nx].v);
			glVertex3fv(wv[k+nx].v);
			glNormal3fv(wn[k+1].v);
			glVertex3fv(wv[k+1].v);

			glNormal3fv(wn[k+nx+1].v);
			glVertex3fv(wv[k+nx+1].v);
			glNormal3fv(wn[k+1].v);
			glVertex3fv(wv[k+1].v);
			glNormal3fv(wn[k+nx].v);
			glVertex3fv(wv[k+nx].v);
		}
	}
	glEnd();
	glDisable(GL_LIGHTING);

}

//Idea: Look at box around circle,traverse grid vertices & test if in interior
void Water::EmbossSphere(float *pos)	
{
	int i,j,k,l, nouter=0;
	float xs,zs,rs=r*r,totvol=0,parvol=0,diff,h; //squared
	
	i = (int)(pos[0] / (int)dx);
	j = (int)(pos[2] / (int)dz);
	memset(outer,-1,osize*sizeof(int));	//init with -1
	
	for (k = max(i-rix, 0); k < min(i+rix+1, nx); k++)
	{
		for (l = max(j-riz, 0); l < min(j+riz+1, nz); l++)
		{
			xs = (pos[0]-k*dx)*(pos[0]-k*dx);
			zs = (pos[2]-l*dz)*(pos[2]-l*dz);
			if (xs + zs <= rs)
			{
				h = pos[1]-(float)sqrt(rs-xs-zs);
				diff = wv[k*nx+l].v[1] - h;
				if (diff<0)
					continue;
				wv[k*nx+l].v[1] = h;
				wvv[k*nx+l] = 0;
				totvol += diff;
			}
			else
			{
				outer[nouter] = k*nx+l;
				nouter++;
			}
		}
	}
	parvol = totvol / nouter;
	for (int k=0;k<nouter;k++)
		wv[outer[k]].v[1] += parvol; 
}