//This code was created by Jordan Isaak
//Code written in March, 2002
//My website is at members.shaw.ca/jordanisaak
//My email is jordanisaak@hotmail.com

#include <GL/glut.h>
#include <math.h>
#include <stdio.h>

/*Constraint structure*/
struct CONSTRAINT
{
	int particle1, particle2;
	float restlength;
};

#define CLOTH_SIZE 5.0f
#define CONSTRAINT_UPDATE 3
#define TEXTURE_SIZE 8
#define CLOTH_RESOLUTION 16
#define TIME_STEP 0.01f
#define TIME_FACTOR TIME_STEP*TIME_STEP

/*Texture object*/
GLuint texName;
/*Texture data array*/
GLubyte texture[3 * TEXTURE_SIZE * TEXTURE_SIZE];
/*Number of particles in the cloth*/
const int numparticles = CLOTH_RESOLUTION * CLOTH_RESOLUTION;
/*Array of forces on the particles in the cloth*/
float forces[numparticles * 3];
/*Array of normals for particles in the cloth*/
float normals[numparticles * 3];
/*First particle array for cloth*/
float cloth1[numparticles * 4];
/*Second particle array for cloth*/
float cloth2[numparticles * 4];
/*Pointer to current cloth array*/
float *currentcloth;
/*Pointer to previous cloth array*/
float *previouscloth;
/*Number of constraints on cloth particle system*/
const int numconstraints = 2 * CLOTH_RESOLUTION * (CLOTH_RESOLUTION - 1)
					+ (CLOTH_RESOLUTION - 1) * (CLOTH_RESOLUTION - 1);
/*Array of constraints on the cloth particle system*/
CONSTRAINT constraints[numconstraints];

/*Position and radius of the sphere*/
float sphere[4] = {0.0f, 0.0f, 0.0f, 1.0f};

/*The gravity force*/
float gravityforce = -9.8f;

/*Variables for controlling the wind*/
float winddirection[3] = {0.0f, 0.0f, -1.0f};
float windspeed = 0.0f;
float windvector[3] = {0.0f, 0.0f, 0.0f};

/*Lighting variables*/
float lightpos[4] = {1.0f, 1.0f, 1.0f, 0.0f};
float lightambient[4] = {0.2f, 0.2f, 0.2f, 1.0f};
float lightdiffuse[4] = {1.0f, 1.0f, 1.0f, 1.0f};
float material[4] = {1.0f, 1.0f, 1.0f, 1.0f};

/*The elapsed time of the program*/
int currenttime = 0;
/*How far into the simulation the program is*/
int simulationtime = 0;

/*Keeps track of where the mouse was last*/
int lastx = 0, lasty = 0;

/*Camera control variables*/
float rotation = 0;
float elevation = 1.0f;
float cameradistance = 8.0f;
float camerapos[3] = {0.0f, 0.0f, 0.0f};

/*Controls what mouse position changes update*/
int controlmode = 0;

/*Stores the current texture used for the flag*/
int currenttexture = 0;

void generatetexture(void)
{
	int i, j;

	/*Depending on which texture is desired, generate different
	coloured checkerboard textures*/
	if(currenttexture == 0)
	{
		for(i = 0; i < TEXTURE_SIZE; i++)
		{
			for(j = 0; j < TEXTURE_SIZE; j++)
			{
				texture[3 * (i * TEXTURE_SIZE + j)] = ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 1] = 255 - ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 2] = 0;
			}
		}
	}
	else if(currenttexture == 1)
	{
		for(i = 0; i < TEXTURE_SIZE; i++)
		{
			for(j = 0; j < TEXTURE_SIZE; j++)
			{
				texture[3 * (i * TEXTURE_SIZE + j)] = 255 - ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 1] = 255 - ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 2] = ((i + j) % 2) * 255;
			}
		}
	}
	else if(currenttexture == 2)
	{
		for(i = 0; i < TEXTURE_SIZE; i++)
		{
			for(j = 0; j < TEXTURE_SIZE; j++)
			{
				texture[3 * (i * TEXTURE_SIZE + j)] = 50 + ((i + j) % 2) * 205;
				texture[3 * (i * TEXTURE_SIZE + j) + 1] = 50 + ((i + j) % 2) * 205;
				texture[3 * (i * TEXTURE_SIZE + j) + 2] = 50 + ((i + j) % 2) * 205;
			}
		}
	}
	else if(currenttexture == 3)
	{
		for(i = 0; i < TEXTURE_SIZE; i++)
		{
			for(j = 0; j < TEXTURE_SIZE; j++)
			{
				texture[3 * (i * TEXTURE_SIZE + j)] = 255 - ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 1] = 255 - ((i + j) % 2) * 255;
				texture[3 * (i * TEXTURE_SIZE + j) + 2] = 255;
			}
		}
	}


	/*Load the texture into OpenGL*/
	glBindTexture(GL_TEXTURE_2D, texName);
	glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, 
					GL_NEAREST);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, 
					GL_NEAREST);
	gluBuild2DMipmaps(GL_TEXTURE_2D, GL_RGB, TEXTURE_SIZE, TEXTURE_SIZE, GL_RGB, GL_UNSIGNED_BYTE, texture);
}

