/*
 *  Programmer:   Yu Zhang
 *  Class:        CAP 5705 Computer Graphics (Spring 2000)
 *  Assignment:   Final project
 *  Date:         April 26, 2000
 *
 *  Title:        Spatial Rational B-Spline Motion
 * USE arrow keys and F1-F8 
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <GL/glut.h>

#define WSIZE 500     /* initial window size */
#define GSIZE 6       /* size of gripper */
#define STEP 1        /* point move step */
#define TIME 0.1      /* time vary step */
#define ASTEP 0.2     /* angle vary step */
#define SPSTEP 2      /* speed step  */
#define DEMOT 12      /* total demo number */

#define G (GLfloat)1.61803398875  /* golden ratio */
#define S (GLfloat) 0.61803398875

#define POINTNUM 50   /* Maximum number of points */
#define BUFSIZE 512
#define PI 3.14159265359
#define NAU 0.000001     /* small float for vertex closeness */
#define INFINITY 999999  /* representation of infinity */

#define TRUE 1
#define FALSE 0
#define MAX 10000 
#define ORDER 3 // degree

/* Prototypes */
GLsizei ww = WSIZE, wh = WSIZE; /* initial window size */
GLfloat vsize = 250.0;  /* half side length of square */

GLfloat ctrlpoints[POINTNUM][3];  /* cube control positions */
GLfloat ctrlaxis[POINTNUM][3];  /* rotate axis */
GLfloat ctrlangle[POINTNUM];  /* rotate angles */

GLfloat intpoints[MAX][3];  /* interpolate points positions */
GLfloat intaxis[MAX][3];  /* interpolate rotate axis */
GLfloat intangle[MAX];  /* interpolate rotate angles */

GLfloat red[] = {1.0, 0.0, 0.0, 1.0};      /* red */
GLfloat green[] = {0.0, 1.0, 0.0, 1.0};    /* green */
GLfloat blue[] = {0.0, 0.0, 1.0, 1.0};     /* blue */
GLfloat magenta[] = {1.0, 0.0, 1.0, 1.0};  /* magenta */
GLfloat black[] = {0.0, 0.0, 0.0, 1.0};    /* black */
GLfloat cyan[] = {0.0, 1.0, 1.0, 1.0};     /* cyan */
GLfloat yellow[] = {1.0, 1.0, 0.0, 1.0};   /* yellow */

GLfloat angleX = 0.0, angleZ = 0.0;
GLfloat xx = 0.0, yy = 0.0, zz = 0.0;

GLfloat atX = 0.0;  /* centering the screen at the center of the coordinate system */
GLfloat atY = 0.0;
GLfloat atZ = 0.0;

GLfloat upX = 0.0;  /* general up direction is */
GLfloat upY = 1.0;
GLfloat upZ = 0.0;

GLfloat cr = 20.0;
GLfloat tupY, ttheta;
GLfloat cam_alpha = 0.3, cam_theta = 0.4;
GLfloat cyl_x = 0.0, cyl_z = 0.0;

int method = 0, surface = 0, ctype = 0, cnum = 5;
int animate = 0, first = 0, selectObj = 0;
int control = 0, curve = 0;
int demoMode = 0;

int ObjNum = 0;   /* picking globals */

int Curx = 0, Cury = 0;
int lastX = 0, lastY = 0;

int totalNum = 0;  
int disNum = 0;  
int interval = 3, tmp_interval = 3;

void drawObjs(GLenum mode);



void normalize(float v[3]) 
{    
    GLfloat d=1.0;

    d = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); 

    if (d == 0.0) {
        printf("\n Zero length vector. \n");    
        return;
    }

    v[0] = v[0]/d; 
    v[1] = v[1]/d; 
    v[2] = v[2]/d; 
}


void normcrossprod(float v1[3], float v2[3], float out[3]) 
{ 
    out[0] = v1[1]*v2[2] - v1[2]*v2[1]; 
    out[1] = v1[2]*v2[0] - v1[0]*v2[2]; 
    out[2] = v1[0]*v2[1] - v1[1]*v2[0]; 
    normalize(out); 
}


void cube(float ratio)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[8][3] = {
        {-1.0, -1.0, -1.0}, {-1.0, -1.0, 1.0}, {-1.0, 1.0, -1.0}, 
        {-1.0, 1.0, 1.0}, {1.0, -1.0, -1.0}, {1.0, -1.0, 1.0}, 
        {1.0, 1.0, -1.0}, {1.0, 1.0, 1.0} };

    GLint tindex[6][4] = { 
        {0, 1, 3, 2}, {4, 6, 7, 5}, {0, 4, 5, 1}, 
        {3, 7, 6, 2}, {0, 2, 6, 4}, {1, 5, 7, 3} };

    for (i=0; i<8; i++)
        for (j=0; j<3; j++)
            vertex[i][j] = vertex[i][j]*ratio;

    for (i=0; i<6; i++) { 
        switch (i) {
            case 0:  glColor3f(1.0, 0.0, 0.0);  /* red */
                     break;

            case 1:  glColor3f(0.0, 1.0, 0.0);  /* green */
                     break;

            case 2:  glColor3f(0.0, 0.0, 1.0);  /* blue */
                     break;

            case 3:  glColor3f(1.0, 0.0, 1.0);  /* magenta */
                     break;

            case 4:  glColor3f(0.0, 1.0, 1.0);  /* cyan */
                     break;

            case 5:  glColor3f(1.0, 1.0, 0.0);  /* yellow */
                     break;
        }

        glBegin(GL_POLYGON);  
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][3]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);

            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
            glVertex3fv(&vertex[tindex[i][3]][0]);    
        glEnd(); 
    }
}


/*  draw cylinder */
void draw_cyl(GLfloat radius, GLfloat height)
{
    GLUquadricObj *cylinder;

    cylinder = gluNewQuadric();
    gluQuadricOrientation(cylinder, GLU_OUTSIDE);
    gluQuadricDrawStyle(cylinder, GLU_FILL);
    gluQuadricNormals(cylinder, GLU_SMOOTH);

    gluCylinder(cylinder, radius, radius, height, 20, 1);
}


void draw_cover(GLfloat radius)
{
    GLUquadricObj *qobj;

    qobj = gluNewQuadric();

    gluQuadricOrientation(qobj, GLU_OUTSIDE);
    gluQuadricDrawStyle(qobj, GLU_FILL);
    gluQuadricNormals(qobj, GLU_SMOOTH);

    gluDisk(qobj, 0.0, radius, 20, 4);
}


