/*
   Mathematical surface, with the ability to move the surface around and
   see all parts.
*/
#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef GLfloat point3[3];

#define SURFACE
#undef  HEAT
#undef  EVEN
#undef  RAINBOW

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

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 200  // needs to be much larger than this for best results
#define ZSIZE XSIZE
#define MINX   0.0
#define MAXX   5.0
#define MINZ   0.0
#define MAXZ   5.0

// define the charge and location of the point charges
int ncharge = 3;

#define MAXCHARGES 3
typedef struct { float charge, xpos, ypos; } charges;

charges Q[MAXCHARGES] = { {  5.0, 3.0, 1.0},
						  { -5.0, 1.0, 4.0},
						  {-10.0 ,2.0, 2.0} };

static float vertices[XSIZE][ZSIZE];

// function prototypes follow
GLfloat XX( int );
GLfloat ZZ( int );
void myinit( void );
void surface( void );
float coulSurf(float,float);
void calcHeat(float);
void calcEven(float);
void calcRainbow(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 ZZ(int j) {
     return (MINZ+((MAXZ-MINZ)/(float)(ZSIZE-1))*(float)(j));
     }

void myinit(void)
{
    int  i, j;
    float x, z;

#ifdef SURFACE
        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
        
        GLfloat mat_specular[]  ={ 0.8, 0.8, 0.8, 1.0 };
        long ii = 1;

        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, &ii ); // two-sided lighting
		
        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_NORMALIZE);  // make normal vectors 1-unit long after transform
#endif

        glClearColor( 0.5, 0.5, 0.5, 1.0 );
        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<ZSIZE; j++ )
       {
          x = XX(i);
          z = ZZ(j);
          vertices[i][j] = coulSurf(x,z);
#ifndef SURFACE
          if (vertices[i][j] > YMAX) YMAX = vertices[i][j];
          if (vertices[i][j] < YMIN) YMIN = vertices[i][j];
#endif
       }
}

float coulSurf(float x, float y)
{
	float retVal, dist;
	int i;
	
	retVal = 0.0;
	for (i=0; i<MAXCHARGES; i++) {
		dist = sqrt ( (x-Q[i].xpos)*(x-Q[i].xpos) + (y-Q[i].ypos)*(y-Q[i].ypos) + 0.1);
		retVal += Q[i].charge/dist;
	}
	retVal = retVal / 6.0;	// scale in vertical direction
	return (retVal);
}

void surface(void)
{
#ifdef SURFACE
    point3 vec1, vec2, triNormal;
#endif
    int i, j;
    float yavg;

#ifdef SURFACE
//  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 };

    glMaterialfv(GL_FRONT, GL_DIFFUSE, white );
    glMaterialfv(GL_BACK, GL_DIFFUSE, yellow );
    glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );
#endif

#ifndef SURFACE
	YRANGE = YMAX - YMIN;
