// Concordia University Graphics Library
// Author:     Peter Grogono
// Thanks to:  Qingzhe Huang, Yiping Le, Mark Kilgard, Ken Shoemake, F.S. Hill, and others.
// Date:       06 June 2003
// Revised:    11 November 2005

// Commenting conventions:
// 1. Specifications given in the header file are not repeated here.
// 2. Comments are given only if the purpose of the code is not obvious.
// 3. Explanations of complex algorithms (e.g., quaternion to matrix)
//    are given in accompanying documentation rather than here.

#include "cugl.h"

#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>

using namespace std;

namespace cugl {

const char *ErrorStrings[] = 
{
   "No error",                                        // NO_ERROR
   "Invalid coordinate index",                        // BAD_INDEX
   "Invalid matrix mode",                             // BAD_MATRIX_MODE
   "Matrix is not a rotation transformation",         // BAD_ROTATION_MATRIX
   "Matrix is singular",                              // SINGULAR_MATRIX
   "Division by zero",                                // ZERO_DIVISOR
   "Normalizing null vector",                         // ZERO_NORM
   "Invalid interpolator argument",                   // BAD_INTERPOLATOR_ARG
   "Failed to open file",                             // OPEN_FAILED
   "File is not in BMP format",                       // NOT_BMP_FILE
   "File is not in 24-bit BMP format",                // NOT_24_BITS
   "File is in compressed BMP format",                // COMPRESSED_BMP_FILE
   "Memory allocation failed",                        // NOT_ENOUGH_MEMORY
   "Pixmap has not been initialized",                 // NO_PIX_MAP
   "Parameters do not define a plane",                // BAD_PLANE
   "Line lies in plane",                              // BAD_LINE
   "Too many materials",                              // TOO_MANY_MATERIALS 
   "Too many points"                                  // TOO_MANY_POINTS 
};

CUGLErrorType cuglError = NO_ERROR;

CUGLErrorType getError()
{
   CUGLErrorType currentValue = cuglError;
   cuglError = NO_ERROR;
   return currentValue;
}

const char *getErrorString(CUGLErrorType err)
{
   return ErrorStrings[err];
}

void checkOpenGLStatus()
{
   GLenum glError = glGetError();
   if (glError != GL_NO_ERROR)
      cerr << 
         "OpenGL error " << glError << ": " << 
         gluErrorString(glError) << endl;
}

void checkCUGLStatus()
{
   CUGLErrorType error = getError();
   if (error != NO_ERROR)
      cerr << 
         "CUGL error " << error << ": " << 
         getErrorString(error) << endl;
}

inline double sqr(double x)
{
   return x * x;
}

// Constructors for class Point.

Point::Point(GLfloat x, GLfloat y, GLfloat z, GLfloat w)
   : x(x), y(y), z(z), w(w) {}


Point::Point (const Quaternion & q)
{
   Vector v = q.vector();
   x = v[0];
   y = v[1];
   z = v[2];
   w = q.scalar();
}

// Member functions for class Point.

// Return a reference to a component of point: i = 0,1,2,3.
GLfloat & Point::operator[](int i)
{
   switch(i)
   {
   default:
      cuglError = BAD_INDEX;
      // return x;
   case 0:
      return x;
   case 1:
      return y;
   case 2:
      return z;
   case 3:
      return w;
   }
}
   
// Return a const reference to a component of point: i = 0,1,2,3.
const GLfloat & Point::operator[](int i) const
{
   switch(i)
   {
   default:
      cuglError = BAD_INDEX;
      // return x;
   case 0:
      return x;
   case 1:
      return y;
   case 2:
      return z;
   case 3:
      return w;
   }
}
   
void Point::normalize()
{
    if (w == 0)
        cuglError = ZERO_DIVISOR;
    else
    {
        x = x / w;
        y = y / w;
        z = z / w;
        w = 1;
    }
}

Point Point::unit() const
{
    if (w == 0)
    {
        cuglError = ZERO_DIVISOR;
        return *this;
    }
    else
        return Point(x/w, y/w, z/w, 1);
}

// Friend functions for class Point

std::ostream & operator<<(std::ostream & os, const Point & p)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   os << setw(0) << '(' << 
      setw(width) << p[0] << ',' << 
      setw(width) << p[1] << ',' << 
      setw(width) << p[2] << ',' << 
      setw(width) << p[3] << ')';
   os.flush();
   return os;
}

// Member functions for class Line

void Line::draw() const
{
   glBegin(GL_LINES);
      s.draw();
      f.draw();
   glEnd();
}

// Friend functions for class Line

Point meet(Line & k, Plane & p)
{
   GLfloat den = p.a*(k.f.x-k.s.x) + p.b*(k.f.y-k.s.y) + p.c*(k.f.z-k.s.z) + p.d*(k.f.w-k.s.w);
   if (den == 0)
   {
      cuglError = BAD_LINE;
      return Point();
   }
   GLfloat x = p.b*(k.s.x*k.f.y-k.s.y*k.f.x) + p.c*(k.s.x*k.f.z-k.s.z*k.f.x) + p.d*(k.s.x*k.f.w-k.s.w*k.f.x); 
   GLfloat y = p.a*(k.s.y*k.f.x-k.s.x*k.f.y) + p.c*(k.s.y*k.f.z-k.s.z*k.f.y) + p.d*(k.s.y*k.f.w-k.s.w*k.f.y);
   GLfloat z = p.a*(k.s.z*k.f.x-k.s.x*k.f.z) + p.b*(k.s.z*k.f.y-k.s.y*k.f.z) + p.d*(k.s.z*k.f.w-k.s.w*k.f.z);
   GLfloat w = p.a*(k.s.w*k.f.x-k.s.x*k.f.w) + p.b*(k.s.y*k.f.w-k.s.w*k.f.y) + p.c*(k.s.z*k.f.w-k.s.w*k.f.z);
   return Point(x, y, z, w);
}

std::ostream & operator<<(std::ostream & os, const Line & k)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   os << 
      setw(width) << k.s << "->" <<
      setw(width) << k.f;
   os.flush();
   return os;
}

// Constructors for class Plane.

// Construct the plane with equation ax + by + cd + dw = 0.
Plane::Plane(GLfloat a, GLfloat b, GLfloat c, GLfloat d)
: a(a), b(b), c(c), d(d)
{
   if (a == 0 && b == 0 && c == 0)
   {
      cuglError = BAD_PLANE;
      a = 0;
      b = 1;
      c = 0;
      d = 0;
   }
}

// Construct the plane that contains points p, q, and r.
Plane::Plane(Point & p, Point & q, Point & r)
{
   a =   p[1] * (q[2]*r[3] - q[3]*r[2]) + q[1] * (p[3]*r[2] - p[2]*r[3]) + r[1] * (p[2]*q[3] - p[3]*q[2]);
   b = - p[2] * (q[3]*r[0] - q[0]*r[3]) + q[2] * (p[0]*r[3] - p[3]*r[0]) + r[2] * (p[3]*q[0] - p[0]*q[3]);
   c =   p[3] * (q[0]*r[1] - q[1]*r[0]) + q[3] * (p[1]*r[0] - p[0]*r[1]) + r[3] * (p[0]*q[1] - p[1]*q[0]);
   d = - p[0] * (q[1]*r[2] - q[2]*r[1]) + q[0] * (p[2]*r[1] - p[1]*r[2]) + r[0] * (p[1]*q[2] - p[2]*q[1]);
   if (a == 0 && b == 0 && c == 0)
   {
      cuglError = BAD_PLANE;
      a = 0;
      b = 1;
      c = 0;
      d = 0;
   }
}

Plane::Plane(Line & k, Point & p)
{
   a =   p[1] * (k.s[2]*k.f[3] - k.s[3]*k.f[2]) + k.s[1] * (p[3]*k.f[2] - p[2]*k.f[3]) + k.f[1] * (p[2]*k.s[3] - p[3]*k.s[2]);
   b = - p[2] * (k.s[3]*k.f[0] - k.s[0]*k.f[3]) + k.s[2] * (p[0]*k.f[3] - p[3]*k.f[0]) + k.f[2] * (p[3]*k.s[0] - p[0]*k.s[3]);
   c =   p[3] * (k.s[0]*k.f[1] - k.s[1]*k.f[0]) + k.s[3] * (p[1]*k.f[0] - p[0]*k.f[1]) + k.f[3] * (p[0]*k.s[1] - p[1]*k.s[0]);
   d = - p[0] * (k.s[1]*k.f[2] - k.s[2]*k.f[1]) + k.s[0] * (p[2]*k.f[1] - p[1]*k.f[2]) + k.f[0] * (p[1]*k.s[2] - p[2]*k.s[1]);
   if (a == 0 && b == 0 && c == 0)
   {
      cuglError = BAD_PLANE;
      a = 0;
      b = 1;
      c = 0;
      d = 0;
   }
}


// Construct the plane orthogonal to v that contains p.
Plane::Plane(Vector & v, Point & p)
{
   a = v[0] * p[3];
   b = v[1] * p[3];
   c = v[2] * p[3];
   d = - (a * p[0] + b * p[1] + c * p[2]);
   if (a == 0 && b == 0 && c == 0)
   {
      cuglError = BAD_PLANE;
      a = 0;
      b = 1;
      c = 0;
      d = 0;
   }
}

// Member functions for class Plane.

// Normalize the plane.
void Plane::normalize()
{
    GLfloat s = float(sqrt(a*a + b*b + c*c));
    if (s == 0)
        cuglError = ZERO_DIVISOR;
    else
    {
        a = a / s;
        b = b / s;
        c = c / s;
        d = d / s;
    }
}

// Return an equivalent normal plane.
Plane Plane::unit() const
{
    GLfloat s = float(sqrt(a*a + b*b + c*c));
    if (s == 0)
    {
        cuglError = ZERO_DIVISOR;
        return *this;
    }
    else
        return Plane(a/s, b/s, c/s, d/s);
}

// Use this plane as an OpenGL clipping plane.
void Plane::clipPlane(GLenum index) const
{
   GLdouble eqn[4];
   eqn[0] = a;
   eqn[1] = b;
   eqn[2] = c;
   eqn[3] = d;
   glClipPlane(index, &eqn[0]);
}

// Friend functions for class Plane

std::ostream & operator<<(std::ostream & os, const Plane & p)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   os << setw(0) << '<' <<
      setw(width) << p.a << ',' <<
      setw(width) << p.b << ',' <<
      setw(width) << p.c << ',' <<
      setw(width) << p.d << '>';
   os.flush();
   return os;
}


// Constructor for class Vector.

// Construct the normal to a polygon defined by a set of points.
Vector::Vector(Point points[], int numPoints)
{
   // Normal to a set of points using Newell's algorithm.
   x = 0;
   y = 0;
   z = 0;
   for (int i = 0; i < numPoints; i++)
   {
      int next = i == numPoints - 1 ? 0 : i + 1;
      x += (points[i][1] - points[next][1]) * (points[i][2] + points[next][2]);
      y += (points[i][2] - points[next][2]) * (points[i][0] + points[next][0]);
      z += (points[i][0] - points[next][0]) * (points[i][1] + points[next][1]);
   }
}

// Member functions for class Vector

// Return a reference to a component of a vector.
GLfloat & Vector::operator[](int i)
{
   switch (i)
   {
   default:
      cuglError = BAD_INDEX;
      // return x;
   case 0:
      return x;
   case 1:
      return y;
   case 2:
      return z;
   }
}

