/*
   Mathematical surface that is everywhere continuous and nowhere differentiable.

   (c) February 2001, Steve Cunningham & Ken Brown
*/

#define OSX

#ifdef MAC
#include "glut.h"
#endif
#ifdef OSX
#include <GLUT/glut.h>
#endif
#ifdef PC
#include <GL/glut.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

typedef GLfloat point3[3];

static char ch; // character from 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

//  Define the parameters of the surface grid
//  For an N x M grid you need XSIZE = N+1 and YSIZE = M+1
#define XSIZE 250  // needs to be 2^N + 1 for some n
#define YSIZE XSIZE
#define MINX   -1.0
#define MAXX    1.0
#define MINY   -1.0
#define MAXY    1.0

// the global array to hold the heights of the points on the surface
static GLfloat vertices[XSIZE][YSIZE];
float ep = 5.0;
int ITER=1;

// function prototypes follow
GLfloat XX( int );
GLfloat YY( int );
void myinit( void );
void surface( void );
float blancmange(float,float,int);
float zigzag(float,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));
     }

// initialize OpenGL and all necessary variables
void myinit(void)
{
    	int i, j;
    	float x, y;

        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
        
//      set up overall light data
        GLfloat mat_specular[]  ={ 0.8, 0.8, 0.8, 1.0 };

        glClearColor( 0.0, 0.0, 1.0, 0.0 ); // blue background
        glShadeModel(GL_SMOOTH); // use gourand shading
        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 );
        
        glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1 ); // two-sided lighting
		
// attributes
        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_DEPTH_TEST); // allow z-buffer display
        glEnable(GL_NORMALIZE);  // make normal vectors 1-unit long after transform

//  Calculate the points of the surface on the grid
    for ( i=0; i<XSIZE; i++ )
       for ( j=0; j<YSIZE; j++ )
       {
          x = XX(i); y = YY(j);
          vertices[i][j] = blancmange(x,y,ITER);
       }
}

// given a 2D location and a number of iterations, this function
// returns the correct height on that point on the blancmange surface
float blancmange(float x, float y, int n)
{
	float retVal, multiplier;
	int i;
	
	retVal = 0.0;
	multiplier = 0.5;
	for (i=0; i<n; i++) {
		multiplier = multiplier * 2.0;
		retVal = retVal + zigzag(x*multiplier,y*multiplier)/multiplier;
	}
	return (retVal);
}

float zigzag(float x, float y)
{
	float smallx, smally;
	int intx, inty;
	
	smallx = fabs(x); intx = (int)smallx;
	smally = fabs(y); inty = (int)smally;
	smallx = smallx - intx;				// move x and y between 0 and 1
	smally = smally - inty;
	smallx = smallx * (1.0 - smallx);	// maximum if value was 1/2
	smally = smally * (1.0 - smally);	// minimum if value was 0 or 1
	return (4.0*smallx*smally);			// scale to height 1
}

// draw the blancmange surface
void surface(void)
{
    point3 vec1, vec2, triNormal;
    int i, j;
    
//  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 };

/*  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, yellow );
		glMaterialfv(GL_BACK, GL_DIFFUSE, white );
		glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess );
		glBegin(GL_POLYGON);
			// calculate the surface normal for correct lighting
			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...

			// now draw the triangle
			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, yellow );
          glMaterialfv(GL_BACK, GL_DIFFUSE, white );
          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();
       }
}

// callback to handle the displaying of the viewing window
void display( void )
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	//	set up view
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    //          eye point   center of view       up
    gluLookAt( ep, 0.0, ep, 0.0, 0.0, 0.0, 0.0, 0.0, 1.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...

    // 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 = ' ';	// remove the control value of the character
    
    //  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,(float)w/(float)h,1.0,100.0);
}

void keyboard(unsigned char key, int xx, int yy)
{
		int i, j;
		float x, y;

		ch = ' ';
        switch (key)
        {
        case 'q' :	//	rotation controls
        case 'w' :
        case 'a' :
        case 's' :
        case 'z' :
        case 'x' :
           ch = key; break;
        case 'i' :
            ITER++;
            //  Calculate the points of the surface on the grid
            for ( i=0; i<XSIZE; i++ )
                for ( j=0; j<YSIZE; j++ ) {
                    x = XX(i); y = YY(j);
                    vertices[i][j] = blancmange(x,y,ITER);
                }
            break;
        case 'f' :
            ep = ep - 0.2; break;
        case 'b' :
            ep = ep + 0.2; break;
        }
        glutPostRedisplay();
}

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("Blancmange function surface");
        glutDisplayFunc(display);
        glutReshapeFunc(reshape);
        glutKeyboardFunc(keyboard);

        myinit();
        glutMainLoop();

		return 0;
}