void tube(float radius, float height)
{
    glPushMatrix();
        glTranslatef(0.0, 0.0, 0.0);

        glRotatef(-90.0, 1.0, 0.0, 0.0);

        draw_cyl(radius, height);  // draw cylinder C2

        glPushMatrix();
            glRotatef(180.0, 1.0, 0.0, 0.0);
            draw_cover(radius);     // draw bottom cover
        glPopMatrix();

        glTranslatef(0.0, 0.0, height);
        draw_cover(radius);     // draw top cover
    glPopMatrix();
}


void cuboid(float length, float width, float height)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[8][3] = {
        {-1.0, -1.0, -1.0}, {-1.0, -1.0, 1.0}, {-1.0, 1.0, -1.0}, 
        {-1.0, 1.0, 1.0}, {1.0, -1.0, -1.0}, {1.0, -1.0, 1.0}, 
        {1.0, 1.0, -1.0}, {1.0, 1.0, 1.0} };

    GLint tindex[6][4] = { 
        {0, 1, 3, 2}, {4, 6, 7, 5}, {0, 4, 5, 1}, 
        {3, 7, 6, 2}, {0, 2, 6, 4}, {1, 5, 7, 3} };

    for (i=0; i<8; i++) {
         vertex[i][0] = vertex[i][0]*length;
         vertex[i][1] = vertex[i][1]*width;
         vertex[i][2] = vertex[i][2]*height;
    }

    for (i=0; i<6; i++) { 
        glBegin(GL_POLYGON);  
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][3]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);

            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
            glVertex3fv(&vertex[tindex[i][3]][0]);    
        glEnd(); 
    }
}


void gripper(float ratio)
{
    float t_radius, t_height;
    float c_length, c_width, c_height;

    t_radius = ratio;
    t_height = ratio*3;

    c_length = ratio*2.5;
    c_width  = ratio*1.2;
    c_height = ratio/2;

    glColor3fv(cyan);    // draw in cyan 
//  glColor3fv(magenta);  // draw in magenta 

    glPushMatrix();
        tube(t_radius, t_height);

        glTranslatef(0.0, t_height, 0.0);
        glRotatef(-90.0, 1.0, 0.0, 0.0);
        cuboid(c_length, c_width, c_height);

        glPushMatrix();
            glTranslatef(0.0, 0.0, c_width);
            glTranslatef(c_length-c_height, 0.0, 0.0);
            glRotatef(90.0, 0.0, 1.0, 0.0);
            cuboid(c_length/2, c_width, c_height);

            glRotatef(-90.0, 1.0, 0.0, 0.0);
            glTranslatef(-c_width/4, 0.0, 0.0);
            tube(t_radius*0.8, t_height/4);
        glPopMatrix();

        glTranslatef(0.0, 0.0, c_width);
        glTranslatef(-c_length+c_height, 0.0, 0.0);
        glRotatef(90.0, 0.0, 1.0, 0.0);
        cuboid(c_length/2, c_width, c_height);

        glRotatef(90.0, 1.0, 0.0, 0.0);
        glTranslatef(-c_width/4, 0.0, 0.0);
        tube(t_radius*0.8, t_height/4);
    glPopMatrix();
}


void Tetrahedron(float ratio)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[4][3] = {
        {1.0, 1.0, 1.0}, {1.0, -1.0, -1.0},
        {-1.0, 1.0, -1.0},{-1.0, -1.0, 1.0} };

    GLint tindex[4][3] = { 
        {0, 1, 2}, {0, 3, 1}, {0, 2, 3}, {1, 3, 2} };

    for (i=0; i<4; i++)
        for (j=0; j<3; j++)
            vertex[i][j] = vertex[i][j]*ratio;  

    for (i=0; i<4; i++) { 
        glBegin(GL_POLYGON); 
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][2]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);
  
            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
        glEnd(); 
    }
}


void Octahedron(float ratio)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[6][3] = {
        {1.0, 0.0, 0.0}, {-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, 
        {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0},{0.0, 0.0, -1.0} };

    GLint tindex[8][3] = { 
        {4, 1, 3}, {4, 3, 0}, {4, 0, 2}, {4, 2, 1},
        {5, 0, 3}, {5, 3, 1}, {5, 1, 2}, {5, 2, 0} };

    for (i=0; i<6; i++)
        for (j=0; j<3; j++)
            vertex[i][j] = vertex[i][j]*ratio;  

    for (i=0; i<8; i++) {    
        glBegin(GL_POLYGON);  
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][2]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);

            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
        glEnd(); 
    }
}


void Icosahedron(float ratio)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[12][3] = {
        {-1.0, 0.0, G}, {1.0, 0.0, G}, {-1.0, 0.0, -G}, {1.0, 0.0, -G},
        {0.0, G, 1.0}, {0.0, G, -1.0}, {0.0, -G, 1.0}, {0.0, -G, -1.0}, 
        {G, 1.0, 0.0}, {-G, 1.0, 0.0}, {G, -1.0, 0.0}, {-G, -1.0, 0.0} };

    GLint tindex[20][3] = { 
        {1, 4, 0}, {4, 9, 0}, {4, 5, 9}, {8, 5, 4}, {1, 8, 4},
        {1, 10, 8}, {10, 3, 8}, {8, 3, 5}, {3, 2, 5}, {3, 7, 2},    
        {3, 10, 7}, {10, 6, 7}, {6, 11, 7}, {6, 0, 11}, {6, 1, 0}, 
        {10, 1, 6}, {11, 0, 9}, {2, 11, 9}, {5, 2, 9}, {11, 2, 7} };

    for (i=0; i<12; i++)
        for (j=0; j<3; j++)
            vertex[i][j] = vertex[i][j]*ratio;  

    for (i=0; i<20; i++) {    
        glBegin(GL_POLYGON);  
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][2]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);

            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
        glEnd(); 
    }
}