// Return a const reference to a component of a vector.
const GLfloat & Vector::operator[](int i) const
{
   switch (i)
   {
   default:
      cuglError = BAD_INDEX;
      // return x;
   case 0:
      return x;
   case 1:
      return y;
   case 2:
      return z;
   }
}

// Normalize this vector in place and return the result.
void Vector::normalize()
{
   GLfloat len = length();
   if (len == 0)
      cuglError = ZERO_NORM;
   else
   {
      x /= len;
      y /= len;
      z /= len;
   }
}

// Return a unit vector with the same direction as this vector.
Vector Vector::unit() const
{
   GLfloat len = length();
   if (len == 0)
   {
      cuglError = ZERO_NORM;
      return Vector(0, 0, 0);
   }
   else
      return Vector(x/len, y/len, z/len);
}

// Divide each component of a vector by a scaling constant.
Vector Vector::operator/(GLfloat scale) const
{
   if (scale == 0)
   {
      cuglError = ZERO_DIVISOR;
      return *this;
   }
   else
      return Vector(x / scale, y / scale, z / scale);
}

// Scale this vector and return the result.
Vector & Vector::operator/=(GLfloat scale)
{
   if (scale == 0)
      cuglError = ZERO_DIVISOR;
   else
   {
      x /= scale;
      y /= scale;
      z /= scale;
   }
   return *this;
}

std::ostream & operator<<(std::ostream & os, const Vector & v)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   os << setw(0) << '(' <<
      setw(width) << v[0] << ',' <<
      setw(width) << v[1] << ',' <<
      setw(width) << v[2] << ')';
   os.flush();
   return os;
}


// Constructors for class Matrix

// Construct the identity matrix.
Matrix::Matrix()
{
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         m[i][j] = i == j ? 1.0f : 0.0f;
}

// Copy matrix r
Matrix::Matrix(const GL_Matrix & r)
{
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         m[i][j] = r[i][j];
}

// Copy GL matrix
Matrix::Matrix(GLenum mode)
{
   if (  mode == GL_PROJECTION_MATRIX || 
         mode == GL_MODELVIEW_MATRIX)
   {
      GL_Matrix g;
      glGetFloatv(mode, &g[0][0]);
      for (int i = 0; i < 4; i++)
         for (int j = 0; j < 4; j++)
            m[i][j] = g[i][j];
   }
   else
   {
      cuglError = BAD_MATRIX_MODE;
      // Construct identity matrix.
      for (int i = 0; i < 4; i++)
         for (int j = 0; j < 4; j++)
            m[i][j] = i == j ? 1.0f : 0.0f;
   }
}

// Gives correct result for glRotatef
Matrix::Matrix(const Vector & axis, double theta)
{
   Vector u = axis.unit();
   double s = sin(theta);
   double c = cos(theta);
   double cc = 1 - c;

   m[0][0] = GLfloat(cc * u[0] * u[0] + c);
   m[1][0] = GLfloat(cc * u[1] * u[0] - s * u[2]);
   m[2][0] = GLfloat(cc * u[2] * u[0] + s * u[1]);
   m[3][0] = 0;

   m[0][1] = GLfloat(cc * u[0] * u[1] + s * u[2]);
   m[1][1] = GLfloat(cc * u[1] * u[1] + c);
   m[2][1] = GLfloat(cc * u[2] * u[1] - s * u[0]);
   m[3][1] = 0;

   m[0][2] = GLfloat(cc * u[0] * u[2] - s * u[1]);
   m[1][2] = GLfloat(cc * u[1] * u[2] + s * u[0]);
   m[2][2] = GLfloat(cc * u[2] * u[2] + c);
   m[3][2] = 0;
   
   m[0][3] = 0;
   m[1][3] = 0;
   m[2][3] = 0;
   m[3][3] = 1;
}
 
Matrix::Matrix(const Vector & u, const Vector & v)
{
   Vector axis = cross(v, u);
   double sina = axis.length();
   double cosa = sqrt(1 - sina * sina);
   axis.normalize();
   double omc = 1 - cosa;

   m[0][0] = GLfloat(omc * axis[0] * axis[0] + cosa);
   m[1][0] = GLfloat(omc * axis[1] * axis[0] - sina * axis[2]);
   m[2][0] = GLfloat(omc * axis[2] * axis[0] + sina * axis[1]);
   m[3][0] = 0;

   m[0][1] = GLfloat(omc * axis[0] * axis[1] + sina * axis[2]);
   m[1][1] = GLfloat(omc * axis[1] * axis[1] + cosa);
   m[2][1] = GLfloat(omc * axis[2] * axis[1] - sina * axis[0]);
   m[3][1] = 0;

   m[0][2] = GLfloat(omc * axis[0] * axis[2] - sina * axis[1]);
   m[1][2] = GLfloat(omc * axis[1] * axis[2] + sina * axis[0]);
   m[2][2] = GLfloat(omc * axis[2] * axis[2] + cosa);
   m[3][2] = 0;
   
   m[0][3] = 0;
   m[1][3] = 0;
   m[2][3] = 0;
   m[3][3] = 1;
}
 

Matrix::Matrix(const Quaternion & q)
{
   m[0][0] = GLfloat(1 - 2 * (q.v.y * q.v.y + q.v.z * q.v.z));
   m[1][0] = GLfloat(2 * (q.v.x * q.v.y - q.v.z * q.s));
   m[2][0] = GLfloat(2 * (q.v.z * q.v.x + q.v.y * q.s));
   m[3][0] = 0;

   m[0][1] = GLfloat(2 * (q.v.x * q.v.y + q.v.z * q.s));
   m[1][1] = GLfloat(1 - 2 * (q.v.z * q.v.z + q.v.x * q.v.x));
   m[2][1] = GLfloat(2 * (q.v.y * q.v.z - q.v.x * q.s));
   m[3][1] = 0;

   m[0][2] = GLfloat(2 * (q.v.z * q.v.x - q.v.y * q.s));
   m[1][2] = GLfloat(2 * (q.v.y * q.v.z + q.v.x * q.s));
   m[2][2] = GLfloat(1 - 2 * (q.v.y * q.v.y + q.v.x * q.v.x));
   m[3][2] = 0;

   m[0][3] = 0;
   m[1][3] = 0;
   m[2][3] = 0;
   m[3][3] = 1;
}

Matrix operator+(const Matrix & lhs, const Matrix & rhs)
{
    Matrix res;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 4; ++j)
            res.m[i][j] = lhs.m[i][j] + rhs.m[i][j];
    return res;
}

Matrix operator-(const Matrix & lhs, const Matrix & rhs)
{
    Matrix res;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 4; ++j)
            res.m[i][j] = lhs.m[i][j] - rhs.m[i][j];
    return res;
}

Matrix operator/(const Matrix & lhs, const Matrix & rhs)
{
    return lhs * rhs.inv();
}

Matrix & Matrix::operator+=(const Matrix & rhs)
{
    Matrix res = *this + rhs;
    return *this = res;
}

Matrix & Matrix::operator-=(const Matrix & rhs)
{
    Matrix res = *this - rhs;
    return *this = res;
}

Matrix & Matrix::operator*=(const Matrix & rhs)
{
    Matrix res = *this * rhs;
    return *this = res;
}

Matrix & Matrix::operator/=(const Matrix & rhs)
{
    Matrix res = *this / rhs;
    return *this = res;
}

Matrix operator*(const Matrix & m, const Matrix & n)
{
   Matrix r;
   r(0,0) = m(0,0) * n(0,0) + m(0,1) * n(1,0) + m(0,2) * n(2,0) + m(0,3) * n(3,0);
   r(0,1) = m(0,0) * n(0,1) + m(0,1) * n(1,1) + m(0,2) * n(2,1) + m(0,3) * n(3,1);
   r(0,2) = m(0,0) * n(0,2) + m(0,1) * n(1,2) + m(0,2) * n(2,2) + m(0,3) * n(3,2);
   r(0,3) = m(0,0) * n(0,3) + m(0,1) * n(1,3) + m(0,2) * n(2,3) + m(0,3) * n(3,3);
   r(1,0) = m(1,0) * n(0,0) + m(1,1) * n(1,0) + m(1,2) * n(2,0) + m(1,3) * n(3,0);
   r(1,1) = m(1,0) * n(0,1) + m(1,1) * n(1,1) + m(1,2) * n(2,1) + m(1,3) * n(3,1);
   r(1,2) = m(1,0) * n(0,2) + m(1,1) * n(1,2) + m(1,2) * n(2,2) + m(1,3) * n(3,2);
   r(1,3) = m(1,0) * n(0,3) + m(1,1) * n(1,3) + m(1,2) * n(2,3) + m(1,3) * n(3,3);
   r(2,0) = m(2,0) * n(0,0) + m(2,1) * n(1,0) + m(2,2) * n(2,0) + m(2,3) * n(3,0);
   r(2,1) = m(2,0) * n(0,1) + m(2,1) * n(1,1) + m(2,2) * n(2,1) + m(2,3) * n(3,1);
   r(2,2) = m(2,0) * n(0,2) + m(2,1) * n(1,2) + m(2,2) * n(2,2) + m(2,3) * n(3,2);
   r(2,3) = m(2,0) * n(0,3) + m(2,1) * n(1,3) + m(2,2) * n(2,3) + m(2,3) * n(3,3);
   r(3,0) = m(3,0) * n(0,0) + m(3,1) * n(1,0) + m(3,2) * n(2,0) + m(3,3) * n(3,0);
   r(3,1) = m(3,0) * n(0,1) + m(3,1) * n(1,1) + m(3,2) * n(2,1) + m(3,3) * n(3,1);
   r(3,2) = m(3,0) * n(0,2) + m(3,1) * n(1,2) + m(3,2) * n(2,2) + m(3,3) * n(3,2);
   r(3,3) = m(3,0) * n(0,3) + m(3,1) * n(1,3) + m(3,2) * n(2,3) + m(3,3) * n(3,3);
   return r;
}

// Matrix to reflect in plane p
void Matrix::reflect(const Plane & p)
{
   GLfloat q = p.a * p.a + p.b * p.b + p.c * p.c;
   m[0][0] = q - 2 * p.a * p.a;
   m[1][0] = - 2 * p.a * p.b;
   m[2][0] = - 2 * p.a * p.c;
   m[3][0] = - 2 * p.a * p.d;

   m[0][1] = - 2 * p.a * p.b;
   m[1][1] = q - 2 * p.b * p.b;
   m[2][1] = - 2 * p.b * p.c;
   m[3][1] = - 2 * p.b * p.d;

   m[0][2] = - 2 * p.a * p.c;
   m[1][2] = - 2 * p.b * p.c;
   m[2][2] = q - 2 * p.c * p.c;
   m[3][2] = - 2 * p.c * p.d;
   
   m[0][3] = 0;
   m[1][3] = 0;
   m[2][3] = 0;
   m[3][3] = q;
}


