cugl.h

Go to the documentation of this file.

#ifndef CUGL_H
#define CUGL_H
 
/** \file cugl.h
 *  Concordia University Graphics Library
 *
 * \authors Peter Grogono
*/
 
/** \mainpage   Introduction
 
   The \c operator<< overloads for classes Point, Quaternion, and Matrix
   use the current settings of format parameters.  \c setw() determines
   the width of each number rather than the entire field width.
 
   Most of the classes do not have copy constructors or \c operator= overloads
   because the default versions generated by the compiler do the right
   thing (for once).  (Class PixelMap is an exception, because it has
   pointer members.)
 
   \c GLfloat is used as the representation type for most floating-point values.
   This is for compatibility with OpenGL.  Although OpenGL provides a choice of
   precision for vertex and normal data, matrices are provided in \c GLfloat form only.
   Single-precision floating-point provides only about six decimal places, but this
   is sufficient for most graphics calculations.  If the number of points, vectors,
   or matrices used is large, substantial amounts of memory may be saved.
 
   Acknowledgements:  Mark Kilgard, Ken Shoemake, F.S. Hill, and others, for
   various publications on mthematical techniques for graphics and OpenGL, and to
   Liping Ye for improving various features, especially bitmaps.
*/
 
/** Extensions
 * \todo Quaternion += and -=
 * \todo Matrix operator*=(Matrix) - done
 * \todo bool isVisible(Point) - harder than it looks
 */
 
#ifdef __APPLE__
#include <OpenGL/gl.h>
#include <OpenGL/glu.h>
#include <GLUT/glut.h>
#else
#ifdef _WIN32
  #include <windows.h>
#endif
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#endif
 
#include <iostream>
#include <cmath>
 
namespace cugl
{
 
/** \defgroup constants Constants */
 
//@{
 
/** A string giving the version of CUGL and the most recent build date. */
const char version[] = "CUGL V2 2009.11.24";
 
/** A well-known mathematical constant. */
const double PI = 4 * atan(1.0);
 
// End of constants
//@}
 
/** The type of OpenGL projection and model view matrices.
 *  It is tempting to write CUGL as a template llibrary, but the fixed type
 *  of OpenGL matrices makes this difficult to do effectively. */
typedef GLfloat GL_Matrix[4][4];
 
/** \defgroup errors Error messages */
 
//@{
 
/** Enumeration values for CUGL function errors. */
enum CUGLErrorType
{
   NO_ERROR = 0,           /**< No error has occurred since the error flag was reset. */
   BAD_INDEX,              /**< Encountered an invalid array index. */
   BAD_MATRIX_MODE,        /**< The matrix mode was not GL_PROJECTION or GL_MODELVIEW. */
   BAD_ROTATION_MATRIX,    /**< The given matrix was not a rotation matrix. */
   SINGULAR_MATRIX,        /**< The matrix is singular. */
   ZERO_DIVISOR,           /**< An attempt was made to divide by zero. */
   ZERO_NORM,              /**< An attempt was made to normalize a null vector. */
   BAD_INTERPOLATOR_ARG,   /**< The angle between the rotations of the interpolator is zero.  */
   OPEN_FAILED,            /**< The program failed to open a BMP file. */
   NOT_BMP_FILE,           /**< The file opened was not a BMP file. */
   NOT_24_BITS,            /**< The BMP file was not in 24-bit format. */
   COMPRESSED_BMP_FILE,    /**< The BMP file was in compressed form. */
   NOT_ENOUGH_MEMORY,      /**< There is not enough memory to store the pixel map. */
   NO_PIX_MAP,             /**< There is no pixel map to draw. */
   BAD_PLANE,              /**< Parameters do not define a plane. */
   BAD_LINE,               /**< The line lies in the plane. */
   TOO_MANY_MATERIALS,     /**< Too many materials. */
   TOO_MANY_POINTS         /**< Too many points. */
};
 
/**
 * Set code to \c NO_ERROR and return last error code.
 * \return the most recent error code.
 */
CUGLErrorType getError();
 
/** Map error codes to descriptive strings.
  * \param err is an error code obtained from \c getError().
  * \return A string describing the eeror.
 */
const char *getErrorString(CUGLErrorType err);
 
/**
 * Display a message if an CUGL error has occurred.
 * This function checks the CUGL error status.
 * If no error has occurred, it has no effect.
 * If CUGL has reported an error, this function writes
 * the CUGL error code and the corresponding message
 * to the stream \c cerr.
 *
 * Since \c getError() clears the CUGL error flag,
 * after calling this function CUGL shows no errors.
 *
 * This function is intended mainly for development and
 * should not normally be included in a production program.
 */
void checkCUGLStatus();
 
/**
 * Display a message if an OpenGL error has occurred.
 * This function checks the OpenGL error status.
 * If no error has occurred, it has no effect.
 * If OpenGL has reported an error, this function writes
 * the OpenGL error code and the corresponding message
 * to the stream \c cerr.
 *
 * Since \c glGetError() clears the OpenGL error flag,
 * after calling this function OpenGL shows no errors.
 *
 * This function is intended mainly for development and
 * should not normally be included in a production program.
 */
void checkOpenGLStatus();
 
// End of error messages
//@}
 
 
 
// Forward declarations.
class Line;
class Plane;
class Vector;
class Quaternion;
 
 
 
/**
 * An instance is a point represented by four homogeneous coordinates.
 * A point is represented by \a (x,y,z,w) in which each component is a \c GLfloat.
 * If \a w=0, the point is `at infinity'.  Points at infinity may be used like other
 * points, but a few operations do not work for them.  CUGL functions do not issue an
 * error message when this happens: they simply return a default result.
 *
 * A Point is in normal form if \a w=1.  A Point may be normalized if \a w!=0.
 */
class Point
{
public:
   friend class Line;
   friend class Vector;
   friend class Matrix;
   friend class Plane;
 
   /**
    * Construct a point with coordinates (x, y, z, w).
    * Default parameter values: (0,0,0,1).
    */
   Point(GLfloat x = 0, GLfloat y = 0, GLfloat z = 0, GLfloat w = 1);
 
   /**
    * Construct a Point from coordinates.
    * \param coordinates is an array containing at least four coordinates.
    * \pre The array must have at least four components.
    */
   Point (GLfloat coordinates[]);
 
   /**
    * Convert a Quaternion to a Point.
    * \note This is a rather peculiar operation and should not normally be used.
    * It is intended for experiments with non-linear transformations.
    * \param q is a quaternion.
    */
   Point (const Quaternion & q);
 
   /**
    * Access coordinates of a point.
    * \param i is the index of the coordinate: p[0] = x, p[1] = y, p[2] = z, p[3] = w.
    * \note Error BAD_INDEX reported if \c i is not one of {0,1,2, 3}.
    */
   GLfloat & operator[](int i);
 
   /**
    * Access coordinates of a point.
    * This form is used with \c const instances (\c p[i] cannot appear on the left of \c =).
    * \param i is the index of the coordinate: p[0] = x, p[1] = y, p[2] = z, p[3] = w.
    * \note Error BAD_INDEX reported if \c i is not one of {0,1,2, 3}.
    */
   const GLfloat & operator[](int i) const;
 
   /**
    * Normalize a point.
    * A Point \a (x,y,z,w) is in normal form if \a w=1.
    * \note Error ZERO_DIVISOR reported if \a w=0.
    * \note The value of this Point is changed by this function.
    */
   void normalize();
 
   /**
    * Return the normalized form of a Point.
    * A Point \a (x,y,z,w) is in normal form if \a w=1.
    * \note Error ZERO_DIVISOR reported if \a w=0.
    */
   Point unit() const;
 
   /**
    * Draw the point using \c glVertex4f().
    */
   void draw() const;
 
   /**
    * Use this point as a light position.
    * If w = 0, the light is directional, otherwise it is positional.
    * \param lightNum specifies which light to use: \a GL_LIGHT0, \a GL_LIGHT1, ...
    */
   void light(GLenum lightNum) const;
 
   /**
    * Translate to the point using \c glTranslatef().
    * \note If the point is at infinity (w = 0), this function has no effect.
    */
   void translate() const;
 
   /**
    * Scale a Point.
    * Scale the coordinates of a Point by dividing by the scalar \a s.
    * This is implemented by multiplying the \a w component of the Point
    * by \a s.  If \a s=0, the result is a Point at infinity.
    * \return the Point at \a (x,y,z,s*w).
    */
   Point operator/(GLfloat s) const;
 
   /**
    * Displace a Point with a Vector.
    * The components of the Vector are added to the corresponding components of the Point.
    * If the Point is not in normal form, the Vector is implicitly scaled.
    * \return Point at (x-w*v.x, y-w*v.y, z-w*v.z, w).
    */
   Point operator+=(const Vector & v);
 
   /**
    * Displace a Point with a Vector
    * The components of the Vector are subtracted from the corresponding components of the Point.
    * If the Point is not in normal form, the Vector is implicitly scaled.
    * \return Point at (x+w*v.x, y+w*v.y, z+w*v.z, w).
    */
   Point operator-=(const Vector & v);
 
   /**
    * Displace a Point with a Vector.
    * The components of the Vector are added to the corresponding components of the Point.
    * If the Point is not in normal form, the Vector is implicitly scaled.
    * \return Point at (p.x+p.w*v.x, p.y+p.w*v.y, p.z+p.w*v.z, p.w).
    */
   friend Point operator+(const Vector & v, const Point & p);
 
   /**
    * Displace a Point with a Vector.
    * The components of the Vector are added to the corresponding components of the Point.
    * If the Point is not in normal form, the Vector is implicitly scaled.
    * \return Point at (p.x+p.w*v.x, p.y+p.w*v.y, p.z+p.w*v.z, p.w).
    */
   friend Point operator+(const Point & p, const Vector & v);
 
   /**
    * Scale a point.
    * \return Point at (p.x + p.w*v.x, p.y + p.w*v.y, p.z + p.w*v.z).
    */
   friend Point operator*(const Point & p, GLfloat s);
 
   /**
    * Scale a point.
    * Multiply each of the components of the Point by s.
    * In homogeneous coordinates, this operation does not change the ``position'' of the Point.
    * However, it may be used, for example, to normalize a Point.
    * \return The scaled Point.
    */
   friend Point operator*(GLfloat s, const Point & p);
 
   /**
    * Return the vector corresponding to the displacement between the two points.
    * This is the vector (\c p.x/p.w - \c q.x/q.w, \c p.y/p.w - \c q.y/q.w, \c p.z/p.w - \c q.z/q.w).
    * If \c p or \c q is a point at infinity, return the zero vector.
    * \return the Vector corresponding to the displacement between the two points.
    */
   friend Vector operator-(const Point & p, const Point & q);
 
   /**
    * Find the point where this line meets the plane p.
    * Sets error BAD_LINE if the line lies within the plane.
    * \return the Point at which Line \c k meets Plane \c p.
    */
   friend Point meet(const Line & k, const Plane & p);
 
   /**
    * Compare two points.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Point & p, const Point & q);
 
   /**
    * Compare two points.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Point & p, const Point & q);
 
   /**
    * Write "(x,y,z,w)" to output stream.
    * Inserts the coordinates of the point into the output stream \c os
    * in the format "(x,y,z,w)".  The current settings of the stream formatting
    * parameters are used.  If a width is specified with \c setw(), it is applied
    * to each coordinate, not to the point as a whole.
    */
   friend std::ostream & operator<<(std::ostream & os, const Point & p);
 
   /**
    * Return the distance between a Point and a Plane.
    * The result has the correct magnitude only if the Point and the Plane
    * are both in normal form.  In particular, the result is incorrect if
    * the Point is at infinity.  However, the sign of the result is correct
    * in all cases, and so it is not necessary to provide normalized
    * arguments if only the sign is important.
    */
   friend GLfloat dist(const Point & p, const Plane & s);
 
   /**
    * Return the distance between a Point and a Plane.
    * The result has the correct magnitude only if the Point and the Plane
    * are both in normal form.  In particular, the result is incorrect if
    * the Point is at infinity.  However, the sign of the result is correct
    * in all cases, and so it is not necessary to provide normalized
    * arguments if only the sign is important.
    */
   friend GLfloat dist(const Plane & s, const Point & p);
 