void Dodecahedron(float ratio)
{
    int i, j;
    GLfloat d1[3], d2[3], norm[3];

    GLfloat vertex[20][3] = {
        {1.0, 1.0, 1.0}, {1.0, 1.0, -1.0}, {1.0, -1.0, 1.0},
        {1.0, -1.0, -1.0}, {-1.0, 1.0, 1.0}, {-1.0, 1.0, -1.0},
        {-1.0, -1.0, 1.0}, {-1.0, -1.0, -1.0}, 
        {S, G, 0.0}, {-S, G, 0.0}, {S, -G, 0.0}, {-S, -G, 0.0},
        {G, 0.0, S}, {G, 0.0, -S}, {-G, 0.0, S}, {-G, 0.0, -S},
        {0.0, S, G}, {0.0, -S, G}, {0.0, S, -G}, {0.0, -S, -G} };

    GLint tindex[12][5] = { 
        {1, 8, 0, 12, 13}, {4, 9, 5, 15, 14}, {2, 10, 3, 13, 12}, 
        {7, 11, 6, 14, 15}, {2, 12, 0, 16, 17}, {1, 13, 3, 19, 18}, 
        {4, 14, 6, 17, 16}, {7, 15, 5, 18, 19}, {4, 16, 0, 8, 9}, 
        {2, 17, 6, 11, 10}, {1, 18, 5, 9, 8}, {7, 19, 3, 10, 11} };

    for (i=0; i<20; i++)
        for (j=0; j<3; j++)
            vertex[i][j] = vertex[i][j]*ratio;  

    for (i=0; i<12; i++) {    
        glBegin(GL_POLYGON);  
            for (j=0; j<3; j++) {    
                d1[j] = vertex[tindex[i][1]][j] - vertex[tindex[i][0]][j];    
                d2[j] = vertex[tindex[i][2]][j] - vertex[tindex[i][0]][j];    
            }
            normcrossprod(d1, d2, norm); 
            glNormal3fv(norm);

            glVertex3fv(&vertex[tindex[i][0]][0]); 
            glVertex3fv(&vertex[tindex[i][1]][0]); 
            glVertex3fv(&vertex[tindex[i][2]][0]);    
            glVertex3fv(&vertex[tindex[i][3]][0]);    
            glVertex3fv(&vertex[tindex[i][4]][0]);    
        glEnd(); 
    }
}


/*  select the object  */
void myObjects(int myselect, float ratio)
{
    switch (myselect) {
        case 0:                    /* draw a cube  */
                cube(ratio);                            
                break;
        case 1:                    /* draw a gripper */
                gripper(ratio/2.0+1.0); 
                break;
        case 2:                    /* draw a tube */
                glColor3f(0.2, 0.7, 0.8);
                tube(ratio, ratio*3);
                break;
        case 3:                    /* draw a tetrahedron */
                glColor3f(0.3, 0.1, 0.5);
                Tetrahedron(ratio); 
                break;
        case 4:                    /* draw a cuboid */
                glColor3f(0.3, 0.1, 0.5);
                cuboid(ratio*2, ratio, ratio/2); 
                break;

        case 5:                    /* draw an octahedron */
                glColor3f(0.4, 0.1, 0.9);
                Octahedron(ratio*2);
                break;
        case 6:                    /* draw a icosahedron */
                glColor3f(0.5, 0.9, 0.1);
                Icosahedron(ratio);
                break;
        case 7:                    /* draw a dodecahedron */
                glColor3f(1.0, 0.6, 0.5);
                Dodecahedron(ratio);
                break;
        case 8:                    /* draw a Utah teapot */
                glColor3f(0.8, 0.1, 0.1);
                glutSolidTeapot(ratio*1.5);
                break;
    }
}


void processHits(GLint hits, GLuint buffer[])
{
    int i, j, names;
    GLuint *ptr;

    printf("hits = %d\n", hits);
    ptr = (GLuint *) buffer;
    for (i=0; i<hits; i++) {  /* for each hit  */
        names = *ptr;
        printf(" number of names for hit = %d\n", names); ptr++;
        printf(" z1 is %g;", (float) *ptr/0x7fffffff); ptr++;
        printf(" z2 is %g\n", (float) *ptr/0x7fffffff); ptr++;
        printf("   the name is ");
        for (j = 0; j < names; j++) {  /* for each name */
            printf("%d ", *ptr); 
            ptr++;
        }
        ptr--;
        ObjNum = (int) *ptr;
        ptr++;
        printf("\n");
    }

    if (hits <= 0)
        ObjNum = 0;
}


void pick(int button, int state, int x, int y)
{
    GLuint selectBuf[BUFSIZE];
    GLint hits;
    GLint viewport[4];
    //control = 0;

    if (button==GLUT_RIGHT_BUTTON && state==GLUT_DOWN) {
        angleX = -45.0;
        angleZ = 45.0;

        if (control == 0) {
            glGetIntegerv(GL_VIEWPORT, viewport);

            glSelectBuffer(BUFSIZE, selectBuf);
            (void) glRenderMode(GL_SELECT);

            glInitNames();
            glPushName(0);
 
            glMatrixMode(GL_PROJECTION);
            glPushMatrix();
                glLoadIdentity();

                /*  create 5x5 pixel picking region near cursor location */
                gluPickMatrix((GLdouble) x, (GLdouble) (viewport[3] - y), 
                              5.0, 5.0, viewport);
                                           /* set the viewing volume */
                glOrtho(-vsize, vsize, -vsize, vsize, -vsize, vsize);  

                glMatrixMode(GL_MODELVIEW);

                glPushMatrix();
                    drawObjs(GL_SELECT);
                glPopMatrix();

                glMatrixMode(GL_PROJECTION);
            glPopMatrix();

            hits = glRenderMode(GL_RENDER);

            processHits(hits, selectBuf);
 
            glMatrixMode(GL_MODELVIEW);
            glLoadIdentity();
            glFlush();
        }
        else {
            ObjNum = 0;
          //angleX = 0.0;
          //angleZ = 0.0;
        }
    }

    control = 1;
    if (ObjNum > 10 || ObjNum < 1) {
        control = 0;
    }
    else if (button==GLUT_LEFT_BUTTON && state==GLUT_DOWN) {
        lastX = x;
        lastY = y;
    }

    glutPostRedisplay();
}


void mouseMotion(int x, int y)
{
    if (ObjNum <= 10 && ObjNum > 0 && control==1) {
        Curx = lastX - x;
        angleZ = angleZ + Curx;
        if (angleZ > 360) {
            angleZ -= 360;
        }
        if (angleZ < -360) {
            angleZ += 360;
        }

        Cury = lastY - y;          
        angleX = angleX - Cury;
        if (angleX > 360) {
            angleX -= 360;
        }
        if (angleX < -360) {
            angleX += 360;
        }
        lastX = x;
        lastY = y;
    }
    else 
        return;
   
    glutPostRedisplay();
}




