/*
   Mathematical surface, with the ability to move the surface around and
   see all parts.  The surface is animated -- it's from a one-parameter
   family of surfaces, and that parameter is animated by the  idle  function
   
   X-rotations  :  q  w
   Y-rotations  :  a  s
   Z-rotations  :  z  x
   
   Sample to illustrate smooth and flat shading for intro graphics class
   copyright © 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
               t=0.0;        // parameter for the family of functions

//  Define the parameters of the surface grid
//  For an N x M grid you need XSIZE = N+1 and ZSIZE = M+1
#define XSIZE 41  // small to allow us to see shading effects
#define YSIZE XSIZE
#define MINX  -3.0 // was -3..3, -3..3 with original equation -- CHOICE = 4
#define MAXX   3.0
#define MINY  -3.0
#define MAXY   3.0

#define SMOOTH
#undef  FLAT

#define CHOICE 4

static GLfloat vertices[XSIZE][YSIZE];

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

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

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

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 specular_color[]={ 0.0, 0.0, 0.0, 1.0 };
        GLfloat mat_specular[]  ={ 0.8, 0.8, 0.8, 1.0 };
        long i = 1;

        glClearColor( 0.0, 0.0, 1.0, 0.0 );
#ifdef SMOOTH
        glShadeModel(GL_SMOOTH);
#else
		glShadeModel(GL_FLAT);
#endif
        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
}

void surface(void)
{
#define f(x,y) 0.3*cos(x*x+y*y+t)		// original function
#define fx(x,y) -0.6*x*sin(x*x+y*y+t)	// partial derivative in x
#define fy(x,y) -0.6*y*sin(x*x+y*y+t)	// partial derivative in y

	point3 vec1, vec2, triNormal;
	int  i, j;
	float x, y;

//  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 };
    
//  Calculate the points of the surface on the grid
	#define EPSILON .001;
	for ( i=0; i<XSIZE; i++ )
		for ( j=0; j<YSIZE; j++ )
		{
			x = XX(i);
			y = YY(j);
			vertices[i][j] = f(x,y);
		}

//  actually draw the surface
	for ( i=0; i<XSIZE-1; i++ )
		for ( j=0; j<YSIZE-1; j++ )
		{
		// first triangle in the quad, front face
		  glMaterialfv(GL_FRONT, GL_DIFFUSE, white );
          glMaterialfv(GL_BACK, GL_DIFFUSE, yellow );
          glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );
          glBegin(GL_POLYGON);
             #ifdef  FLAT
             vec1[0] = XX(i+1)-XX(i  );
             vec1[1] = YY(j+1)-YY(j  );
             vec1[2] = vertices[i+1][j+1]-vertices[i  ][j  ];
             vec2[0] = XX(i+1)-XX(i  );
             vec2[1] = YY(j+1)-YY(j+1);
             vec2[2] = vertices[i+1][j+1]-vertices[i  ][j+1];
             triNormal[0] = -vec1[1] * vec2[2] + vec1[2] * vec2[1];
             triNormal[1] = -vec1[2] * vec2[0] + vec1[0] * vec2[2];
             triNormal[2] = -vec1[0] * vec2[1] + vec1[1] * vec2[0];
             glNormal3fv(triNormal); // hack together the normal vector...
             #endif
             #ifdef  SMOOTH
		     x = XX(i);
		     y = YY(j);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             glVertex3f(XX(i  ),YY(j  ),vertices[i  ][j  ]);
             #ifdef  SMOOTH
		     x = XX(i+1);
		     y = YY(j+1);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             glVertex3f(XX(i+1),YY(j+1),vertices[i+1][j+1]);
             #ifdef  SMOOTH
		     x = XX(i);
		     y = YY(j+1);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             glVertex3f(XX(i  ),YY(j+1),vertices[i  ][j+1]);
          glEnd();

        // second triangle in the quad, front face  *****************
          glMaterialfv(GL_FRONT, GL_DIFFUSE, white );
          glMaterialfv(GL_BACK, GL_DIFFUSE, yellow );
          glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );
          glBegin(GL_POLYGON);
             #ifdef  FLAT
             vec1[0] = XX(i  )-XX(i+1);
             vec1[1] = YY(j  )-YY(j+1);
             vec1[2] = vertices[i  ][j  ]-vertices[i+1][j+1];
             vec2[0] = XX(i  )-XX(i+1);
             vec2[1] = YY(j  )-YY(j  );
             vec2[2] = vertices[i  ][j  ]-vertices[i+1][j  ];
             triNormal[0] = -vec1[1] * vec2[2] + vec1[2] * vec2[1];
             triNormal[1] = -vec1[2] * vec2[0] + vec1[0] * vec2[2];
             triNormal[2] = -vec1[0] * vec2[1] + vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             #ifdef  SMOOTH
		     x = XX(i);
		     y = YY(j);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             glVertex3f(XX(i  ),YY(j  ),vertices[i  ][j  ]);
             #ifdef  SMOOTH
		     x = XX(i+1);
		     y = YY(j);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
             glVertex3f(XX(i+1),YY(j  ),vertices[i+1][j  ]);
             #ifdef  SMOOTH
		     x = XX(i+1);
		     y = YY(j+1);
			 vec1[0] = 1.0;
             vec1[1] = 0.0;
             vec1[2] = fx(x,y);
             vec2[0] = 0.0;
             vec2[1] = 1.0;
             vec2[2] = fy(x,y);
             triNormal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1];
             triNormal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2];
             triNormal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
             glNormal3fv(triNormal);
             #endif
            glVertex3f(XX(i+1),YY(j+1),vertices[i+1][j+1]);
          glEnd();
       }
}

void display( void )
{

//  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.  Sounds a little like
//  Humpty-Dumpty, but...

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    // 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 = ' ';
    
    //  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 );

    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(30.0,1.0,1.0,100.0);
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();
        /*           eye point        center of view      up   */
        gluLookAt( 10.0, 0.0, 10.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; i = positive, j = negative
           ch = key; break;
        case 'w' :
           ch = key; break;
        case 'a' :    // rotate around Y; k = positive, m = negative
           ch = key; break;
        case 's' :
           ch = key; break;
        case 'z' :    // rotate around Z; a = positive, s = negative
           ch = key; break;
        case 'x' :
           ch = key; break;
        }
//        glutPostRedisplay(); /* not needed because this is covered by idle */
}

void animate(void)
{
	switch (CHOICE) {
		case 1: t += 0.01; break;
		case 2: t += 0.1; break;
		case 3: t += 0.01; break;
		case 4:	t += 0.3;
				if (t > 2*3.14159) t-= 2*3.14159; 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("moving waves of math");
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutIdleFunc(animate);      // an empty function this time
        glutKeyboardFunc(keyboard); // enable keyboard callback

        myinit();

        glutMainLoop();
}