// Shadow matrix for point and plane.
// OpenGL does not seem to like matrices with m[3][3] < 0.
// Consequently, we ensure that m[3][3] >= 0.
void Matrix::shadow(const Point & s, const Plane & p)
{
   GLfloat k = p.a * s[0] + p.b * s[1] + p.c * s[2] + p.d * s[3];
   GLfloat m33 = p.d * s[3] - k;
   if (m33 >= 0)
   {
      m[0][0] = p.a * s[0] - k;
      m[1][0] = p.b * s[0];
      m[2][0] = p.c * s[0];
      m[3][0] = p.d * s[0];
   
      m[0][1] = p.a * s[1];
      m[1][1] = p.b * s[1] - k;
      m[2][1] = p.c * s[1];
      m[3][1] = p.d * s[1];

      m[0][2] = p.a * s[2];
      m[1][2] = p.b * s[2];
      m[2][2] = p.c * s[2] - k;
      m[3][2] = p.d * s[2];
   
      m[0][3] = p.a * s[3];
      m[1][3] = p.b * s[3];
      m[2][3] = p.c * s[3];
      m[3][3] = m33;
   }
   else
   {
      m[0][0] = k - p.a * s[0];
      m[1][0] = - p.b * s[0];
      m[2][0] = - p.c * s[0];
      m[3][0] = - p.d * s[0];
   
      m[0][1] = - p.a * s[1];
      m[1][1] = k - p.b * s[1];
      m[2][1] = - p.c * s[1];
      m[3][1] = - p.d * s[1];

      m[0][2] = - p.a * s[2];
      m[1][2] = - p.b * s[2];
      m[2][2] = k - p.c * s[2];
      m[3][2] = - p.d * s[2];
   
      m[0][3] = - p.a * s[3];
      m[1][3] = - p.b * s[3];
      m[2][3] = - p.c * s[3];
      m[3][3] = - m33;
   }
}

// Member functions for class Matrix

Matrix Matrix::transpose() const
{
   Matrix t;
   for (int r = 0; r < 4; r++)
      for (int c = 0; c < 4; c++)
         t(r,c) = m[c][r];
   return t;
}

GLfloat Matrix::det() const
{
    return
         - m[0][0] * m[1][1] * m[2][2] * m[3][3]
         + m[0][0] * m[1][1] * m[2][3] * m[3][2]
         + m[0][0] * m[2][1] * m[1][2] * m[3][3]
         - m[0][0] * m[2][1] * m[1][3] * m[3][2]
         - m[0][0] * m[3][1] * m[1][2] * m[2][3]
         + m[0][0] * m[3][1] * m[1][3] * m[2][2]
         + m[1][0] * m[0][1] * m[2][2] * m[3][3]
         - m[1][0] * m[0][1] * m[2][3] * m[3][2]
         - m[1][0] * m[2][1] * m[0][2] * m[3][3]
         + m[1][0] * m[2][1] * m[0][3] * m[3][2]
         + m[1][0] * m[3][1] * m[0][2] * m[2][3]
         - m[1][0] * m[3][1] * m[0][3] * m[2][2]
         - m[2][0] * m[0][1] * m[1][2] * m[3][3]
         + m[2][0] * m[0][1] * m[1][3] * m[3][2]
         + m[2][0] * m[1][1] * m[0][2] * m[3][3]
         - m[2][0] * m[1][1] * m[0][3] * m[3][2]
         - m[2][0] * m[3][1] * m[0][2] * m[1][3]
         + m[2][0] * m[3][1] * m[0][3] * m[1][2]
         + m[3][0] * m[0][1] * m[1][2] * m[2][3]
         - m[3][0] * m[0][1] * m[1][3] * m[2][2]
         - m[3][0] * m[1][1] * m[0][2] * m[2][3]
         + m[3][0] * m[1][1] * m[0][3] * m[2][2]
         + m[3][0] * m[2][1] * m[0][2] * m[1][3]
         - m[3][0] * m[2][1] * m[0][3] * m[1][2];
}

Matrix Matrix::inv() const
{
   GLfloat d =  det();
   Matrix inv;
   if (d == 0)
      cuglError = SINGULAR_MATRIX;
   else
   {
      inv.m[0][0] =  (-m[1][1] * m[2][2] * m[3][3] + m[1][1] * m[2][3] * m[3][2]
                     + m[2][1] * m[1][2] * m[3][3] - m[2][1] * m[1][3] * m[3][2]
                     - m[3][1] * m[1][2] * m[2][3] + m[3][1] * m[1][3] * m[2][2]) / d;
      inv.m[0][1] = -(-m[0][1] * m[2][2] * m[3][3] + m[0][1] * m[2][3] * m[3][2]
                     + m[2][1] * m[0][2] * m[3][3] - m[2][1] * m[0][3] * m[3][2]
                     - m[3][1] * m[0][2] * m[2][3] + m[3][1] * m[0][3] * m[2][2]) / d;
      inv.m[0][2] = - (m[0][1] * m[1][2] * m[3][3] - m[0][1] * m[1][3] * m[3][2]
                     - m[1][1] * m[0][2] * m[3][3] + m[1][1] * m[0][3] * m[3][2]
                     + m[3][1] * m[0][2] * m[1][3] - m[3][1] * m[0][3] * m[1][2]) / d;
      inv.m[0][3] =- (-m[0][1] * m[1][2] * m[2][3] + m[0][1] * m[1][3] * m[2][2]
                     + m[1][1] * m[0][2] * m[2][3] - m[1][1] * m[0][3] * m[2][2]
                     - m[2][1] * m[0][2] * m[1][3] + m[2][1] * m[0][3] * m[1][2]) / d;
      inv.m[1][0] = -(-m[1][0] * m[2][2] * m[3][3] + m[1][0] * m[2][3] * m[3][2]
                     + m[2][0] * m[1][2] * m[3][3] - m[2][0] * m[1][3] * m[3][2]
                     - m[3][0] * m[1][2] * m[2][3] + m[3][0] * m[1][3] * m[2][2]) / d;
      inv.m[1][1] =  (-m[0][0] * m[2][2] * m[3][3] + m[0][0] * m[2][3] * m[3][2]
                     + m[2][0] * m[0][2] * m[3][3] - m[2][0] * m[0][3] * m[3][2]
                     - m[3][0] * m[0][2] * m[2][3] + m[3][0] * m[0][3] * m[2][2]) / d;
      inv.m[1][2] = -(-m[0][0] * m[1][2] * m[3][3] + m[0][0] * m[1][3] * m[3][2]
                     + m[1][0] * m[0][2] * m[3][3] - m[1][0] * m[0][3] * m[3][2]
                     - m[3][0] * m[0][2] * m[1][3] + m[3][0] * m[0][3] * m[1][2]) / d;
      inv.m[1][3] =  (-m[0][0] * m[1][2] * m[2][3] + m[0][0] * m[1][3] * m[2][2]
                     + m[1][0] * m[0][2] * m[2][3] - m[1][0] * m[0][3] * m[2][2]
                     - m[2][0] * m[0][2] * m[1][3] + m[2][0] * m[0][3] * m[1][2]) / d;
      inv.m[2][0] =  (-m[1][0] * m[2][1] * m[3][3] + m[1][0] * m[2][3] * m[3][1]
                     + m[2][0] * m[1][1] * m[3][3] - m[2][0] * m[1][3] * m[3][1]
                     - m[3][0] * m[1][1] * m[2][3] + m[3][0] * m[1][3] * m[2][1]) / d;
      inv.m[2][1] = -(-m[0][0] * m[2][1] * m[3][3] + m[0][0] * m[2][3] * m[3][1]
                     + m[2][0] * m[0][1] * m[3][3] - m[2][0] * m[0][3] * m[3][1]
                     - m[3][0] * m[0][1] * m[2][3] + m[3][0] * m[0][3] * m[2][1]) / d;
      inv.m[2][2] =  -(m[0][0] * m[1][1] * m[3][3] - m[0][0] * m[1][3] * m[3][1]
                     - m[1][0] * m[0][1] * m[3][3] + m[1][0] * m[0][3] * m[3][1]
                     + m[3][0] * m[0][1] * m[1][3] - m[3][0] * m[0][3] * m[1][1]) / d;
      inv.m[2][3] =   (m[0][0] * m[1][1] * m[2][3] - m[0][0] * m[1][3] * m[2][1]
                     - m[1][0] * m[0][1] * m[2][3] + m[1][0] * m[0][3] * m[2][1]
                     + m[2][0] * m[0][1] * m[1][3] - m[2][0] * m[0][3] * m[1][1]) / d;
      inv.m[3][0] = -(-m[1][0] * m[2][1] * m[3][2] + m[1][0] * m[2][2] * m[3][1]
                     + m[2][0] * m[1][1] * m[3][2] - m[2][0] * m[1][2] * m[3][1]
                     - m[3][0] * m[1][1] * m[2][2] + m[3][0] * m[1][2] * m[2][1]) / d;
      inv.m[3][1] =  (-m[0][0] * m[2][1] * m[3][2] + m[0][0] * m[2][2] * m[3][1]
                     + m[2][0] * m[0][1] * m[3][2] - m[2][0] * m[0][2] * m[3][1]
                     - m[3][0] * m[0][1] * m[2][2] + m[3][0] * m[0][2] * m[2][1]) / d;
      inv.m[3][2] =   (m[0][0] * m[1][1] * m[3][2] - m[0][0] * m[1][2] * m[3][1]
                     - m[1][0] * m[0][1] * m[3][2] + m[1][0] * m[0][2] * m[3][1]
                     + m[3][0] * m[0][1] * m[1][2] - m[3][0] * m[0][2] * m[1][1]) / d;
      inv.m[3][3] =  -(m[0][0] * m[1][1] * m[2][2] - m[0][0] * m[1][2] * m[2][1]
                     - m[1][0] * m[0][1] * m[2][2] + m[1][0] * m[0][2] * m[2][1]
                     + m[2][0] * m[0][1] * m[1][2] - m[2][0] * m[0][2] * m[1][1]) / d;
   }
return inv;
}

double Matrix::angle() const
{
   double cosAngle = 0.5 * (m[0][0] + m[1][1] + m[2][2] - 1);
   if (fabs(cosAngle) > 1)
   {
      cuglError = BAD_ROTATION_MATRIX;
      return 0;
   }
   else
      return acos(cosAngle);
}

Vector Matrix::axis() const
{
   GLfloat twoSine = GLfloat(2 * sin(angle()));
   if (twoSine == 0)
   {
      cuglError = BAD_ROTATION_MATRIX;
      return Vector();
   }
   else
      return Vector(
         (m[1][2] - m[2][1]) / twoSine,
         (m[2][0] - m[0][2]) / twoSine,
         (m[0][1] - m[1][0]) / twoSine );
}

std::ostream & operator<<(std::ostream & os, const Matrix & m)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   for (int i = 0; i < 4; i++)
   {
      os << setw(0) << '(';
      for (int j = 0; j < 4; j++)
         os << setw(width) << m(j,i);
      os << ')' << endl;
   }
   os.flush();
   return os;
}

Quaternion Matrix::quaternion() const
{
   // Gives consistent results.
   double sqs = m[0][0] + m[1][1] + m[2][2] + m[3][3];
   if (sqs < 0)
   {
      cuglError = BAD_ROTATION_MATRIX;
      return Quaternion();
   }
   GLfloat scale = GLfloat(sqrt(sqs));
   return Quaternion(scale / 2, 
      Vector
      (
         (m[1][2] - m[2][1]) / (2 * scale),
         (m[2][0] - m[0][2]) / (2 * scale),
         (m[0][1] - m[1][0]) / (2 * scale) 
      ));
}

// Constructors for class Quaternion

Quaternion::Quaternion(Matrix m)
{
   double sqs = m(0,0) + m(1,1) + m(2,2) + m(3,3);
   if (sqs < 0)
   {
      cuglError = BAD_ROTATION_MATRIX;
      s = 1;
      v = Vector(0,0,0);
   }
   GLfloat scale = GLfloat(sqrt(sqs));
   s = scale / 2;
   v = Vector
   (
      (m(1,2) - m(2,1)) / (2 * scale),
      (m(2,0) - m(0,2)) / (2 * scale),
      (m(0,1) - m(1,0)) / (2 * scale) 
   );
}

