/*
   Plane surface with standard form Ax+By+Cz+D=0 in a volume, with a function of
   three variables whose domain includes the volume and with the value of the function
   in the plane region shown in color, with the ability to move the surface around and
   see all parts
   
   X-rotations  :  q  w
   Y-rotations  :  a  s
   Z-rotations  :  z  x
   
   Sample for intro computer graphics course.  © 2000, Steve Cunningham
*/
#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

static char ch; // the character from the keyboard that controls the rotation
static GLfloat saveState[16] = {1.0,0.0,0.0,0.0,
                                0.0,1.0,0.0,0.0,
                                0.0,0.0,1.0,0.0,
                                0.0,0.0,0.0,1.0},
               viewProj[16]; // hold transforms

typedef GLfloat point3[3];

// define coefficients of plane and the function of three variables we will explore
float A=0.0, B=0.0, D=1.0;
float C=1.0;	//	cannot be zero since we compute the z-values of the plane from here

//  Define the parameters of the surface grid
//  For an N x M grid you need XSIZE = N+1 and ZSIZE = M+1
#define XSIZE 161
#define YSIZE XSIZE
#define ZSIZE XSIZE

#define MINX  -4.0
#define MAXX   4.0
#define MINY  -4.0
#define MAXY   4.0
#define MINZ  -4.0
#define MAXZ   4.0

point3 vertices[XSIZE][YSIZE];
float fnvalmax = -1000.0, fnvalmin = 1000.0;

void box(void);
float f(float,float,float);
float function(float,float);
void getColor(float, float*, float*, float*);
void myinit(void);
void surface(void);
void display(void);
void reshape(int,int);
void keyboard(unsigned char, int, int);
void animate(void);

void box(void) {
	glLineWidth(3.0);
	glColor3f(1.0,1.0,1.0);
	glBegin(GL_LINE_STRIP);
		glVertex3f(MINX,MINY,MINZ);
		glVertex3f(MINX,MINY,MAXZ);
		glVertex3f(MINX,MAXY,MAXZ);
		glVertex3f(MINX,MAXY,MINZ);
		glVertex3f(MINX,MINY,MINZ);
	glEnd();
	glBegin(GL_LINE_STRIP);
		glVertex3f(MAXX,MINY,MINZ);
		glVertex3f(MAXX,MINY,MAXZ);
		glVertex3f(MAXX,MAXY,MAXZ);
		glVertex3f(MAXX,MAXY,MINZ);
		glVertex3f(MAXX,MINY,MINZ);
	glEnd;
	glBegin(GL_LINES);
		glVertex3f(MINX,MINY,MINZ);
		glVertex3f(MAXX,MINY,MINZ);
	glEnd();
	glBegin(GL_LINES);
		glVertex3f(MINX,MINY,MAXZ);
		glVertex3f(MAXX,MINY,MAXZ);
	glEnd();
	glBegin(GL_LINES);
		glVertex3f(MINX,MAXY,MAXZ);
		glVertex3f(MAXX,MAXY,MAXZ);
	glEnd();
	glBegin(GL_LINES);
		glVertex3f(MINX,MAXY,MINZ);
		glVertex3f(MAXX,MAXY,MINZ);
	glEnd;
}

float f(float x, float y, float z) {
	float result;
	result = sin(x*y*z);
	return result;
}

void getColor( float c, float *r, float *g, float *b )
{
	float yval;
	
	yval = (c - fnvalmin)/(fnvalmax - fnvalmin);
	if (yval < 0.30) { *r=yval/0.3; *g=0.0; *b=0.0; return; }
	if ((yval >= 0.30) && (yval < 0.89))
		{ *r=1.0; *g=(yval-0.3)/0.59; *b=0.0; return; }
	if (yval >= 0.89) { *r=1.0; *g=1.0; *b=(yval-0.89)/0.11; }
	return;
}

float function( float x, float y)
{
	return (-D-A*x-B*y)/C;
	}
	
void myinit(void)
{
int i, j, k;
float x,y,z,c;

//	set up initial data
	glClearColor( 0.0, 0.0, 1.0, 0.0 );
	glEnable(GL_DEPTH_TEST);
	for (i=0; i<XSIZE; i++)
		for (j=0; j<YSIZE; j++)
			for (k=0; k<ZSIZE; k++) {
				x=MINX + (MAXX-MINX)*( (float)i / (float) (XSIZE-1) );
				y=MINY + (MAXY-MINY)*( (float)j / (float) (YSIZE-1) );
				z=MINZ + (MAXZ-MINZ)*( (float)j / (float) (ZSIZE-1) );
				c = f(x,y,z);
				if (c<fnvalmin) fnvalmin = c;
				if (c>fnvalmax) fnvalmax = c;
			}
}