void accumulateforces(void)
{
	int i;
	float windforce;

	//Determine the wind force vector
	windvector[0] = winddirection[0] * windspeed;
	windvector[1] = winddirection[1] * windspeed;
	windvector[2] = winddirection[2] * windspeed;

	for(i = 0; i < numparticles; i++)
	{
		//Initialize forces
		forces[3 * i] = 0.0f;
		forces[3 * i + 1] = 0.0f;
		forces[3 * i + 2] = 0.0f;

		//Gravity force
		forces[3 * i + 1] += gravityforce;

		//Wind force
		windforce = windvector[0] * normals[3 * i]
					+ windvector[1] * normals[3 * i + 1]
					+ windvector[2] * normals[3 * i + 2];
		forces[3 * i] += normals[3 * i] * windforce;
		forces[3 * i + 1] += normals[3 * i + 1] * windforce;
		forces[3 * i + 2] += normals[3 * i + 2] * windforce;

		//Dampening force
		windforce = (previouscloth[4 * i] - currentcloth[4 * i]) * normals[3 * i]
					+ (previouscloth[4 * i + 1] - currentcloth[4 * i + 1]) * normals[3 * i + 1]
					+ (previouscloth[4 * i + 2] - currentcloth[4 * i + 2]) * normals[3 * i + 2];
		forces[3 * i] += 80.0f * normals[3 * i] * windforce;
		forces[3 * i + 1] += 80.0f * normals[3 * i + 1] * windforce;
		forces[3 * i + 2] += 80.0f * normals[3 * i + 2] * windforce;
	}
}

void verlet(void)
{
	int i;
	float *tempfloatptr;

	//Update each particle via verlet integration
	for(i = 0; i < numparticles; i++)
	{
		previouscloth[4 * i] = 2.0 * currentcloth[4 * i]
							-  previouscloth[4 * i]
							+ forces[3 * i] * currentcloth[4 * i + 3] * TIME_FACTOR;
		previouscloth[4 * i + 1] = 2.0 * currentcloth[4 * i + 1]
							-  previouscloth[4 * i + 1]
							+ forces[3 * i + 1] * currentcloth[4 * i + 3] * TIME_FACTOR;
		previouscloth[4 * i + 2] = 2.0 * currentcloth[4 * i + 2]
							- previouscloth[4 * i + 2]
							+ forces[3 * i + 2] * currentcloth[4 * i + 3] * TIME_FACTOR;
	}

	//Update the current and previous cloth pointers
	tempfloatptr = previouscloth;
	previouscloth = currentcloth;
	currentcloth = tempfloatptr;
}