/*  Find knots for a given set of 3D points.  */
void yield_knots(float pts_x[], float pts_y[], float pts_z[], 
                 int k, float knot[])
{
    int i, k1, k2;
    float delta;

    k1 = k-1; 
    k2 = k+2;
    knot[0] = 0.0;   /* initialization */

    delta = sqrt((pts_x[2]-pts_x[0])*(pts_x[2]-pts_x[0]) +
                 (pts_y[2]-pts_y[0])*(pts_y[2]-pts_y[0]) +
                 (pts_z[2]-pts_z[0])*(pts_z[2]-pts_z[0]));

    knot[1]= sqrt(delta); 

    for (i=2; i<k; i++) {
        delta = sqrt((pts_x[i+1] - pts_x[i])*(pts_x[i+1] - pts_x[i]) +
                     (pts_y[i+1] - pts_y[i])*(pts_y[i+1] - pts_y[i]) +
                     (pts_z[i+1] - pts_z[i])*(pts_z[i+1] - pts_z[i]));

        knot[i] = knot[i-1] + sqrt(delta); 
    }

    delta = sqrt((pts_x[k2] - pts_x[k])*(pts_x[k2] - pts_x[k]) +
                 (pts_y[k2] - pts_y[k])*(pts_y[k2] - pts_y[k]) +
                 (pts_z[k2] - pts_z[k])*(pts_z[k2] - pts_z[k]));

    knot[k]= knot[k-1] + sqrt(delta);  
}


/*  set up the linear system for B-spline interpolation
    alpha,beta,gamma: 1-D arrays that constitute the elements of 
    the interpolation matrix.
*/
void setupsys(float knot[], int k, float alpha[],
              float beta[], float gamma[])
{
    int i, k1;
    float delta_im2, delta_im1, delta_i, delta_ip1, sum;

    k1 = k - 1;

    alpha[0] = 0.0; 
    beta[0]  = 1.0; 
    gamma[0] = 0.0; 

    if (k == 1) {  // special case 
        alpha[1] = 0.0;
        beta[1] = 1.0;
        gamma[1] = 0.0;
        return;
    }

    if (k == 2) {  // special case 
        delta_im1 = knot[1] - knot[0];
        delta_i = knot[2] - knot[1];
        sum = delta_im1 + delta_i;

        alpha[1] = delta_i*delta_i/sum;
        beta[1]  = (delta_i*delta_im1)/sum + delta_im1*delta_i/ sum;
        gamma[1] = delta_im1*delta_im1/sum;

        alpha[1] = alpha[1]/sum;
        beta[1]  = beta[1]/sum;
        gamma[1] = gamma[1]/sum;

        beta[2]  = 1.0;
        alpha[2] = 0.0;
        gamma[2] = 0.0;
        return;
    }

    /* do rest when k>2 */

    delta_im1 = knot[1] - knot[0];
    delta_i   = knot[2] - knot[1];
    delta_ip1 = knot[3] - knot[2];
    sum = delta_im1 + delta_i;

    alpha[1] = delta_i*delta_i/sum;
    beta[1]  = (delta_i*delta_im1)/sum + 
                delta_im1*(delta_i + delta_ip1)/(sum + delta_ip1);
    gamma[1] = delta_im1*delta_im1/(sum + delta_ip1);

    alpha[1] = alpha[1]/sum;
    beta[1]  = beta[1]/sum;
    gamma[1] = gamma[1]/sum;

    /* Now for the main loop: */
    for (i=2; i<k1; i++) {
        /* compute delta_i_minus_2,...  */
        delta_im2 = knot[i-1] - knot[i-2];
        delta_im1 = knot[i]   - knot[i-1];
        delta_i   = knot[i+1] - knot[i];
        delta_ip1 = knot[i+2] - knot[i+1];
        
        sum = delta_im1 + delta_i;

        alpha[i] = delta_i*delta_i/(delta_im2 + sum);
        beta[i]  = delta_i*(delta_im2 + delta_im1)/(delta_im2 + sum) +
                   delta_im1*(delta_i + delta_ip1)/(sum + delta_ip1);
        gamma[i] = delta_im1*delta_im1/(sum + delta_ip1);
                
        alpha[i] = alpha[i]/sum;
        beta[i]  = beta[i]/sum;
        gamma[i] = gamma[i]/sum;
    }

    /*  special care at the end:  */
    delta_im2 = knot[k-2] - knot[k-3];
    delta_im1 = knot[k1] - knot[k-2];
    delta_i   = knot[k] - knot[k1];
    sum = delta_im1 + delta_i;

    alpha[k1] = delta_i*delta_i/(delta_im2 + sum);
    beta[k1]  = delta_i*(delta_im2 + delta_im1)/(delta_im2 + sum) +
                delta_im1*delta_i/sum;
    gamma[k1] = delta_im1*delta_im1/sum;

    alpha[k1] = alpha[k1]/sum;
    beta[k1]  = beta[k1]/sum;
    gamma[k1] = gamma[k1]/sum;

    alpha[k]=0.0; 
    beta[k]=1.0; 
    gamma[k]=0.0;
}


/*  perform LU decomposition of tridiagonal system with
    lower diagonal alpha, diagonal beta, upper diagonal gamma.

    low: L-matrix entries
    up:  U-matrix entries
*/
void lu_decompse(float alpha[], float beta[], float gamma[], 
                 int k, float up[], float low[])
{
    int i;

    up[0]  = beta[0];
    low[0] = 0;

    for (i=1; i<=k; i++) {
        low[i] = alpha[i]/up[i-1];
        up[i]  = beta[i] - low[i]*gamma[i-1];
    }
}


/*  solve  tridiagonal linear system of size (k+1)(k+1) whose 
    LU decompostion has entries up and low, and whose right 
    hand side is rhs, and whose original matrix had upper 
    diagonal gamma. Solution is d[0],...,d[k+2];

    k:    size of system: k+1 eqs in k+1 unknowns.
    rhs:  right hand side
    d:  solution vector.
*/
void solve(float up[], float low[], float gamma[], int k,
           float rhs[], float d[])
{
    int i;
    float aux[500];

    d[0] = rhs[0];
    d[1] = rhs[1];

    /* forward substitution:  */
    aux[0] = rhs[1];

    for (i=1; i<=k; i++) 
        aux[i] = rhs[i+1] - low[i]*aux[i-1];
        
    /* backward substitution:  */
    d[k+1] = aux[k]/up[k];

    for (i=k-1; i>0; i--) 
        d[i+1] = (aux[i] - gamma[i]*d[i+2])/up[i];

    d[k+2] = rhs[k+2];
}


