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

#include <windows.h>
#include <mmsystem.h>


struct Point
{
	GLfloat v[3];
	GLfloat height;
	GLfloat velocity;
	int neighbors;
	int nbr[12];		// Usually six, sometimes five.  Can be higher if MERGE_NODES is false.
	float ndst[12];		// Weights for the weighted average.
	float xTex, yTex;	// Texture values.
};

const int n=3;					// facet rate
const int MERGE_NODES=1;
const int FIX_VOLUME=1;
const float PICK_RADIUS_FACTOR=0.01f;
const float turnSpeed=0.1f;
const float damping=1;
const int USE_KLUDGEY_DAMPING_ALG=1;		// This is the only way I've found to make sure
											// that the waves don't (a) peter out too soon or
											// (b) go on forever.

GLubyte myTexels[512][256][3];

const int maxTriangles=5120;     // 5*(4**(n+1)), much more if MERGE_NODES==0.

int triangles[maxTriangles][3];  
GLfloat temp[maxTriangles][3];

Point pts[maxTriangles/2+2];

GLfloat eyex=0,eyey=0,eyez=1,upx=0,upy=1,upz=0;
int midX=250, midY=250, sizeX=125, sizeY=125;


float frameSpeed=0.1f;
int wireFrame=0;
int outputFramerate=0;
int outerTexture=0;
int innerTexture=0;
int loadedTexture=1;
int activePoint=-1;
GLfloat ambientLight=0.1f;
GLfloat background=1;

//////////////////////////////// forward declaration

void showpoint(int);
float xTexCoord(Point& p);
float yTexCoord(Point& p);


//////////////////////////////// Initialization code

static int numTriangles=0;
int nextTriangle()
{
	assert(numTriangles<maxTriangles);
	return numTriangles++;
}

static int numPoints=0;
int nextPoint()
{
	assert(numPoints<maxTriangles/2+2);
	pts[numPoints].neighbors=0;
	pts[numPoints].height=1;
	pts[numPoints].velocity=0;
	return numPoints++;
}

void makeNeighbors(int a, int b)
{
	int i;
	for(i=0;i<pts[a].neighbors;i++)
		if(pts[a].nbr[i]==b)
			return;
	pts[a].nbr[pts[a].neighbors]=b;
	pts[b].nbr[pts[b].neighbors]=a;
	pts[a].neighbors++;
	pts[b].neighbors++;
}

void triangle(int a, int b, int c)
{
	int d=nextTriangle();
	triangles[d][0]=a;
	triangles[d][1]=b;
	triangles[d][2]=c;
	makeNeighbors(a,b);
	makeNeighbors(b,c);
	makeNeighbors(c,a);
}

void normal(Point &p)
{
    float d=1/sqrt(p.v[0]*p.v[0]+p.v[1]*p.v[1]+p.v[2]*p.v[2]);
	assert(d<100.0);
	p.v[0]*=d;
	p.v[1]*=d;
	p.v[2]*=d;
}


int addPoints(int a, int b)
{
	Point p;
	int i;
	p.v[0]=pts[a].v[0]+pts[b].v[0];
	p.v[1]=pts[a].v[1]+pts[b].v[1];
	p.v[2]=pts[a].v[2]+pts[b].v[2];
	normal(p);
	for(i=0;MERGE_NODES && i<numPoints;i++)
		if(p.v[0]==pts[i].v[0] && p.v[1]==pts[i].v[1] && p.v[2]==pts[i].v[2])
			return i;
	i=nextPoint();
	pts[i].v[0]=p.v[0];
	pts[i].v[1]=p.v[1];
	pts[i].v[2]=p.v[2];
	return i;
}

void divideTriangle(int a, int b, int c, int m)
{
    if(m>0)
    {
		int d=addPoints(a,b);
		int e=addPoints(b,c);
		int f=addPoints(a,c);
        divideTriangle(a, d, f, m-1);
        divideTriangle(b, d, e, m-1);
        divideTriangle(c, e, f, m-1);
        divideTriangle(d, e, f, m-1);
    }
    else
		triangle(a,b,c);
}

