/*
   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
   
   (c) February 2001, Steve Cunningham & Ken Brown
*/
#define OSX

#ifdef MAC
#include "glut.h"
#endif
#ifdef OSX
#include <GLUT/glut.h>
#endif
#ifdef UNIX
#include <GL/glut.h>
#endif
#ifdef PC
#include <GL/glut.h>
#endif

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define TRUE 1
#define FALSE 0

// variables to keep track of the current rotation on each axis
GLfloat camPitch = 1.0, camYaw = 0.0, camRadius = 16.0;
int mx,my; // used to calculate mouse displacement
int mdown; // is a button down?
int buttonDown; // if so, which button?

typedef GLfloat point3[3];

GLfloat 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 125  // needs to be much larger than this for best results
#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

int CHOICE = 4; // which surface to display (1-4)

// global array to hold the height of each point on the surface
static GLfloat vertices[XSIZE][YSIZE];
float ep = 12.0;

//	data for texture map
float D1, D2;
float texParms[4];
static GLuint texName;
float ramp[256][3];

// function prototypes follow
GLfloat XX( int );
GLfloat YY( int );
void makeRamp(void);
void hsv2rgb(float, float, 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 );

// two functions to return X and Y 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 makeRamp(void)
{
	int i;
	float h, s, v, r, g, b;

	// Make color ramp for 1D texture: starts at 0, ends at 240, 256 steps
	for (i=0; i<256; i++) {
		h = (float)i*240.0/255.0;
		s = 1.0; v = 1.0;
		hsv2rgb( h, s, v, &r, &g, &b );
		ramp[i][0] = r; ramp[i][1] = g; ramp[i][2] = b;
		}
}

void hsv2rgb( float h, float s, float v, float *r, float *g, float *b)
{
//	conversion from Foley et.al., fig. 13.4
//	don't worry about all possible cases because we know s = v = 1.0
//	and 0 <= h <=240

	float f, p, q, t;
	int   k;

	h = h/60.0;
	k = (int)h;
	f = h - (float)k;
	p = v * (1.0 - s);
	q = v * (1.0 - (s * f));
	t = v * (1.0 - (s * (1.0 - f)));
	switch (k) {
		case 0: *r = v; *g = t; *b = p; break;
		case 1: *r = q; *g = v; *b = p; break;
		case 2: *r = p; *g = v; *b = t; break;
		case 3: *r = p; *g = q; *b = v; break;
		case 4: *r = t; *g = p; *b = v; break;
		case 5: *r = v; *g = p; *b = q; break;
	}
}

// initialize OpenGL and all necessary global variables
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 };

        glClearColor( 0.1, 0.1, 0.1, 0.0 );
        glShadeModel(GL_SMOOTH); // use gourand shading
        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 );

        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1 ); // 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

//	set up for texture map
	makeRamp();
	glPixelStorei( GL_UNPACK_ALIGNMENT, 1 );
	glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL );
	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP );
 	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR );
	glTexParameterf( GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR );
	glTexImage1D( GL_TEXTURE_1D, 0, 3, 256, 0, GL_RGB, GL_FLOAT, ramp );
	glEnable( GL_TEXTURE_GEN_S );
	glEnable( GL_TEXTURE_1D );
}

// draw the surface
void surface(void)
{
    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);
          switch (CHOICE) {
          case 1: vertices[i][j] = x*(y+t)*(x*x-(y+t)*(y*+t)); break;
          case 2: vertices[i][j] = (x*x+2.0*y*y)/exp(x*x+2.0*y*y-t); break;
          case 3: vertices[i][j] = x*x + y*y + t*x*x*x + t*y*y*y -1; break;
          case 4: vertices[i][j] = 0.3*cos(x*x+y*y+t); break;
          }
       }