/*  Computes B-spline points pts[1] and pts[k+1] according to 
    Bessel end condition.
*/
void bessel_ends(float pts[], float knot[], int k)
{
    float alpha, beta; 
        
    if (k == 1) {  // when only have one interval
        pts[1] = (2.0*pts[0] + pts[3])/3.0;
        pts[2] = (2.0*pts[3] + pts[0])/3.0;
    } 
    else if (k == 2) {
        /* beginning: */
        alpha = (knot[2] - knot[1])/(knot[2] - knot[0]);
        beta  = 1.0 - alpha;

        pts[1] = (pts[2] - alpha*alpha*pts[0] - 
                  beta*beta*pts[4])/(2.0*alpha*beta);

        pts[1] = 2.0*(alpha*pts[0] + beta*pts[1])/3.0 + pts[0]/3.0;

        /* end: */
        alpha = (knot[2] - knot[1])/(knot[2] - knot[0]);
        beta  = 1.0 - alpha;

        pts[3] = (pts[2] - alpha*alpha*pts[0] - 
                  beta*beta*pts[4])/(2.0*alpha*beta);

        pts[3] = 2.0*(alpha*pts[3] + beta*pts[4])/3.0 + pts[4]/3.0;     
    }
    else {
        /* beginning: */
        alpha = (knot[2] - knot[1])/(knot[2] - knot[0]);
        beta  = 1.0 - alpha;

        pts[1] = (pts[2] - alpha*alpha*pts[0] - 
                  beta*beta*pts[3])/(2.0*alpha*beta);

        pts[1] = 2.0*(alpha*pts[0] + beta*pts[1])/3.0 + pts[0]/3.0;

        /* end: */
        alpha = (knot[k] - knot[k-1])/(knot[k] - knot[k-2]);
        beta  = 1.0 - alpha;

        pts[k+1] = (pts[k] - alpha*alpha*pts[k-1] - 
                    beta*beta*pts[k+2])/(2.0*alpha*beta);

        pts[k+1] = 2.0*(alpha*pts[k+1] + beta*pts[k+2])/3.0 + pts[k+2]/3.0;
    }
}


/*  uses de Boor algorithm to compute one coordinate on B-spline 
    curve for param. value u in interval i. 

        degree: polynomial degree of each piece of curve
        coeff:  B-spline control points
        knot:   knot sequence
        u:      evaluation abscissa
        i:      u's interval: u[i]<= u < u[i+1]
*/
float deboor(int degree, float coeff[], float knot[], float u, int i)
{
    int k, j;
    float t1, t2;
    float coeffa[300];  /* might need adjustment! */
        
    for (j=i-degree+1; j<=i+1; j++)
        coeffa[j]=coeff[j];

    for (k=1; k<=degree; k++) {
        for (j=i+1; j>=i-degree+k+1; j--){
            t1 = (knot[j+degree-k] - u)/(knot[j+degree-k] - knot[j-1]);
            t2 = 1.0 - t1;

            coeffa[j] = t1*coeffa[j-1] + t2*coeffa[j];
        }       
    }

    return (coeffa[i+1]);
}


/*  generates points on B-spline curve (one coordinate).
        degree: polynomial degree of each piece of curve
        k:      number of intervals
        coeff:  B-spline control points
        knot:   knot sequence: knot[0]...knot[k+2*degree-2]
        dense:  how many points per segment
        points: output array.
        point_num: how many points are generated.
*/
void sppoints(int degree, int k, float coeff_x[],       float knot[], 
              int dense, float points_x[], int *point_num)
{
    int i, ii, k2;
    float u;

    k2 = k + degree - 1;

    *point_num=0;
    for (i=degree-1; i<k+degree-1; i++) {
        //printf("\n i= %d", i); 

        if (knot[i+1] > knot[i]) {

            for (ii=0; ii<=dense; ii++) {       

                u = knot[i] + ii*(knot[i+1] - knot[i])/dense;
                //printf("\n u= %d", u); 

                points_x[*point_num] = deboor(degree, coeff_x, knot, u, i);

                *point_num= (*point_num) + 1;
            }   
        }
    }
}


/* Finds the spline interpolant to the points in 
   pts_x, pts_y, pts_z.
*/
void int_spline(float knot[], int k, float pts_x[], float pts_y[], 
                float pts_z[], float spline_x[], float spline_y[], 
                float spline_z[], int tt)
{
    float alpha[100], beta[100], gamma[100], up[100], low[100];
    float tmp_x[100], tmp_y[100], tmp_z[100];
    int i;

    for (i=0; i<=k; i++) {
        tmp_x[i] = pts_x[i];
        tmp_y[i] = pts_y[i];
        tmp_z[i] = pts_z[i];
    }

    pts_x[0] = tmp_x[0];
    pts_y[0] = tmp_y[0];
    pts_z[0] = tmp_z[0];

    pts_x[1] = 0.0;
    pts_y[1] = 0.0;
    pts_z[1] = 0.0;

    for (i=2; i<=k; i++) {
        pts_x[i] = tmp_x[i-1];
        pts_y[i] = tmp_y[i-1];
        pts_z[i] = tmp_z[i-1];
    }

    pts_x[k+1] = 0.0;
    pts_y[k+1] = 0.0;
    pts_z[k+1] = 0.0;

    pts_x[k+2] = tmp_x[k];
    pts_y[k+2] = tmp_y[k];
    pts_z[k+2] = tmp_z[k];

    if (tt == 1)
        yield_knots(pts_x, pts_y, pts_z, k, knot);

    setupsys(knot, k, alpha, beta, gamma);

    lu_decompse(alpha, beta, gamma, k, up, low);

    bessel_ends(pts_x, knot, k);
    bessel_ends(pts_y, knot, k);
    bessel_ends(pts_z, knot, k);

    solve(up, low, gamma, k, pts_x, spline_x);
    solve(up, low, gamma, k, pts_y, spline_y);
    solve(up, low, gamma, k, pts_z, spline_z);
}