void tetrahedron( int m)
{
	int i,j;
    double a=.8506508;
    double b=sqrt(1-a*a);
    double PI=3.1415927;

    for(i=0;i<12;i++)
        assert(nextPoint()==i);

    pts[0].v[0]=0;       pts[0].v[1]=0;       pts[0].v[2]=1;
    pts[11].v[0]=0;      pts[11].v[1]=0;      pts[11].v[2]=-1;

    for(i=0;i<5;i++)
    {
        double t1=2*PI*i/5, t2=2*PI*(i-.5)/5;
        pts[i+1].v[0]=a*cos(t1);
        pts[i+1].v[1]=a*sin(t1);
        pts[i+1].v[2]=b;

        pts[i+6].v[0]=a*cos(t2);
        pts[i+6].v[1]=a*sin(t2);
        pts[i+6].v[2]=-b;
    }

    for(i=0;i<4;i++)
    {
        divideTriangle(0,i+1,i+2,m);
        divideTriangle(i+1,i+7,i+2,m);
        divideTriangle(i+1,i+6,i+7,m);
        divideTriangle(i+6,11,i+7,m);
    }

    divideTriangle(0,5,1,m);
    divideTriangle(5,6,1,m);
    divideTriangle(5,10,6,m);
    divideTriangle(10,11,6,m);


	for(i=0;i<numPoints;i++)
	{
		float total=0;
		for(j=0;j<pts[i].neighbors;j++)
		{
			int k=pts[i].nbr[j];
			pts[i].ndst[j]=1/sqrt(
					(pts[i].v[0]-pts[k].v[0])*(pts[i].v[0]-pts[k].v[0])
					+(pts[i].v[1]-pts[k].v[1])*(pts[i].v[1]-pts[k].v[1])
					+(pts[i].v[2]-pts[k].v[2])*(pts[i].v[2]-pts[k].v[2]));
			total+=pts[i].ndst[j];
		}
		total=1/total;
		for(j=0;j<pts[i].neighbors;j++)
			pts[i].ndst[j]*=total;
		pts[i].xTex=xTexCoord(pts[i]);
		pts[i].yTex=yTexCoord(pts[i]);

	}
}

void myinit()
{
    GLfloat mat_specular[]={1.0, 1.0, 1.0, 1.0};
    GLfloat mat_diffuse[]={1.0, 1.0, 1.0, 1.0};
    GLfloat mat_ambient[]={1.0, 1.0, 1.0, 1.0};
    GLfloat mat_shininess={100.0};
    GLfloat light_ambient[]={ambientLight, ambientLight, ambientLight, 1.0};
    GLfloat light_diffuse[]={0.9f, 0.9f, 0.9f, 1.0};
    GLfloat light_specular[]={0, 0, 0, 0};
	GLfloat light_pos[]={10.0,10.0,10.0,0.0};

/* set up ambient, diffuse, and specular components for light 0 */

    glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
    glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
    glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
	glLightfv(GL_LIGHT0, GL_POSITION, light_pos);

/* define material proerties for front face of all polygons */

    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_specular);
    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, mat_ambient);
    glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, mat_diffuse);
    glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess);

    glShadeModel(GL_SMOOTH); /*enable smooth shading */
    glEnable(GL_LIGHTING); /* enable lighting */
    glEnable(GL_LIGHT0);  /* enable light 0 */
    glEnable(GL_DEPTH_TEST); /* enable z buffer */
	glEnable(GL_COLOR_MATERIAL);

    glClearColor (background, background, background, 1.0);
	tetrahedron(n);

/////////////////// Texture code

	int i,j,k;
	FILE* f = fopen("texture1.rgb","r");
    if(f==NULL)
    {
		cout << "texture1.rgb not found!" << endl;
		loadedTexture=0;
		cout << "Texturing disabled." << endl;
		return;
    }		
	for(i=0;i<54;i++)
		assert(fscanf(f,"%c",&myTexels[0][0][0])==1);
	for(i=0;i<512;i++)
		for(j=0;j<256;j++)
			for(k=2;k>-1;k--)
				assert(fscanf(f,"%c",&myTexels[i][j][k])==1);

	fclose(f);

	glTexImage2D(GL_TEXTURE_2D, 0, 3, 512, 256, 0, 
			GL_RGB, GL_UNSIGNED_BYTE, myTexels);
    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_MIN_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

	cout << "triangles:" << numTriangles << " Points: " << numPoints << endl;
}

////////////////////////////// Initialization code ends here