   /**
    * Return the disgtance bewteen two points.
    * \note Does not check that the points are normalized.
    */
   friend double dist(const Point & p, const Point & q)
   {
      return sqrt((p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y) + (p.z - q.z) * (p.z - q.z));
   }
 
private:
 
   /** X coordinate of point. */
   GLfloat x;
 
   /** Y coordinate of point. */
   GLfloat y;
 
   /** Z coordinate of point. */
   GLfloat z;
 
   /** W coordinate of point. */
   GLfloat w;
 
};
 
 
/**
 * An instance is a line defined by two Points.
 * Lines are occasionally useful.  For example, it is simple to
 * display surface normals using this class.
 * A Line is represented by the two Points which are its end-points.
 */
 
class Line
{
public:
 
   friend class Plane;
 
   /**
    * Construct the Line from Point \c p to Point \c q.
    */
   Line(const Point & p, const Point & q);
 
   /**
    * Construct the Line from Point \c p to Point \c p+v.
    */
   Line(const Point & p, const Vector & v);
 
   /**
    * Find the point where this line meets the plane p.
    */
   friend Point meet(const Line & k, const Plane & p);
 
   /**
    * Draw the line using \c glBegin(GL_LINE) ....
    */
   void draw() const;
 
   /**
    * Compare two lines.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Line & x, const Line & y);
 
   /**
    * Compare two lines.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Line & x, const Line & y);
 
   /**
    * Write a representation of the line to the output stream.
    * The format is \a pt->pt where "\a pt" is the format used for points.
    * If \c setw is used to set the width, it is passed to the inserter
    * for Point.
    */
   friend std::ostream & operator<<(std::ostream & os, const Line & k);
 
private:
   /** Start point of line. */
   Point s;
 
   /** Finish point of line. */
   Point f; // finish point
};
 
 
/**
 * An instance is a plane defined by its equation \a Ax+By+Cz+D=0.
 *
 * Planes are used for clipping, shadows, and reflections
 * (see class Matrix).
 *
 * Homogeneous points (x,y,z,w) lie on the plane if Ax+By+Cz+Dw=0.
 * The notation for a plane is \a <A,B,C,D>.  The plane \a <0,0,0,D> is undefined.
 * An attempt to construct such a plane sets the error flag to BAD_PLANE and the
 * plane is set to <0,1,0,0> (in the conventional OpenGL frame, this often corresponds
 * to the ground, \a y=0).
 *
 * A Plane is in normal form if the Vector \a (A,B,C) is a unit vector.
 * A Plane can be normalized by scaling \a A,B,C,D.
 */
 
class Plane
{
   friend class Matrix;
 
public:
   /**
    * Construct a plane given the coefficients of its equation \a Ax+By+Cx+D=0.
    * If no arguments are provided, construct the plane y = 0.
    */
   Plane(GLfloat a = 0, GLfloat b = 1, GLfloat c = 0, GLfloat d = 0);
 
   /**
    * Construct a plane containing the three given points.
    */
   Plane(const Point & p, const Point & q, const Point & r);
 
   /**
    * Construct a plane containing the line \c s and the point \c p.
    */
   Plane(const Line & s, const Point & p);
 
   /**
    * Construct the plane orthogonal to \c v and passing through \c p.
    */
   Plane(const Vector & v, const Point & p);
 
   /**
    * Normalize this plane.
    * The Plane \a (A,B,C,D) is in normal form when \a (A,B,C) is a unit vector.
    * Error ZERO_DIVISOR reported if \a |(A,B,C|=0 (which should not be the case for a well-formed plane).
    * \note The value of the Plane is changed by this function.
    */
   void normalize();
 
   /**
    * Return a normalized Plane equivalent to this plane.
    * The Plane \a (A,B,C,D) is in normal form when \a (A,B,C) is a unit vector.
    * Error ZERO_DIVISOR reported if \a |(A,B,C|=0 (which should not be the case for a well-formed plane).
    */
   Plane unit() const;
 
   /**
    * Return a vector of unit length that is normal to the plane.
    */
   Vector normal() const;
 
   /**
    * Use this plane as a clipping plane.
    * This function calls \c glClipPlane
    * to suppress rendering of objects on one side of the plane.
    * \param index must be one of GL_CLIP_PLANE0, GL_CLIP_PLANE1, ...  A
    * An implementation of OpenGL is supposed to provide at least six
    * clipping planes, numbered 0,1,...,5.
    */
   void clipPlane(GLenum index) const;
 
   /**
    * Find the point where this line meets the plane p.
    */
   friend Point meet(const Line & k, const Plane & p);
 
   /**
    * Compare two planes.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Plane & x, const Plane & y);
 
   /**
    * Compare two planes.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Plane & x, const Plane & y);
 
   /**
    * Write a description of the plane to the output stream as \a <a,b,c,d>.
    * Inserts the components of the plane into the output
    * stream \c os in the format \a <a,b,c,d>.  The current settings
    * of the stream formatting parameters are used.  If a width
    * is specified with \c setw(), it is applied to each component
    * rather than to the plane as a whole.
    */
   friend std::ostream & operator<<(std::ostream & os, const Plane & p);
 
   /**
    * Return the distance between a Point and a Plane.
    * The result has the correct magnitude only if the Point and the Plane
    * are both in normal form.  In particular, the result is incorrect if
    * the Point is at infinity.  However, the sign of the result is correct
    * in all cases, and so it is not necessary to provide normalized
    * arguments if only the sign is important.
    */
   friend GLfloat dist(const Point & p, const Plane & s);
 
   /**
    * Return the distance between a Point and a Plane.
    * The result has the correct magnitude only if the Point and the Plane
    * are both in normal form.  In particular, the result is incorrect if
    * the Point is at infinity.  However, the sign of the result is correct
    * in all cases, and so it is not necessary to provide normalized
    * arguments if only the sign is important.
    */
   friend GLfloat dist(const Plane & s, const Point & p);
 
   /** Return the A component of the plane defined by Ax+By+Cz+d=0. */
   GLfloat getA() { return a; }
 
   /** Return the B component of the plane defined by Ax+By+Cz+d=0. */
   GLfloat getB() { return b; }
 
   /** Return the C component of the plane defined by Ax+By+Cz+d=0. */
   GLfloat getC() { return c; }
 
   /** Return the D component of the plane defined by Ax+By+Cz+d=0. */
   GLfloat getD() { return d; }
 
 
private:
   /** Coefficient of \a x in the plane equation. */
   GLfloat a;
 
   /** Coefficient of \a y in the plane equation. */
   GLfloat b;
 
   /** Coefficient of \a z in the plane equation. */
   GLfloat c;
 
   /** Coefficient of \a w in the plane equation. */
   GLfloat d;
};
 
class Matrix;
 
/**
 * An instance is a vector with 3 real components.
 * A vector v is a 3-vector represented by three orthogonal components
 * \a vx, \a vy, and \a vz.  The norm of the vector \a v is vx*vx+vy*vy+vz*vz.
 * The length of \c v is \c sqrt(norm(v)).
 */
class Vector
{
   friend class Point;
   friend class Matrix;
   friend class Quaternion;
 
public:
   /**
    * Construct the zero vector: (0,0,0).
    */
   Vector() : x(0), y(0), z(0)
   {}
 
   /**
    * Construct the vector (x,y,z).
    */
   Vector(GLfloat x, GLfloat y, GLfloat z)
         : x(x), y(y), z(z)
   {}
 
   /**
    * Construct a vector from the array \c coordinates.
    * \pre The array must have at least three components.
    */
   Vector(GLfloat coordinates[]);
 
   /**
    * Construct a vector normal to the polygon defined
    * by the given points using Martin Newell's algorithm.
    * The normal vector will be exact if the points lie in a plane,
    * otherwise it will be a sort of average value.  As with OpenGL,
    * the vector will point in the direction from which the points
    * are enumerated in a counter-clockwise direction.
    *
    * Unlike other functions, this function does \b not use homogeneous coordinates.
    * The points are assumed to have (x,y,z) coordinates; the w component is ignored.
    * \param points is an array of points.
    * \param numPoints is the number of points in the array.
    * \return the vector normal to the plane defined by \a points.
    * \note The vector is \b not a unit vector because it will probably
    * be averaged with other vectors.
    */
   Vector(Point points[], int numPoints);
 
   /**
    * Construct a vector from two points.
    * Vector(p, q) is equivalent to p - q.
    */
   Vector(const Point & p, const Point & q);
 
   /**
    * Construct a vector from a quaternion by ignoring the scalar component of the quaternion.
    * Vector(q) constructs the vector returned by \c q.vector().
    * \param q is the \a Quaternion whose vector component is to be used.
    */
   explicit Vector(const Quaternion & q);
 
   /**
    * Add vector \c v to this vector.
    */
   Vector operator+=(const Vector & v);
 
   /**
    * Subtract vector \c v from this vector.
    */
   Vector operator-=(const Vector & v);
 
   /**
    * Return Vector (-x,-y,-z).
    */
   Vector operator-() const;
 
   /**
    * Multiply each component of the vector by \c scale.
    */
   Vector operator*=(GLfloat scale);
 
   /**
    * Return the vector (x/scale,y/scale,z/scale).
    * Error ZERO_DIVISOR reported if \c scale is zero.
    * \return the Vector (x/scale,y/scale,z/scale).
    */
   Vector operator/(GLfloat scale) const;
 
   /**
    * Divide each component of this vector by \c scale and return the new value.
    * Error ZERO_DIVISOR reported if \c scale is zero.
    */
   Vector operator/=(GLfloat scale);
 
   /**
    * Normalize this vector.
    * See also \a Vector::unit().
    * Reports error ZERO_NORM if the vector is zero.
    * \note  The value of this vector is changed by this operation.
    */
   void normalize();
 
   /**
    * Return a unit vector with the same direction as this vector.
    * See also \a Vector::normalize().
    * Reports error ZERO_NORM if the vector is zero.
    * \note The value of this vector is not changed by this operation.
    */
   Vector unit() const;
 
   /**
    * Return the norm of this vector.
    * The norm of a vector is the sum of its squared components
    * and is also the square of the length.
    * \return the norm of this vector.
    */
   GLfloat norm() const;
 
   /**
    * Return the length of this vector.
    * The length of a vector is the square root of its norm.
    * \return the length of this vector.
    */
   GLfloat length() const;
 
   /**
    * Use this vector to translate an object.
    * This is implemented by passing the vector to \c glTranslate().
    */
   void translate() const;
 
   /**
    * Use this vector as a normal vector when drawing a surface.
    * This is implemented by passing the vector to \c glNormal().
    */
   void drawNormal() const;
 
   /**
    * Construct the skew-symmetric matrix corresponding to this vector.
    */
   Matrix skew();
 
   /**
    * Draw this vector as a line in the graphics window.
    * The line joins \a P and \a P+V, where \a P is the point provided,
    * or the origin if no point is given.
    * \param p is the start point for the vector.
    */
   void draw(const Point & p = Point()) const;
 
   /**
    * Get or set a reference to a component of this vector.
    * \c v[0] is the \a x component; \c v[1] is the \a y component; \c v[2] is the \a z component.
    * Error BAD_INDEX reported if \c i is not one of 0, 1, 2.
    */
   GLfloat & operator[](int i);
 
   /**
    * Get a reference to a component of this vector.
    * This form is used for \c const instances (\c v[i] cannot appear on the left of \c =).
    * \c v[0] is the \a x component; \c v[1] is the \a y component; \c v[2] is the \a z component.
    * Error BAD_INDEX reported if \c i is not one of 0, 1, 2.
    */
   const GLfloat & operator[](int i) const;
 
   /**
    * Return cross product of vectors \c u and \c v.
    */
   friend Vector cross(const Vector & u, const Vector & v);
 
   /**
    * Return cross product of vectors \c u and \c v.
    */
   friend Vector operator*(const Vector & u, const Vector & v);
 
   /**
    * Return dot product of vectors \c u and \c v.
    */
   friend GLfloat dot(const Vector & u, const Vector & v);
 
   /**
    * Return Vector \a s*v.
    */
   friend Vector operator*(const Vector & v, GLfloat s);
 