void satisfyconstraints(void)
{
	int i;
	float sphereradius;
	float deltalength;
	float delta[3];
	float diff;
	float invmass1, invmass2;

	sphereradius = sphere[3] + 0.05;

	//Satisfy stick constraints within cloth
	for(i = 0; i < numconstraints; i++)
	{
		delta[0] = currentcloth[4 * constraints[i].particle1]
					- currentcloth[4 * constraints[i].particle2];
		delta[1] = currentcloth[4 * constraints[i].particle1 + 1]
					- currentcloth[4 * constraints[i].particle2 + 1];
		delta[2] = currentcloth[4 * constraints[i].particle1 + 2]
					- currentcloth[4 * constraints[i].particle2 + 2];

		invmass1 = currentcloth[4 * constraints[i].particle1 + 3];
		invmass2 = currentcloth[4 * constraints[i].particle2 + 3];

		deltalength = (float)sqrt(delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2]);

		diff = (deltalength - constraints[i].restlength)
			/ (deltalength * (invmass1 + invmass2));

		currentcloth[4 * constraints[i].particle1] -= invmass1 * delta[0] * diff;
		currentcloth[4 * constraints[i].particle1 + 1] -= invmass1 * delta[1] * diff;
		currentcloth[4 * constraints[i].particle1 + 2] -= invmass1 * delta[2] * diff;

		currentcloth[4 * constraints[i].particle2] += invmass2 * delta[0] * diff;
		currentcloth[4 * constraints[i].particle2 + 1] += invmass2 * delta[1] * diff;
		currentcloth[4 * constraints[i].particle2 + 2] += invmass2 * delta[2] * diff;
	}

	//Satisfy constraint that cloth points can't enter the sphere
	for(i = 0; i < numparticles; i++)
	{
		delta[0] = currentcloth[4 * i] - sphere[0];
		delta[1] = currentcloth[4 * i + 1] - sphere[1];
		delta[2] = currentcloth[4 * i + 2] - sphere[2];

		deltalength = delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2];

		//If the particle is in the sphere, move it to outside the sphere
		if(deltalength < sphereradius * sphereradius)
		{
			deltalength = (float)sqrt(deltalength);
			currentcloth[4 * i] += delta[0] * (sphereradius - deltalength) / sphereradius;
			currentcloth[4 * i + 1] += delta[1] * (sphereradius - deltalength) / sphereradius;
			currentcloth[4 * i + 2] += delta[2] * (sphereradius - deltalength) / sphereradius;
		}
	}
}

void calculatenormals(void)
{
	int i, j;
	float vectorlength;
	float trianglenormal[3];
	float v1[3], v2[3];

	//Set all normals to 0
	for(i = 0; i < numparticles; i++)
	{
		normals[3 * i] = 0.0f;
		normals[3 * i + 1] = 0.0f;
		normals[3 * i + 2] = 0.0f;
	}

	//Calculate the normal of each triangle
	//in the mesh and factor it into all the vertices
	//that are a part of it
	for(i = 0; i < CLOTH_RESOLUTION - 1; i++)
	{
		for(j = 0; j < CLOTH_RESOLUTION - 1; j++)
		{
			v1[0] = currentcloth[4 * (i * CLOTH_RESOLUTION + j)]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j)];
			v1[1] = currentcloth[4 * (i * CLOTH_RESOLUTION + j) + 1]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j) + 1];
			v1[2] = currentcloth[4 * (i * CLOTH_RESOLUTION + j) + 2]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j) + 2];

			v2[0] = currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1)]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j)];
			v2[1] = currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1) + 1]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j) + 1];
			v2[2] = currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1) + 2]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j) + 2];

			trianglenormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
			trianglenormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
			trianglenormal[2] = v1[0] * v2[1] - v1[1] * v2[0];

			normals[3 * (i * CLOTH_RESOLUTION + j)] += trianglenormal[0];
			normals[3 * (i * CLOTH_RESOLUTION + j) + 1] += trianglenormal[1];
			normals[3 * (i * CLOTH_RESOLUTION + j) + 2] += trianglenormal[2];

			normals[3 * (i * CLOTH_RESOLUTION + j + 1)] += trianglenormal[0];
			normals[3 * (i * CLOTH_RESOLUTION + j + 1) + 1] += trianglenormal[1];
			normals[3 * (i * CLOTH_RESOLUTION + j + 1) + 2] += trianglenormal[2];

			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j)] += trianglenormal[0];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j) + 1] += trianglenormal[1];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j) + 2] += trianglenormal[2];

			v1[0] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1)]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1)];
			v1[1] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 1]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1) + 1];
			v1[2] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 2]
					- currentcloth[4 * (i * CLOTH_RESOLUTION + j + 1) + 2];

			v2[0] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j)]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1)];
			v2[1] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j) + 1]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 1];
			v2[2] = currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j) + 2]
					- currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 2];

			trianglenormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
			trianglenormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
			trianglenormal[2] = v1[0] * v2[1] - v1[1] * v2[0];

			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j + 1)] += trianglenormal[0];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 1] += trianglenormal[1];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j + 1) + 2] += trianglenormal[2];

			normals[3 * (i * CLOTH_RESOLUTION + j + 1)] += trianglenormal[0];
			normals[3 * (i * CLOTH_RESOLUTION + j + 1) + 1] += trianglenormal[1];
			normals[3 * (i * CLOTH_RESOLUTION + j + 1) + 2] += trianglenormal[2];

			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j)] += trianglenormal[0];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j) + 1] += trianglenormal[1];
			normals[3 * ((i + 1) * CLOTH_RESOLUTION + j) + 2] += trianglenormal[2];
		}
	}


	//Normalize vertex normals
	for(i = 0; i < numparticles; i++)
	{
		vectorlength = 1.0f / (float)sqrt(normals[3 * i] * normals[3 * i]
							+ normals[3 * i + 1] * normals[3 * i + 1]
							+ normals[3 * i + 2] * normals[3 * i + 2]);
		normals[3 * i] *= vectorlength;
		normals[3 * i + 1] *= vectorlength;
		normals[3 * i + 2] *= vectorlength;
	}
}

