/*
   Mathematical visualization of a function from RxR to RxR -- differential equation in two
   variables, complex-valued function, etc...
*/
#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef GLfloat point3[3], vector2[2];

#undef  HEAT
#define EVEN
#undef  RAINBOW
#undef  GRAY
#undef  MULTI

#define BIGNUMBER 10000.0
point3 myColor;
float   YMIN = BIGNUMBER, YMAX = -BIGNUMBER, YRANGE;

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 20  // needs to be much larger than this for best results
#define YSIZE XSIZE
#define MINX  -5.0
#define MAXX   5.0
#define MINY  -5.0
#define MAXY   5.0
#define EPSILON 0.1

static float length[XSIZE][YSIZE];
static vector2 vectors[XSIZE][YSIZE];

// function prototypes follow
GLfloat XX( int );
GLfloat YY( int );
void myinit( void );
void surface( void );
float getLength(float, float);
void getVector(float, float, float *, float *);
void calcHeat(float);
void calcEven(float);
void calcRainbow(float);
void calcGray(float);
void calcMulti(float);
void display( void );
void reshape( int ,int );
void keyboard(unsigned char, int, int );

// 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)
{
    int  i, j;
    float x, y, u, v;

#ifndef GRAY
        glClearColor( 0.5, 0.5, 0.5, 1.0 );
#endif
#ifdef GRAY
		glClearColor( 0.0, 0.0, 1.0, 1.0 );
#endif
        glEnable(GL_DEPTH_TEST); // allow z-buffer display
        glShadeModel(GL_FLAT);

//  Calculate the points of the surface on the grid and log min and max
    for ( i=0; i<XSIZE; i++ )
       for ( j=0; j<YSIZE; j++ )
       {
          x = 0.5*(XX(i)+XX(i+1));
          y = 0.5*(YY(j)+YY(j+1));;
          getVector(x, y, &u, &v);
          vectors[i][j][0] = u; vectors[i][j][1] = v;
          length[i][j] = getLength(u,v);
          if (length[i][j] > YMAX) YMAX = length[i][j];
          if (length[i][j] < YMIN) YMIN = length[i][j];
       }
}

float getLength(float a, float b)
{
	return (sqrt(a*a+b*b));
}

//	Function that can either represent differential equations or complex functions
void getVector(float x, float y, float *u, float *v)
{
	*u = y*y - 1.0;	// Buchanan et al., Visualization in Teaching & Learning Math, p. 140
	*v = x*x - 1.0;	// x' = y^2-1, y' = x^2-1; value is vector <x',y'>
//	*u = x*x*x - 3.0*x*y*y +12.0*x + 12.0;	// Braden, op cit, p. 192
//	*v = 3.0*x*x*y - y*y*y +12.0*y;	// w = z^3+12z+12; value is complex number w
	return;
}

void surface(void)
{
    int i, j;
    float yavg, len, x, y;

	YRANGE = YMAX - YMIN;
//	draw the vector and the surface
	for ( i=0; i<XSIZE-1; i++ )
		for ( j=0; j<YSIZE-1; j++ )
		{
//	draw the direction vector
			if ((i%10 == 5) && (j%10 == 5)) {
			glColor4f(0.0, 1.0, 1.0, 1.0);
			x = 0.5*(XX(i)+XX(i+1));
			y = 0.5*(YY(j)+YY(j+1));
			glBegin(GL_LINE_STRIP);
				glVertex3f(x,y,0.0);
				glVertex3f(x,y,EPSILON);
				len = 5.0 * sqrt(vectors[i][j][0]*vectors[i][j][0]+vectors[i][j][1]*vectors[i][j][1]);
				glVertex3f(x+vectors[i][j][0]/len,y+vectors[i][j][1]/len, EPSILON);
			glEnd();
			}
// first triangle in the quad
			glBegin(GL_POLYGON);
				yavg = (length[i][j]+length[i+1][j]+length[i+1][j+1])/3.0;
#ifdef  HEAT
				calcHeat((yavg-YMIN)/YRANGE);
#endif
#ifdef  EVEN
				calcEven((yavg-YMIN)/YRANGE);
#endif
#ifdef  RAINBOW
				calcRainbow((yavg-YMIN)/YRANGE);
#endif
#ifdef  GRAY
				calcGray((yavg-YMIN)/YRANGE);
#endif
#ifdef  MULTI
				calcMulti((yavg-YMIN)/YRANGE);
#endif
				glColor3f(myColor[0],myColor[1],myColor[2]);
            	glVertex3f(XX(i),YY(j), 0.0);	// colors in the plane, not a surface
            	glVertex3f(XX(i+1),YY(j), 0.0);
            	glVertex3f(XX(i+1),YY(j+1), 0.0);
          glEnd();

        // second triangle in the quad  
          glBegin(GL_POLYGON);
			 yavg = (length[i][j]+length[i][j+1]+length[i+1][j+1])/3.0;
#ifdef  HEAT
			 calcHeat((yavg-YMIN)/YRANGE);
#endif
#ifdef  EVEN
			 calcEven((yavg-YMIN)/YRANGE);
#endif
#ifdef  RAINBOW
			 calcRainbow((yavg-YMIN)/YRANGE);
#endif
#ifdef  GRAY
			 calcGray((yavg-YMIN)/YRANGE);
#endif
#ifdef  MULTI
			 calcMulti((yavg-YMIN)/YRANGE);
#endif
			 glColor3f(myColor[0],myColor[1],myColor[2]);
             glVertex3f(XX(i),YY(j), 0.0);	// colors in the plane, not a surface
             glVertex3f(XX(i+1),YY(j+1), 0.0);
             glVertex3f(XX(i),YY(j+1), 0.0);
          glEnd();
       }
}