   /**
    * Return Vector \a s*v.
    */
   friend Vector operator*(GLfloat s, const Vector & v);
 
   /**
    * Return Vector \a u+v.
    */
   friend Vector operator+(const Vector & u, const Vector & v);
 
   /**
    * Return Vector \a u-v.
    */
   friend Vector operator-(const Vector & u, const Vector & v);
 
   /**
    * Displace a Point with a Vector.
    * \return Point at (p.x+p.w*v.x, p.y+p.w*v.y, p.z+p.w*v.z, p.w).
    */
   friend Point operator+(const Vector & v, const Point & p);
 
   /**
    * Return Point \c p displaced by vector \c v.
    */
   friend Point operator+(const Point & p, const Vector & v);
 
   /**
    * Compare two vectors.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Vector & x, const Vector & y);
 
   /**
    * Compare two vectors.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Vector & x, const Vector & y);
 
   /**
    * Write vector to output stream as \a (x,y,z).
    * Inserts the components of the vector into the output
    * stream \c os in the format \a (x,y,z).  The current settings
    * of the stream formatting parameters are used.  If a width
    * is specified with \c setw(), it is applied to each coordinate,
    * not to the vector as a whole.
    */
   friend std::ostream & operator<<(std::ostream & os, const Vector & v);
 
private:
 
   /** X component of vector. */
   GLfloat x;
 
   /** Y component of vector. */
   GLfloat y;
 
   /** Z component of vector. */
   GLfloat z;
};
 
/**
 * Unit vector parallel to \a X axis.
 */
const Vector I(1, 0, 0);
 
/**
 * Unit vector parallel to \a Y axis.
 */
const Vector J(0, 1, 0);
 
 
/**
 * Unit vector parallel to \a Z axis.
 */
const Vector K(0, 0, 1);
 
/**
 * An instance is a matrix compatible with an OpenGL transformation matrix.
 * An instance of class Matrix is a 4 by 4 matrix with
 * components of type \c GLfloat.  The components are ordered
 * in the same way as an OpenGL matrix (column-row order,
 * for compatibility with FORTRAN).
 
 * Note that OpenGL performs matrix calculations very efficiently.
 * As far as possible, construct transformations using sequences of
 * OpenGL matrix operations.  Construct and use the matrices here
 * only if OpenGL does not provide the required facilities.
 */
class Matrix
{
public:
 
   /**
    * Construct the identity matrix.
    */
   Matrix();
 
   /**
    * Construct a copy of an arbitrary OpenGL matrix.
    */
   Matrix(GL_Matrix r);
 
   /**
    * Construct a copy of an OpenGL projection or model view matrix.
    \param mode should be \c GL_PROJECTION_MATRIX or \c GL_MODELVIEW_MATRIX.
 
    Report error BAD_MATRIX_MODE and construct identity matrix
    for other values of mode.
    */
   explicit Matrix(GLenum mode);
 
   /**
    * Construct a rotation matrix.
    \param axis is the axis of rotation.
    \param theta is the magnitude of the rotation.
    \note The angle \c theta should be in radians (unlike OpenGL, which uses degrees).
    */
   Matrix(const Vector & axis, double theta);
 
   /**
    * Construct a matrix that reflects a point in the given plane.
    * The matrix can be used to simulate a mirror in an OpenGL program.
    * See also Matrix::reflect().
    \param refl is the plane of reflection.
    */
   explicit Matrix(const Plane & refl);
 
   /**
    * Construct a shadow matrix from a point light source and a plane.
    * The matrix transforms an object into its shadow on the given plane.
    * It is your job to ensure that the shadow has the appropriate colour, transparency, etc.
    * The light source may be a local point (w > 0) or a point at infinity (w = 0).
    * See also \c Matrix::shadow().
      \param lightPos is the position of the light source.
      \param plane is the plane onto which the shadow is projected.
 
    */
   Matrix(const Point & lightPos, const Plane & plane);
 
   /**
    * Construct a rotation matrix from a quaternion.
    * \pre The quaternion must be a unit quaternion.
    */
   explicit Matrix(const Quaternion & q);
 
   /**
    * Construct the matrix that rotates one vector to another.
    * \param u is a vector representing an initial orientation.
    * \param v is a vector representing the final orientation.
    * The matrix, applied to \a u, will yield \a v.
    * \pre The vectors \a u and \a v must be unit vectors.
    */
   Matrix(const Vector & u, const Vector & v);
 
   /**
    * Construct a matrix from 16 floats.
    */
   Matrix(
      const GLfloat m00, const GLfloat m10, const GLfloat m20, const GLfloat m30,
      const GLfloat m01, const GLfloat m11, const GLfloat m21, const GLfloat m31,
      const GLfloat m02, const GLfloat m12, const GLfloat m22, const GLfloat m32,
      const GLfloat m03, const GLfloat m13, const GLfloat m23, const GLfloat m33 );
 
   /**
    * Return the quaternion corresponding to this matrix.
    * This function may report BAD_ROTATION_MATRIX.
    * \note The result may be imprecise if the rotation angle is close to 180 degrees.
    * \pre The matrix must be a rotation matrix.
    */
   Quaternion quaternion() const;
 
   /**
    * Return a pointer to the first element of the matrix.
    * This function may be used in conjunction with
    * \c glMultMatrix() to apply this matrix to the current OpenGL matrix.
    * \return a pointer to the first element of the matrix.
    */
   GLfloat *get()
   ;
 
   /**
    * Return the transpose of this matrix.
    */
   Matrix transpose() const;
 
   /**
    * Return the trace of this matrix.
    * The trace will include the bottom right element, m[3][3].
    */
   GLfloat trace()
   {
      return m[0][0] + m[1][1] + m[2][2] + m[3][3];
   }
 
   /**
    * Compute the inverse of this matrix using Gauss-Jordan elimination.
    * If the matrix is singular, report \c SINGULAR_MATRIX and return the zero matrix.
    * Calculations are performed in double precision because this gives
    * slightly better results.
    * The prefix operator ~ has the same effect.
    * \return the inverse of this matrix.
    */
   Matrix inv() const;
 
   /**
    * Return the inverse of this matrix.
    * This provides an alternative syntax for inv().
    * \return the inverse of this matrix.
    */
   Matrix operator~() const;
 
   /**
    * Multiply the current OpenGL matrix by this matrix.
    */
   void apply() const;
 
   /**
    * Apply this matrix to a point.
    */
   Point apply(const Point & p) const;
 
   /**
    * Apply this matrix to a vector.
    */
   Vector apply(const Vector & v) const;
 
   /**
    * Return the axis of a rotation matrix.
    * Report error \c BAD_ROTATION_MATRIX and leave result undefined,
    * if the matrix is not a rotation.
    * \return the axis of a rotation matrix.
    */
   Vector axis() const;
 
   /**
    * Return the angle (in radians) of a rotation matrix.
    * Report error \c BAD_ROTATION_MATRIX and leave result undefined,
    * if the matrix is not a rotation.
    * \return the angle (in radians) of a rotation matrix.
    */
   double angle() const;
 
   /**
    * Return a reference to the element \c m[i][j] of the matrix.
    * \pre The values of \c i and \c j must be in the range [0,3].
    * \return the element \c m[i][j] of the matrix.
    */
   GLfloat & operator()(int i, int j);
 
   /**
    * Return a reference to the element \c m[i][j] of the \c const matrix.
    * This version is used for \c const instances (\c m(i,j) cannot appear on the left of \c =).
    * \pre The values of \c i and \c j must be in the range [0,3].
    * \return the element \c m[i][j] of the matrix.
    */
   const GLfloat & operator()(int i, int j) const;
 
   /**
    * Set the matrix so that it reflects a point in the plane \c p.
    * The matrix can be used to simulate a mirror in an OpenGL program.
    * See also constructors for class Matrix.
      \param p is the plane of reflection.
    */
   void reflect(const Plane & p);
 
   /**
    * Set the matrix so that it creates a shadow on the plane \c p from a light source \c lightPos.
    * The matrix transforms an object into its shadow on the given plane.
    * It is your job to ensure that the shadow has the appropriate colour, transparency, etc.
    * The point \c p may be either a local point (w > 0) or a point at infinity (w = 0).
    * See also constructors for class Matrix.
    \param lightPos is the position of the light source causing the shadow.
    \param plane is the plane upon which the shadow is cast.
    */
   void shadow(const Point & lightPos, const Plane & plane);
 
   /** Add and assign matrices. */
   Matrix operator+=(const Matrix & rhs);
 
   /** Add two matrices. */
   friend Matrix operator+(const Matrix & m, const Matrix & n);
 
   /** Subtract and assign matrices. */
   Matrix operator-=(const Matrix & rhs);
 
   /** Subtract two matrices. */
   friend Matrix operator-(const Matrix & m, const Matrix & n);
 
   /** Multiply and assign matrices. */
   Matrix operator*=(const Matrix & rhs);
 
   /** Multiply two matrices. */
   friend Matrix operator*(const Matrix & m, const Matrix & n);
 
   /** Multiply scalar and matrix. */
   friend Matrix operator*(GLfloat s, const Matrix & m);
 
   /** Multiply matrix and scalar. */
   friend Matrix operator*(const Matrix & m, GLfloat s);
 
   /** Multiply by scalar and assign. */
   Matrix & operator*=(GLfloat s);
 
   /** Divide by scalar and assign. */
   Matrix & operator/=(GLfloat s);
 
   /** Divide by scalar. */
   friend Matrix operator/(const Matrix & m, GLfloat s);
 
   /**
    * Compare two matrices.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Matrix & x, const Matrix & y);
 
   /**
    * Compare two matrices.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Matrix & x, const Matrix & y);
 
   /**
    * Write a four-line image of the matrix to the output stream.
    * The current settings of the stream formatting parameters are used.
    * If a width is specified with \c setw(), it is applied to each
    * element of the matrix, not the matrix as a whole.
    */
   friend std::ostream & operator<<(std::ostream & os, const Matrix & m);
 
private:
 
   /** The elements of the matrix. */
   GL_Matrix m;
};
 
 
 
/**
 * An instance is a quaternion represented as a (scalar, vector) pair.
 * We use the notation \a (s,v) for a quaternion with scalar component
 *   \a s and vector component \a v.  Although arbitrary quaternions can be
 *   constructed, all of the applications of quaternions provided by
 *   the class assume that the quaternion is a unit quaternion
 *   representing a rotation.
 */
class Quaternion
{
public:
 
   friend class Vector;
   friend class Matrix;
 
   /**
    * Construct the quaternion (1,(0,0,0)) (the null rotation).
    */
   Quaternion() : s(1)
   {}
 
   /**
    * Construct the quaternion (s, (x,y,z)).
    */
   Quaternion(GLfloat s, GLfloat x, GLfloat y, GLfloat z)
         : s(s), v(Vector(x,y,z))
   {}
 
   /**
    * Construct the quaternion \a (0,v).
    * The result will be a unit quaternion if \a v is a unit vector,
    * in which case the quaternion represents a rotation through
    * 90 degrees about the axis \a v.
    */
   Quaternion(const Vector & v) : s(0), v(v)
   {}
 
   /**
    * Construct the quaternion \a (s,v).
    * \note Don't confuse this constructor with Quaternion(axis,angle).
    * \param s is the scalar component of the quaternion.
    * \param v is the vector component of the quaternion.
    */
   Quaternion(GLfloat s, const Vector & v) : s(s), v(v)
   {}
 
   /**
    * Construct a quaternion with the given axis of rotation and angle.
    * \note Don't confuse this constructor with Quaternion(s,v).
    * \pre The axis must be a unit vector.
    * \param axis gives the axis of rotation.
    * \param angle gives the amount of the rotation.
    */
   Quaternion(Vector axis, double angle)
         : s(GLfloat(cos(angle/2))), v(GLfloat(sin(angle/2)) * axis)
   {}
 
   /**
    * Construct the quaternion corresponding to an OpenGL rotation matrix.
    * The result may be imprecise if the rotation angle is close to 180 degrees.
    * This function may report BAD_ROTATION_MATRIX.
    * \pre The matrix must be a rotation matrix.
    */
   Quaternion(Matrix m);
 
