/*
   Surface by identifying points (x,y,z) with f(x,y,z)=C for a single-valued function
   of three variables.  Useful in many places in science.
   
   Rough working start on a project that might be of interest
   
   Sample for introductory computer graphics course -- copyright (c) 2000, Steve Cunningham
*/
#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef GLfloat point3[3];

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

//  Define the parameters of the volume grid
#define XSIZE 40		// larger values overrun working memory space
#define YSIZE XSIZE
#define ZSIZE XSIZE
#define MINX	-3.0
#define MAXX	 3.0
#define MINY	-3.0
#define MAXY	 3.0
#define MINZ	-3.0
#define MAXZ	 3.0

#define CHOICE 2

static GLfloat vertices[XSIZE][ZSIZE];
GLfloat ep = 20.0;
GLfloat C = 4.0;		// f(x,y,z)=C is our surface

// function prototypes follow
GLfloat XX( int );
GLfloat YY( int );
GLfloat ZZ( int );
void myinit( void );
void surface( float );
void cube(float,float,float,float);
void display( void );
void reshape( int ,int );
void keyboard(unsigned char, int, int );

// functions to return X, Y, and Z values for array indices
GLfloat XX(int i) {
     return (MINX+((MAXX-MINX)/(float)(XSIZE-1))*(float)(i));
     }

GLfloat YY(int i) {
     return (MINY+((MAXY-MINY)/(float)(YSIZE-1))*(float)(i));
     }

GLfloat ZZ(int i) {
     return (MINZ+((MAXZ-MINZ)/(float)(ZSIZE-1))*(float)(i));
     }

void myinit(void)
{
        GLfloat light_pos0[]={  0.0, 10.0, 10.0,  1.0 }; // first light over z-axis
        GLfloat light_col0[]={  1.0,  0.0,  0.0,  1.0 }; // and red
        GLfloat amb_color0[]={  0.3,  0.0,  0.0,  1.0 }; // even ambiently
        
        GLfloat light_pos1[]={  5.0, 10.0, -7.26, 1.0 }; // second light back/right
        GLfloat light_col1[]={  0.0,  1.0,  0.0,  1.0 }; // and green
        GLfloat amb_color1[]={  0.0,  0.3,  0.0,  1.0 }; // even ambiently
        
        GLfloat light_pos2[]={ -5.0, 10.0, -7.26, 1.0 }; // third light back/left
        GLfloat light_col2[]={  0.0,  0.0,  1.0,  1.0 }; // and blue
        GLfloat amb_color2[]={  0.0,  0.0,  0.3,  1.0 }; // even ambiently
        
//      set up overall light data
        GLfloat mat_specular[]  ={ 0.8, 0.8, 0.8, 1.0 };
        long i = 1;

        glClearColor( 0.5, 0.5, 0.5, 0.0 );
        glShadeModel(GL_SMOOTH);
        glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_specular );
        
        glLightfv(GL_LIGHT0, GL_POSITION, light_pos0 ); // light 0
        glLightfv(GL_LIGHT0, GL_AMBIENT, amb_color0 );
        glLightfv(GL_LIGHT0, GL_SPECULAR, light_col0 );
        glLightfv(GL_LIGHT0, GL_DIFFUSE, light_col0 );
        
        glLightfv(GL_LIGHT1, GL_POSITION, light_pos1 ); // light 1
        glLightfv(GL_LIGHT1, GL_AMBIENT, amb_color1 );
        glLightfv(GL_LIGHT1, GL_SPECULAR, light_col1 );
        glLightfv(GL_LIGHT1, GL_DIFFUSE, light_col1 );
        
        glLightfv(GL_LIGHT2, GL_POSITION, light_pos2 ); // light 2
        glLightfv(GL_LIGHT2, GL_AMBIENT, amb_color2 );
        glLightfv(GL_LIGHT2, GL_SPECULAR, light_col2 );
        glLightfv(GL_LIGHT2, GL_DIFFUSE, light_col2 );
        
        glLightModeliv(GL_LIGHT_MODEL_TWO_SIDE, &i ); // two-sided lighting

/* attributes */

        glEnable(GL_LIGHTING);   // so lighting models are used
        glEnable(GL_LIGHT0);     // we'll use LIGHT0
        glEnable(GL_LIGHT1);     // ... and LIGHT1
        glEnable(GL_LIGHT2);     // ... and LIGHT2
        glEnable(GL_DEPTH_TEST); // allow z-buffer display
        glEnable(GL_NORMALIZE);  // make normal vectors 1-unit long after transform

	glNewList(1, GL_COMPILE);	// speed up the surface process and avoid recomputation
		surface(C);
	glEndList();
}