#ifdef HEAT
void calcHeat(float yval)
{
	if (yval < 0.30) { myColor[0]=yval/0.3; myColor[1]=0.0; myColor[2] = 0.0; return; }
	if ((yval >= 0.30) && (yval < 0.89))
		{ myColor[0]=1.0; myColor[1]=(yval-0.3)/0.59; myColor[2] = 0.0; return; }
	if (yval >= 0.89) { myColor[0]=1.0; myColor[1]=1.0; myColor[2]=(yval-0.89)/0.11; }
	return;
}
#endif

#ifdef EVEN
void calcEven(float yval)
{
	if (yval < 0.33) { myColor[0]=yval/0.3; myColor[1]=0.0; myColor[2] = 0.0; return; }
	if ((yval >= 0.33) && (yval < 0.66))
		{ myColor[0]=1.0; myColor[1]=(yval-0.33)/0.33; myColor[2] = 0.0; return; }
	if (yval >= 0.66) { myColor[0]=1.0; myColor[1]=1.0; myColor[2]=(yval-0.66)/0.33; }
	return;
}
#endif

#ifdef RAINBOW
void calcRainbow(float yval)
{
	if (yval < 0.2)	// purple to blue ramp
		{ myColor[0]=0.5*(1.0-yval/0.2);myColor[1]=0.0;myColor[2]=0.5+(0.5*yval/0.2);return;}
	if ((yval >= 0.2) && (yval < 0.40))	// blue to cyan ramp
		{ myColor[0]= 0.0; myColor[1]= (yval-0.2)*5.0; myColor[2] = 1.0; return; }
	if ((yval >= 0.40) && (yval < 0.6))	// cyan to green ramp
		{ myColor[0]= 0.0; myColor[1]= 1.0; myColor[2] = (0.6-yval)*5.0; return; }
	if ((yval >= 0.6) && (yval < 0.8))	// green to yellow ramp
		{ myColor[0]= (yval-0.6)*5.0; myColor[1]= 1.0; myColor[2] = 0.0; return; }
	if (yval >= 0.8)	// yellow to red ramp
		{ myColor[0]= 1.0; myColor[1]= (1.0-yval)*5.0; myColor[2]= 0.0; }
	return;
}
#endif

#ifdef GRAY
void calcGray(float yval)
{
	myColor[0] = myColor[1] = myColor[2] = yval;
	return;
}
#endif

#ifdef MULTI
void calcMulti(float yval)
{
	if (yval < 0.2)	// purple to blue ramp, starting with black
		{ myColor[0]=5.0*yval*0.5*(1.0-yval/0.2);
		  myColor[1]=0.0;myColor[2]=5.0*yval*(0.5+(0.5*yval/0.2));return;}
	if ((yval >= 0.2) && (yval < 0.40))	// blue to cyan ramp, starting with black
		{ myColor[0]= 0.0; myColor[1]= 5.0*(yval-0.2)*(yval-0.2)*5.0;
		  myColor[2] = 5.0*yval*1.0; return; }
	if ((yval >= 0.40) && (yval < 0.6))	// cyan to green ramp, starting with black
		{ myColor[0]= 0.0; myColor[1]= 5.0*(yval-0.4)*1.0;
		  myColor[2] = 5.0*(yval-0.4)*(0.6-yval)*5.0; return; }
	if ((yval >= 0.6) && (yval < 0.8))	// green to yellow ramp, starting with black
		{ myColor[0]= 5.0*(yval-0.6)*(yval-0.6)*5.0; 
		  myColor[1]= 5.0*(yval-0.6)*1.0; myColor[2] = 0.0; return; }
	if (yval >= 0.8)	// yellow to red ramp, starting with black
		{ myColor[0]= 5.0*(yval-0.8)*1.0;
		  myColor[1]= 5.0*(yval-0.8)*(1.0-yval)*5.0; myColor[2]= 0.0; }
	return;
}
#endif

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();
        gluLookAt( 0.0, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.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; q = positive, s = negative
		case 's' :
		case 'z' :	// rotate around Z; z = positive, x = negative
		case 'x' :
			ch = key; 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("Complex function");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);

	myinit();
	glutMainLoop();
}