   /**
    * Construct a quaternion from Euler angles.
    */
   Quaternion(double xr, double yr, double zr);
 
   /**
    * Construct a Quaternion from a Point.
    * The quaternion consists of the vector \a (p[0],p[1],p[2]) and the scalar \a p[3].
    * \note This is a unusual operation and should not normally be used.
    * It is intended for experiments with non-linear transformations.
    */
   Quaternion(const Point & p);
 
   /**
    * Construct the quaternion that rotates one vector to another.
    * \param u is a vector representing an initial orientation.
    * \param v is a vector representing the final orientation.
    * The quaternion, applied to \a u, will yield \a v.
    * \pre The vectors \a u and \a v must be unit vectors.
    */
   Quaternion(const Vector & u, const Vector & v);
 
   /**
    * Add the quaternion q to this quaternion.
    */
   Quaternion operator+=(const Quaternion & q);
 
   /**
    * Subtract the quaternion q from this quaternion.
    */
   Quaternion operator-=(const Quaternion & q);
 
   /**
    * Return the quaternion \a q+r.
    * \note The sum of two unit quaternions is not in general a unit quaternion.
    * However, it occasionally appears in expressions such as k q1 + (1 - k) q2,
    * which does yield a unit quaternion.
    * \return the Quaternion \a q+r.
    */
   friend Quaternion operator+(const Quaternion & q, const Quaternion & r);
 
   /**
    * Return the quaternion \a q-r.
    * \note The difference of two unit quaternions is not in general a unit quaternion.
    * \return the Quaternion \a q-r.
    */
   friend Quaternion operator-(const Quaternion & q, const Quaternion & r);
 
   /**
    * Return the quaternion product \a q*r.
    * If \c q and \c r represent
    * rotations, then \a q*r represents the rotation \a q followed by the rotation \a r.
    * \note Quaternion multiplication is not commutative: \a q*r is not equal to \a r*q.
    * \return the Quaternion product \a q*r.
    */
   friend Quaternion operator*(const Quaternion & q, const Quaternion & r);
 
   /**
    * Promote the vector \a v to a quaternion \a qv and
    * return the quaternion product \a qv*q.
    */
   friend Quaternion operator*(const Vector & v, const Quaternion & q);
 
   /**
    * Promote the vector \a v to a quaternion \a qv and
    * return the quaternion product \a q*qv.
    */
   friend Quaternion operator*(const Quaternion & q, const Vector & v);
 
   /** REMOVED: ambiguous operattion - need left and right quotients.
    * Return the quaternion ratio \a q/r.
    * If q and r represent
    * rotations, then \a q/r represents the rotation \a q followed by
    * the rotation that is the inverse of \a r.
    * \return the Quaternion ratio \a q/r.
    */
   // friend Quaternion operator/(const Quaternion & q, const Quaternion & r);
 
   /**
    * Multiply this quaternion by \c q and return the result.
    */
   Quaternion & operator*=(const Quaternion & q);
 
   /**
    * Divide this quaternion by \c q and return the result.
    */
   Quaternion & operator/=(const Quaternion & q);
 
   /**
    * Return the quaternion \a a*q, where \a a is a scalar.
    * \a a is a scalar and a*(s,v) = (a*s, a*v).
    * \return the Quaternion \a a*q.
    */
   friend Quaternion operator*(const Quaternion & q, GLfloat a);
 
   /**
    * Multiply by scalar and assign.
    */
   Quaternion & operator*=(GLfloat s);
 
   /**
    * Divide by scalar and assign.
    */
   Quaternion & operator/=(GLfloat s);
 
   /**
    * Return the quaternion \a a*q, where \a a is a scalar.
    * \a a is a scalar and a*(s,v) = (a*s, a*v).
    * \return the Quaternion \a a*q.
    */
   friend Quaternion operator*(GLfloat a, const Quaternion & q);
 
   /**
    * Return the quaternion \a q/a, where \a a is a scalar.
    * \a a is a scalar and (s,v)/a = (s/a, v/a).
    * Report error ZERO_DIVISOR if \a a = 0.
    * \return the Quaternion \a q/a.
    */
   Quaternion operator/(GLfloat scale) const;
 
   /**
    * Normalize this quaternion.
    * See also Quaternion::unit().
    * Report error ZERO_DIVISOR if q = (0,(0,0,0)).
    * \note The value of this quaternion is changed by this operation.
    */
   void normalize();
 
   /**
    * Return a unit quaternion corresponding to this quaternion.
    * See also Quaternion::normalize().
    * Report error ZERO_DIVISOR if q = (0,(0,0,0)).
    * \note The value of this quaternion is not changed by this operation.
    */
   Quaternion unit() const;
 
   /**
    * Return the conjugate of this quaternion.
    * The conjugate of (s,v) is (s,-v).
    * The inverse and conjugate of a unit quaternion are equal.
    * \return the conjugate of this quaternion.
    */
   Quaternion conj() const;
 
   /**
    * Return Inverse of this quaternion.
    * q.inv() = q.conj()/q.norm().
    * The inverse and conjugate of a unit quaternion are equal.
    * Report error ZERO_DIVISOR if q.norm() = 0.
    * If \c q represents a rotation then \c q.inv() represents
    * the opposite, or inverse, rotation.
    * The prefix operator ~ has the same effect.
    * \return the inverse of this quaternion.
    */
   Quaternion inv() const;
 
   /**
    * Return Inverse of this quaternion.
    * This provides alternative syntax for \c inv().
    * \return the inverse of this quaternion.
    */
   Quaternion operator~() const;
 
   /**
    * Apply this quaternion to the vector \c w.
    * The vector \c w is not changed.
    * \return the rotated vector \c q.inv()*w*q.
    */
   Vector apply(const Vector & w) const;
 
   /**
    * Return Vector component \a v of the quaternion \a q = \a (s,v).
    * The same effect can be achieved with the constructor \c Vector::Vector(q).
    */
   Vector vector() const;
 
   /**
    * Return Scalar component \a s of the quaternion \a q = \a (s,v).
    */
   GLfloat scalar() const;
 
   /**
    * Return the norm of this quaternion.
    * The norm is the sum of the squared components
    * and is also the square of the magnitude.
    * \return the norm of this quaternion.
    */
   GLfloat norm() const;
 
   /**
    * Return the magnitude of this quaternion.
    * The magnitude is the square root of the norm.
    * \return the magnitude of this quaternion.
    */
   GLfloat magnitude() const;
 
   /**
    * Compute the rotation matrix corresponding to this quaternion
    * and store it in the Matrix \c m.
    */
   void matrix(Matrix & m) const;
 
   /**
    * Compute the rotation matrix corresponding to this quaternion
    * and store it in the OpenGL matrix \c m.
    */
   void matrix(GL_Matrix m) const;
 
   /**
    * Apply the current quaternion to the current OpenGL matrix.
    */
   void apply() const;
 
   /**
    * Return the axis of rotation of this quaternion.
    * \pre The quaternon must be a unit quaternion.
    * \return a unit vector giving the axis of rotation of the quaternion.
    */
   Vector axis() const;
 
   /**
    * Return the amount of rotation of this quaternion.
    * \pre The quaternon must be a unit quaternion.
    * \return the amount of rotation provided by this quaternion in radians.
    */
   double angle() const;
 
   /**
    * Integrate an angular velocity vector.
    * The rotation represented by this quaternion will be updated by
    * applying the angular velocity \c omega for a short interval of
    * time \c dt.
    * \param omega is an angular velocity vector.
    * \param dt is the time increment for integration.
    */
   void integrate(const Vector & omega, double dt);
 
   /**
    * Return Euler angles for this quaternion.
    * The quaternion gives the rotation that would be obtained by calling
    * \c glRotatef three times, with arguments (\c xr,1,0,0), (\c yr,0,1,0), and
    * (\c zr,0,0,1).
    * \pre The angle transformations are applied in \a x, \a y, \a z order.
    * \return the Euler angles corresponding to this quaternion.
    */
   void euler(double & xr, double & yr, double & zr) const;
 
   /**
    * Use this quaternion to simulate a trackball.
    * To use this function in an OpenGL program, use the default constructor
    * to construct the quaternion \a q = (1,(0,0,0)).  In the display function,
    * use \c q.rotate() to modify the OpenGL model view matrix.
    *
    * Call \c trackball(x1,y1,x2,y2) to record a small mouse movement
    * from \a (x1,y1) to \a (x2,y2).  The coordinates should be normalized
    * so that both \a x and \a y are in [-1,1].
    *
    * This function will compute a rotation will apply this rotation to
    * the current quaternion.  The rotation simulates a trackball.
    */
   void trackball(GLfloat x1, GLfloat y1, GLfloat x2, GLfloat y2);
 
   /**
    * log(q) of a quaternion is a pure quaternion (scalar part 0).
    * log(q) is an element of the Lie algebra of the quaternion group.
    * log(1;0) = (0;0).
    */
   friend Quaternion log(const Quaternion & q)
   {
      return (q.s == 1) ?
             Quaternion(0.0, Vector()) :
             Quaternion(0.0, (acos(q.s) / sqrt(1 - q. s * q.s))* q.v);
   }
 
   /**
    * exp(q) is defined for pure quaternions (scalar part zero) only.
    * (that is, members of the Lie algebra of the quaternion group).
    * This function ignores the scalar part - it does not check for zero.
    */
   friend Quaternion exp(const Quaternion & q)
   {
      Vector a = q.v;
      GLfloat angle = a.length();
      return (angle == 0) ?
             Quaternion() :
             Quaternion(cos(angle), a * sin(angle)/angle);
   }
 
   /**
    * exp(v) is a quaternion.  The vector is treated as a "pure" quaternion
    * (that is, as a member of the Lie algebra of the quaternion group).
    */
   friend Quaternion exp(const Vector & v)
   {
      GLfloat angle = v.length();
      return (angle == 0) ?
             Quaternion() :
             Quaternion(cos(angle), v * sin(angle)/angle);
   }
 
//   /**
//    * The 'logarithm' of a unit quaternion is a vector.
//    * This function is defined because it appears in descriptions of quaternions.
//    * Although the name 'ln' is appropriate in some ways, this function must be used with
//    * caution because it may not have the properties you expect of a logarithm.
//    * For example, \a ln(q1)+ln(q2) makes sense only if \a q1 and \a q2 have the
//    * same axis.
//    * \pre The quaternion must be a unit quaternion.
//    */
//   friend Vector ln(const Quaternion & q);
//
//   /**
//    * The 'exponent' of a vector is a unit quaternion.
//    * Although the name 'exp' is appropriate in some ways, this function must be used with
//    * caution because it may not have the properties you expect of an exponent.
//    * For example, \a exp(v1+v2)=exp(v1)*exp(v2) only if \a v1 and \a v2 have the same
//    * direction.
//    */
//   friend Quaternion exp(const Vector & v);
 
   /**
    * Return the dot product of quaternions \c q and \c r.
    * The dot product of \a (s,u) and \a (t,v) is \a st+dot(u,v)
    * where \a dot(u,v) denotes the vector dot product of \a u and \a v.
    * \return the dot product of quaternions \c q and \c r.
    */
   friend GLfloat dot(const Quaternion & q, const Quaternion & r);
 
   /**
    * Compare two quaternions.
    * This function returns \c true only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator==(const Quaternion & x, const Quaternion & y);
 
   /**
    * Compare two quaternions.
    * This function returns \c false only if corresponding components are exactly equal.
    * Values that are theoretically equal but computed in different ways are likely
    * to be unequal according to this function.
    */
   friend bool operator!=(const Quaternion & x, const Quaternion & y);
 
   /**
    * Write the quaternion to the output stream as \a s \a (x,y,z).
    */
   friend std::ostream & operator<<(std::ostream & os, const Quaternion & q);
 
private:
 
   /** Scalar component of quaternion. */
   GLfloat s;
 