Quaternion::Quaternion(double xr, double yr, double zr)
{
   double cos_z = cos(zr);
   double sin_z = sin(zr);
   double cos_y = cos(yr);
   double sin_y = sin(yr);
   double cos_x = cos(xr);
   double sin_x = sin(xr);
   double ds = sqrt(1 + cos_z * cos_y + cos_z * cos_x + sin_z * sin_y * sin_x + cos_y * cos_x) / 2;
   double s4 = 4 * ds;
   v.x = GLfloat((cos_y * sin_x + cos_z * sin_x - sin_z * sin_y * cos_x) / s4);
   v.y = GLfloat((sin_z * sin_x + cos_z * sin_y * cos_x + sin_y) / s4);
   v.z = GLfloat((sin_z * cos_y + sin_z * cos_x - cos_z * sin_y * sin_x) / s4);
   s = GLfloat(ds);
}

Quaternion::Quaternion (const Point & p)
{
   v[0] = p[0];
   v[1] = p[1];
   v[2] = p[2];
   s    = p[3];
}

Quaternion::Quaternion(const Vector & u, const Vector & w)
{
    Vector axis = cross(w, u);
    double angle = asin(axis.length());
    axis.normalize();
    v = GLfloat(sin(angle/2)) * axis;
    s = GLfloat(cos(angle/2));
}

// Member functions for class Quaternion

Quaternion Quaternion::operator/(GLfloat scale) const
{
   if (scale == 0)
   {
      cuglError = ZERO_DIVISOR;
      return *this;
   }
   else
      return Quaternion(s / scale, v / scale);
}

void Quaternion::normalize()
{
   GLfloat m = magnitude();
   if (m == 0)
      cuglError = ZERO_DIVISOR;
   else
   {
      s /= m;
      v /= m;
   }
}

Quaternion Quaternion::unit() const
{
   GLfloat m = magnitude();
   if (m == 0)
   {
      cuglError = ZERO_DIVISOR;
      return *this;
   }
   else
      return Quaternion(s/m, v/m);
}

Quaternion Quaternion::inv() const
{
   GLfloat n = norm();
   if (n == 0)
   {
      cuglError = ZERO_DIVISOR;
      return *this;
   }
   else
      return (*this).conj() / n;
}

void Quaternion::matrix(Matrix & m) const
{
   m(0,0) = GLfloat(1 - 2 * (v.y * v.y + v.z * v.z));
   m(1,0) = GLfloat(2 * (v.x * v.y - v.z * s));
   m(2,0) = GLfloat(2 * (v.z * v.x + v.y * s));
   m(3,0) = 0;

   m(0,1) = GLfloat(2 * (v.x * v.y + v.z * s));
   m(1,1) = GLfloat(1 - 2 * (v.z * v.z + v.x * v.x));
   m(2,1) = GLfloat(2 * (v.y * v.z - v.x * s));
   m(3,1) = 0;

   m(0,2) = GLfloat(2 * (v.z * v.x - v.y * s));
   m(1,2) = GLfloat(2 * (v.y * v.z + v.x * s));
   m(2,2) = GLfloat(1 - 2 * (v.y * v.y + v.x * v.x));
   m(3,2) = 0;

   m(0,3) = 0;
   m(1,3) = 0;
   m(2,3) = 0;
   m(3,3) = 1;
}

void Quaternion::matrix(GL_Matrix m) const
{
   // Gives same direction as glRotatef
   m[0][0] = GLfloat(1 - 2 * (v.y * v.y + v.z * v.z));
   m[1][0] = GLfloat(2 * (v.x * v.y - v.z * s));
   m[2][0] = GLfloat(2 * (v.z * v.x + v.y * s));
   m[3][0] = 0;

   m[0][1] = GLfloat(2 * (v.x * v.y + v.z * s));
   m[1][1] = GLfloat(1 - 2 * (v.z * v.z + v.x * v.x));
   m[2][1] = GLfloat(2 * (v.y * v.z - v.x * s));
   m[3][1] = 0;

   m[0][2] = GLfloat(2 * (v.z * v.x - v.y * s));
   m[1][2] = GLfloat(2 * (v.y * v.z + v.x * s));
   m[2][2] = GLfloat(1 - 2 * (v.y * v.y + v.x * v.x));
   m[3][2] = 0;

   m[0][3] = 0;
   m[1][3] = 0;
   m[2][3] = 0;
   m[3][3] = 1;
}

void Quaternion::apply() const
{
   GL_Matrix m;
   matrix(m);
   glMultMatrixf(&m[0][0]);
}

void Quaternion::euler(double & xr, double & yr, double & zr) const
{
   double sqs = s * s;
   double sqx = v.x * v.x;
   double sqy = v.y * v.y;
   double sqz = v.z * v.z;

   xr = yr = zr = 0;
   xr = atan2(static_cast<double>(2 * (v.y * v.z + v.x * s)), - sqx - sqy + sqz + sqs);
   yr = asin(- 2 * (v.x * v.z - v.y * s));
   zr = atan2(static_cast<double>(2 * (v.x * v.y + v.z * s)),sqx - sqy - sqz + sqs);
}

void Quaternion::integrate(Vector & omega, double dt)
{
   *this *= Quaternion(omega.unit(), omega.length() * dt);
}

const double BALLRADIUS = 0.8f;
const double BRADSQ = BALLRADIUS * BALLRADIUS;

// Helper function for trackball().
GLfloat project(GLfloat x, GLfloat y)
   // Expect: |x| <= 1 and |y| <= 1.
   // Project the point (x,y) onto either a sphere or a 
   // hyperboloid, depending on its distance from the origin. 
{
   double dsq = x * x + y * y;
   double d = ::sqrt(dsq);
   if (d < BALLRADIUS * 0.5 * sqrt(2.0))
      // sphere
      return GLfloat(sqrt(BRADSQ - dsq));
   else
      // hyperbola
      return GLfloat(BRADSQ / (2 * d));
}

void Quaternion::trackball(GLfloat x1, GLfloat y1, GLfloat x2, GLfloat y2)
{
   // Project the mouse points onto a sphere or hyperboloid.
   Vector v1(x1, y1, project(x1, y1));
   Vector v2(x2, y2, project(x2, y2));

   // Construct the quaternion from these vectors.
   Vector a = v2 * v1;
   Vector d = v1 - v2;
   double t = d.length() / (2 * BALLRADIUS);
   if (t > 1) t = 1;
   if (t < -1) t = -1;
   double theta = asin(t);
   (*this) *= Quaternion(GLfloat(cos(theta)), a.unit() * GLfloat(sin(theta)));
};

// Friend functions for class Quaternion

Vector ln(const Quaternion & q)
{
   double theta = q.s;
   double scale = theta == 0 ? 1 : theta / sin(theta);
   return (GLfloat(theta == 0 ? 1 : theta / sin(theta))) * q.vector();
}

Quaternion exp(const Vector & v)
{
   double theta = v.length();
   return Quaternion(GLfloat(cos(theta)), GLfloat(theta == 0 ? 1 : sin(theta) / theta) * v);
}

std::ostream & operator<<(std::ostream & os, const Quaternion & q)
{
   if ( ! os.good() )
      return os;
   int width = os.width();
   os << 
      setw(width) << q.s << ' ' << 
      setw(width) << q.v;
   os.flush();
   return os;
}

// Implementation of class Rotation

Rotation::Rotation(unsigned int maxSteps, bool cycle)
: Interpolator(maxSteps, cycle)
{
    initialize(I, J);
}

Rotation::Rotation(const Vector & first, const Vector & last, unsigned int maxSteps, bool cycle)
: Interpolator(maxSteps, cycle)
{
    initialize(first, last);
}

void Rotation::set(const Vector & first, const Vector & last)
{
    initialize(first, last);
    begin();
}

void Rotation::set(const Vector & next)
{
    initialize(Vector(get() * vFirst), next);
    begin();
}

Quaternion Rotation::get() const
{
    double t = getParam();
    return sinOmega == 0 ? first :
        first * GLfloat((sin(omega * (1 - t)) / sinOmega)) +
        last * GLfloat((sin(omega * t) / sinOmega));
}

void Rotation::apply()
{
    Quaternion q = get();
    q.apply();
}

void Rotation::initialize(const Vector & u, const Vector & v)
{
    double cosOmega = dot(u, v);
    if (cosOmega < 0)
        cosOmega = - cosOmega;
    omega = acos(cosOmega);
    sinOmega = sin(omega);
    if (sinOmega == 0)
        cuglError = BAD_INTERPOLATOR_ARG;
    Vector axis = cross(u, v);
    first = Quaternion(axis, 0);
    last = Quaternion(axis, omega);
}

// Implementation of class Translation

Translation::Translation(int maxSteps, bool repeat)
    : Interpolator(maxSteps, repeat) {}

Translation::Translation(const Point & first, const Point & last, int maxSteps, bool repeat)
    : Interpolator(maxSteps, repeat), first(first), last(last) {}

void Translation::set(const Point & pFirst, const Point & pLast)
{
    first = pFirst;
    last = pLast;
    begin();
}

void Translation::set(const Point & pNew)
{
    first = get();
    last = pNew;
    begin();
}

void Translation::apply()
{
    Point p = get();
    p.translate();
}

Point Translation::get() const
{
    double t = getParam();
    return Point(
        GLfloat((1 - t) * first[0] + t * last[0]),
        GLfloat((1 - t) * first[1] + t * last[1]),
        GLfloat((1 - t) * first[2] + t * last[2]),
        GLfloat((1 - t) * first[3] + t * last[3]) );
}

// Implementation of class Scaling

Scaling::Scaling(int maxSteps, bool repeat)
    : Interpolator(maxSteps, repeat), first(0), last(1) {}

Scaling::Scaling(double first, double last, int maxSteps, bool repeat)
    : Interpolator(maxSteps, repeat), first(first), last(last) {}

void Scaling::set(double pFirst, double pLast)
{
    first = pFirst;
    last = pLast;
    begin();
}

void Scaling::set(double vNew)
{
    first = get();
    last = vNew;
    begin();
}

double Scaling::get() const
{
    double t = getParam();
    return (1 - t) * first + t * last;
}

// Implementation of class Camera

Point Camera::eyeGet() const
{
    double t = getParam();
    return Point(
        GLfloat(eyeFirst[0] * (1 - t) + eyeLast[0] * t),
        GLfloat(eyeFirst[1] * (1 - t) + eyeLast[1] * t),
        GLfloat(eyeFirst[2] * (1 - t) + eyeLast[2] * t),
        GLfloat(eyeFirst[3] * (1 - t) + eyeLast[3] * t) );
}

Vector Camera::dirGet() const
{
    double t = getParam();
    Quaternion q;
    if (sinOmega == 0)
        q = qFirst;
    else
        q = qFirst * GLfloat((sin(omega * (1 - t)) / sinOmega)) +
        qLast * GLfloat((sin(omega * t) / sinOmega));
    return Vector(q * dirFirst);
}

void Camera::apply()
{
    Point eye;
    Point mod;
    switch (mode)
    {
    case TRANSLATE:
        eye = eyeGet();
        mod = eye + distance * dirLast;
        break;
    case ROTATE:
        eye = eyeLast;
        mod = eye + distance * dirGet();
        break;
    }
    lookAt(eye, mod, up);
}

void Camera::set(Point p)
{
    cout << "Set point 1 " << *this;
    if (mode == ROTATE)
        dirFirst = dirLast = dirGet();
    eyeFirst = eyeGet();
    eyeLast = p;
    mode = TRANSLATE;
    begin();
    cout << "Set point 2 " << *this;
}