void showPoint(int number, Point &p)
{
	int i;
	cout << "Point number "<< number << endl;
	cout << "(" << p.v[0] << "," << p.v[1] << "," << p.v[2] << ")" << endl;
	cout << "h:" << p.height << "   v: " << p.velocity << endl;
	cout << "neighbors: " << p.neighbors << endl;
	for(i=0;i<p.neighbors;i++)
		cout << p.nbr[i] << ':' << p.ndst[i] << endl;
}

void showPoint(int number)
{
	showPoint(number, pts[number]);
}

////////////////////////////////Debug code ends here

const float INVPI=1/3.1415927f;
float xTexCoord(Point& p)
{
	float r;
	if(p.v[0]==0)
	{
		if(p.v[1]>0)
			r=.5;
		else
			r=1;
	}
	else
	{
		r=atan(p.v[1]/p.v[0])*INVPI*.5+.25;
		if(p.v[0]<0)
			r+=.5;
	}
	return r;
}

float yTexCoord(Point& p)
{
	float r;
	if(p.v[2]==1)
		r=1;
	if(p.v[2]==-1)
		r=0;
	else
		r=atan(p.v[2]/sqrt(1-p.v[2]*p.v[2]))*INVPI+.5;
	return r;
}

int isFaceUp(int tri[])
{
  int num=0;
  if(pts[tri[0]].height>0)
    ++num;
  if(pts[tri[1]].height>0)
    ++num;
  if(pts[tri[2]].height>0)
    ++num;
  return (num>1);
}

///////////////////////////////Display code starts here