   /** Vector component of quaternion. */
   Vector v;
};
 
 
 
/**
 * An instance represents a camera.
 * An instance of class \c Camera is used to control the view.
 * The function \c Camera::apply() calls \c gluLookAt() with
 * appropriate parameters.  The other functions of this class
 * set up these parameters.
 *
 * Camera movements may be smooth (interlpolated between end
 * points) or abrupt.  The function Camera::setMode() determines
 * the kind of movement.
 * Smooth transitions are obtained by calling \c Camera::idle()
 * in the GLUT idle callback function, together with the various
 * camera movement functions.
 */
 
 
 
//typedef void (BaseCamera::*IdleMemFunc)();
 
//IdleMemFunc gIdleMemFunc;
 
/**
* All derived class must implement at least two
* virtual functions \c BaseCamera::idle() which will
* typically be called in the callback function for glutIdleFunc,
* i.e. within callback function "idle()" DerivedClass::idle() will
* be called and later in "main" callback "idle()" will be called
* as parameter to "glutIdleFunc".
* void idle() { DerivedClass::idle();}
* int main() { ... glutIdleFunc(idle); ...}
* \c BaseCamera::apply() const. Similarly inside callback function
* for "glutDisplayFunc", say "display()", DerivedClass::apply() will be
* called.
* i.e. void display() { DerivedClass::apply();..}
* int main(){...glutDisplayFunc(display); ...}
*
*/
 
class BaseCamera
{
public:
   /** Construct an instance of a \a BaseCamera.
    * This constructor is usually called from a subclass.
    */
   BaseCamera() : steps(0), maxSteps(50)
   {}
 
   /**
   * Set the number of steps required for a camera movement.
   * \param s is the number of steps executed during a movement.
   * \a s = 10 will give a rapid, possibly jerky, movement.
   * \a s = 100 will give a slow, smooth movement.
   */
   void setSteps(int s)
   {
      maxSteps=s;
   }
 
   /**
   * This function should be called from the GLUT idle() function
   * or its equivalent.  It performs one step of the motion until
   * the motion is complete, after which it has no effect.
   */
   virtual void idle()=0;
 
   /**
    * This function should be called from the GLUT display() function
    * or its equivalent.
    */
   virtual void apply() const =0;
 
   /**
   * This function is just a suggested wrapper for all derived class
   * for initializations, such as set up first and last position for
   * motions in /c Interpolator, or set up new eye and view parameter
   * for /c Camera. I don't want to force new class to wrap their
   * initialization actions within this function, so, I don't make this
   * method as abstract one.
   * (added by nick)
   */
   virtual void update()
   {}
 
protected:
 
   /** Current value of step counter. */
   int steps;
 
   /** Maximum number of steps for a smooth movement. */
   int maxSteps;
 
};
 
 
/**
 * An instance represents a camera.
 * An instance of class \c Camera is used to control the view.
 * The function \c Camera::apply() calls \c gluLookAt() with
 * appropriate parameters.  The other functions of this class
 * set up these parameters.
 *
 * Camera movements may be smooth (interlpolated between end
 * points) or abrupt.  The function Camera::setMode() determines
 * the kind of movement.
 * Smooth transitions are obtained by calling \c Camera::idle()
 * in the GLUT idle callback function, together with the various
 * camera movement functions.
 */
class Camera : public BaseCamera
{
   //friend void myIdle();
public:
 
   /**
    * Construct a default camera.  The view is set as follows:
    * Eye = (0,0,1).  Model=(0,0,0).  Up=(0,1,0).
    */
   Camera();
 
   /**
    * Construct a camera for a particular view.
    * \param eye sets the eye coordinates for \c gluLookAt().
    * \param model sets the model coordinates for \c gluLookAt().
    * \param up sets the 'up' vector for \c gluLookAt().
    */
   Camera(const Point & eye, const Point & model, const Vector & up);
 
   /**
    * Construct a camera for a particular view.
    * The 'up' vector is set to (0,1,0).
    * \param eye sets the eye coordinates for \c gluLookAt().
    * \param model sets the model coordinates for \c gluLookAt().
    */
   Camera(const Point & eye, const Point & model);
 
   /**
    * Construct a camera for a particular view.
    * The 'up' vector is set to (0,1,0).
    * The model point is set to (0,0,0).
    * \param eye sets the eye coordinates for \c gluLookAt().
    */
   explicit Camera(const Point & eye);
 
   /**
    * Set the camera for a particular view.
    * \param eye sets the eye coordinates for \c gluLookAt().
    * \param model sets the model coordinates for \c gluLookAt().
    * \param up sets the 'up' vector for \c gluLookAt().
    */
   void set(const Point & eye, const Point & model, const Vector & up)
   ;
 
   /**
    * Set the camera for a particular view.
    * The 'up' vector is set to (0,1,0).
    * \param eye sets the eye coordinates for \c gluLookAt().
    * \param model sets the model coordinates for \c gluLookAt().
    */
   void set(const Point & eye, const Point & model)
   ;
 
   /**
    * Set the camera for a particular view.
    * The model point is set to (0,0,0).
    * The 'up' vector is set to (0,1,0).
    * \param eye sets the eye coordinates for \c gluLookAt().
    */
   void set(const Point & eye)
   ;
 
 
   /**
    * Update the camera position.
    * This function should be called from the GLUT idle() function
    * or its equivalent.  It performs one step of the motion until
    * the motion is complete, after which it has no effect.
    */
   virtual void idle();
 
   /**
    * Apply the camera position to the model-view matrix.
    * This function calls \c gluLookAt() with appropriate parameters.
    * It should be called in the GLUT display function or its equivalent
    * after the model-view matrix has been initialized.
    */
   void apply() const;
 
   /**
    * Move the camera upwards (+Y direction).
    * \param distance gives the amount of movement.  Negative values move
    * the camera downwards.
    */
   void moveUp(GLfloat distance);
 
   /**
    * Move the camera forwards.
    * The camera moves towards the model.
    * \param distance gives the amount of movement.  Negative values move
    * the camera away from the model..
    */
   void moveForward(GLfloat distance);
 
   /**
    * Move the camera to its left.
    * \param distance gives the amount of movement.  Negative values move
    * the camera to its right.
    */
   void moveLeft(GLfloat distance);
 
   /**
    * Tilt the camera upwards.
    * \param angle gives the angle through which the camera should be moved.
    * Negative values give a downward tilt.
    */
   void tiltUp(double angle);
 
   /**
    * Pan the camera to the left.
    * \param angle gives the angle through which the camera should be panned.
    * Negative values give panning to the right.
    */
   void panLeft(double angle);
 
   /**
   * Set the kind of motion required.
   * \param mode should be set to \c true for smooth movements
   * or to \c false for abrupt movements.
   */
   void setMode(bool mode)
   {
      smooth = mode;
   }
 
   /**
    * Write a representation of the camera to a stream.
    * This function displays the current settings of the eye point,
    * the model point, and the up vector.
    * \param os is a reference to an output stream.
    * \param c is a reference to a camera.
    */
 
   friend std::ostream & operator<<(std::ostream & os, const Camera & c);
 
protected:
 
   /** Update the camera position. */
   void update(const Point & e, const Point & m);
 
   /** Eye coordinates for gluLookAt(). */
   Point eye;
 
   /** Previous eye coordinates. */
   Point eyeOld;
 
   /** Updated eye coordinates. */
   Point eyeNew;
 
 
   /** Model coordinates for gluLookAt() */
   Point model;
 
   /** Previous model coordinates */
   Point modelOld;
 
   /** Updated model coordinates */
   Point modelNew;
 
   /** Up vector for gluLookAt() */
   Vector up;
 
   /** Mode for motion */
   bool smooth;
};
 
/**
 * An instance is an interpolator for rotations.
 * An instance of Interpolator is used to interpolate between two rotations.
 * In other words, given two orientations of an object, an interpolator can
 * be used to move the object smoothly from one orientation to the other.
 * An interpolator is represented by two quaternions representing the
 * extreme values of the rotations.  When a real number \a t in the interval
 * [0,1] is passed to the interpolator, the interpolator constructs a third
 * quaternion representing an intermediate orientation.  This quaternion can
 * be used to generate a rotation matrix or transform the OpenGL model view.
 *
 * The simplest way to use this class is as follows:
 *       -# Construct two quaternions, \a q and \a r, corresponding to the
 *          first and last positions of the desired movement.
 *       -# Construct an interpolator \a i(q,r).
 *       -# In the OpenGL display function, call \c i.apply(t), giving
 *          \a t a sequence of values ranging from 0 to 1.
 *
 *
 */
class Interpolator: public BaseCamera
{
public:
 
   /**
    * Construct an interpolator with dummy values for the quaternions.
    * The interpolator may fail if it is used without first setting the
    * start and finish quaternions.  Use \c Quaternion::set() to change
    * the first and last quaternions to useful values.
    */
   Interpolator();
 
   /**
    * Construct an interpolator to rotate from \c qFirst to \c qLast.
    * Report error \c BAD_INTERPOLATOR_ARG if \c qFirst = \c qLast.
    */
   Interpolator(const Quaternion & qFirst, const Quaternion & qLast);
 
   /**
    * Change the first and last quaternions of an existing interpolator.
    * Report error \c BAD_INTERPOLATOR_ARG if \c qFirst = \c qLast.
    */
   void set(const Quaternion & qFirst, const Quaternion & qLast)
   ;
 
   /**
    * Perform the rotation at position \a t, 0 <= \a t <= 1
    * by modifying the current OpenGL transformation matrix.
    */
   void apply(double t) const;
 
   /**
    * Compute the rotation matrix for position \a t, 0 <= \a t <= 1.
    */
   void getMatrix(double t, GL_Matrix m) const;
 
   /**
    * Return the quaternion for position \a t, 0 <= \a t <= 1.
    */
   Quaternion getQuaternion(double t) const;
 
   /**
   * Update the Interpolator position.
   * This function should be called from the GLUT idle() function
   * or its equivalent.  It performs one step of the motion until
   * the motion is complete, after which it has no effect.
   */
   virtual void idle();
 
   /**
      * Apply the Interpolator rotation to the model-view matrix.
      * This function calls apply(t) with parameters t which is actually steps/maxSteps.
      * It should be called in the GLUT display function or its equivalent
      * after the model-view matrix has been initialized.
      */
   virtual void apply() const;
 
   /**
   * Update Interpolator first and last Quaternion.
      */
   virtual void update()
   {
      ;
   }
 
protected:
 
   /** Initialize the interpolator. */
   void initialize();
 
   /** The quaternion currently in use. */
   Quaternion current;
 
   /** The start quaternion for the interpolation. */
   Quaternion first;
 
   /** The final quaternion for the interpolation. */
   Quaternion last;
 
   /** The angle between the first and last positions. */
   double omega;
 
   /** The sine of the angle between the first and last positions. */
   double sinOmega;
 
   /** The cosine of the angle between the first and last positions. */
   double cosOmega;
};
 
 
 
/**
  * An instance is an RGB pixel array.
  * An instance of PixelMap holds an array that OpenGL can use
  * as a texture or to display directly.
  *
  * This class is implemented so that it does not depend on Windows.
  * In principle, it should even work with Linux and other operating
  * systems, although this has not been tested thoroughly.
  *
  * If you want to use a \a BMP image as a texture, the width and height
  * must be powers of two.  Thus a 256 x 512 \a BMP image is acceptable
  * but a 640 x 400 image is not.  This restriction does not apply to
  * images that are displayed as bitmaps.
  *
  * Since an instance of \c PixelMap contains pointers, there are memory
  * management issues.  The destructor deletes pixel and other data
  * associated with a pixel map.  The copy constructor is declared privately
  * and not implemented: this prevents instances being passed by value,
  * although they may be passed by reference.  Since there is a default
  * constructor, you can construct arrays of PixelMap's or pointers to
  * PixelMap.
  */
class PixelMap
{
public:
 
   /**
    * Construct an empty pixel map.
    */
   PixelMap();
 
   /**
    * Read a pixel map from a file.
    * \param bmpFileName is a path to a BMP file.
    */
   PixelMap(const char *bmpFileName);
 
   /**
    * Construct a pixel map from a region of the pixel buffer.
    * \param x is the \a X coordinate of the lower left corner of the region.
    * \param y is the \a Y coordinate of the lower left corner of the region.
    * \param w is the width of the region in pixels.
    * \param h is the height of the region in pixels.
    * \note See also PixelMap::read(x, y, w, h).
    */
   PixelMap(GLint x, GLint y, GLsizei w, GLsizei h);
 