void Camera::set(Quaternion q)
{
    cout << "Set quat 1 " << *this;
    if (mode == TRANSLATE)
        eyeFirst = eyeLast = eyeGet();
    dirFirst = dirGet();
    dirLast = Vector(q * dirFirst);
    double cosOmega = dot(dirFirst, dirLast);
    if (cosOmega < 0)
        cosOmega = - cosOmega;
    omega = acos(cosOmega);
    sinOmega = sin(omega);
    qFirst = Quaternion(I, 0);
    qLast = q;
    mode = ROTATE;
    begin();
    cout << "Set quat 2 " << *this;
}

void Camera::moveForward(GLfloat distance)
{
    set(eyeLast + distance * dirLast);
}

void Camera::moveBackward(GLfloat distance)
{
    set(eyeLast - distance * dirLast);
}

void Camera::moveUp(GLfloat distance)
{
    set(eyeLast + distance * up);
}

void Camera::moveDown(GLfloat distance)
{
    set(eyeLast - distance * up);
}

void Camera::moveLeft(GLfloat distance)
{
    set(eyeLast + distance * cross(up, dirLast));
}

void Camera::moveRight(GLfloat distance)
{
    set(eyeLast - distance * cross(up, dirLast));
}

void Camera::panLeft(double angle)
{
    set(Quaternion(up, -angle));
}

void Camera::panRight(double angle)
{
    set(Quaternion(up, angle));
}

void Camera::tiltUp(double angle)
{
    set(Quaternion(cross(dirFirst, up), angle));
}

void Camera::tiltDown(double angle)
{
    set(Quaternion(cross(dirFirst, up), -angle));
}

ostream & operator<<(ostream & os, const Camera & cam)
{
    Point eye;
    Point mod;
    switch (cam.mode)
    {
    case Camera::TRANSLATE:
        eye = cam.eyeGet();
        mod = eye + cam.distance * cam.dirLast;
        break;
    case Camera::ROTATE:
        eye = cam.eyeLast;
        mod = eye + cam.distance * (cam.dirGet());
        break;
    }
    return os << "Eye: " << eye << "  Model: " << mod << "  Up: " << cam.up << endl;
}


// Constructor and destructor for class Revolute.

Revolute::Revolute(int numSteps, GLfloat profile[][2]) : 
   numSteps(numSteps), 
   numSlices(45), 
   eccentricity(0),
   ready(false)
{
   coor = new GLfloatArray[numSteps];
   for (int n = 0; n < numSteps; n++)
   {
      coor[n] = new GLfloat[2];
      for (int c = 0; c < 2; c++)
         coor[n][c] = profile[n][c];
   }
}

Revolute::~Revolute()
{
   for (int n = 0; n < numSteps; n++)
      delete [] coor[n];
   delete [] coor;
   delete [] texCoor;
   delete [] points;
   delete [] normals;
   delete [] faceNormals;
}

// Member functions for class Revolute.

void Revolute::setSlices(int slices)
{
   if (slices % 2 == 0)
      slices++;
   if (slices < 3)
      slices = 3;
   numSlices = slices;
   ready = false;
}

void Revolute::setEccentricity(double ecc)
{
   if (ecc < 0 || ecc >= 1)
      ecc = 0;
   eccentricity = ecc;
   ready = false;
}

void Revolute::process()
{
   // Allocate space for working arrays.
   texCoor = new GLfloat[numSteps];
   points = new Point[numSteps * numSlices];
   normals = new Vector[numSteps * numSlices];
   faceNormals = new Vector[numSteps * numSlices];
   
   // Indexes
   int dCurr, dNext, pCurr, pNext, sCurr, sNext;

   // Compute texture coordinates along axis and then normalize them.
   GLfloat texmax = 0;
   texCoor[0] = 0;
   for (pCurr = 1; pCurr < numSteps; pCurr++)
   {
      GLfloat dist = GLfloat(sqrt(sqr(coor[pCurr][0] - coor[pCurr-1][0]) + sqr(coor[pCurr][1] - coor[pCurr-1][1])));
      texCoor[pCurr] = texCoor[pCurr-1] + dist;
      texmax += dist;
   }
   for (pCurr = 1; pCurr < numSteps; pCurr++)
      texCoor[pCurr] = texCoor[pCurr] / texmax;
   
   // Compute 3D points on the surface of the solid.
   for (sCurr = 0; sCurr < numSlices; sCurr++)
   {
      double theta = (2 * PI * sCurr) / numSlices;
      GLfloat sint = GLfloat(sin(theta));
      GLfloat cost = GLfloat(cos(theta));
      if (eccentricity != 0)
      {
         double e = 1 - sqr(eccentricity);
         sint *= GLfloat(e);
         cost /= GLfloat(sqrt(e));
      }
      for (pCurr = 0; pCurr < numSteps; pCurr++)
      {
         points[numSteps * sCurr + pCurr] = 
            Point(
            coor[pCurr][0] * cost, 
            coor[pCurr][0] * sint, 
            coor[pCurr][1] );
      }
   }
   
   // Find the normal for each quadrilateral face 
   // of the solid, using Newell's method.
   for (pCurr = 0; pCurr < numSteps - 1; pCurr++)
   {
      pNext = pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr++)
      {
         sNext = sCurr == numSlices - 1 ? 0 : sCurr + 1;
         
         // Get the vertexes.
         Point quad[4];
         quad[0] = points[numSteps * sCurr + pCurr];
         quad[1] = points[numSteps * sNext + pCurr];
         quad[2] = points[numSteps * sNext + pNext];
         quad[3] = points[numSteps * sCurr + pNext];
         
         // Find the normal and store it.
         GLfloat xnorm = 0;
         GLfloat ynorm = 0;
         GLfloat znorm = 0;
         for (dCurr = 0; dCurr < 3; dCurr++)
         {
            dNext = dCurr == 2 ? 0 : dCurr + 1;
            xnorm += (quad[dCurr][1] - quad[dNext][1]) * (quad[dCurr][2] + quad[dNext][2]);
            ynorm += (quad[dCurr][2] - quad[dNext][2]) * (quad[dCurr][0] + quad[dNext][0]);
            znorm += (quad[dCurr][0] - quad[dNext][0]) * (quad[dCurr][1] + quad[dNext][1]);
         }
         faceNormals[numSteps * sCurr + pCurr] = Vector(xnorm, ynorm, znorm);
      }
   }
   
   // Find normals at each vertex by averaging face normals.
   for (pCurr = 0; pCurr < numSteps; pCurr++)
   {
      pNext = pCurr == numSteps - 1 ? 0 : pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr++)
      {
         sNext = sCurr == 0 ? numSlices - 1 : sCurr - 1;
         Vector norm;
         if (pCurr == 0)
            norm = 
            faceNormals[numSteps * sCurr + pCurr] +
            faceNormals[numSteps * sNext + pCurr];
         else if (pCurr == numSteps - 1)
            norm = 
            faceNormals[numSteps * sCurr + pCurr - 1] +
            faceNormals[numSteps * sNext + pCurr - 1];
         else
            norm = 
            faceNormals[numSteps * sCurr + pCurr - 1] +
            faceNormals[numSteps * sNext + pCurr - 1] +
            faceNormals[numSteps * sCurr + pCurr] +
            faceNormals[numSteps * sNext + pCurr];
         normals[numSteps * sCurr + pCurr] = norm.unit();
      }
   }
   ready = true;
}

void Revolute::draw(bool showNormals)
{
   if (!ready)
      return;

   int pCurr, pNext, sCurr, sNext;
   GLfloat slices = GLfloat(numSlices);
   glBegin(GL_QUAD_STRIP);
   for (pCurr = 0; pCurr < numSteps - 1; pCurr++)
   {
      pNext = pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr += 2)
      {
         sNext = sCurr == numSlices - 1 ? 0 : sCurr + 1;

         glTexCoord2f(texCoor[pCurr], sCurr / slices);
         normals[numSteps * sCurr + pCurr].drawNormal(); 
         points[numSteps * sCurr + pCurr].draw();

         glTexCoord2f(texCoor[pNext], sCurr / slices);
         normals[numSteps * sCurr + pNext].drawNormal(); 
         points[numSteps * sCurr + pNext].draw();

         glTexCoord2f(texCoor[pCurr], sNext / slices);
         normals[numSteps * sNext + pCurr].drawNormal(); 
         points[numSteps * sNext + pCurr].draw();

         glTexCoord2f(texCoor[pNext], sNext / slices);
         normals[numSteps * sNext + pNext].drawNormal(); 
         points[numSteps * sNext + pNext].draw();
      }
   }
   glEnd();
   if (showNormals)
   {
      GLboolean lighting;
      glGetBooleanv(GL_LIGHTING, &lighting);
      if (lighting)
         glDisable(GL_LIGHTING);
      glColor3f(1, 1, 1);
      glBegin(GL_LINES);
      for (pCurr = 0; pCurr < numSteps; pCurr++)
         for (sCurr = 0; sCurr < numSlices; sCurr++)
         {
            int i = numSteps * sCurr + pCurr; 
            points[i].draw();
            (points[i] + normals[i]).draw();
         }
      glEnd();
      if (lighting)
         glEnable(GL_LIGHTING);
   }
}

// Construct a solid of revolution
    
void revolve(int numSteps, GLfloat coor[][2], int numSlices, bool drawNormals)
{
   // The number of slices must be odd.
   if (numSlices % 2 == 0)
      numSlices++;

   // Allocate space for working arrays.
   Point *points = new Point[numSteps * numSlices];
   Vector *normals = new Vector[numSteps * numSlices];
   Vector *faceNormals = new Vector[numSteps * numSlices];

   // Indexes
   int dCurr, dNext, pCurr, pNext, sCurr, sNext;

   // Compute 3D points on the surface of the solid.
   for (sCurr = 0; sCurr < numSlices; sCurr++)
   {
      double theta = (2 * PI * sCurr) / numSlices;
      GLfloat sint = GLfloat(sin(theta));
      GLfloat cost = GLfloat(cos(theta));
      for (pCurr = 0; pCurr < numSteps; pCurr++)
      {
         points[numSteps * sCurr + pCurr] = 
            Point(
               coor[pCurr][0] * cost, 
               coor[pCurr][0] * sint, 
               coor[pCurr][1] );
      }
   }

   // Find the normal for each quadrilateral face 
   // of the solid, using Newell's method.
   for (pCurr = 0; pCurr < numSteps - 1; pCurr++)
   {
      pNext = pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr++)
      {
         sNext = sCurr == numSlices - 1 ? 0 : sCurr + 1;

         // Get the vertexes.
         Point quad[4];
         quad[0] = points[numSteps * sCurr + pCurr];
         quad[1] = points[numSteps * sNext + pCurr];
         quad[2] = points[numSteps * sNext + pNext];
         quad[3] = points[numSteps * sCurr + pNext];

         // Find the normal and store it.
         GLfloat xnorm = 0;
         GLfloat ynorm = 0;
         GLfloat znorm = 0;
         for (dCurr = 0; dCurr < 3; dCurr++)
         {
            dNext = dCurr == 2 ? 0 : dCurr + 1;
            xnorm += (quad[dCurr][1] - quad[dNext][1]) * (quad[dCurr][2] + quad[dNext][2]);
            ynorm += (quad[dCurr][2] - quad[dNext][2]) * (quad[dCurr][0] + quad[dNext][0]);
            znorm += (quad[dCurr][0] - quad[dNext][0]) * (quad[dCurr][1] + quad[dNext][1]);
         }
         faceNormals[numSteps * sCurr + pCurr] = Vector(xnorm, ynorm, znorm);
      }
   }

   // Find normals at each vertex by averaging face normals.
   for (pCurr = 0; pCurr < numSteps; pCurr++)
   {
      pNext = pCurr == numSteps - 1 ? 0 : pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr++)
      {
         sNext = sCurr == 0 ? numSlices - 1 : sCurr - 1;
         Vector norm;
         if (pCurr == 0)
            norm = 
               faceNormals[numSteps * sCurr + pCurr] +
               faceNormals[numSteps * sNext + pCurr];
         else if (pCurr == numSteps - 1)
            norm = 
               faceNormals[numSteps * sCurr + pCurr - 1] +
               faceNormals[numSteps * sNext + pCurr - 1];
         else
            norm = 
               faceNormals[numSteps * sCurr + pCurr - 1] +
               faceNormals[numSteps * sNext + pCurr - 1] +
               faceNormals[numSteps * sCurr + pCurr] +
               faceNormals[numSteps * sNext + pCurr];
         normals[numSteps * sCurr + pCurr] = norm.unit();

         // Draw normals if requested.
         if (drawNormals)
         {
            glBegin(GL_LINES);
               points[numSteps * sCurr + pCurr].draw();
               (points[numSteps * sCurr + pCurr] + norm).draw();
            glEnd();
         }
      }
   }
   
   // Tell OpenGL about each vertex and its normal.
   glBegin(GL_QUAD_STRIP);
   for (pCurr = 0; pCurr < numSteps - 1; pCurr++)
   {
      pNext = pCurr + 1;
      for (sCurr = 0; sCurr < numSlices; sCurr += 2)
      {
         sNext = sCurr == numSlices - 1 ? 0 : sCurr + 1;
         normals[numSteps * sCurr + pCurr].drawNormal();
         points[numSteps * sCurr + pCurr].draw();
         normals[numSteps * sCurr + pNext].drawNormal();
         points[numSteps * sCurr + pNext].draw();
         normals[numSteps * sNext + pCurr].drawNormal();
         points[numSteps * sNext + pCurr].draw();
         normals[numSteps * sNext + pNext].drawNormal();
         points[numSteps * sNext + pNext].draw();
      }
   }
   glEnd();

   delete [] points;
   delete [] normals;
   delete [] faceNormals;
}