void updatecloth(void)
{
	int i;

	//Accumulate forces on the particles
	accumulateforces();
	//Do verlet integration
	verlet();
	//Attempt to satisfy all the constraints on the system
	for(i = 0; i < CONSTRAINT_UPDATE; i++)
		satisfyconstraints();
	//Calculate the triangle normals
	calculatenormals();
}

void initializecloth(int setuptype)
{
	int i, j, currentconstraint = 0;
	float restlength;
	float diagonallength;

	//Initialize cloth position and point weights
	if(setuptype == 0)
	{
		for(i = 0; i < CLOTH_RESOLUTION; i++)
		{
			for(j = 0; j < CLOTH_RESOLUTION; j++)
			{
				cloth1[4 * (CLOTH_RESOLUTION * i + j)] = -CLOTH_SIZE / 2.0f + j * CLOTH_SIZE / CLOTH_RESOLUTION;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 1] = 2.0f;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 2] = -i * CLOTH_SIZE / CLOTH_RESOLUTION;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 3] = 1.0f;
			}
		}
	}
	else if(setuptype == 1)
	{
		for(i = 0; i < CLOTH_RESOLUTION; i++)
		{
			for(j = 0; j < CLOTH_RESOLUTION; j++)
			{
				cloth1[4 * (CLOTH_RESOLUTION * i + j)] = 0.0f;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 1] = -CLOTH_SIZE / 2.0f + j * CLOTH_SIZE / CLOTH_RESOLUTION;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 2] = -i * CLOTH_SIZE / CLOTH_RESOLUTION;
				cloth1[4 * (CLOTH_RESOLUTION * i + j) + 3] = 1.0f;
			}
		}
	}

	//Initialize second cloth array to match first
	for(i = 0; i < 4 * numparticles; i++)
		cloth2[i] = cloth1[i];

	//Initialize 2 corner points to infinite weight
	cloth1[3] = 0.0f;
	cloth1[4 * (CLOTH_RESOLUTION - 1) + 3] = 0.0f;

	for(i = 0; i < 4 * CLOTH_RESOLUTION * CLOTH_RESOLUTION; i++)
		cloth2[i] = cloth1[i];

	//Set up pointers for the previous frame and the current frame for the cloth
	currentcloth = cloth1;
	previouscloth = cloth2;

	restlength = CLOTH_SIZE / CLOTH_RESOLUTION;

	//Initialize cloth constraints
	//Constraints in the horizontal direction along the cloth
	for(i = 0; i < CLOTH_RESOLUTION; i++)
	{
		for(j = 0; j < CLOTH_RESOLUTION - 1; j++)
		{
			constraints[currentconstraint].particle1 = i * CLOTH_RESOLUTION + j;
			constraints[currentconstraint].particle2 = i * CLOTH_RESOLUTION + j + 1;
			constraints[currentconstraint].restlength = restlength;
			currentconstraint++;
		}
	}
	//Constraints in the vertical direction along the cloth
	for(i = 0; i < CLOTH_RESOLUTION - 1; i++)
	{
		for(j = 0; j < CLOTH_RESOLUTION; j++)
		{
			constraints[currentconstraint].particle1 = i * CLOTH_RESOLUTION + j;
			constraints[currentconstraint].particle2 = (i + 1) * CLOTH_RESOLUTION + j;
			constraints[currentconstraint].restlength = restlength;
			currentconstraint++;
		}
	}

	//A single diagonal constraint across every square in the cloth grid
	diagonallength = (float) sqrt(2.0f * restlength * restlength);
	for(i = 0; i < CLOTH_RESOLUTION - 1; i++)
	{
		for(j = 0; j < CLOTH_RESOLUTION - 1; j++)
		{
			constraints[currentconstraint].particle1 = i * CLOTH_RESOLUTION + j;
			constraints[currentconstraint].particle2 = (i + 1) * CLOTH_RESOLUTION + j + 1;
			constraints[currentconstraint].restlength = diagonallength;
			currentconstraint++;
		}
	}
}