   /**
    * Delete a pixel map.
    */
   ~PixelMap();
 
   /**
    * Read a BMP file and convert it to a pixel map.
    * Previous data associated with this object will be deleted.
    * The BMP file must satisfy several criteria if it is to
    * be converted successfully.  Conversion failure is
    * indicated by the following CUGL errors.
    * \li \c OPEN_FAILED: the file could not be opened
    * \li \c NOT_BMP_FILE: the file is not a BMP file
    * \li \c NOT_24_BITS: the format is not 24 bits/pixel
    * \li \c COMPRESSED_BMP_FILE: the file is compressed
    * \li \c NOT_ENOUGH_MEMORY: there is not enough memory to store the data
    *
    * OpenGL cannot use a bitmap as a texture if its dimensions are not powers of 2.
    * To check whether the bitmap has acceptable dimensions, call \c PixelMap::badSize().
    * To convert the bitmap dimensions to acceptable values, call \c PixelMap::rescale().
    * \param bmpFileName is the name of the file to be read, with the extension '.bmp'.
    */
   void read(const char *bmpFileName);
 
   /**
    * Read a pixel map from a region of the framebuffer.
    * This function is similar to the constructor with the same parameters,
    * but allocates new memory only if necessary.
    * \param x is the \a X coordinate of the lower left corner of the region.
    * \param y is the \a Y coordinate of the lower left corner of the region.
    * \param w is the width of the region in pixels.
    * \param h is the height of the region in pixels.
    * \param mode is passed to glReadBuffer before reading the pixels.
    * \note See also the constructor PixelMap(x, y, w, h).
    */
   void read(GLint x, GLint y, GLsizei w, GLsizei h, GLenum mode = GL_FRONT);
 
   /**
    * Write a pixel map to an output stream as a \a BMP file.
    * \param bmpFileName is the name of the file to be written, with the extension '.bmp'.
    */
   void write(const char *bmpFileName);
 
   /**
    * Check the dimensions of the bit map.
    * If the dimensions are not powers of 2, return \c true.
    * If the dimensions of the bitmap are not powers of 2,
    * OpenGL cannot use it as a texture.
    * You should call \c PixelMap::rescale() to resize the bit map.
    */
   bool badSize();
 
   /**
    * Rescale a bit map whose dimensions are not powers of 2.
    * The new image will be distorted; the amount of distortion depends
    * on how much the dimensions have to be altered.
    * Use \c PixelMap::badSize() to determine whether the dimensions are powers of 2.
    */
   void rescale();
 
   /**
    * Draw the pixel map at the current raster position.
    * Error \c NO_PIX_MAP if there is no pixel map avaialble.
    */
   void draw();
 
   /**
    * Set texture parameters for the pixel map.
    * \param name is an OpenGL index for the texture parameters
    * provided by the caller.
    * \note This function sets \c GL_TEXTURE_MAG_FILTER and
    * \c GL_TEXTURE_MIN_FILTER to \c GL_NEAREST.  Call glTexParameter()
    * to change these settings.
    * \note When this function returns, OpenGL has copied the pixel map
    * into its own memory space.  It is therefore safe to delete the
    * PixelMap instance after calling setTexture().
    */
   void setTexture(GLuint name);
 
   /**
    * Set texture parameters for the pixel mipmaps.
    * Construct a family of mipmaps for texturing.
    * \param name is an OpenGL index for the texture parameters
    * provided by the caller.
    * \note This function sets \c GL_TEXTURE_MIN_FILTER to
    * \c GL_LINEAR_MIPMAP_NEAREST.  Call glTexParameter()
    * to change this setting.
    * \note When this function returns, OpenGL has copied the pixel map
    * into its own memory space.  It is therefore safe to delete the
    * PixelMap instance after calling setMipmaps().
    * \note Call this function once only for each texture you need in
    * your program.  Use \c glBindTexture to select textures in the
    * display function.
    */
   void setMipmaps(GLuint name);
 
   /** Check that two pixel maps are compatible for combining. */
   friend bool compatible(const PixelMap & m1, const PixelMap & m2);
 
   /** Combine two pixel maps.
    * \param m1 is the first map to be combined.
    * \param m2 is the second map to be combined.
    * \param res is the resulting pixel map.
             The caller is responsible for constructing this map;
             it is not constructed by the function.
      \param prop is the mixing proportion: 0 gives m1,
             0.5 gives half of each, and 1 gives m2.
    */
   friend void mix(const PixelMap & m1, const PixelMap & m2, PixelMap & res, double prop);
 
   /** Select a region from a pixel map.
    * \param src is the pixel map from which the data is extracted.
    * \param xp defines the left side of the selected region.
    * \param yp defines the right side of the selected region.
    * \param width is the width of the selected region.
    * \param height is the height of the selected region.
    */
   void select(const PixelMap & src, int xp, int yp, int width, int height);
 
   /**
    * Return number of rows in pixel map.
    */
   unsigned long getRows() const
   {
      return numRows;
   }
 
   /**
    * Return number of columns in pixel map.
    */
   unsigned long getColumns() const
   {
      return numCols;
   }
 
   /**
    * Return bytes of memory used by pixel map.
    */
   unsigned long getSize() const
   {
      return size;
   }
 
   /**
    * Return name of BMP file.
    */
   char *getName() const
   {
      return fileName;
   }
 
   /**
    * Return \c true if a pixel map has been read successfully.
    */
   bool ready() const
   {
      return pixels != 0;
   }
 
   /**
    * Write a description of the pixel map to the output stream.
    */
   friend std::ostream & operator<<(std::ostream & os, const PixelMap & pm);
 
private:
 
   /**
    * The copy constructor is private and unimplemented.
    * This prevents pixel maps from being copied by mistake.
    */
   PixelMap(const PixelMap & pm);
 
   /** Allocate memory for a new pixel map if necessary. */
   bool allocate(unsigned long newSize);
 
   /** Number of rows in the pixel map. */
   unsigned long numRows;
 
   /** Number of columns in the pixel map. */
   unsigned long numCols;
 
   /** Size, in bytes, of the pixel map. */
   unsigned long size;
 
   /** Name of the file used to store the pixel map. */
   char *fileName;
 
   /** Pointer to the pixels of the map. */
   unsigned char *pixels;
};
 
 
/** The type of an array of GLfloats. */
typedef GLfloat* GLfloatArray;
 
/**
 * An instance is a surface of revolution represented by its profile.
 
 * The shape of the surface is determined by a 'profile'.
 *
 * The profile is defined by an array.  Each component of the
 * array is a pair \a (h,r) in which \a h is measured along
 * the axis of revolution and \a r is the distance from the
 * axis.  (If the axis is assumed to be vertical, then \a h
 * suggests 'height' and \a r suggests 'radius'.)
 * The surface is generated by rotating the profile
 * about the \a Z-axis.
 *
 * The function \c Revolute::draw() generates
 * vertex coordinates, normals, and texture coordinates for the
 * surface.
 *
 * It is up to the caller to set the OpenGL state for rendering:
 * this includes setting material properties and defining rules for polygon shading,
 * or binding an appropriate texture.
 *
 * The normals generated by \c Revolute::process() are obtained by averaging
 * the normals of the polygons that meet at each vertex.  Consequently,
 * if \c GL_SMOOTH shading is used and sufficient points and slices are specified,
 * the object should look fairly smooth.
 */
class Revolute
{
public:
   /**
    * Construct a surface of revolution.
    *
    * \param numSteps is the number of coordinates in the profile.
    * \param profile is the address of a 2D array giving the points of the profile.
    */
   Revolute(int numSteps, GLfloat profile[][2]);
 
   /**
    * Delete a surface of revolution.
    */
   ~Revolute();
 
   /**
    * Set the number of "slices".
    * The value determines the visual quality of the object.
    * The number of slices should be an odd number greater than 2.
    * If it is not, Revolute::process() will change it.
    * By default, the number of slices is 45, corresponding to
    * 8 degrees per slice.
    * \note After changing the number of slices, you must call
    * \c Revolute::process() again before calling \c Revolute::draw().
    */
   void setSlices(int slices);
 
   /**
    * Set the eccentricity of the surface.
    * By default, the eccentricity is 0, giving a surface
    * with a circular cross-section.  Setting the eccentricity
    * to a non-zero value gives a surface with an elliptical
    * cross-section.  As the eccentricity approaches 1,
    * the major axis of the ellipse increases without limit
    * as the cross-section approaches a parabola.
    * \param ecc must be greater than or equal to zero and
    * less than 1.
    * \note After changing the eccentricity, you must call
    * \c Revolute::process() again before calling \c Revolute::draw().
    */
   void setEccentricity(double ecc);
 
   /**
    * Compute vertexes, normals, and texture coordinates for the surface of revolution.
    * This function must be called \b before calling \c Revolute::draw().
    */
   void process();
 
   /**
    * Draw the surface of revolution.
    * \param showNormals determines whether surface normals will be drawn.
    * Note that surface normals are computed for lighting purposes anyway:
    * \c showNormals is provided mainly as a debugging aid and should not
    * normally be needed.
    * \note Revolute::process() must be called before this function.
    */
   void draw(bool showNormals = false);
 
private:
 
   /**
    * The copy constructor is private and not implemented:
    * instances cannot be copied implicitly.
    */
   Revolute(const Revolute &);
 
   /**
    * The assignment operator is private and not implemented:
    * instances cannot be assigned.
    */
   const Revolute& operator=(const Revolute &);
 
   /** */
   int numSteps;
 
   /** */
   int numSlices;
 
   /** */
   double eccentricity;
 
   /** */
   bool ready;
 
   /** */
   GLfloatArray *coor;
 
   /** */
   GLfloat *texCoor;
 
   /** */
   Point *points;
 
   /** */
   Vector *normals;
 