// Utility functions for reading a BMP file.

// Read low-order byte then high-order byte
unsigned short getShort(ifstream & is)
{
   unsigned char c1, c2;
   c1=is.get();
   c2=is.get();
   return ( (c2 << 8) | c1 );
}

// Read bytes from low-order to high-order
unsigned long getULong(ifstream & is)
{
   unsigned char c1, c2, c3, c4;
   c1=is.get();
   c2=is.get();
   c3=is.get();
   c4=is.get();
   return ( (c4 << 24) | (c3 << 16) | (c2 << 8) | c1);
}

// Read bytes from low-order to high-order
long getLong(ifstream & is)
{
   unsigned char c1, c2, c3, c4;
   c1=is.get();
   c2=is.get();
   c3=is.get();  
   c4=is.get();
   return (long(c4) << 24) | (long(c3) << 16) | (long(c2) << 8) | long(c1);
}

// Write low-order byte then high-order byte
void putShort(ofstream & os, unsigned short us)
{
   os.put((unsigned char)(us & 0xFF));
   os.put((unsigned char)(us >> 8));
};

// Write bytes from low-order to high-order
void putULong(ofstream & os, unsigned long ul)
{
   os.put((unsigned char)( ul        & 0xFF));
   os.put((unsigned char)((ul >>  8) & 0xFF));
   os.put((unsigned char)((ul >> 16) & 0xFF));
   os.put((unsigned char)((ul >> 24) & 0xFF));
};

// Write bytes from low-order to high-order
void putLong(ofstream & os, long n)
{
   os.put((unsigned char)( n        & 0xFF));
   os.put((unsigned char)((n >>  8) & 0xFF));
   os.put((unsigned char)((n >> 16) & 0xFF));
   os.put((unsigned char)((n >> 24) & 0xFF));
};

int padding(unsigned long cols)
{
   int rowLength = ((3 * cols + 3) / 4) * 4;
   return rowLength - 3 * cols;
}

bool isPowerOfTwo(unsigned long ul)
{
   int c = 0;
   while (ul)
   {
      c += ul & 1;
      ul >>= 1;
   }
   return c == 1;
}

unsigned long highBit(unsigned long ul)
{
   unsigned long m = 1 << 31;
   while ((m & ul) == 0)
      m >>= 1;
   return m << 1;
}

// Constructors and destructor for class PixelMap

PixelMap::PixelMap() : 
   numRows(0), 
   numCols(0), 
   pixels(NULL), 
   size(0)
   {}

PixelMap::PixelMap(const string & bmpFileName) : 
   numRows(0), 
   numCols(0), 
   pixels(NULL), 
   size(0)
   {
      read(bmpFileName);
   }

PixelMap::PixelMap(GLint x, GLint y, GLsizei w, GLsizei h) :
   numCols(w - w % 4),
   numRows(h - h % 4),
   pixels(NULL),
   size(0)
{
   if (allocate(3 * numCols * numRows))
      glReadPixels(x, y, numCols, numRows, GL_RGB, GL_UNSIGNED_BYTE, pixels);
}

PixelMap::~PixelMap()
{
   delete pixels;
}

// Member functions of class PixelMap

bool PixelMap::allocate(unsigned long newSize)
{
   // If the new size is not greater than the space we have
   // already allocated, there is no need to reallocate.
   if (newSize == size)
      return true;

   // Attempt to get the memory we need, return if failure.
   unsigned char *mem = new unsigned char[newSize];
   if (mem == NULL)
   {
      cuglError = NOT_ENOUGH_MEMORY;
      return false;
   }

   // Update data members.
   delete [] pixels;
   pixels = mem;
   size = newSize;
   return true;
}

void PixelMap::read(const string & bmpFileName)
{
   //ifstream inf(bmpFileName, ios::in | ios::binary | ios::nocreate);
   //in linux, there is no mode of "nocreate"
   ifstream inf(bmpFileName.c_str(), ios::in | ios::binary );
   if (!inf)
   {
      cuglError = OPEN_FAILED;
      return;
   }

   char b, m;
   inf.get(b);
   inf.get(m);
   if (b != 'B' || m != 'M')
   {
      cuglError = NOT_BMP_FILE;
      inf.close();
      return;
   }

   unsigned long  fileSize   = getULong(inf);   
   unsigned short res1       = getShort(inf);      // should be 0
   unsigned short res2       = getShort(inf);      // should be 0
   unsigned long  offset     = getULong(inf);      // unreliable offset to image
   unsigned long  headerSize = getULong(inf);      // should be 40
   long           cols       = getLong(inf);	   	// # columns
   long           rows       = getLong(inf);       // # rows
   unsigned short planes     = getShort(inf);      // usually 1
   unsigned short bpx        = getShort(inf);      // bits per pixel
   unsigned long  compressed = getULong(inf);      // 0 for uncompressed
   unsigned long  imageSize  = getULong(inf);      // bytes in image (unreliable - may be zero)
   long           xPels      = getLong(inf);       // may be 0
   long           yPels      = getLong(inf);       // may be 0
   unsigned long  entries    = getULong(inf);      // should be 0 or 256   
   unsigned long  impColours = getULong(inf);      // should be 0

   if (bpx != 24)
   {
      cuglError = NOT_24_BITS;
      inf.close();
      return;
   }

   if (compressed != 0)
   {
      cuglError = COMPRESSED_BMP_FILE;
      inf.close();
      return;
   }

   // Deal with the variant of BMP file: if height is negative
   if (rows <= 0)
	   rows = -rows;

   // Allocate memory
   if ( ! allocate(3 * rows * cols))
   {
      cuglError = NOT_ENOUGH_MEMORY;
      return;
   }

   // Ensure each row has a multiple of 4 bytes
   int padLength = padding(cols);

   fileName = bmpFileName;
   numRows = rows;
   numCols = cols;

   // Transfer data from BMP file to memory.
   long count = 0;
   for (long row = 0; row < rows; row++)
   {
      for (long col = 0; col < cols; col++)		
      {  
         unsigned char red, green, blue;
         
         blue=inf.get(); 
         green=inf.get();
         red=inf.get();
	 //inf.get(blue);
         //inf.get(green);
         //inf.get(red);
         pixels[count++] = red;
         pixels[count++] = green;
         pixels[count++] = blue;
      }
      for (int k = 0; k < padLength; k++)
      {
         char c;
         inf >> c;
      }
   }
   inf.close();
}

void PixelMap::read(GLint x, GLint y, GLsizei w, GLsizei h) 
{
    if (allocate(3 * (w - w % 4) * (h - h % 4)))
    {
        numCols = w - w % 4;
        numRows = h - h % 4;
        glReadPixels(x, y, numCols, numRows, GL_RGB, GL_UNSIGNED_BYTE, pixels);
        fileName = "";
    }
}

void PixelMap::write(const string & bmpFileName) const
{
   const unsigned long HEADER_SIZE = 40;
   const unsigned long FILE_HEADER_SIZE = 14;
   const unsigned short PLANES = 1;
   const unsigned short BITS_PER_PIXEL = 24;
   const unsigned long UNCOMPRESSED = 0;
   const long XPELS = 0;
   const long YPELS = 0;
   const unsigned long ENTRIES = 0;
   const unsigned long COLOURS = 0;

   ofstream os(bmpFileName.c_str(), ios::binary);
   if (!os)
   {
      cuglError = OPEN_FAILED;
      return;
   }

   // Ensure each row has a multiple of 4 bytes
   int padLength = padding(numCols);
   unsigned long imageSize = numRows * (3 * numCols + padLength);
   unsigned long fileSize = imageSize + HEADER_SIZE + FILE_HEADER_SIZE;

   // Write file header
   os.put('B');
   os.put('M');
   putULong(os, fileSize);
   putShort(os, 0);
   putShort(os, 0);
   putULong(os, HEADER_SIZE + FILE_HEADER_SIZE);

   // Write header
   putULong(os, HEADER_SIZE);
   putLong (os, numCols);
   putLong (os, numRows);
   putShort(os, PLANES);
   putShort(os, BITS_PER_PIXEL);
   putULong(os, UNCOMPRESSED);
   putULong(os, imageSize);
   putLong (os, XPELS);
   putLong (os, YPELS);
   putULong(os, ENTRIES);
   putULong(os, COLOURS);

   // Write data
   long count = 0;
   for (unsigned long row = 0; row < numRows; row++)
   {
      for (unsigned long col = 0; col < numCols; col++)
      {
         unsigned char red   = pixels[count++];
         unsigned char green = pixels[count++];
         unsigned char blue  = pixels[count++];
         os.put(blue);
         os.put(green);
         os.put(red);
      }
      for (int k = 0; k < padLength; k++)
         os.put(' ');
   }
   os.close();
}

void PixelMap::draw() const
{
   if (pixels == NULL)
   {
      cuglError = NO_PIX_MAP;
      return;
   }
   glDrawPixels(GLsizei(numCols), GLsizei(numRows), GL_RGB, GL_UNSIGNED_BYTE, pixels);
}

void PixelMap::setTexture(GLuint name)
{
   if (pixels == NULL)
   {
      cuglError = NO_PIX_MAP;
      return;
   }
   glBindTexture(GL_TEXTURE_2D, name);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
   glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, numCols, numRows, 0, GL_RGB, GL_UNSIGNED_BYTE, pixels);
}

