/*
   Mathematical surface defined by three parametric functions of two variables:
       X(u,v)
       Y(u,v)
       Z(u,v)
   
   Example of a mathematical visualization for the introductory computer graphics
   course, adapted to allow the creation of an STL file to create physical hardcopy
   of the surface if desired.  This has been done via the LOM and Z-Corp systems.
   
   Ref:  Thomas Banchoff et al., Student-Generated Software for Differential
         Geometry, in Visualization in Teaching and Learning Mathematics,
         Walter Zimmermann and Steve Cunningham, eds., MAA Notes No. 19,
         Mathematical Association of America, 1991
*/

#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 the parameters of the surface grid
#define MAXALONG  100
//	one more than number of sides in torus cross-section so 4 -> triangular section
#define MAXAROUND 4
//	maximum number of pieces in a strip
#define MAXSTRIP 25

//	torus
#define RNGRAD 4.0	// radius of torus ring
#define TUBRAD 2.0	// radius of torus tube
#define X(u,v) (RNGRAD+TUBRAD*cos(1.3333333*u+v+t))*cos(u)	// includes t-animation
#define Y(u,v) (RNGRAD+TUBRAD*cos(1.3333333*u+v+t))*sin(u)	// that works for SCREEN
#define Z(u,v) TUBRAD*sin(1.3333333*u+v+t)
#define AXISLENGTH 3.0

float umin=-3.14159, umax=3.14159,
	  vmin=-3.14159, vmax=3.14159;
int usteps=MAXALONG;
int vsteps=MAXAROUND;

typedef GLfloat point3[3];

static char ch; // character from the keyboard that controls rotations
GLfloat ep = 25.0;
float t = 0.0;
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

typedef struct surfpoint {
	GLfloat x, y, z;
	} surfpoint;
struct surfpoint surface[MAXALONG][MAXAROUND], theStrip[MAXSTRIP][2];

FILE	*fp;

// function prototypes follow
void myinit( void );
void doSurface( void );
void doQuad(int, int, surfpoint, surfpoint, surfpoint, surfpoint);
void doStrip(int, surfpoint, surfpoint, surfpoint, surfpoint);
void emit(surfpoint *);
void display( void );
void reshape( int ,int );
void rotate(void);
void keyboard(unsigned char, int, int );
void animate(void);
GLfloat UU(int);
GLfloat VV(int);

// two functions to return u and v values for array indices  i  and  j  respectively
GLfloat UU(int i) {
     return (umin+((umax-umin)/(float)(usteps-1))*(float)(i));
     }
GLfloat VV(int j) {
     return (vmin+((vmax-vmin)/(float)(vsteps-1))*(float)(j));
     }

// initialize OpenGL system and global variables
void myinit(void)
{
        GLfloat light_pos0[]={  0.0, 20.0, 20.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[]={ 10.0, 20.0,-14.52, 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[]={-10.0, 20.0,-14.52, 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 };
        long i = 1;

        glClearColor( 0.0, 0.0, 1.0, 0.0 );
        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);
		glEnable(GL_LIGHT0);
		glEnable(GL_LIGHT1);
		glEnable(GL_LIGHT2);
		glEnable(GL_DEPTH_TEST);
		glEnable(GL_NORMALIZE);
//		glShadeModel(GL_SMOOTH);
        
//	Create a call list to display the surface
		glNewList(1, GL_COMPILE);
			doSurface();
		glEndList();
}

void doSurface(void)
{
    int  i, j;
    float u, v;

//  Set up the surface color
    GLfloat yellow[]= {1.0, 1.0, 0.0, 1.0};
    GLfloat mat_shininess[]={ 30.0 };

//  Calculate the points of the surface boundaries on the grid
	for ( i=0; i<usteps; i++ )
		for ( j=0; j<vsteps; j++ )
		{
			u = UU(i);
			v = VV(j);
			surface[i][j].x = X(u,v);
			surface[i][j].y = Y(u,v);
			surface[i][j].z = Z(u,v);
		}

//  actually draw the surface
    for ( i=0; i<usteps-1; i++ )	// along torus -- larger
       for ( j=0; j<vsteps-1; j++ )	// around torus -- smaller
       {	// front face
                glMaterialfv(GL_FRONT, GL_AMBIENT, yellow );
                glMaterialfv(GL_FRONT, GL_DIFFUSE, yellow );
                glMaterialfv(GL_BACK, GL_AMBIENT, yellow );
                glMaterialfv(GL_BACK, GL_DIFFUSE, yellow );
                glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );

                // new approach -- define separate quad function to deal with the geometry
                // sequence (i,j), (i+1,j), (i,j+1), (i+1,j+1)
                doQuad(20, 1, surface[i][j], surface[i+1][j], surface[i][j+1], surface[i+1][j+1]);
       }
}