#endif
//	draw the surface
    for ( i=0; i<XSIZE-1; i++ )
       for ( j=0; j<ZSIZE-1; j++ )
       {  // first triangle in the quad
          glBegin(GL_POLYGON);
#ifdef SURFACE
             vec1[0] = XX(i+1)-XX(i);
             vec1[1] = vertices[i+1][j]-vertices[i][j];
             vec1[2] = ZZ(j)-ZZ(j);
             vec2[0] = XX(i+1)-XX(i+1);
             vec2[1] = vertices[i+1][j+1]-vertices[i+1][j];
             vec2[2] = ZZ(j+1)-ZZ(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...
#endif
#ifdef  HEAT
			 yavg = (vertices[i][j]+vertices[i+1][j]+vertices[i+1][j+1])/3.0;
			 calcHeat((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
#ifdef  EVEN
			 yavg = (vertices[i][j]+vertices[i+1][j]+vertices[i+1][j+1])/3.0;
			 calcEven((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
#ifdef  RAINBOW
			 yavg = (vertices[i][j]+vertices[i+1][j]+vertices[i+1][j+1])/3.0;
			 calcRainbow((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
             glVertex3f(XX(i),vertices[i  ][j  ],ZZ(j));
             glVertex3f(XX(i+1),vertices[i+1][j  ],ZZ(j));
             glVertex3f(XX(i+1),vertices[i+1][j+1],ZZ(j+1));
          glEnd();

        // second triangle in the quad  
          glBegin(GL_POLYGON);
#ifdef SURFACE
             vec1[0] = XX(i+1)-ZZ(i);
             vec1[1] = vertices[i+1][j+1]-vertices[i][j];
             vec1[2] = ZZ(j+1)-ZZ(j);
             vec2[0] = XX(i)-XX(i+1);
             vec2[1] = vertices[i][j+1]-vertices[i+1][j+1];
             vec2[2] = ZZ(j+1)-ZZ(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);
#endif
#ifdef  HEAT
			 yavg = (vertices[i][j]+vertices[i][j+1]+vertices[i+1][j+1])/3.0;
			 calcHeat((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
#ifdef  EVEN
			 yavg = (vertices[i][j]+vertices[i][j+1]+vertices[i+1][j+1])/3.0;
			 calcEven((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
#ifdef  RAINBOW
			 yavg = (vertices[i][j]+vertices[i][j+1]+vertices[i+1][j+1])/3.0;
			 calcRainbow((yavg-YMIN)/YRANGE);
			 glColor3f(myColor[0],myColor[1],myColor[2]);
#endif
             glVertex3f(XX(i),vertices[i  ][j  ],ZZ(j));
             glVertex3f(XX(i+1),vertices[i+1][j+1],ZZ(j+1));
             glVertex3f(XX(i),vertices[i  ][j+1],ZZ(j+1));
          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)
//	commented out line creates a green "bathtub ring" in the surface
{
//	if ((yval > .49) && (yval < .51 )) {myColor[0] = myColor[2] = 0.0;  
//			myColor[1] = 1.0; return; }
	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)	// question: red at both ends?  Seems wrong...
{
	if (yval < 0.167)	// red to magenta ramp
		{ myColor[0]= 1.0 ; myColor[1]= 0.0 ; myColor[2] = yval * 6.0 ; return; }
	if ((yval >= 0.167) && (yval < 0.33))	// magenta to blue ramp
		{ myColor[0]= (0.33-yval)*6.0 ; myColor[1]= 0.0 ; myColor[2] = 1.0 ; return; }
	if ((yval >= 0.33) && (yval < 0.50))	// blue to cyan ramp
		{ myColor[0]= 0.0; myColor[1]= (yval-0.33)*6.0; myColor[2] = 1.0; return; }
	if ((yval >= 0.50) && (yval < 0.667))	// cyan to green ramp
		{ myColor[0]= 0.0; myColor[1]= 1.0; myColor[2] = (0.667-yval)*6.0; return; }
	if ((yval >= 0.667) && (yval < 0.833))	// green to yellow ramp
		{ myColor[0]= (yval-0.667)*6.0; myColor[1]= 1.0; myColor[2] = 0.0; return; }
	if (yval >= 0.833)	// yellow to red ramp
		{ myColor[0]= 1.0; myColor[1]= (1.0-yval)*6.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();
        /*           eye point    center of view       up   */
        gluLookAt( 15.0, 5.0, 0.0, 2.5, -2.0, 2.5, 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
           ch = key; break;
        case 'w' :
           ch = key; break;
        case 'a' :    // rotate around Y; q = positive, s = negative
           ch = key; break;
        case 's' :
           ch = key; break;
        case 'z' :    // rotate around Z; z = positive, x = negative
           ch = key; break;
        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("Coulombic surface");
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutKeyboardFunc(keyboard); // enable keyboard callback

        myinit();
        glutMainLoop();
}