void PixelMap::setMipmaps(GLuint name)
{
   if (pixels == NULL)
   {
      cuglError = NO_PIX_MAP;
      return;
   }
   glBindTexture(GL_TEXTURE_2D, name);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST );
   gluBuild2DMipmaps(GL_TEXTURE_2D, GL_RGB, numCols, numRows, GL_RGB, GL_UNSIGNED_BYTE, pixels);
}

bool PixelMap::badSize() const
{
   return ( ! isPowerOfTwo(numCols) || ! isPowerOfTwo(numRows) );
}

void PixelMap::rescale()
{
   unsigned long cols = isPowerOfTwo(numCols) ? numCols : highBit(numCols);
   unsigned long rows = isPowerOfTwo(numRows) ? numRows : highBit(numRows);
   unsigned char *mem = new unsigned char[3 * rows * cols];
   if (mem == NULL)
   {
      cuglError = NOT_ENOUGH_MEMORY;
      return;
   }
	gluScaleImage(GL_RGB, numCols, numRows, GL_UNSIGNED_BYTE, pixels, 
					  cols, rows, GL_UNSIGNED_BYTE, mem);            
   delete [] pixels;
   pixels = mem;
	numCols = cols;             
	numRows = rows;   
   size = 3 * rows * cols;
}

std::ostream & operator<<(std::ostream & os, const PixelMap & pm)
{
   if (pm.fileName.length() > 0)
      os << pm.fileName;
   os << " (" <<
      pm.numRows << " rows, " <<
      pm.numCols << " columns, " <<
      pm.size << " bytes)";
   return os;
}

// Implementation of normals for triangle strips
void triStripNormals(Point points[], Vector normals[], int numPoints)
{
    const int MAXPOINTS = 100;
    Vector sides[MAXPOINTS];
    Vector triangleNormals[MAXPOINTS];
    int i;

    if (numPoints > MAXPOINTS)
    {
        cuglError = TOO_MANY_POINTS;
        return;
    }

    for (i = 0; i < numPoints - 1; i++)
        sides[i] = points[i + 1] - points[i];

    bool even = true;
    for (i = 0; i < numPoints - 2; i++)
    {
        if (even)
            triangleNormals[i] = cross(sides[i], sides[i + 1]);
        else
            triangleNormals[i] = cross(sides[i + 1], sides[i]);
        even = ! even;
    }

    for (i = 0; i < numPoints; i++)
    {
        if (i - 2 >= 0) normals[i] += triangleNormals[i-2];
        if (i - 1 >= 0 && i < numPoints - 1) normals[i] += triangleNormals[i-1];
        if (i < numPoints - 2) normals[i] += triangleNormals[i];
        normals[i].normalize();
    }
}

// Implementation of model functions

void axes(GLfloat size)
{
   const int STACKS = 12;
   const int SLICES = 20;
   GLfloat POINTS[STACKS][2] = 
   {
      { 0,     -0.1f  },
      { 0.04f,  0     },
      { 0.04f,  0.2f  },
      { 0.04f,  0.4f  },
      { 0.044f, 0.6f  },
      { 0.056f, 0.75f },
      { 0.1f,   0.75f },
      { 0.064f, 0.8f  },
      { 0.036f, 0.85f },
      { 0.016f, 0.9f  },
      { 0.004f, 0.95f },
      { 0,       1    }
   };

   GLfloat pts[STACKS][2];
   for (int x = 0; x < STACKS; x++)
      for (int y = 0; y < 2; y++)
         pts[x][y] = size * POINTS[x][y];

   // Store one axis in a call list.
   GLuint list = glGenLists(1);
   glNewList(list, GL_COMPILE);
   revolve(STACKS, pts, STACKS);
   glEndList();

   // X axis
   setMaterial(RED);
   glPushMatrix();
   glRotatef(90, 0, 1, 0);
   glCallList(list);
   glPopMatrix();

   // Y axis
   setMaterial(GREEN);
   glPushMatrix();
   glRotatef(-90, 1, 0, 0);
   glCallList(list);
   glPopMatrix();

   // Z axis
   setMaterial(BLUE);
   glPushMatrix();
   glCallList(list);
   glPopMatrix();

   glDeleteLists(list, 1);
}


GLuint makePlaneList(bool shadow)
{
   GLuint index = glGenLists(1);
   glNewList(index, GL_COMPILE);
   buildPlane(shadow);
   glEndList();
   return index;
}

void buildPlane(bool shadow)
{
   int p;
   // Set up for Bezier surface generation.
   glEnable(GL_MAP2_VERTEX_3);
   glEnable(GL_AUTO_NORMAL);
   if (shadow)
      setMaterial(BLACK);
   else
      setMaterial(METAL);

   // Fuselage
   const int fuWidth = 4;
   const int fuLength = 6;
   const int fuLoops = 20;
   const int fuSlices = 20;
   const GLfloat fuShapeFactor = 0.9f;
   GLfloat fuPoints[fuLength][fuWidth][3];

   struct { GLfloat len; GLfloat size; } fuParameters[fuLength] =
   {
      {   -10, 0 },
      { -9.6f, 1.4f },
      {    -9, 1.6f },
      {     8, 1.4f },
      {  9.9f, 1 },
      {    10, 0 }
   };
   
   for (p = 0; p < fuLength; p++)
   {
      for (int y = 0; y < fuWidth; y++)
         fuPoints[p][y][2] = fuParameters[p].len;
      fuPoints[p][0][0] = 0;
      fuPoints[p][1][0] = fuParameters[p].size;
      fuPoints[p][2][0] = fuParameters[p].size;
      fuPoints[p][3][0] = 0;

      fuPoints[p][0][1] = - fuShapeFactor * fuParameters[p].size;
      fuPoints[p][1][1] = - fuShapeFactor * fuParameters[p].size;
      fuPoints[p][2][1] = fuShapeFactor * fuParameters[p].size;
      fuPoints[p][3][1] = fuShapeFactor * fuParameters[p].size;
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1, 3,           fuWidth,
      0, 1, 3 * fuWidth, fuLength, 
      &fuPoints[0][0][0]);
   glMapGrid2f(fuLoops, 0, 1, fuSlices, 0, 1);

   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);
   glScalef(-1, 1, 1);
   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);

   // Wings
   const int wiLength = 6;
   const GLfloat wiAspectRatio = 0.2f;
   const GLfloat wiHeight = -0.6f;
   const GLfloat wiDihedral = 0.1f;
   GLfloat xVals[wiLength] = { 0.5f, 5, 8, 11, 14, 15 };
   GLfloat wiPoints[wiLength][4][3];

   for (p = 0; p < wiLength; p++)
   {
      GLfloat x = xVals[p];
      GLfloat y = wiHeight + x * wiDihedral;
      GLfloat zTrailingEdge = (x - 1) / 7 + 2;
      GLfloat zGuide = 8 * (x - 1) / 14 - 7;
      GLfloat yGuide = wiAspectRatio * (zGuide - zTrailingEdge);

      wiPoints[p][0][0] = x;
      wiPoints[p][0][1] = y;
      wiPoints[p][0][2] = zTrailingEdge;

      wiPoints[p][1][0] = x;
      wiPoints[p][1][1] = p == 5 ? y : y + yGuide;
      wiPoints[p][1][2] = zGuide;

      wiPoints[p][2][0] = x;
      wiPoints[p][2][1] = p == 5 ? y: y - yGuide;
      wiPoints[p][2][2] = zGuide;

      wiPoints[p][3][0] = x;
      wiPoints[p][3][1] = y;
      wiPoints[p][3][2] = zTrailingEdge;
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1,  3, 4,
      0, 1, 12, wiLength, 
      &wiPoints[0][0][0]);
   glMapGrid2f(fuLoops, 0, 1, fuSlices, 0, 1);

   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);
   glScalef(-1, 1, 1);
   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);

   // Tailplane
   const int taWidth = 4;
   const int taLength = 6;
   const GLfloat taAspectRatio = 0.3f;
   const GLfloat taHeight = 0.1f;
   const GLfloat taDihedral = 0.3f;
   GLfloat xTa[taLength] = { 0.3f, 1.7f, 2.6f, 3.3f, 4.6f, 5 };
   GLfloat taPoints[taLength][taWidth][3];

   for (p = 0; p < taLength; p++)
   {
      GLfloat x = xTa[p];
      GLfloat y = taHeight + x * taDihedral;
      GLfloat zTrailingEdge = (x - 1) / 4 + 9.5f;
      GLfloat zGuide = (x - 1) / 2 + 7.5f;
      GLfloat yGuide = taAspectRatio * (zGuide - zTrailingEdge);

      taPoints[p][0][0] = x;
      taPoints[p][0][1] = y;
      taPoints[p][0][2] = zTrailingEdge;

      taPoints[p][1][0] = x;
      taPoints[p][1][1] = p == 5 ? y : y + yGuide;
      taPoints[p][1][2] = zGuide;

      taPoints[p][2][0] = x;
      taPoints[p][2][1] = p == 5 ? y: y - yGuide;
      taPoints[p][2][2] = zGuide;

      taPoints[p][3][0] = x;
      taPoints[p][3][1] = y;
      taPoints[p][3][2] = zTrailingEdge;
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1, 3,           taWidth,
      0, 1, 3 * taWidth, taLength, 
      &taPoints[0][0][0]);
   glMapGrid2f(fuLoops, 0, 1, fuSlices, 0, 1);

   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);
   glScalef(-1, 1, 1);
   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);

   // Rudder
   const int ruWidth = 4;
   const int ruLength = 6;
   const GLfloat ruAspectRatio = 0.2f;
   GLfloat ruPoints[ruLength][ruWidth][3];

   for (p = 0; p < ruLength; p++)
   {
      GLfloat y = 0.7f * p;
      GLfloat zTrailingEdge = y / 3 + 9.5f;
      GLfloat zGuide = y + 6;
      GLfloat xGuide = ruAspectRatio * (zGuide - zTrailingEdge);

      ruPoints[p][0][0] = 0;
      ruPoints[p][0][1] = y;
      ruPoints[p][0][2] = zTrailingEdge;

      ruPoints[p][1][0] = p == 5 ? 0 : - xGuide;
      ruPoints[p][1][1] = y;
      ruPoints[p][1][2] = zGuide;

      ruPoints[p][2][0] = p == 5 ? 0 : xGuide;
      ruPoints[p][2][1] = y;
      ruPoints[p][2][2] = zGuide;

      ruPoints[p][3][0] = 0;
      ruPoints[p][3][1] = y;
      ruPoints[p][3][2] = zTrailingEdge;
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1, 3,           ruWidth,
      0, 1, 3 * ruWidth, ruLength, 
      &ruPoints[0][0][0]);
   glMapGrid2f(fuLoops, 0, 1, fuSlices, 0, 1);

   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);

   const GLfloat xMotor = 3;
   const GLfloat yMotor = -1.8f;
   const GLfloat zMotor = -3;

   // Motor support
   const int suWidth = 4;
   const int suLength = 4;
   const GLfloat suAspectRatio = 0.2f;
   GLfloat suPoints[ruLength][ruWidth][3];

   for (p = 0; p < suLength; p++)
   {
      GLfloat y = 0.2f * p + 0.7f;
      GLfloat zTrailingEdge = y / 2 + 2.5f;
      GLfloat zGuide = y / 2 + 0.5f;
      GLfloat xGuide = suAspectRatio * (zGuide - zTrailingEdge);

      suPoints[p][0][0] = 0;
      suPoints[p][0][1] = y;
      suPoints[p][0][2] = zTrailingEdge;

      suPoints[p][1][0] = p == 3 ? 0 : - xGuide;
      suPoints[p][1][1] = y;
      suPoints[p][1][2] = zGuide;

      suPoints[p][2][0] = p == 3 ? 0 : xGuide;
      suPoints[p][2][1] = y;
      suPoints[p][2][2] = zGuide;

      suPoints[p][3][0] = 0;
      suPoints[p][3][1] = y;
      suPoints[p][3][2] = zTrailingEdge;
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1, 3, suWidth,
      0, 1, 3 * suWidth, suLength, 
      &suPoints[0][0][0]);
   glMapGrid2f(fuLoops, 0, 1, fuSlices, 0, 1);

   glPushMatrix();
   glTranslatef(xMotor, yMotor, zMotor);
   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);
   glPopMatrix();

   glPushMatrix();
   glTranslatef(- xMotor, yMotor, zMotor);
   glEvalMesh2(GL_FILL, 0, fuLoops, 0, fuSlices);
   glPopMatrix();

   // Motor
   const double PI = 4 * atan(1.0);
   const int moLength = 5;
   const int moRound = 8;
   const GLfloat zVal[moLength] = { 1, 0, 1, 4, 3 };
   const GLfloat rad[moLength] = { 0.5, 0.5, 1.5, 0.5, 0.5 };
   GLfloat moPoints[moLength][moRound][3];

   for (int a = 0; a < moRound; a++)
   {
      const double theta = a * 2 * PI / (moRound - 1);
      const GLfloat sint = GLfloat(sin(theta));
      const GLfloat cost = GLfloat(cos(theta));
      for (int z = 0; z < moLength; z++)
      {
         moPoints[z][a][0] = - rad[z] * sint;
         moPoints[z][a][1] = rad[z] * cost;
         moPoints[z][a][2] = zVal[z];
      }
   }

   glMap2f(GL_MAP2_VERTEX_3, 
      0, 1, 3, moRound,
      0, 1, 3 * moRound, moLength, 
      &moPoints[0][0][0]);
   glMapGrid2f(20, 0, 1, 20, 0, 1);

   glPushMatrix();
   glTranslatef(xMotor, yMotor, zMotor);
   glEvalMesh2(GL_FILL, 0, 20, 0, 20);
   glPopMatrix();

   glPushMatrix();
   glTranslatef(- xMotor, yMotor, zMotor);
   glEvalMesh2(GL_FILL, 0, 20, 0, 20);
   glPopMatrix();

   GLUquadricObj *quad = gluNewQuadric();
   gluQuadricDrawStyle(quad, GLU_FILL);
   gluQuadricNormals(quad, GLU_SMOOTH);
   setMaterial(BLACK);

   glPushMatrix();
   glTranslatef(xMotor, yMotor + 0.1f, zMotor + 1);
   // Cover front of engine
   gluDisk(quad, 0, 0.45f, 10, 5);   
   glTranslatef(0, 0, 2);
   // Cover rear of engine
   gluDisk(quad, 0, 0.35f, 10, 5);   
   glPopMatrix();

   glPushMatrix();
   glTranslatef(- xMotor, yMotor + 0.1f, zMotor + 1);
   // Cover front of engine
   gluDisk(quad, 0, 0.45f, 10, 5);   
   glTranslatef(0, 0, 2);
   // Cover rear of engine
   gluDisk(quad, 0, 0.35f, 10, 5);   
   glPopMatrix();

   gluDeleteQuadric(quad);
}