/*
	Function to take a quad whose points are given as in the following figure:
	
		P0 -------m------- P1
		 |                 |
		 |                 |
               n |                 |
		 |                 |
		 |                 |
		P2 --------------- P3
		
	and divide the quad into  m  strips vertically, then pass each of those
	strips to another function that images the strip.  Eventually the single-
	strip piece will also emit STL information to get the 4-3 torus done right
	on the LOM machine.
*/
void doQuad(int n, int m, surfpoint p0, surfpoint p1, surfpoint p2, surfpoint p3)
{
	int i;
//	float ;
	surfpoint A, B, C, D;	//  A and B are the top points of each separate
							//	strip; C and D are the bottom points.

	for (i=0; i<m; i++) {
		A.x = (p0.x*(float)(m-i)   + p1.x*(float)i)/(float)m;
		A.y = (p0.y*(float)(m-i)   + p1.y*(float)i)/(float)m;
		A.z = (p0.z*(float)(m-i)   + p1.z*(float)i)/(float)m;
		B.x = (p0.x*(float)(m-i-1) + p1.x*(float)(i+1))/(float)m;
		B.y = (p0.y*(float)(m-i-1) + p1.y*(float)(i+1))/(float)m;
		B.z = (p0.z*(float)(m-i-1) + p1.z*(float)(i+1))/(float)m;

		C.x = (p2.x*(float)(m-i)   + p3.x*(float)i)/(float)m;
		C.y = (p2.y*(float)(m-i)   + p3.y*(float)i)/(float)m;
		C.z = (p2.z*(float)(m-i)   + p3.z*(float)i)/(float)m;
		D.x = (p2.x*(float)(m-i-1) + p3.x*(float)(i+1))/(float)m;
		D.y = (p2.y*(float)(m-i-1) + p3.y*(float)(i+1))/(float)m;
		D.z = (p2.z*(float)(m-i-1) + p3.z*(float)(i+1))/(float)m;
		doStrip(n, A, B, C, D);
	}
}

/*
	Now we have one vertical strip, and we must subdivide that into n pieces,
	each of which will be two triangles.  We actually create a nx2 array of
	surfpoints, and then divide each quad going down the array into two
	triangles, each of which is "emitted" with its own function call.  That
	function call is where we decide whether to draw the triangle or to write
	the triangle out in an STL file.
		P0 --1-- P1
		 |       |
		 |       |
               n |       |
		 |       |
		 |       |
		P2 --1-- P3
*/

void doStrip(int n, surfpoint p0, surfpoint p1, surfpoint p2, surfpoint p3)
{
	int i, j;
	surfpoint A, B, buffer[3];
	
	for (i=0; i<=n; i++) {
		A.x = (p0.x*(float)(n-i) + p2.x*(float)i)/(float)n;
		A.y = (p0.y*(float)(n-i) + p2.y*(float)i)/(float)n;
		A.z = (p0.z*(float)(n-i) + p2.z*(float)i)/(float)n;
		B.x = (p1.x*(float)(n-i) + p3.x*(float)i)/(float)n;
		B.y = (p1.y*(float)(n-i) + p3.y*(float)i)/(float)n;
		B.z = (p1.z*(float)(n-i) + p3.z*(float)i)/(float)n;
		theStrip[i][0] = A;
		theStrip[i][1] = B;
	}
	//	now manipulate the strip to send out the triangles one at a time to the
	//	actual output function.  Note that we use a rolling buffer to manage the
	//	triangle sequence.
	buffer[0] = theStrip[0][0];
	buffer[1] = theStrip[0][1];
	for (i=1; i<=n; i++)
		for (j=0; j<2; j++) {
			buffer[2] = theStrip[i][j];
			emit(buffer);
			buffer[0] = buffer[1];
			buffer[1] = buffer[2];
		}
}