   /** */
   Vector *faceNormals;
};
 
 
/**
 * \defgroup freefun Miscellaneous functions
 */
 
//@{
 
/**
 * Convert degrees to radians.
 */
inline double radians(double angle)
{
   return angle * PI / 180;
}
 
/**
 * Convert radians to degrees.
 */
inline double degrees(double angle)
{
   return angle * 180 / PI;
}
 
/** Return the square of the argument. */
inline double sqr(double x)
{
   return x * x;
}
 
/** Return a random integer in [0, max). */
inline unsigned int randInt(unsigned int max)
{
   static unsigned int seed = 12345678;
   unsigned int result;
   const unsigned int bucket = static_cast<unsigned int>(4294967296.0 / max);
   do
   {
      seed = 1664525 * seed + 1013904223;
      result = seed / bucket;
   }
   while (result >= max);
   return result;
}
 
/** Return a random integer in [-max, max]. */
inline int randSym(unsigned int max)
{
   return randInt(2 * max + 1) - max;
}
 
/** Return a random double in [0, 1). */
inline double randReal()
{
   int BIG = 10000000;
   return double(randInt(BIG)) / double(BIG);
}
 
/**
 * Construct normals for OpenGL triangle strips.
 * Given a set of vertexes definining a triangle strip in OpenGL format,
 * this functions constructs a normal corresponding to each vertex.
 * \param points is an array of points giving the vertex coordinates.
 * \param normals is an array of vectors in which the vertex normals will be stored.
 * \param numPoints is the number of vertexes provided and the number of normals
 * that will be calculated.
 * \param neg specifies negated normals, if \a true.
          The default is \a false.
 * \note To avoid allocation and deallocation overhead, this function uses a
 * a fixed amount of workspace that allows up to 100 vertexes to be processed.
 * If numPoints > 100, the function will have no effect.
 * \note For efficiency, it is better to compute the normals during initialization
 * rather than each time the model is displayed.
 */
void triStripNormals(Point points[], Vector normals[], int numPoints, bool neg = false);
 
/**
 * Constructs a surface of revolution.
 * Constructs and renders a surface of revolution obtained by rotating
 * a curve about the \a Y axis.
 
 * \param numSteps is the number of points in the array \c coor
 * \param coor is the 2D coordinates of points of the profile
 * \param numSlices is the number of pie-shaped slices used to render the surface.
 * \param drawNormals determines whether normals are generated.  By default, normals
 *        are not generated.
 
 * For example, if \c numSlices = 20, points will be constructed at 360/20 = 18
 * degree intervals.
 
 * This function constructs an array of points in 3D space and then issues
 * calls to \c glVertex().  If \c drawNormals is true, it also issues calls to
 * \c glNormal().  The effect of these calls is to define a 3D mesh.
 * It is up to the caller to set the OpenGL state for rendering:
 * this includes setting material properties and defining rules for polygon shading.
 
 * The normals generated by \c revolve() are obtained by averaging the normals of the
 * polygons that meet at each vertex.  Consequently, if \c GL_SMOOTH shading is used
 * and enough points are specified, the object should look fairly smooth.
 
 * \note This function has been replaced by class \c Revolute, which provides more
 * functionality and is more efficient.
 */
void revolve(int numSteps, GLfloat coor[][2], int numSlices, bool drawNormals = false);
 
 
 
// Additional inline functions for Point and Vector classes.
 
/**
 * Call gluLookAt() looking at the origin of the model (0,0,0)
 * with 'up' vector (0,1,0).
 * \param eye is the position of the viewer's eye.
 */
inline void lookAt(Point eye)
{
   gluLookAt
   (
      eye[0], eye[1], eye[2],
      0, 0, 0,
      0, 1, 0
   );
}
 
/**
 * Call gluLookAt() with 'up' vector (0,1,0).
 * \param eye is the position of the viewer's eye in model coordinates.
 * \param model is the point at the centre of the view in model coordinates.
 */
inline void lookAt(Point eye, Point model)
{
   gluLookAt
   (
      eye[0], eye[1], eye[2],
      model[0], model[1], model[2],
      0, 1, 0
   );
}
 
/**
 * Call gluLookAt().
 * \param eye is the position of the viewer's eye in model coordinates.
 * \param model is the point at the centre of the view in model coordinates.
 * \param up is a vector giving the upwards direction in the model.
 */
inline void lookAt(Point eye, Point model, Vector up)
{
   gluLookAt
   (
      eye[0], eye[1], eye[2],
      model[0], model[1], model[2],
      up[0], up[1], up[2]
   );
}
 
// End of free functions
//@}
 
 
/**
 * \defgroup model Materials and Models
 */
 
//@{
 
/**
 * \todo Models should probably be put into a separate file.
 */
 
// Axes
 
/**
 * Draw coordinate axes.
 * This function draws the three principal axes in the current
 * position, using \c size to determine the length of each axis.
 * The axes are colour-coded: \a X = red, \a Y = green, \a Z = blue.
 */
void axes(GLfloat size = 1);
 
/**
 * Draw an aircraft.
 * The argument determines whether the plane is drawn with a metallic
 * colour (\c shadow = \c false, the default) or black, to act as a
 * shadow (\c shadow = \c true).
 * The aircraft is built using Bezier surfaces. If you use glScalef
 * to change its size, you must enable \c GL_NORMALIZE to correct
 * the normals. Otherwise, the lighting will be wrong.
 */
void buildPlane(bool shadow = false);
 
/**
 * Construct a GL call list for the aircraft and return its index.
 */
GLuint makePlaneList(bool shadow = false);
 
 
/**
 * Enumeration for materials.
 */
enum MATERIAL
{
   BLACK,            /**< Black. */
   WHITE,            /**< White. */
   RED,              /**< Red. */
   GREEN,            /**< Green. */
   BLUE,             /**< Blue. */
   METAL,            /**< General metallic colour.*/
   GLASS,            /**< General transparent marterial.*/
   BRASS,            /**< Brass. */
   BRONZE,           /**< Matt bronze. */
   POLISHED_BRONZE,  /**< Polished bronze. */
   CHROME,           /**< Chrome. */
   COPPER,           /**< Matt copper. */
   POLISHED_COPPER,  /**< Polished copper. */
   GOLD,             /**< Matt gold. */
   POLISHED_GOLD,    /**< Polished gold. */
   PEWTER,           /**< Pewter. */
   SILVER,           /**< Matt silver. */
   POLISHED_SILVER,  /**< Polished silver. */
   EMERALD,          /**< Emerald. */
   JADE,             /**< Jade. */
   OBSIDIAN,         /**< Obsidian. */
   PEARL,            /**< Pearl. */
   RUBY,             /**< Ruby. */
   TURQUOISE,        /**< Turquoise. */
   BLACK_PLASTIC,    /**< Black plastic. */
   BLACK_RUBBER,      /**< Black rubber. */
   LAST_VALUE
};
 
/**
 * Set material values from the enumeration.
 \param m is chosen from the enumeration \c MATERIAL
 \param face should be one of \c GL_FRONT (the default),
 \c GL_BACK, or \c GL_FRONT_AND_BACK.
 */
void setMaterial(const int m, GLenum face = GL_FRONT);
 
// The following array describes material properties.
//    ambiance  [4],
//    diffuse   [4],
//    specular  [4],
//    shininess [1].
// 13 values are required to define a material.
 
/**
 * Add a material to the set of built-in materials.
 * The array of material has a fixed size of 100.
 * An attempt to create more than 100 materials will fail.
 * The parameters specify:
 *
 * - RGBA values for ambient light
 * - RGBA values for diffuse light
 * - RGBA values for specular light
 * - Value for shininess exponent
 * \return the index of the new material.
 */
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);
 
/**
 * Add a material to the set of built-in materials.
 * The array of material has a fixed size of 100.
 * An attempt to create more than 100 materials will fail.
 * \param params specifies:
 * - RGBA values for ambient light
 * - RGBA values for diffuse light
 * - RGBA values for specular light
 * - Value for shininess exponent
 * \return the index of the new material.
 */
int addMaterial(GLfloat params[]);
 
// End of free functions
//@}
 
//=============================== End of declarations ==============================//
 
// Code below this point consists of inline function definitions.
// Functions that are not defined here are defined in cugl.cpp.
 
// class Point: inlined constructors and member functions.
 
inline Point::Point (GLfloat coordinates[])
{
   x = coordinates[0];
   y = coordinates[1];
   z = coordinates[2];
   w = coordinates[3];
}
 
inline Point Point::operator+=(const Vector & v)
{
   x += w * v.x;
   y += w * v.y;
   z += w * v.z;
   return *this;
}
 
inline Point Point::operator-=(const Vector & v)
{
   x -= w * v.x;
   y -= w * v.y;
   z -= w * v.z;
   return *this;
}
 
inline Point Point::operator/(GLfloat s) const
{
   return Point(x, y, z, s * w);
}
 
inline void Point::draw() const
{
   if (w == 0)
      return;
   glVertex4f(x, y, z, w);
}
 
inline void Point::light(GLenum lightNum) const
{
   GLfloat p[4];
   p[0] = x;
   p[1] = y;
   p[2] = z;
   p[3] = w;
   glLightfv(lightNum, GL_POSITION, p);
}
 
inline void Point::translate() const
{
   glTranslatef(x/w, y/w, z/w);
}
 
// class Point: inlined friend functions.
 
inline Point operator+(const Point & p, const Vector & v)
{
   return Point(p.x + p.w * v.x, p.y + p.w * v.y, p.z + p.w * v.z, p.w);
}
 
inline Point operator+(const Vector & v, const Point & p)
{
   return Point(p.x + p.w * v.x, p.y + p.w * v.y, p.z + p.w * v.z, p.w);
}
 
inline Point operator*(const Point & p, GLfloat s)
{
   return Point(s * p.x, s * p.y, s * p.z, s * p.w);
}
 
inline Point operator*(GLfloat s, const Point & p)
{
   return Point(s * p.x, s * p.y, s * p.z, s * p.w);
}
 
inline bool operator==(const Point & p, const Point & q)
{
   return p.x == q.x && p.y == q.y && p.z == q.z && p.w == q.w;
}
 
inline bool operator!=(const Point & p, const Point & q)
{
   return !(p == q);
}
 
inline GLfloat dist(const Point & p, const Plane & s)
{
   return p.x * s.a + p.y * s.b + p.z * s.c + p.w * s.d;
}
 
inline GLfloat dist(const Plane & s, const Point & p)
{
   return s.a * p.x + s.b * p.y + s.c * p.z + s.d * p.w;
}
 
// Class Line: inlined member functions.
 
inline Line::Line(const Point & p, const Point & q)
      : s(p), f(q)
{ }
 
inline Line::Line(const Point & p, const Vector & v)
      : s(p), f(p + v)
{ }
 
 
// Class Line: inlined friend functions.
 
inline bool operator==(const Line & m, const Line & n)
{
   return m.s == n.s && m.f == n.f;
}
 
inline bool operator!=(const Line & m, const Line & n)
{
   return !(m == n);
}
 
 
 
// Class Plane: inlined member functions.
 
inline Vector Plane::normal() const
{
   return Vector(a,b,c).unit();
}
 
// Class Plane: inlined friend functions.
 
inline bool operator==(const Plane & p, const Plane & q)
{
   return p.a == q.a && p.b == q.b && p.c == q.c && p.d == q.d;
}
 
inline bool operator!=(const Plane & p, const Plane & q)
{
   return !(p == q);
}
 
 
 
// class Vector: inlined constructors
 
inline Vector::Vector (GLfloat coordinates[])
{
   x = coordinates[0];
   y = coordinates[1];
   z = coordinates[2];
}
 
inline Vector::Vector(const Point & p, const Point & q)
{
   if (p.w == 0 || q.w == 0)
   {
      x = 0;
      y = 0;
      z = 0;
   }
   else
   {
      x = p.x/p.w - q.x/q.w;
      y = p.y/p.w - q.y/q.w;
      z = p.z/p.w - q.z/q.w;
   }
}
 
inline Vector::Vector(const Quaternion & q)
{
   x = q.v.x;
   y = q.v.y;
   z = q.v.z;
}
 
// class Vector: inlined member functions.
 
inline Vector Vector::operator+=(const Vector & v)
{
   x += v.x;
   y += v.y;
   z += v.z;
   return *this;
}
 
inline Vector Vector::operator-=(const Vector & v)
{
   x -= v.x;
   y -= v.y;
   z -= v.z;
   return *this;
}
 
inline Vector Vector::operator-() const
{
   return Vector(-x, -y, -z);
}
 
inline Vector Vector::operator*=(GLfloat scale)
{
   x *= scale;
   y *= scale;
   z *= scale;
   return *this;
}
 
inline GLfloat Vector::norm() const
{
   return x * x + y * y + z * z;
}
 
inline GLfloat Vector::length() const
{
   return GLfloat(sqrt(norm()));
}
 
inline void Vector::translate() const
{
   glTranslatef(x, y, z);
}
 
inline Matrix Vector::skew()
{
   return Matrix(
             0.0f,    z,    -y, 0.0f,
             -z, 0.0f,     x, 0.0f,
             y,   -x, -0.0f, 0.0f,
             0.0f, 0.0f,  0.0f, 1.0f );
}
 
inline void Vector::drawNormal() const
{
   glNormal3f(x, y, z);
}
 
inline void Vector::draw(const Point & p) const
{
   glBegin(GL_LINES);
   p.draw();
   (p + (*this)).draw();
   glEnd();
}
 
// class Vector: inlined friend functions.
 
inline Vector operator-(const Point & p, const Point & q)
{
   if (p.w == 0 || q.w == 0)
      return Vector();
   else
      return Vector(p.x/p.w - q.x/q.w, p.y/p.w - q.y/q.w, p.z/p.w - q.z/q.w);
}
 
inline Vector operator+(const Vector & u, const Vector & v)
{
   return Vector(u.x + v.x, u.y + v.y, u.z + v.z);
}
 
inline Vector operator-(const Vector & u, const Vector & v)
{
   return Vector(u.x - v.x, u.y - v.y, u.z - v.z);
}
 
inline Vector operator*(const Vector & v, GLfloat s)
{
   return Vector(s * v.x, s * v.y, s * v.z);
}
 
inline Vector operator*(GLfloat s, const Vector & v)
{
   return Vector(s * v.x, s * v.y, s * v.z);
}
 
inline Vector cross(const Vector & u, const Vector & v)
{
   return Vector(
             u.y * v.z - u.z * v.y,
             u.z * v.x - u.x * v.z,
             u.x * v.y - u.y * v.x );
}
 