/* function to produce spline */
void spline(void)
{
    float knot[POINTNUM], knot1[POINTNUM];
    float tmp;
    //float knot[] = {0, 1, 2, 3, 4, 5};
    //float knot1[] = {0, 0, 0, 1, 2, 3, 4, 5, 5, 5};
    float pts_x[POINTNUM], pts_y[POINTNUM], pts_z[POINTNUM];
    float spline_x[POINTNUM], spline_y[POINTNUM], spline_z[POINTNUM];
    float pts_ax[POINTNUM], pts_ay[POINTNUM], pts_az[POINTNUM];
    float spline_ax[POINTNUM], spline_ay[POINTNUM], spline_az[POINTNUM];

    int i, k, k2, dense, degree, points_num;
    float points_x[MAX], points_y[MAX], points_z[MAX];
    float points_ax[MAX], points_ay[MAX], points_az[MAX];

    degree = 3;
    dense = interval;

    k  = cnum - 1;
    k2 = k + degree - 1;

    for (i=0; i<cnum; i++) {
         pts_x[i] = ctrlpoints[i][0];
         pts_y[i] = ctrlpoints[i][1];
         pts_z[i] = ctrlpoints[i][2];

         pts_ax[i] = ctrlaxis[i][0];
         pts_ay[i] = ctrlaxis[i][1];
         pts_az[i] = ctrlaxis[i][2];
    }

    int_spline(knot, k, pts_x, pts_y, pts_z, spline_x, spline_y, spline_z, 1);
    int_spline(knot, k, pts_ax, pts_ay, pts_az, spline_ax, spline_ay, spline_az, 2);

/*
    for (i=0; i<=k+2; i++) {
        printf("\n knot= %f", knot[i]);
        printf("\n x= %f,  y= %f,  z= %f", ctrlpoints[i][0], 
               ctrlpoints[i][1], ctrlpoints[i][2]);
      //printf("\n cx= %f,  cy= %f,  cz= %f", spline_x[i], 
      //       spline_y[i], spline_z[i]);               

        printf("\n ax= %f,  ay= %f,  az= %f", ctrlaxis[i][0], 
               ctrlaxis[i][1], ctrlaxis[i][2]);
      //printf("\n cax= %f,  cay= %f,  caz= %f\n", spline_ax[i], 
      //       spline_ay[i], spline_az[i]);             
    }
*/

    tmp = knot[0];
    for (i=0; i<=2; i++) {
        knot1[i] = tmp;
    }

    for (i=3; i<=k+1; i++) {
        knot1[i] = knot[i-2];
    }

    tmp = knot[k];
    for (i=k+2; i<=k+4; i++) {
        knot1[i] = tmp;
      //printf("\n knot= %f", knot1[i]);
    }

    sppoints(degree, k, spline_x, knot1, dense, points_x, &points_num);
    sppoints(degree, k, spline_y, knot1, dense, points_y, &points_num);
    sppoints(degree, k, spline_z, knot1, dense, points_z, &points_num);

    sppoints(degree, k, spline_ax, knot1, dense, points_ax, &points_num);
    sppoints(degree, k, spline_ay, knot1, dense, points_ay, &points_num);
    sppoints(degree, k, spline_az, knot1, dense, points_az, &points_num);

    totalNum = points_num;

    for (i=0; i<points_num; i++) {
         intpoints[i][0] = points_x[i];
         intpoints[i][1] = points_y[i];
         intpoints[i][2] = points_z[i];

         intaxis[i][0] = points_ax[i];
         intaxis[i][1] = points_ay[i];
         intaxis[i][2] = points_az[i];
    }

  //printf("\n points_num = %d", points_num);
/*  
    for (i=0; i<points_num; i++) {
        printf("\n px= %f,  py= %f,  pz= %f", intpoints[i][0], 
               intpoints[i][1], intpoints[i][2]);
        printf("\n pax= %f,  pay= %f,  paz= %f", intaxis[i][0], 
               intaxis[i][1], intaxis[i][2]);
    }
*/
}


/*  create the initial object positions  */
void create_pos(int pos_num, int pos_mode)
{
    int i;
    GLfloat demopoints[DEMOT][3] = { 
        {-180.0, 80.0, 90.0}, {-100.0, 50.0, 50.0}, {0.0, 0.0, 10.0}, 
        {100.0, 40.0, 0.0}, {180.0, 90.0, -30.0}, {150.0, 80.0, -80.0}, 
        {80.0, 50.0, -20.0}, {0.0, 30.0, -30.0}, {-30.0, 30.0, 10.0}, 
        {-10.0, -20.0, 10.0}, {-100.0, 0.0, 50.0}, {-50.0, -60.0, 80.0} };

    GLfloat demoaxis[DEMOT][3] = { 
        {20.0, 30.0, 20.0}, {10.0, 80.0, 60.0}, { 0.0, 43.0, 30.0}, 
        {60.0, 50.0, 40.0}, {90.0, 60.0, 60.0}, { 70.0, 83.0, 72.0}, 
        {80.0, 60.0, 60.0}, {70.0, 40.0, 30.0}, { 60.0, 43.0, 30.0}, 
        {20.0, 30.0, 20.0}, {40.0, 40.0, 30.0}, { 50.0, 3.0, 30.0} };

  //GLfloat demoangle[DEMOT];

    if (pos_mode==0) {
        for (i=0; i<pos_num; i++) {                                
            ctrlpoints[i][0] = demopoints[i][0];
            ctrlpoints[i][1] = demopoints[i][1];
            ctrlpoints[i][2] = demopoints[i][2];

            ctrlaxis[i][0] = demoaxis[i][0];
            ctrlaxis[i][1] = demoaxis[i][1];
            ctrlaxis[i][2] = demoaxis[i][2];

            //intangle[i] = demoangle[i];
        }
    }

    if (pos_mode == 1) {
        for (i=0; i<pos_num; i++) {
                            /* yield pos_num random control points */
            ctrlpoints[i][0] = (GLfloat)(rand()%4001)/10.0 - 200.0;
            ctrlpoints[i][1] = (GLfloat)(rand()%2001)/10.0 - 100.0;
            ctrlpoints[i][2] = (GLfloat)(rand()%2001)/10.0 - 100.0;

            ctrlaxis[i][0] = (GLfloat)(rand()%901)/10.0;
            ctrlaxis[i][1] = (GLfloat)(rand()%901)/10.0;
            ctrlaxis[i][2] = (GLfloat)(rand()%901)/10.0;

            //intangle[i] = (GLfloat)(rand()%1810)/10.0;
        }
    }

    if (pos_mode==2) {
        ctrlpoints[cnum-1][0] = demopoints[cnum-1][0];
        ctrlpoints[cnum-1][1] = demopoints[cnum-1][1];
        ctrlpoints[cnum-1][2] = demopoints[cnum-1][2];

        ctrlaxis[cnum-1][0] = demoaxis[cnum-1][0];
        ctrlaxis[cnum-1][1] = demoaxis[cnum-1][1];
        ctrlaxis[cnum-1][2] = demoaxis[cnum-1][2];

        //intangle[cnum-1] = demoangle[cnum-1];
    }

    if (pos_mode == 3) {
                            /* yield one more random control point */
        ctrlpoints[cnum-1][0] = (GLfloat)(rand()%4001)/10.0 - 200.0;
        ctrlpoints[cnum-1][1] = (GLfloat)(rand()%2001)/10.0 - 100.0;
        ctrlpoints[cnum-1][2] = (GLfloat)(rand()%2001)/10.0 - 100.0;

        ctrlaxis[cnum-1][0] = (GLfloat)(rand()%901)/10.0;
        ctrlaxis[cnum-1][1] = (GLfloat)(rand()%901)/10.0;
        ctrlaxis[cnum-1][2] = (GLfloat)(rand()%901)/10.0;

        //intangle[cnum-1] = (GLfloat)(rand()%1810)/10.0;
    }
}