//	Handle one triangle, sent as an array of three surfpoints.  First calculate the
//	normal to the triangle, then either call OpenGL output or write the points to the
//	STL file.
void emit( surfpoint *buffer ) {
	surfpoint Normal, diff1, diff2;

	diff1.x = buffer[1].x - buffer[0].x;
	diff1.y = buffer[1].y - buffer[0].y;
	diff1.z = buffer[1].z - buffer[0].z;
	diff2.x = buffer[2].x - buffer[1].x;
	diff2.y = buffer[2].y - buffer[1].y;
	diff2.z = buffer[2].z - buffer[1].z;
	Normal.x = diff1.y*diff2.z - diff2.y*diff1.z;
	Normal.y = diff1.z*diff2.x - diff1.x*diff2.z;
	Normal.z = diff1.x*diff2.y - diff1.y*diff2.x;

	glBegin(GL_POLYGON);
		glNormal3f(Normal.x,Normal.y,Normal.z);
		glVertex3f(buffer[0].x,buffer[0].y,buffer[0].z);
		glVertex3f(buffer[1].x,buffer[1].y,buffer[1].z);
		glVertex3f(buffer[2].x,buffer[2].y,buffer[2].z);
	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);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    //          eye point   center of view      up
    gluLookAt(0.0, 0.0, ep, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
	//	save the viewing projection for later use
	glGetFloatv( GL_MODELVIEW_MATRIX, viewProj );
    //	do the rotation and draw the surface by calling the display list
	rotate();
	doSurface();
//    glCallList(1);	// then the 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);
}

//	Function to perform rotation around coordinate axes, not around object axes
//	rotation control uses two keypads embedded in standard keyboard
//    q  w                     o  p  -- rotate around x-axis + -
//     a  s                   k  l   -- rotate around y-axis + -
//      z  x                 m  ,    -- rotate around z-axis + -
void rotate(void )
{
	#define ANGLE 2.0	// angle in degrees
	glLoadIdentity();
	if (ch != ' ') {	// keyboard indicated a rotation...
	//	Put identity onto modelview stack and start with the initial rotation
		switch(ch) {
			case 'q':
			case 'o':
				glRotatef( ANGLE, 1.0, 0.0, 0.0); break;
			case 'w':
			case 'p':
				glRotatef(-ANGLE, 1.0, 0.0, 0.0); break;
        	case 'a':
        	case 'k':
            	glRotatef( ANGLE, 0.0, 1.0, 0.0); break;
        	case 's':
        	case 'l':
            	glRotatef(-ANGLE, 0.0, 1.0, 0.0); break;
        	case 'z':
        	case 'm':
            	glRotatef( ANGLE, 0.0, 0.0, 1.0); break;
        	case 'x':
        	case ',':
            	glRotatef(-ANGLE, 0.0, 0.0, 1.0); break;
	    }
	}
    //  NOW 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.  Rebuild the overall modelview matrix by starting with the
    //  viewing transformation, saved from the reshape() function, and then
    //	multiplying by the modeling transformation created in the first two
    //	statements.
    	glMultMatrixf( saveState );
    	glGetFloatv( GL_MODELVIEW_MATRIX, saveState );
    	glLoadIdentity();
    	glMultMatrixf( viewProj );
    	glMultMatrixf( saveState );
}

void animate(void)
{
    t = t + 0.01;
    if (t > 2*3.14159) t-= 2*3.14159;
    ch = ' ';
    glutPostRedisplay();
}

void keyboard(unsigned char key, int x, int y)
{
    ch = ' ';
    switch (key) {
        case 'q' :	// rotate around X; q,o = positive, w,p = negative
        case 'o' :
        case 'w' :
        case 'p' :
        case 'a' :	// rotate around Y; a,k = positive, s,l = negative
        case 'k' :
        case 's' :
        case 'l' :
        case 'z' :	// rotate around Z; z,m = positive, x,, = negative
        case 'm' :
        case 'x' :
        case ',' :
           ch = key; break; // save the key
        case 'f' :	// move forward 'f' or back 'b'
        	ep -= 0.1; break;
        case 'b' :
        	ep += 0.1; break;
    }
    glutPostRedisplay(); // perform display again
}

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("Parametric Surface");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard); // enable keyboard callback
	glutIdleFunc(animate);

	myinit();
	glutMainLoop();
}