void surface(float C)
{
	// define your function of three variables here -- this is a simple one
//	#define f(a,b,c) a*a+b*b+c*c
	#define f(a,b,c) a*b*c

    int  i, j, k;
    float x, x1, y, y1, z, z1;
    float p1, p2, p3, p4, p5, p6, p7, p8;
    float rad;

//  Set up the surface color
    GLfloat white[] = {1.0, 1.0, 1.0, 1.0};
    GLfloat yellow[]= {1.0, 1.0, 0.0, 1.0};
    GLfloat mat_shininess[]={ 30.0 };

//	Draw a bounding cube for the space where the surface lies
	glBegin(GL_LINE_LOOP);
		glVertex3f(MINX,MINY,MINZ); glVertex3f(MINX,MAXY,MINZ);
		glVertex3f(MAXX,MAXY,MINZ); glVertex3f(MAXX,MINY,MINZ);
	glEnd();
	glBegin(GL_LINE_LOOP);
		glVertex3f(MINX,MINY,MAXZ); glVertex3f(MINX,MAXY,MAXZ);
		glVertex3f(MAXX,MAXY,MAXZ); glVertex3f(MAXX,MINY,MAXZ);
	glEnd();
	glBegin(GL_LINES);
		glVertex3f(MINX,MINY,MINZ); glVertex3f(MINX,MINY,MAXZ);
		glVertex3f(MINX,MAXY,MINZ); glVertex3f(MINX,MAXY,MAXZ);
		glVertex3f(MAXX,MINY,MINZ); glVertex3f(MAXX,MINY,MAXZ);
		glVertex3f(MAXX,MAXY,MINZ); glVertex3f(MAXX,MAXY,MAXZ);
	glEnd();
    
//  Identify where the value of the function on the grid passes through the constant
	rad = 0.7*(MAXX-MINX)/(float)XSIZE;
    for ( i=0; i<XSIZE; i++ )
    	for (j=0; j<YSIZE; j++)
			for ( k=0; k<ZSIZE; k++ ) {
				x = XX(i); x1 = XX(i+1);
				y = YY(j); y1 = YY(j+1);
				z = ZZ(k); z1 = ZZ(k+1);
				p1 = f( x, y, z);
				p2 = f( x, y,z1);
				p3 = f(x1, y,z1);
				p4 = f(x1, y, z);
				p5 = f( x,y1, z);
				p6 = f( x,y1,z1);
				p7 = f(x1,y1,z1);
				p8 = f(x1,y1, z);
				if( ((p1-C)*(p2-C)<0.0) || ((p2-C)*(p3-C)<0.0) ||
					((p3-C)*(p4-C)<0.0) || ((p1-C)*(p4-C)<0.0) ||
					((p1-C)*(p5-C)<0.0) || ((p2-C)*(p6-C)<0.0) ||
					((p3-C)*(p7-C)<0.0) || ((p4-C)*(p8-C)<0.0) ||
					((p5-C)*(p6-C)<0.0) || ((p6-C)*(p7-C)<0.0) ||
					((p7-C)*(p8-C)<0.0) || ((p5-C)*(p8-C)<0.0) ) {
						cube(x,y,z,rad);
					}
			}

}

//	function  cube  takes the lowest point (minimal x,y,z) and side as parameters
void cube(float x, float y, float z, float d)
{
	GLUquadric *myQuad;
	
	myQuad = gluNewQuadric();
	glTranslatef(x,y,z);
	gluSphere(myQuad,d,3,4);
	glTranslatef(-x,-y,-z);
}

void display( void )
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	//           eye point        center of view      up
	gluLookAt( 0.0, ep, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);

	// We save the original viewing projection in the  viewProj  array
	glPushMatrix();
	glGetFloatv( GL_MODELVIEW_MATRIX, viewProj );

	//  Put the identity onto the modelview stack to start the viewing transform
	glLoadIdentity();

#define Angle 5.0

    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;
    }
    ch = ' ';
    
	glMultMatrixf( saveState );
	glGetFloatv( GL_MODELVIEW_MATRIX, saveState );
	glLoadIdentity();
    
    glMultMatrixf( viewProj );
    glMultMatrixf( saveState );

	glCallList(1);

    glPopMatrix();
    glutSwapBuffers();
 }

void reshape(int w,int h)
{
	glViewport(0,0,(GLsizei)w,(GLsizei)h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(30.0,1.0,1.0,100.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 '+' :
			C = C + 0.2;
			glNewList(1, GL_COMPILE);
				surface(C);
			glEndList(); break;
		case '-' :
			C = C - 0.2;
			glNewList(1, GL_COMPILE);
				surface(C);
			glEndList(); break;
		case 'f' :
			ep = ep - 0.2; break;
        case 'b' :
			ep = ep + 0.2; break;
	}
	glutPostRedisplay();
}

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("implicit surface");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);

	myinit();
	glutMainLoop();
}