/*  initialization  */
void init (void)  
{
    GLfloat diffuseMaterial[4] = {0.5, 0.1, 0.1, 1.0};
    GLfloat mat_ambient[] = {0.5, 0.1, 0.2, 1.0};
    GLfloat mat_specular[] = {0.9, 0.8, 0.5, 1.0};
    GLfloat mat_shininess[] = {20.0};

    GLfloat light0_position[] = {1.0, 1.0, 1.0, 0.0};
    GLfloat model_ambient[] = {0.5, 0.5, 0.5, 1.0};

    glClearColor(1.0, 1.0, 1.0, 1.0);  // white background 
    glShadeModel(GL_SMOOTH);
    //glShadeModel(GL_FLAT);

    glMaterialfv(GL_FRONT, GL_DIFFUSE, diffuseMaterial);
    glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
    glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
    glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess); 

    glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, model_ambient);
    //glLightModelf(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);

    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);

    glDepthFunc(GL_LEQUAL);
    glEnable(GL_DEPTH_TEST);  // attributes 

    glEnable(GL_AUTO_NORMAL);
    glEnable(GL_NORMALIZE);

    glColorMaterial(GL_FRONT, GL_DIFFUSE);
    glEnable(GL_COLOR_MATERIAL);

    glViewport(0, 0, ww, wh);  // set the viewport to the size of the window 
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity(); 
                                       // set the viewing volume
    glOrtho(-vsize, vsize, -vsize, vsize, -vsize, vsize);  

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    create_pos(cnum, 0);  /* yield init 5 random initial control points */
}


/*  reshape when window size is changed  */
void reshape(int w, int h)
{
    glViewport(0, 0, w, h);  /* set the viewport to the size of the window */
    ww = w;
    wh = h;
    glFlush();
}


/* increase frame */
void moveObj(void)
{
    //t = t + TIME;
    disNum++;
    glutPostRedisplay();        
}


/*  draw coordinate axes  */
void draw_coord(float len)
{
    glBegin(GL_LINES);
        glLineWidth(2.0);
        glColor3fv(red);     /* draw x axis in red */
        glVertex3f(-len, 0.0, 0.0);
        glVertex3f(2*len, 0.0, 0.0);
                
        glColor3fv(green);   /* draw y axis in green */
        glVertex3f(0.0, -len, 0.0);
        glVertex3f(0.0, 2*len, 0.0);
                
        glColor3fv(blue);    /* draw z axis in blue */
        glVertex3f(0.0, 0.0, -len);
        glVertex3f(0.0, 0.0, 2*len);
    glEnd();
}


/* demo mode */
void largeObj(int p)
{
    float ratio = 66.0;

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glPushMatrix();
        glRotatef(angleX, 1, 0, 0);
        glRotatef(angleZ, 0, 0, 1);

        myObjects(selectObj, ratio);
    glPopMatrix();     
}


/* select draw mode */
void drawObjs(GLenum mode)
{
    int i;

    if (control == 0) {
        for (i=0; i<cnum; i++) {
            if (mode == GL_SELECT) 
                glLoadName(i+1);

            glPushMatrix();
                glTranslatef(ctrlpoints[i][0], ctrlpoints[i][1], ctrlpoints[i][2]);

                glRotatef(ctrlaxis[i][0], 1, 0, 0);
                glRotatef(ctrlaxis[i][1], 0, 1, 0);
                glRotatef(ctrlaxis[i][2], 0, 0, 1);

              //glRotatef(intangle[i] , ctrlaxis[i][0], ctrlaxis[i][1], ctrlaxis[i][2]);
                myObjects(selectObj, 10.0);
            glPopMatrix();
        }
    }
    else 
        largeObj(ObjNum);
}


/*  display graphics  */
void display(void)
{
    int i;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);  /*clear the window */

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    gluLookAt(atX + cr*cos(cam_theta)*sin(cam_alpha),   /* center the camera at */
              atY + cr*sin(cam_theta),
              atZ + cr*cos(cam_theta)*cos(cam_alpha),
              atX, atY, atZ,     /* Focus */
              upX, upY, upZ);    /* Up */

    draw_coord(20.0);

    drawObjs(GL_RENDER);

    if (method == 1) {  // display intermidate objects
        spline();
        for (i=0; i<totalNum; i++) {
            glPushMatrix();
                glTranslatef(intpoints[i][0], intpoints[i][1], intpoints[i][2]);

                glRotatef(intaxis[i][0], 1, 0, 0);
                glRotatef(intaxis[i][1], 0, 1, 0);
                glRotatef(intaxis[i][2], 0, 0, 1);

                //glRotatef(intangle[i] , ctrlaxis[i][0], ctrlaxis[i][1], ctrlaxis[i][2]);

                myObjects(selectObj, 10.0);
            glPopMatrix();
        }
    }

    if (animate == 1) {  // animation
        spline();

        if (disNum < totalNum) {
            glPushMatrix();
                glTranslatef(intpoints[disNum][0], intpoints[disNum][1], intpoints[disNum][2]);

                glRotatef(intaxis[disNum][0], 1, 0, 0);
                glRotatef(intaxis[disNum][1], 0, 1, 0);
                glRotatef(intaxis[disNum][2], 0, 0, 1);

                //glRotatef(intangle[i] , ctrlaxis[i][0], ctrlaxis[i][1], ctrlaxis[i][2]);

                myObjects(selectObj, 10.0);
            glPopMatrix();
        }
        else {
            glutIdleFunc(NULL); 
            animate = 0;
        }
    }

    if (curve == 1) {  // draw the curve
        spline();

        glColor3fv(magenta);  // draw in magenta 

        glBegin(GL_LINE_STRIP);
            for (i=0; i<totalNum; i++) {

                glPushMatrix();
                    glTranslatef(intpoints[i][0], intpoints[i][1], intpoints[i][2]);

                    glVertex3fv(intpoints[i]);
                    //glRotatef(intangle[i] , ctrlaxis[i][0], ctrlaxis[i][1], ctrlaxis[i][2]);
                glPopMatrix();
            }
        glEnd();
    }

    glFlush();
    glutSwapBuffers();
}