// The following array describes material properties.
//    ambiance  [4],
//    diffuse   [4],
//    specular  [4],
//    shininess [1].
// 13 values are required to define a material.

// The order of the definitions must correspond with
// the enumeration MATERIAL in the header file.

const int MAX_MATERIALS = 100;

GLfloat materialProperties[][MAX_MATERIALS] =
{
   { // Black
        0, 0, 0, 1,
        0, 0, 0, 1,
        0, 0, 0, 1, 20
   },
   { // White
        1, 1, 1, 1,
        1, 1, 1, 1,
        1, 1, 1, 1, 75
   },
   { // Red
        1, 0, 0, 1,
        1, 0, 0, 1,
        1, 1, 1, 1, 75
   },
   { // Green
        0, 1, 0, 1,
        0, 1, 0, 1,
        1, 1, 1, 1, 75
   },
   { // Blue
        0, 0, 1, 1,
        0, 0, 1, 1,
        1, 1, 1, 1, 75
   },
   { // Metal
        0.2f, 0.2f, 0.2f, 1,
        0.65f, 0.65f, 0.65f, 1,
        1, 1, 1, 1, 75
   },
   { // Glass
        0.2f, 0.2f, 0.2f, 0.2f,
        0.2f, 0.2f, 0.2f, 0.2f,
        1, 1, 1, 1, 120
   },
   { // Brass
        0.329412f, 0.223529f, 0.027451f, 1.000000f,
        0.780392f, 0.568627f, 0.113725f, 1.000000f,
        0.992157f, 0.941176f, 0.807843f, 1.000000f,
        27.897400f 
   }, 
   { // Bronze
        0.212500f, 0.127500f, 0.054000f, 1.000000f,
        0.714000f, 0.428400f, 0.181440f, 1.000000f,
        0.393548f, 0.271906f, 0.166721f, 1.000000f,
        25.600000f
   },
   { // Polished Bronze
        0.250000f, 0.148000f, 0.064750f, 1.000000f,
        0.400000f, 0.236800f, 0.103600f, 1.000000f,
        0.774597f, 0.458561f, 0.200621f, 1.000000f,
        76.800003f
   },
   { // Chrome
        0.250000f, 0.250000f, 0.250000f, 1.000000f,
        0.400000f, 0.400000f, 0.400000f, 1.000000f,
        0.774597f, 0.774597f, 0.774597f, 1.000000f,
        76.800003f
   },
   { // Copper
        0.191250f, 0.073500f, 0.022500f, 1.000000f,
        0.703800f, 0.270480f, 0.082800f, 1.000000f,
        0.256777f, 0.137622f, 0.086014f, 1.000000f,
        12.800000f
   },
   { // Polished Copper
        0.229500f, 0.088250f, 0.027500f, 1.000000f,
        0.550800f, 0.211800f, 0.066000f, 1.000000f,
        0.580594f, 0.223257f, 0.069570f, 1.000000f,
        51.200001f
   },
   { // Gold
        0.247250f, 0.199500f, 0.074500f, 1.000000f,
        0.751640f, 0.606480f, 0.226480f, 1.000000f,
        0.628281f, 0.555802f, 0.366065f, 1.000000f,
        51.200001f
   },
   { // Polished Gold
        0.247250f, 0.224500f, 0.064500f, 1.000000f,
        0.346150f, 0.314300f, 0.090300f, 1.000000f,
        0.797357f, 0.723991f, 0.208006f, 1.000000f,
        83.199997f
   },
   { // Pewter
        0.105882f, 0.058824f, 0.113725f, 1.000000f,
        0.427451f, 0.470588f, 0.541176f, 1.000000f,
        0.333333f, 0.333333f, 0.521569f, 1.000000f,
        9.846150f
   },
   { // Silver
        0.192250f, 0.192250f, 0.192250f, 1.000000f,
        0.507540f, 0.507540f, 0.507540f, 1.000000f,
        0.508273f, 0.508273f, 0.508273f, 1.000000f,
        51.200001f
   },
   { // Polished Silver
        0.231250f, 0.231250f, 0.231250f, 1.000000f,
        0.277500f, 0.277500f, 0.277500f, 1.000000f,
        0.773911f, 0.773911f, 0.773911f, 1.000000f,
        89.599998f
   },
   { // Emerals
        0.021500f, 0.174500f, 0.021500f, 0.550000f,
        0.075680f, 0.614240f, 0.075680f, 0.550000f,
        0.633000f, 0.727811f, 0.633000f, 0.550000f,
        76.800003f
   },
   { // Jade
        0.135000f, 0.222500f, 0.157500f, 0.950000f,
        0.540000f, 0.890000f, 0.630000f, 0.950000f,
        0.316228f, 0.316228f, 0.316228f, 0.950000f,
        12.800000f
   },
   { // Obsidian
        0.053750f, 0.050000f, 0.066250f, 0.820000f,
        0.182750f, 0.170000f, 0.225250f, 0.820000f,
        0.332741f, 0.328634f, 0.346435f, 0.820000f,
        38.400002f
   },
   { // Pearl
        0.250000f, 0.207250f, 0.207250f, 0.922000f,
        1.000000f, 0.829000f, 0.829000f, 0.922000f,
        0.296648f, 0.296648f, 0.296648f, 0.922000f,
        11.264000f
   },
   { // Ruby
        0.174500f, 0.011750f, 0.011750f, 0.550000f,
        0.614240f, 0.041360f, 0.041360f, 0.550000f,
        0.727811f, 0.626959f, 0.626959f, 0.550000f,
        76.800003f
   },
   { // Turquoise
        0.100000f, 0.187250f, 0.174500f, 0.800000f,
        0.396000f, 0.741510f, 0.691020f, 0.800000f,
        0.297254f, 0.308290f, 0.306678f, 0.800000f,
        12.800000f
   },
   { // Black Plastic
        0.000000f, 0.000000f, 0.000000f, 1.000000f,
        0.010000f, 0.010000f, 0.010000f, 1.000000f,
        0.500000f, 0.500000f, 0.500000f, 1.000000f,
        32.000000f
   },
   { // Black Rubber
        0.020000f, 0.020000f, 0.020000f, 1.000000f,
        0.010000f, 0.010000f, 0.010000f, 1.000000f,
        0.400000f, 0.400000f, 0.400000f, 1.000000f,
        10.000000f
   },
};

void setMaterial(const int m, GLenum face)
{
  glMaterialfv(face, GL_AMBIENT,  &(materialProperties[m][0]));
  glMaterialfv(face, GL_DIFFUSE,  &(materialProperties[m][4]));
  glMaterialfv(face, GL_SPECULAR, &(materialProperties[m][8]));
  glMaterialf(face, GL_SHININESS, materialProperties[m][12]);
}

int matIndex = LAST_VALUE;

int addMaterial(
   GLfloat ambR, GLfloat ambG, GLfloat ambB, GLfloat ambA,
   GLfloat difR, GLfloat difG, GLfloat difB, GLfloat difA,
   GLfloat speR, GLfloat speG, GLfloat speB, GLfloat speA,
   GLfloat shine)
{
    if (matIndex < MAX_MATERIALS)
    {
        materialProperties[matIndex][ 0] = ambR;
        materialProperties[matIndex][ 1] = ambG;
        materialProperties[matIndex][ 2] = ambB;
        materialProperties[matIndex][ 3] = ambA;
        materialProperties[matIndex][ 4] = difR;
        materialProperties[matIndex][ 5] = difG;
        materialProperties[matIndex][ 6] = difB;
        materialProperties[matIndex][ 7] = difA;
        materialProperties[matIndex][ 8] = speR;
        materialProperties[matIndex][ 9] = speG;
        materialProperties[matIndex][10] = speB;
        materialProperties[matIndex][11] = speA;
        materialProperties[matIndex][12] = shine;
        return matIndex++;
    }
    else
    {
        cuglError = TOO_MANY_MATERIALS;
        return 0;
    }
}

int addMaterial(GLfloat params[])
{
    if (matIndex < MAX_MATERIALS)
    {
        for (int i = 0; i < 13; i++)
            materialProperties[matIndex][i] = params[i];
        return matIndex++;
    }
    else
    {
        cuglError = TOO_MANY_MATERIALS;
        return 0;
    }
}

};