DWORD lastTime=0;
int frameCount=0;
void display(void)
{
	int i,j;

	++frameCount;
	if(frameCount>=20)
	{
		DWORD time=timeGetTime();
		if(outputFramerate)
			cout << "Framerate:"<<20000.0/(time-lastTime)<<endl;
		frameCount=0;
		lastTime=time;
	}

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	gluLookAt(eyex,eyey,eyez,0,0,0,upx,upy,upz);
	for(i=0;i<numPoints;i++)
	{
		float dst=0;
		for(j=0;j<pts[i].neighbors;j++)
			dst+=(pts[pts[i].nbr[j]].height-pts[i].height)*pts[i].ndst[j];
		pts[i].velocity+=frameSpeed*dst;
		if(USE_KLUDGEY_DAMPING_ALG)
			pts[i].height+=pts[i].velocity*(frameSpeed<0?-frameSpeed:frameSpeed);
	}
	if(!USE_KLUDGEY_DAMPING_ALG)
		for(i=0;i<numPoints;i++)
			pts[i].height+=pts[i].velocity*(frameSpeed<0?-frameSpeed:frameSpeed);
	if(FIX_VOLUME && frameSpeed!=0)
	{
		float totalVolume=0;
		for(i=0;i<numPoints;i++)
		{
			double blah=pts[i].height*pts[i].height*pts[i].height;
			if(blah<0)
				totalVolume-=blah;
			else
				totalVolume+=blah;
		}
		totalVolume=(float)pow(numPoints/totalVolume,0.33333333333);
		float velocityMult=totalVolume;
		if(velocityMult>.995f)
			velocityMult=.995f;
		for(i=0;i<numPoints;i++)
		{
			pts[i].height*=totalVolume;
			pts[i].velocity*=velocityMult;
		}
	}

	for(i=0;i<numPoints;i++)
	{
		temp[i][0]=pts[i].v[0]*pts[i].height;
		temp[i][1]=pts[i].v[1]*pts[i].height;
		temp[i][2]=pts[i].v[2]*pts[i].height;
	}
	if(wireFrame)
		glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	else
		glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
	for(i=0;i<numTriangles;i++)
	{
		if(isFaceUp(triangles[i]))
		{
			if(loadedTexture && outerTexture)
			{
				glEnable(GL_TEXTURE_2D);
				glBegin(GL_TRIANGLES);
				glColor3f(1,1,1);
				for(j=0;j<3;j++)
				{
					glNormal3fv(pts[triangles[i][j]].v);
					if(pts[triangles[i][j]].xTex-pts[triangles[i][(j+1)%3]].xTex>.5
								|| pts[triangles[i][j]].xTex-pts[triangles[i][(j+2)%3]].xTex>.5)
						glTexCoord2f(pts[triangles[i][j]].xTex-1,pts[triangles[i][j]].yTex);
					else
						glTexCoord2f(pts[triangles[i][j]].xTex,pts[triangles[i][j]].yTex);
					glVertex3fv(temp[triangles[i][j]]);
				}
				glEnd();
			}
			else
			{
				glDisable(GL_TEXTURE_2D);
				glBegin(GL_TRIANGLES);
				glColor3f(.5,.5,1);
				for(j=0;j<3;j++)
				{
					if(pts[triangles[i][j]].height>0)
					{
						glColor3f(.5,.5,1);
						glNormal3fv(pts[triangles[i][j]].v);
					}
					else
					{
						glColor3f(.5,1,.5);
						GLfloat* v=pts[triangles[i][j]].v;
						glNormal3f(-v[0],-v[1],-v[2]);
					}
					glVertex3fv(temp[triangles[i][j]]);
				}
				glEnd();
			}
		}
		else
		{
			if(loadedTexture && innerTexture)
			{
				glEnable(GL_TEXTURE_2D);
				glBegin(GL_TRIANGLES);
				glColor3f(1,1,1);
				for(j=0;j<3;j++)
				{
					GLfloat* v=pts[triangles[i][j]].v;
					glNormal3f(-v[0],-v[1],-v[2]);
					if(pts[triangles[i][j]].xTex-pts[triangles[i][(j+1)%3]].xTex>.5
								|| pts[triangles[i][j]].xTex-pts[triangles[i][(j+2)%3]].xTex>.5)
						glTexCoord2f(pts[triangles[i][j]].xTex-1,pts[triangles[i][j]].yTex);
					else
						glTexCoord2f(pts[triangles[i][j]].xTex,pts[triangles[i][j]].yTex);
					glVertex3fv(temp[triangles[i][j]]);
				}
				glEnd();
			}
			else
			{
				glDisable(GL_TEXTURE_2D);
				glBegin(GL_TRIANGLES);
				glColor3f(.5,1,.5);
				for(j=0;j<3;j++)
				{
					if(pts[triangles[i][j]].height>0)
					{
						glColor3f(.5,.5,1);
						glNormal3fv(pts[triangles[i][j]].v);
					}
					else
					{
						glColor3f(.5,1,.5);
						GLfloat* v=pts[triangles[i][j]].v;
						glNormal3f(-v[0],-v[1],-v[2]);
					}
					glVertex3fv(temp[triangles[i][j]]);
				}
				glEnd();
			}
		}
	}

	if(activePoint>-1)
	{
		glTranslatef(temp[activePoint][0],temp[activePoint][1],temp[activePoint][2]);
		glColor3f(1,0,0);
		glutSolidSphere(.1,6,6);
	}

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

int dir=1;
void keys(unsigned char key, int x, int y)
{
	if(key == '0') pts[0].velocity+=dir;
	if(key == '1') pts[1].velocity+=dir;
	if(key == '2') pts[2].velocity+=dir;
	if(key == '-') dir=-dir;
	if(key == 'w') wireFrame=!wireFrame;
	if(key == 'z') frameSpeed=0.5f;
	if(key == 'x') frameSpeed=0.25f;
	if(key == 'c') frameSpeed=0.1f;
	if(key == 'v') frameSpeed=0.01f;
	if(key == 'b') frameSpeed=0;
	if(key == 'n') frameSpeed=-0.05f;
	if(key == 'm') frameSpeed=-0.1f;
	if(key == 't') outerTexture=!outerTexture;
	if(key == 'y') innerTexture=!innerTexture;
	if(key == 'f') outputFramerate=!outputFramerate;
	if(key == 'd') 
	{
		background=1-background;
		glClearColor(background,background,background,1);
	}
}


void arrowKeys(int key, int x, int y)
{
   GLfloat factor=sqrt(1/(1+turnSpeed*turnSpeed));
   if(key == GLUT_KEY_UP) 
   {
	   GLfloat tempx=-eyex, tempy=-eyey, tempz=-eyez;
	   eyex=(eyex+turnSpeed*upx)*factor;
	   eyey=(eyey+turnSpeed*upy)*factor;
	   eyez=(eyez+turnSpeed*upz)*factor;
	   upx=(upx+turnSpeed*tempx)*factor;
	   upy=(upy+turnSpeed*tempy)*factor;
	   upz=(upz+turnSpeed*tempz)*factor;
   }
   if(key == GLUT_KEY_DOWN) 
   {
	   GLfloat tempx=-eyex, tempy=-eyey, tempz=-eyez;
	   eyex=(eyex-turnSpeed*upx)*factor;
	   eyey=(eyey-turnSpeed*upy)*factor;
	   eyez=(eyez-turnSpeed*upz)*factor;
	   upx=(upx-turnSpeed*tempx)*factor;
	   upy=(upy-turnSpeed*tempy)*factor;
	   upz=(upz-turnSpeed*tempz)*factor;
   }
   GLfloat rightx=eyey*upz-eyez*upy;
   GLfloat righty=eyez*upx-eyex*upz;
   GLfloat rightz=eyex*upy-eyey*upx;
   if(key == GLUT_KEY_LEFT) 
   {
	   eyex=(eyex+turnSpeed*rightx)*factor;
	   eyey=(eyey+turnSpeed*righty)*factor;
	   eyez=(eyez+turnSpeed*rightz)*factor;
   }
   if(key == GLUT_KEY_RIGHT)
   {
	   eyex=(eyex-turnSpeed*rightx)*factor;
	   eyey=(eyey-turnSpeed*righty)*factor;
	   eyez=(eyez-turnSpeed*rightz)*factor;
   }
   if(key == GLUT_KEY_HOME) 
   {
	   upx=(upx+turnSpeed*rightx)*factor;
	   upy=(upy+turnSpeed*righty)*factor;
	   upz=(upz+turnSpeed*rightz)*factor;
   }
   if(key == GLUT_KEY_PAGE_UP)
   {
	   upx=(upx-turnSpeed*rightx)*factor;
	   upy=(upy-turnSpeed*righty)*factor;
	   upz=(upz-turnSpeed*rightz)*factor;
   }
}

void mouse(int button, int state, int x, int y)
{
	int i;
    float rightx=eyey*upz-eyez*upy;
    float righty=eyez*upx-eyex*upz;
    float rightz=eyex*upy-eyey*upx;
	int best;
	float bestVal;
	float bestZ;
	for(i=0;i<numPoints;i++)
	{
		float xComp=temp[i][0]*rightx+temp[i][1]*righty+temp[i][2]*rightz;
		float zComp=temp[i][0]*eyex+temp[i][1]*eyey+temp[i][2]*eyez;
		float yComp=temp[i][0]*upx+temp[i][1]*upy+temp[i][2]*upz;
		float screenX=midX-sizeX*xComp;
		float screenY=midY-sizeY*yComp;
		float value=((screenX-x)*(screenX-x)+(screenY-y)*(screenY-y))*(2.5-zComp);
                if(value<sizeX*sizeY*PICK_RADIUS_FACTOR)
                    value=-zComp;
		if(i == 0 || value<bestVal)
		{
			best=i;
			bestVal=value;
			bestZ=zComp;
		}
	}
	if(button==GLUT_LEFT_BUTTON && state==GLUT_DOWN)
		pts[best].velocity+=dir;
	if(button==GLUT_MIDDLE_BUTTON && state==GLUT_DOWN)
		for(i=0;i<numPoints;i++)
		{
			float ratio=pts[i].v[0]*pts[best].v[0]+pts[i].v[1]*pts[best].v[1]
									+pts[i].v[2]*pts[best].v[2];
			if(ratio>.8)
				pts[i].velocity+=ratio*dir*.25;
		}
	if(button==GLUT_RIGHT_BUTTON && state==GLUT_DOWN)
		for(i=0;i<numPoints;i++)
		{
			float ratio=pts[i].v[0]*pts[best].v[0]+pts[i].v[1]*pts[best].v[1]
									+pts[i].v[2]*pts[best].v[2];
			if(ratio>0)
				pts[i].velocity+=ratio*dir*.125;
		}
}

void myReshape(int w, int h)
{
	midX=w/2;
	midY=h/2;
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
	if (w <= h)
	{
        glOrtho(-2,2,-2.0*h/w,2.0*h/w,-10,10);
		sizeX=w/4;
		sizeY=w/4;
	}
    else
	{
        glOrtho(-2.0*w/h,2.0*w/h,-2,2,-10,10);
		sizeX=h/4;
		sizeY=h/4;
	}
    glMatrixMode(GL_MODELVIEW);
    display();
}


void
main(int argc, char **argv)
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowSize(500, 500);
    glutCreateWindow("sphere");
	glutKeyboardFunc(keys);
    glutSpecialFunc(arrowKeys);
	glutMouseFunc(mouse);
    myinit();
    glutReshapeFunc(myReshape);
    glutDisplayFunc(display);
    glutMainLoop();
}