/*
	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.
	However, the STL file is not perfect because of roundoffs at the "seams" in the
	parametric space; some detailed work is needed to ensure that these seams have
	exactly equal coordinates.
	
	The surface is animated but does not include any interaction.
   
	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

	Source file to be used with
	Cunningham, Computer Graphics: Programming in OpenGL for Visual Communication, Prentice-Hall, 2007

	Intended for class use only
*/

#include <GLUT/glut.h>
#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
#define ANGLE 2.0	// angle in degrees

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

GLfloat ep = 25.0;
float t = 0.0;
float umin=-3.14159, umax=3.14159,
	  vmin=-3.14159, vmax=3.14159;
int usteps=MAXALONG;
int vsteps=MAXAROUND;

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 animate(void);
GLfloat UU(int);
GLfloat VV(int);

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

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

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 ii = 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, &ii ); // two-sided lighting

//	attributes
		glEnable(GL_LIGHTING);
		glEnable(GL_LIGHT0);
		glEnable(GL_LIGHT1);
		glEnable(GL_LIGHT2);
		glEnable(GL_DEPTH_TEST);
		glEnable(GL_NORMALIZE);
}

void doSurface(void)
{
    int  i, j;

//  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 );
			// To avoid problems with long thin quads, we define a 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 )
{
    int  i, j;
    float u, v;

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

void animate(void)
{
    t = t + 0.01;
    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(800, 800);
	glutInitWindowPosition(0, 0);
	glutCreateWindow("Parametric Surface");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutIdleFunc(animate);

	myinit();
	glutMainLoop();
}