void drawcloth(void)
{
	int i, j;

	//Draw a series of triangle strips for the cloth
	for(i = 0; i < CLOTH_RESOLUTION - 1; i++)
	{
		glBegin(GL_TRIANGLE_STRIP);
		for(j = 0; j < CLOTH_RESOLUTION; j++)
		{
			glNormal3fv(&(normals[3 * ((i + 1) * CLOTH_RESOLUTION + j)]));
			glTexCoord2f(((float)j) / (CLOTH_RESOLUTION - 1), ((float)(i + 1)) / (CLOTH_RESOLUTION - 1));
			glVertex3fv(&(currentcloth[4 * ((i + 1) * CLOTH_RESOLUTION + j)]));
			glNormal3fv(&(normals[3 * (i * CLOTH_RESOLUTION + j)]));
			glTexCoord2f(((float)j) / (CLOTH_RESOLUTION - 1), ((float)i) / (CLOTH_RESOLUTION - 1));
			glVertex3fv(&(currentcloth[4 * (i * CLOTH_RESOLUTION + j)]));
		}
		glEnd();
	}
}

void menu(int selection)
{
	controlmode = selection;
}

void idle(void)
{
	currenttime = glutGet(GLUT_ELAPSED_TIME);

	//Update the simulation to match the current time
	while(simulationtime < currenttime)
	{
		updatecloth();
		simulationtime += TIME_STEP * 1000;
	}
	glutPostRedisplay();
}

void init(void)
{
	initializecloth(0);

	glShadeModel(GL_SMOOTH);
	glDisable(GL_CULL_FACE);
	glEnable(GL_DEPTH_TEST);

	//Setup texturing
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
	glGenTextures(1, &texName);
	generatetexture();
	glEnable(GL_TEXTURE_2D);

	//Setup lighting
	glLightfv(GL_LIGHT0, GL_AMBIENT, lightambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, lightdiffuse);
	glLightfv(GL_LIGHT0, GL_POSITION, lightpos);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, material);
}

void display(void)
{
	//Set the camera position
	camerapos[0] = sin(rotation) * sin(elevation) * cameradistance;
	camerapos[1] = cos(elevation) * cameradistance;
	camerapos[2] = cos(rotation) * sin(elevation) * cameradistance;

	//Set up for rendering
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	gluLookAt(camerapos[0], camerapos[1], camerapos[2],
				0.0f, 0.0f, 0.0f,
				0.0f, 1.0f, 0.0f);

	//Draw the cloth
	glColor3f(1.0f, 1.0f, 1.0f);
	drawcloth();

	glDisable(GL_TEXTURE_2D);

	//Draw the wind direction
	glColor3f(0.0f, 1.0f, 0.0f);
	glDisable(GL_LIGHTING);
	glBegin(GL_LINES);
	glVertex3fv(windvector);
	glVertex3f(0.0f, 0.0f, 0.0f);
	glEnd();
	glEnable(GL_LIGHTING);
	glColor3f(1.0f, 1.0f, 1.0f);

	//Draw the sphere
	glTranslatef(sphere[0], sphere[1], sphere[2]);
	glutSolidSphere(sphere[3], 16, 16);

	glEnable(GL_TEXTURE_2D);

	glutSwapBuffers();
}