void surface(void)
{
//  define a point data type and general variables
    point3 myColor;
    int  i, j;
    float x, y, z, c;

//  Calculate the points of the surface on the grid
    for ( i=0; i<XSIZE; i++ )
       for ( j=0; j<YSIZE; j++ )
       {
          vertices[i][j][0]=x=MINX + (MAXX-MINX)*( (float)i / (float) (XSIZE-1) );
          vertices[i][j][1]=y=MINY + (MAXY-MINY)*( (float)j / (float) (YSIZE-1) );
          vertices[i][j][2]=z=function(x,y);
          // if ((i<2)&&(j<2)) printf("i %d, j %d -> %f\n",i,j,z);
       }

//  draw the plane surface -- note that are drawing quads instead of triangles
//	because the surface is a plane
	for ( i=0; i<XSIZE-1; i++ )
		for ( j=0; j<YSIZE-1; j++ )	// for each quad in the domain
		{	//	compute "average" point in the quad
			x = (vertices[i][j][0] + vertices[i+1][j][0] + 
				vertices[i][j+1][0] + vertices[i+1][j+1][0])/4.0;
			y = (vertices[i][j][1] + vertices[i+1][j][1] + 
				vertices[i][j+1][1] + vertices[i+1][j+1][1])/4.0;
			z = (vertices[i][j][2] + vertices[i+1][j][2] + 
				vertices[i][j+1][2] + vertices[i+1][j+1][2])/4.0;
			c = f(x,y,z);	//	compute function at the "average" point
			getColor(c, &myColor[0], &myColor[1], &myColor[2]);
			glColor3fv(myColor);	// color determined by the value of the function
			if ((i==0)&&(j==0)) glColor3f(1.0,0.0,0.0);	// GLITCH IN THIS CASE !!!!!
			glBegin(GL_POLYGON);
				glVertex3fv(vertices[i  ][j  ]);
				glVertex3fv(vertices[i+1][j  ]);
				glVertex3fv(vertices[i+1][j+1]);
				glVertex3fv(vertices[i  ][j+1]);
			glEnd();
		}
}

void display( void )
{
#define Angle 5.0
//  The GL_MODELVIEW_MATRIX starts out including the viewing transformation,
//  so we'll first take that off, then construct the modeling transformation,
//  then take *THAT* off, then put them back together.
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glPushMatrix();
    glGetFloatv( GL_MODELVIEW_MATRIX, viewProj );
    glLoadIdentity();
    
    switch(ch)
    {
        case 'q':
            glRotatef( Angle, 1.0, 0.0, 0.0); break;
        case 'w':
            glRotatef(-Angle, 1.0, 0.0, 0.0); break;
        case 'a':
            glRotatef( Angle, 0.0, 1.0, 0.0); break;
        case 's':
            glRotatef(-Angle, 0.0, 1.0, 0.0); break;
        case 'z':
            glRotatef( Angle, 0.0, 0.0, 1.0); break;
        case 'x':
            glRotatef(-Angle, 0.0, 0.0, 1.0); break;
    }
    
    //  NOW we apply the rest of the modeling transformation by POST-multiplying
    //  by the saved matrix, and then save that and put the identity back on the
    //  stack.  *whew*
    glMultMatrixf( saveState );
    glGetFloatv( GL_MODELVIEW_MATRIX, saveState );
    glLoadIdentity();
    
    //  And finally rebuild the overall modelview matrix by multiplying by the
    //  viewing transformation and then by the modeling transformation
    glMultMatrixf( viewProj );
    glMultMatrixf( saveState );

    box();
    surface();

    //  And put ourselves back into the original GL_MODELVIEW_MATRIX by popping
    //  off the new work.  The stack again has one thing on it, and that is the
    //  viewing transformation.
    glPopMatrix();
    glutSwapBuffers();
 }

void reshape(int w,int h)
{
	glViewport(0,0,(GLsizei)w,(GLsizei)h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(60.0,1.0,1.0,40.0);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//             eye point     center of view      up
	gluLookAt( 10.0,  10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
}

void keyboard(unsigned char key, int x, int y)
{
	ch = ' ';
	switch (key)
	{
        case 'q' :    // rotate around X; q = positive, w = negative
        case 'w' :
        case 'a' :    // rotate around Y; a = positive, s = negative
        case 's' :
        case 'z' :    // rotate around Z; z = positive, x = negative
        case 'x' :
           ch = key; break;
        case '+':
        	D += 0.2; break;
        case '-':
        	D -= 0.2; break;
	}
	glutPostRedisplay(); /* perform display again */
}

void animate(void)
{
}

void main(int argc, char** argv)
{
// Standard GLUT initialization
	glutInit(&argc,argv);
	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(70,70);
	glutCreateWindow("function values on plane in space");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutIdleFunc(animate);
	glutKeyboardFunc(keyboard);

	myinit();
	glutMainLoop();
}