/*
   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.
   This program includes the use of a 1D texture map to provide Chromadepthª
   3D viewing.
   
   X-rotations  :  q  w
   Y-rotations  :  a  s
   Z-rotations  :  z  x
   
   Sample for graphics classes -- Dr. 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 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

#define CHOICE 4

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;
	}
}

void myinit(void)
{
	long i = 1;

	GLfloat light_pos0[]={  0.0, 10.0, 10.0,  1.0 };
	GLfloat light_col0[]={  1.0,  1.0,  1.0,  1.0 };
	GLfloat amb_color0[]={  0.3,  0.3,  0.3,  1.0 };
//	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);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_specular );
	glLightfv(GL_LIGHT0, GL_POSITION, light_pos0 );
	glLightfv(GL_LIGHT0, GL_AMBIENT, amb_color0 );
	glLightfv(GL_LIGHT0, GL_SPECULAR, light_col0 );
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light_col0 );
	glLightModeliv(GL_LIGHT_MODEL_TWO_SIDE, &i );

//	attributes
	glEnable(GL_LIGHTING);   // so lighting models are used
	glEnable(GL_LIGHT0);
	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_MODULATE );
	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 );
}

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);
             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 )
{
#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.  Sounds a little like
//  Humpty-Dumpty, but...
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(30.0,1.0,1.0,100.0);
    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);

	gluLookAt( ep, ep, ep/4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0);
    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;
    }
	ch = ' ';
    
    //  NOW we apply the rest of the modeling transformation by POST-multiplying
	//  by the saved matrix, then save that and put the identity back on the stack.
	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);
}

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;
	}
}

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();
}

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);
        glutKeyboardFunc(keyboard);

        myinit();
        glutMainLoop();
}