/*  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);
			 // build the normal for the triangle
             vec1[0] = XX(i+1)-XX(i);
             vec1[1] = YY(j)-YY(j);
             vec1[2] = vertices[i+1][j]-vertices[i][j];
             vec2[0] = XX(i+1)-XX(i+1);
             vec2[1] = YY(j+1)-YY(j);
             vec2[2] = vertices[i+1][j+1]-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); // hack together the normal vector...

             glVertex3f(XX(i),  YY(j),  vertices[i  ][j  ]);
             glVertex3f(XX(i+1),YY(j),  vertices[i+1][j  ]);
             glVertex3f(XX(i+1),YY(j+1),vertices[i+1][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);
             vec1[0] = XX(i+1)-YY(i);
             vec1[1] = YY(j+1)-YY(j);
             vec1[2] = vertices[i+1][j+1]-vertices[i][j];
             vec2[0] = XX(i)-XX(i+1);
             vec2[1] = YY(j+1)-YY(j+1);
             vec2[2] = vertices[i][j+1]-vertices[i+1][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);

             glVertex3f(XX(i),  YY(j),  vertices[i  ][j  ]);
             glVertex3f(XX(i+1),YY(j+1),vertices[i+1][j+1]);
             glVertex3f(XX(i),  YY(j+1),vertices[i  ][j+1]);
          glEnd();
       }
}

void display( void )
{
    static GLfloat camX, camY, camZ;

    // convert from spherical coords to cartesian
    camX = cos(camYaw) * cos(camPitch) * camRadius;
    camY = sin(camYaw) * cos(camPitch) * camRadius;
    camZ = sin(camPitch) * camRadius;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
  
//  Now build the texture information before the transformations are made
    glEnable(GL_TEXTURE_1D);
    glEnable( GL_TEXTURE_GEN_S );
    glTexGeni( GL_S, GL_TEXTURE_GEN_MODE, GL_EYE_LINEAR );
    D1 = ep + 1.0; D2 = ep + 10.0;
    texParms[0] = texParms[1] = 0.0;
    texParms[2] = -1.0/(D2-D1);
    texParms[3] = -D1/(D2-D1);
    glTexGenfv( GL_S, GL_EYE_PLANE, texParms);
    glBindTexture(GL_TEXTURE_1D, texName);

    //          eye point     center of view      up
    gluLookAt( camX, camY, camZ, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);

    surface();

    glutSwapBuffers();
 }

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

void keyboard(unsigned char key, int x, int y)
{
        switch (key)
        {
			case 'w' :
				camYaw += .08;
				if(camYaw > 6.28)
					camYaw = camYaw-6.28;
				break;
    	    case 'q' :
				camYaw -= .08;
				if(camYaw < 0.0)
					camYaw = 6.28+camYaw;
				break;
    	    case 's' :
				camPitch += .08;
				if(camPitch > 1.57)
					camPitch = 1.57;
				break;
    	    case 'a' :
				camPitch -= .08;
				if(camPitch < -1.57)
					camPitch = -1.57;
				break;
    	    case 'x' :
				camRadius += 0.5; break;
    	    case 'z' :
    			camRadius -= 0.5;
				if(camRadius < 1.0)
					camRadius = 1.0;
				break;
			case '1' :
				CHOICE = 1; break;
			case '2' :
				CHOICE = 2; break;
			case '3' :
				CHOICE = 3; break;
			case '4' :
				CHOICE = 4; break;
        }
}

void mouse(int button, int state, int x, int y)
{
	// save the mouse position for the camera rotation
	mx = x;
	my = y;
	if(state == GLUT_DOWN)
		mdown = TRUE;
	else
		mdown = FALSE;

	buttonDown = button; // remember which button is down
}

void motion(int x, int y)
{
	if(buttonDown == GLUT_LEFT_BUTTON)
	{
		// adjust the camera's yaw
		camYaw += 0.02*(GLfloat)(mx-x);

		// wrap the yaw values at 0 and 2*pi
		if(camYaw < 0.0)
				camYaw = 6.28+camYaw;
		else if(camYaw > 6.28)
				camYaw = camYaw-6.28;

		// adjust the camera's pitch
		camPitch += 0.02*(GLfloat)(y-my);

		// make sure we don't pitch too far
		if(camPitch > 1.57)
			camPitch = 1.57;
		else if(camPitch < -1.57)
			camPitch = -1.57;
	}
	else
	{
		// zoom in and out if right mouse button is down
		camRadius += 0.5*(GLfloat)(y-my);
		if(camRadius < 1.0)
			camRadius = 1.0;
	}

	mx = x;
	my = y;
}

// update the surface
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; break;
    }
    if (t > 2*3.14159) t-= 2*3.14159;
    glutPostRedisplay();
}

int 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
        glutMouseFunc(mouse);
        glutMotionFunc(motion);

        myinit();
        glutMainLoop();

        return 0;
}