inline Vector operator*(const Vector & u, const Vector & v)
{
   return Vector(
             u.y * v.z - u.z * v.y,
             u.z * v.x - u.x * v.z,
             u.x * v.y - u.y * v.x );
}
 
inline GLfloat dot(const Vector & u, const Vector & v)
{
   return u.x * v.x + u.y * v.y + u.z * v.z;
}
 
inline bool operator==(const Vector & u, const Vector & v)
{
   return u.x == v.x && u.y == v.y && u.z == v.z;
}
 
inline bool operator!=(const Vector & u, const Vector & v)
{
   return !(u == v);
}
 
 
 
// Class Matrix: inlined constructors
 
inline Matrix::Matrix(const Plane & p)
{
   reflect(p);
}
 
inline Matrix::Matrix(const Point & lightPos, const Plane & plane)
{
   shadow(lightPos, plane);
}
 
inline Matrix::Matrix(
   GLfloat m00, GLfloat m10, GLfloat m20, GLfloat m30,
   GLfloat m01, GLfloat m11, GLfloat m21, GLfloat m31,
   GLfloat m02, GLfloat m12, GLfloat m22, GLfloat m32,
   GLfloat m03, GLfloat m13, GLfloat m23, GLfloat m33 )
{
   m[0][0] = m00;
   m[1][0] = m10;
   m[2][0] = m20;
   m[3][0] = m30;
   m[0][1] = m01;
   m[1][1] = m11;
   m[2][1] = m21;
   m[3][1] = m31;
   m[0][2] = m02;
   m[1][2] = m12;
   m[2][2] = m22;
   m[3][2] = m32;
   m[0][3] = m03;
   m[1][3] = m13;
   m[2][3] = m23;
   m[3][3] = m33;
}
 
// Class Matrix: inlined member functions
 
inline Point Matrix::apply(const Point & p) const
{
   return Point
          (
             m[0][0] * p.x + m[0][1] * p.y + m[0][2] * p.z + m[0][3] * p.w,
             m[1][0] * p.x + m[1][1] * p.y + m[1][2] * p.z + m[1][3] * p.w,
             m[2][0] * p.x + m[2][1] * p.y + m[2][2] * p.z + m[2][3] * p.w,
             m[3][0] * p.x + m[3][1] * p.y + m[3][2] * p.z + m[3][3] * p.w
          );
}
 
inline void Matrix::apply() const
{
   glMultMatrixf(&m[0][0]);
}
 
inline GLfloat *Matrix::get()
{
   return &m[0][0];
}
 
inline GLfloat & Matrix::operator()(int i, int j)
{
   return m[i][j];
}
 
inline const GLfloat & Matrix::operator()(int i, int j) const
{
   return m[i][j];
}
 
inline Matrix Matrix::operator~() const
{
   return inv();
}
 
inline Vector Matrix::apply(const Vector & v) const
{
   // Gives same result as quaternion.
   return Vector
          (
             m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2],
             m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2],
             m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2]
          );
}
 
inline Matrix & Matrix::operator/=(GLfloat s)
{
   m[0][0] /= s;
   m[0][1] /= s;
   m[0][2] /= s;
   m[0][3] /= s;
   m[1][0] /= s;
   m[1][1] /= s;
   m[1][2] /= s;
   m[1][3] /= s;
   m[2][0] /= s;
   m[2][1] /= s;
   m[2][2] /= s;
   m[2][3] /= s;
   m[3][0] /= s;
   m[3][1] /= s;
   m[3][2] /= s;
   m[3][3] /= s;
   return *this;
}
 
inline Matrix & Matrix::operator*=(GLfloat s)
{
   m[0][0] *= s;
   m[0][1] *= s;
   m[0][2] *= s;
   m[0][3] *= s;
   m[1][0] *= s;
   m[1][1] *= s;
   m[1][2] *= s;
   m[1][3] *= s;
   m[2][0] *= s;
   m[2][1] *= s;
   m[2][2] *= s;
   m[2][3] *= s;
   m[3][0] *= s;
   m[3][1] *= s;
   m[3][2] *= s;
   m[3][3] *= s;
   return *this;
}
 
// Class Matrix: inlined friend functions
 
inline Matrix Matrix::operator+=(const Matrix & rhs)
{
   for (int c = 0; c < 4; ++c)
      for (int r = 0; r < 4; ++r)
         m[c][r] += rhs(c,r);
   return *this;
}
 
inline Matrix operator+(const Matrix & m, const Matrix & n)
{
   Matrix sum(m);
   sum += n;
   return sum;
}
 
inline Matrix Matrix::operator-=(const Matrix & rhs)
{
   for (int c = 0; c < 4; ++c)
      for (int r = 0; r < 4; ++r)
         m[c][r] -= rhs(c,r);
   return *this;
}
 
inline Matrix operator-(const Matrix & m, const Matrix & n)
{
   Matrix sum(m);
   sum -= n;
   return sum;
}
 
inline Matrix Matrix::operator*=(const Matrix & rhs)
{
   *this = (*this) * rhs;
   return *this;
}
 
inline 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;
}
 
inline bool operator==(const Matrix & m, const Matrix & n)
{
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         if (m(i,j) != n(i,j))
            return false;
   return true;
}
 
inline bool operator!=(const Matrix & m, const Matrix & n)
{
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         if (m(i,j) != n(i,j))
            return true;
   return false;
}
 
inline Matrix operator*(GLfloat s, const Matrix & m)
{
   Matrix sm;
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         sm(i, j) = s * m(i, j);
   return sm;
}
 
inline Matrix operator*(const Matrix & m, GLfloat s)
{
   Matrix ms;
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         ms(i, j) = s * m(i, j);
   return ms;
}
 
inline Matrix operator/(const Matrix & m, GLfloat s)
{
   Matrix md = m;
   return md /= s;
}
 
 
// Class Quaternion: inlined member functions.
 
inline Quaternion Quaternion::operator+=(const Quaternion & q)
{
   s += q.s;
   v += q.v;
   return *this;
}
 
inline Quaternion Quaternion::operator-=(const Quaternion & q)
{
   s -= q.s;
   v -= q.v;
   return *this;
}
 
inline Quaternion & Quaternion::operator*=(GLfloat x)
{
   v *= x;
   s *= x;
   return *this;
}
 
inline Quaternion & Quaternion::operator/=(GLfloat x)
{
   v /= x;
   s /= x;
   return *this;
}
 
inline Vector Quaternion::apply(const Vector & w) const
{
   return Vector
          (
             -(-w[0] * v[0] - w[1] * v[1] - w[2] * v[2]) * v[0] + s * (s * w[0] + w[1] * v[2] - w[2] * v[1])
             - v[1] * (s * w[2] + w[0] * v[1] - w[1] * v[0]) + v[2] * (s * w[1] + w[2] * v[0] - w[0] * v[2]),
             -(-w[0] * v[0] - w[1] * v[1] - w[2] * v[2]) * v[1] + s * (s * w[1] + w[2] * v[0] - w[0] * v[2])
             - v[2] * (s * w[0] + w[1] * v[2] - w[2] * v[1]) + v[0] * (s * w[2] + w[0] * v[1] - w[1] * v[0]),
             -(-w[0] * v[0] - w[1] * v[1] - w[2] * v[2]) * v[2] + s * (s * w[2] + w[0] * v[1] - w[1] * v[0])
             - v[0] * (s * w[1] + w[2] * v[0] - w[0] * v[2]) + v[1] * (s * w[0] + w[1] * v[2] - w[2] * v[1])
          );
}
 
inline Vector Quaternion::vector() const
{
   return v;
}
 
inline GLfloat Quaternion::scalar() const
{
   return s;
}
 
inline GLfloat Quaternion::norm() const
{
   return s * s + v.norm();
}
 
inline GLfloat Quaternion::magnitude() const
{
   return GLfloat(sqrt(norm()));
}
 
inline Quaternion Quaternion::conj() const
{
   return Quaternion(s, -v);
}
 
inline Vector Quaternion::axis() const
{
   return v.unit();
}
 
inline double Quaternion::angle() const
{
   return 2 * acos(s);
}
 
inline Quaternion Quaternion::operator~() const
{
   return inv();
}
 
inline Quaternion & Quaternion::operator*=(const Quaternion & q)
{
   GLfloat ns = s * q.s - v[0] * q.v[0] - v[1] * q.v[1] - v[2] * q.v[2];
   v = Vector(
          q.s * v[0] + s * q.v[0] + v[1] * q.v[2] - v[2] * q.v[1],
          q.s * v[1] + s * q.v[1] + v[2] * q.v[0] - v[0] * q.v[2],
          q.s * v[2] + s * q.v[2] + v[0] * q.v[1] - v[1] * q.v[0] );
   s = ns;
   return *this;
 
}
 
// class Quaternion: inlined friend functions.
 
inline GLfloat dot(const Quaternion & q, const Quaternion & r)
{
   return q.s * r.s + dot(q.v, r.v);
}
 
inline Quaternion operator+(const Quaternion & q, const Quaternion & r)
{
   return Quaternion(q.s + r.s, q.v + r.v);
}
 
inline Quaternion operator-(const Quaternion & q, const Quaternion & r)
{
   return Quaternion(q.s - r.s, q.v - r.v);
}
 
inline Quaternion operator*(const Quaternion & q, const Quaternion & r)
{
   return Quaternion
          (
             q.s * r.s - q.v[0] * r.v[0] - q.v[1] * r.v[1] - q.v[2] * r.v[2],
             r.s * q.v[0] + q.s * r.v[0] + q.v[1] * r.v[2] - q.v[2] * r.v[1],
             r.s * q.v[1] + q.s * r.v[1] + q.v[2] * r.v[0] - q.v[0] * r.v[2],
             r.s * q.v[2] + q.s * r.v[2] + q.v[0] * r.v[1] - q.v[1] * r.v[0]
          );
}
 
inline Quaternion operator*(const Vector & v, const Quaternion & q)
{
   return Quaternion
          (
             - v[0] * q.v[0] - v[1] * q.v[1] - v[2] * q.v[2],
             q.s * v[0] +   v[1] * q.v[2] - v[2] * q.v[1],
             q.s * v[1] +   v[2] * q.v[0] - v[0] * q.v[2],
             q.s * v[2] +   v[0] * q.v[1] - v[1] * q.v[0]
          );
}
 
inline Quaternion operator*(const Quaternion & q, const Vector & v)
{
   return Quaternion
          (
             - q.v[0] * v[0] - q.v[1] * v[1] - q.v[2] * v[2],
             q.s * v[0] + q.v[1] * v[2] - q.v[2] * v[1],
             q.s * v[1] + q.v[2] * v[0] - q.v[0] * v[2],
             q.s * v[2] + q.v[0] * v[1] - q.v[1] * v[0]
          );
}
 
inline Quaternion & Quaternion::operator/=(const Quaternion & q)
{
   GLfloat den = q.norm();
   GLfloat ns = (s * q.s + v[0] * q.v[0] + v[1] * q.v[1] + v[2] * q.v[2]) / den;
   v = Vector
       (
          (q.s * v[0] - s * q.v[0] - v[1] * q.v[2] + v[2] * q.v[1]) / den,
          (q.s * v[1] - s * q.v[1] - v[2] * q.v[0] + v[0] * q.v[2]) / den,
          (q.s * v[2] - s * q.v[2] - v[0] * q.v[1] + v[1] * q.v[0]) / den
       );
   s = ns;
   return *this;
}
 
inline Quaternion operator*(const Quaternion & q, GLfloat a)
{
   return Quaternion(a * q.s, a * q.v);
}
 
inline Quaternion operator*(GLfloat a, const Quaternion & q)
{
   return Quaternion(a * q.s, a * q.v);
}
 
inline bool operator==(const Quaternion & q, const Quaternion & r)
{
   return q.s == r.s && q.v == r.v;
}
 
inline bool operator!=(const Quaternion & q, const Quaternion & r)
{
   return !(q == r);
}
 
 
 
 
 
// class Interpolator: inline member functions
 
inline void Interpolator::getMatrix(double t, GL_Matrix m) const
{
   getQuaternion(t).matrix(m);
}
 
}
; // end of namespace
 
#endif