/*
	Project example - heat flow in a thin rectangular body
	
	Assumptions are simple -- body is modeled as a set of small squares, and heat flows
	between squares based on the relative heats of the two squares and a diffusion
	function that includes a constant that models the heat conductivity of the bodies.
	
	The diffusion function mimics the error-diffusion approach of some standard halftoning
	techniques
*/
#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define ROWS 10			//	body is assumed to be ROWSxCOLS (unitless) squares
#define COLS 30
#define AMBIENT 25.0;	//	initial ambient temperature in degrees Celsius
#define HOT	50.0		//	hot temperature of heat-source cell
#define COLD 0.0		//	cold temperature of cold-sink cell
#define NHOTS 4
#define NCOLDS 5

#define BOTH
#undef  GEOM
#undef  COLOR

GLfloat angle = 0.0;
GLfloat temps[ROWS][COLS], back[ROWS+2][COLS+2];
GLfloat theta = 0.0, vp = 24.0;
int hotspots[NHOTS][2] = { {ROWS/2,0},{ROWS/2-1,0},{ROWS/2-2,0},{0,3*COLS/4} };
int coldspots[NCOLDS][2] = { {ROWS-1,COLS/3}, {ROWS-1,1+COLS/3}, {ROWS-1,2+COLS/3}, {ROWS-1,3+COLS/3}, {ROWS-1,4+COLS/3} };

int onceonly = 0;

void myinit(void);
void cube(void);
void display(void);
void reshape(int, int);
void animate(void);

void myinit(void)
{
	int i,j;
	
	glEnable (GL_DEPTH_TEST);

	glClearColor(0.0, 1.0, 1.0, 1.0); // cyan background

//	set up graph
	for (i=0; i<ROWS; i++) {
		for (j=0; j < COLS; j++) {
			temps[i][j] = AMBIENT;
		}
	}
	for (i=0; i<NHOTS; i++) {
//		fprintf(stderr,"HOT  -> %d, %d\n",hotspots[i][0],hotspots[i][1]);
		temps[hotspots[i][0]][hotspots[i][1]]=HOT; }
	for (i=0; i<NCOLDS; i++) {
//		fprintf(stderr,"COLD -> %d, %d\n",coldspots[i][0],coldspots[i][1]);
		temps[coldspots[i][0]][coldspots[i][1]]=COLD; }
}

//	Unit cube in first octant; used to model the behavior of individual squares
//	in the materials
void cube (void)
{
	typedef GLfloat point [3];

	point v[8] = {
		{0.0, 0.0, 0.0},
		{0.0, 0.0, 1.0},
		{0.0, 1.0, 0.0},
		{0.0, 1.0, 1.0},
		{1.0, 0.0, 0.0},
		{1.0, 0.0, 1.0},
		{1.0, 1.0, 0.0},
		{1.0, 1.0, 1.0} };

	glBegin (GL_QUAD_STRIP);
		glVertex3fv(v[4]);
		glVertex3fv(v[5]);
		glVertex3fv(v[0]);
		glVertex3fv(v[1]);
		glVertex3fv(v[2]);
		glVertex3fv(v[3]);
		glVertex3fv(v[6]);
		glVertex3fv(v[7]);
	glEnd();

	glBegin (GL_QUAD_STRIP);
		glVertex3fv(v[1]);
		glVertex3fv(v[3]);
		glVertex3fv(v[5]);
		glVertex3fv(v[7]);
		glVertex3fv(v[4]);
		glVertex3fv(v[6]);
		glVertex3fv(v[0]);
		glVertex3fv(v[2]);
	glEnd();
	
#ifdef GEOM
	glColor3f(0.0,0.0,0.0);
	glBegin (GL_LINE_STRIP);
		glVertex3fv(v[1]), glVertex3fv(v[3]);
		glVertex3fv(v[7]), glVertex3fv(v[5]);
		glVertex3fv(v[1]);
	glEnd();
#endif
}

void display( void )
{
	#define SCALE 10.0
	
	int i,j;
	
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);  

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//         eye poin      center of view       up
	gluLookAt(vp, vp, vp/2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
    
	glPushMatrix();
    glRotatef(theta, 0.0, 0.0, 1.0);

	//	Draw the bars
	for (i = 0; i < ROWS; i++) {
		for (j = 0; j < COLS; j++) {
#ifndef GEOM
			glColor3f( temps[i][j]/HOT, 0.0, 1.0-temps[i][j]/HOT );	// hotter is redder; colder bluer
#endif
#ifdef GEOM
			glColor3f( 1.0, 0.3, 0.3 );
#endif
			glPushMatrix();
			glTranslatef((float)i-(float)ROWS/2.0,(float)j-(float)COLS/2.0,0.0);
#ifndef COLOR
			glScalef(1.0, 1.0, 0.1+4.0*temps[i][j]/HOT);		// height between 0 and HOT
#endif
			cube();
			glPopMatrix();
			}
		}
    glPopMatrix();  // undo last rotation you set
	glutSwapBuffers();
}

void reshape(int w,int h)
{
	glViewport(0,0,(GLsizei)w,(GLsizei)h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(60.0, (float)w/(float)h, 1.0, 300.0);
}

void animate(void)
{
	int i, j, m, n;
	float filter[3][3]={{ 0.0625, 0.125, 0.0625 },
						{ 0.125,  0.25,  0.125  },
						{ 0.0625, 0.125, 0.0625 } };
	
	//	update angle for rotating display
	theta += 2.0;
	if (theta > 360.0) theta = theta - 360.0;
	
	//	increment temperatures throughout the material
	for (i=0; i<ROWS; i++)	// first back temps up so you can recreate temps
		for (j=0; j<COLS; j++)
			back[i+1][j+1] = temps[i][j];
	for (i=1; i<ROWS+2; i++) {
		back[i][0]=back[i][1];
		back[i][COLS+1]=back[i][COLS];
		}
	for (j=0; j<COLS+2; j++) {
		back[0][j] = back[1][j];
		back[ROWS+1][j]=back[ROWS][j];
		}
	for (i=0; i<ROWS; i++)	// use diffusion based on back values to compute temp
		for (j=0; j<COLS; j++) {
			temps[i][j]=0.0;
			for (m=-1; m<=1; m++)
				for (n=-1; n<=1; n++)
					temps[i][j]+=back[i+1+m][j+1+n]*filter[m+1][n+1];
		}
	for (i=0; i<NHOTS; i++) {
		temps[hotspots[i][0]][hotspots[i][1]]=HOT; }
	for (i=0; i<NCOLDS; i++) {
		temps[coldspots[i][0]][coldspots[i][1]]=COLD; }
	
	//	finished updating temps; now do the display
	glutPostRedisplay(); /* perform display again */
}

void main(int argc, char** argv)
{
/* Standard GLUT initialization */
	    glutInit(&argc,argv);
	    glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH); 
	    glutInitWindowSize(500,500);
	    glutInitWindowPosition(50,50);
		glutCreateWindow("Temperature in bar");
	    glutDisplayFunc(display);
	    glutReshapeFunc(reshape);
	    glutIdleFunc(animate);
		
        myinit();
	    glutMainLoop(); /* enter event loop */
}