void reshape(int w, int h)
{
	glViewport(0, 0, (GLsizei) w, (GLsizei) h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(60.0, (GLfloat) w/(GLfloat) h, 0.1, 100.0);
	glMatrixMode(GL_MODELVIEW);
}

void mousemove(int x, int y)
{
	double modelviewmatrix[16], projectionmatrix[16];
	int viewport[4];
	double projectedpoint[3];
	float u;
	float veclength;

	//Move sphere
	if(controlmode == 1)
	{
		glLoadIdentity();
		gluLookAt(camerapos[0], camerapos[1], camerapos[2],
				0.0f, -0.5f, 0.0f,
				0.0f, 1.0f, 0.0f);
		glGetDoublev(GL_MODELVIEW_MATRIX, modelviewmatrix);
		glGetDoublev(GL_PROJECTION_MATRIX, projectionmatrix);
		glGetIntegerv(GL_VIEWPORT, viewport);

		gluUnProject(x, (viewport[3] - y), 1.0, modelviewmatrix, projectionmatrix, viewport,
				&(projectedpoint[0]), &(projectedpoint[1]), &(projectedpoint[2]));

		u = camerapos[1] / (camerapos[1] - projectedpoint[1]);

		sphere[0] = camerapos[0] + u * (projectedpoint[0] - camerapos[0]);
		sphere[2] = camerapos[2] + u * (projectedpoint[2] - camerapos[2]);
	}
	//Move camera
	else if(controlmode == 0)
	{
		rotation += (x - lastx) / 100.0f;
		elevation += (y - lasty) / 100.0f;

		if(elevation < 0.1f)
			elevation = 0.1f;
		if(elevation > 3.0f)
			elevation = 3.0f;

		camerapos[0] = sin(rotation) * sin(elevation) * cameradistance;
		camerapos[1] = cos(elevation) * cameradistance;
		camerapos[2] = cos(rotation) * sin(elevation) * cameradistance;
	}
	//Change wind direction
	else if(controlmode == 2)
	{
		glLoadIdentity();
		gluLookAt(camerapos[0], camerapos[1], camerapos[2],
				0.0f, -0.5f, 0.0f,
				0.0f, 1.0f, 0.0f);
		glGetDoublev(GL_MODELVIEW_MATRIX, modelviewmatrix);
		glGetDoublev(GL_PROJECTION_MATRIX, projectionmatrix);
		glGetIntegerv(GL_VIEWPORT, viewport);

		gluUnProject(x, (viewport[3] - y), 1.0, modelviewmatrix, projectionmatrix, viewport,
				&(projectedpoint[0]), &(projectedpoint[1]), &(projectedpoint[2]));

		//Calculate wind direction
		u = camerapos[1] / (camerapos[1] - projectedpoint[1]);
		winddirection[0] = camerapos[0] + u * (projectedpoint[0] - camerapos[0]);
		winddirection[2] = camerapos[2] + u * (projectedpoint[2] - camerapos[2]);

		//Normalize wind direction
		veclength = (float) sqrt(winddirection[0] * winddirection[0]
								+ winddirection[2] * winddirection[2]);
		winddirection[0] /= veclength;
		winddirection[2] /= veclength;
	}

	lastx = x;
	lasty = y;
	idle();
}

void mousedown(int button, int state, int x, int y)
{
	lastx = x;
	lasty = y;
}

void keyboard (unsigned char key, int x, int y)
{
	switch (key)
	{
		case 27:
			exit(0);
			break;
		case '+':
		case '=':
			windspeed += 1.0f;
			break;
		case '-':
			windspeed -= 1.0f;
			break;
		case 'f':
			initializecloth(1);
			break;
		case 'c':
			initializecloth(0);
			break;
		case 'z':
			cameradistance += 0.5f;
			break;
		case 'x':
			cameradistance -= 0.5f;
			if(cameradistance < 4.0f)
				cameradistance = 4.0f;
			break;
		case 't':
			currenttexture++;
			currenttexture %= 4;
			generatetexture();
			break;
		default:
			break;
	}
}

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

	//Setup menu
	glutCreateMenu(menu);
	glutAddMenuEntry("Mouse changes view direction", 0);
	glutAddMenuEntry("Mouse changes sphere position", 1);
	glutAddMenuEntry("Mouse changes wind direction", 2);
	glutAttachMenu(GLUT_RIGHT_BUTTON);

	//Set up control functions
	glutIdleFunc(idle);
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);
	glutMotionFunc(mousemove);
	glutMouseFunc(mousedown);
	glutMainLoop();
	return 0; 
}