/*  function deals with the general keys  */
void keyboard(unsigned char key, int x, int y)
{
    switch (key) {
        case 32:           /* when SPACE bar is pushed */
                if (demoMode==0 && cnum<DEMOT)
                    create_pos(cnum, 0);
                else
                    create_pos(cnum, 1);

                first = 0;

                glutPostRedisplay();
                break;

        case 'i': 
        case 'I':          /* when 'i' is pushed */
                cnum += STEP;

                if (cnum < POINTNUM) {
                    if (demoMode==0 && cnum<DEMOT)
                        create_pos(cnum, 2);
                    else
                        create_pos(cnum, 3);
                }
                else 
                    cnum = POINTNUM;

                glutPostRedisplay();
                break;

        case 'd': 
        case 'D':          /* when 'd' is pushed */
                cnum -= STEP;
                if (cnum <= 2) 
                    cnum = 2;

                glutPostRedisplay();
                break;

        case 'a': 
        case 'A':          /* when 'a' is pushed */
                interval -= SPSTEP;
                if (interval <= 0) 
                    interval = 0;

                glutPostRedisplay();
                break;

        case 's': 
        case 'S':          /* when 's' is pushed */
                interval += SPSTEP;
                if (interval >= BUFSIZE) 
                    interval = BUFSIZE;

                glutPostRedisplay();
                break;

        case 27:          /* when ESC is pushed */
                exit(0);
    }
}


/*  function deals with the special keys  */
void specialKey(int skey, int x, int y)
{
    switch (skey) {
        case GLUT_KEY_UP:              /* when UP arrow is pushed */
                if (control == 0) {
                    cam_theta = cam_theta + ASTEP;
                    if (cam_theta > 2*PI) 
                        cam_theta = cam_theta - 2*PI;
                    if (cam_theta>PI/2 && cam_theta<3*PI/2)
                        upY = -1.0;
                    else
                        upY = 1.0;
                }
                else
                    cyl_z -= STEP;

                glutPostRedisplay();
                break;

        case GLUT_KEY_DOWN:            /* when DOWN arrow is pushed */
                if (control == 0) {
                    cam_theta = cam_theta - ASTEP;
                    if (cam_theta < 0) 
                        cam_theta = cam_theta + 2*PI;
                    if (cam_theta > PI/2 && cam_theta < 3*PI/2)
                        upY = -1.0;
                    else
                        upY = 1.0;
                }
                else
                    cyl_z += STEP;

                glutPostRedisplay();
                break;

        case GLUT_KEY_LEFT:            /* when LEFT arrow is pushed */
                if (control == 0) {
                    cam_alpha = cam_alpha - ASTEP;
                }
                else
                    cyl_x -= STEP;

                glutPostRedisplay();
                break;

        case GLUT_KEY_RIGHT:           /* when RIGHT arrow is pushed */
                if (control == 0) {
                    cam_alpha = cam_alpha + ASTEP;
                }
                else
                    cyl_x += STEP;

                glutPostRedisplay();
                break;

        case GLUT_KEY_F1:              /* when 'F1' is pushed */
                method = !method;
                interval = 3;
                first = 0;
                surface = 0;
                curve = 0;

                glutPostRedisplay();
                break;

        case GLUT_KEY_F2:              /* when 'F2' is pushed */
                animate = 1;

              //tmp_interval = interval;
                method = 0;
                surface = 0;
                disNum = 0;
                if (first == 0) {
                   interval = 60;
                   first = 1;
                }
                if (interval <= 2) 
                    interval = 2;

                glutIdleFunc(moveObj);

                glutPostRedisplay();
                break;

        case GLUT_KEY_F3:              /* when 'F3' is pushed */
                surface = !surface;
                first = 0;
                if (surface == 1) {
                    tmp_interval = interval;
                    interval = 200;
                    method = 1;
                }
                else {
                    interval= tmp_interval;
                    method = 0;
                }

                glutPostRedisplay();
                break;

        case GLUT_KEY_F4:              /* when 'F4' is pushed */
                selectObj = !selectObj;

                glutPostRedisplay();
                break;

        case GLUT_KEY_F5:              /* when 'F5' is pushed */
                curve = !curve;
                surface = 0;

                if (curve == 1) {
                    tmp_interval = interval;
                    interval = 200;
                }
                else {
                    interval= tmp_interval;
                }

                glutPostRedisplay();
                break;

        case GLUT_KEY_F6:              /* when 'F6' is pushed */
                if (selectObj >= 8) {
                    selectObj = 0;
                }
                else {
                    selectObj = selectObj + 1;
                }

                glutPostRedisplay();
                break;

        case GLUT_KEY_F8:              /* when 'F8' is pushed */
                demoMode = !demoMode;

                glutPostRedisplay();
                break;
    }
}


int main(int argc, char** argv)
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
    glutInitWindowSize(ww, wh);        /* width x height pixel window */
    glutInitWindowPosition(0, 0);      /* place window top left on display */
    glutCreateWindow("Final Project -- Yu Zhang"); /* window title */
    init();                            /* set attributes */
    glutDisplayFunc(display);          /* display function */
    glutMouseFunc(pick);               /* mouse function */ 
    glutMotionFunc(mouseMotion);
    glutKeyboardFunc(keyboard);        /* general keys */
    glutSpecialFunc(specialKey);       /* special keys  */
    glutReshapeFunc(reshape);          /* resize the window */
    glutMainLoop();
    return 0;
}