#ifndef CUGL_H
#define CUGL_H

/** \file cugl.h
 *  Concordia University Graphics Library
 *
 * \authors Peter Grogono
 * \date November 2005
 */

/** \mainpage   Introduction

CUGL is a library intended to be used with OpenGL.  Unlike most other such libraries,
CUGL is not intended to provide an abstraction layer over OpenGL, with classes, for example.
Instead, it complements OpenGL by providing capabilities that are often useful in graphics
programs.

CUGL defines classes for mathematical objects such as points, lines, planes, vectors, 
matrices, and quaternions.  Each class has its own reasonably complete set of operations,
as well as operations that relate it to other classes.  For example, two points define a
line and a plane defines a vector normal to it.

\section rep Representation

\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.

\section eq Equality

Classes \c Point, \c Line, \c Plane, \c Vector, \c Matrix, and \c Quaternion provide
two tests for equality.  \c operator== tests for \b exact equality.  If the coordinates
of a point, for example, are calculated in two different ways, giving \c Points \c p and \c q,
then \c p==q will probably yield \c false.

Each of these classes provides another function, \c eq, which tests for \b approximate
equality: values are considered equal if they agree to about six decimal places.
The \c eq test is slower than \c operator==, but may be more useful in some situations.

The functions \c operator!= and \c neq are the negations of \c operator== and \c eq, respectively.

Structured objects are comparaed by comparing their components.  For example, planes
\a (a,b,c,d) and \a (a',b',c'd') are considered equal if \a a=a', \a b=b', \a c=c', and \a d=d'.
This method may not be the most precise in a geometric sense.

\section insert Insertion operators

\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.

\section cons Copy constructors and assignment operators

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.)

\section templates Why doesn't CUGL use templates?

It is tempting to convert the CUGL classes Point, Line, Plane, Vector, Matrix, and Quaternion
to template classes.  The reason that I have not done this is that OpenGL does not use C++
capabilities and, in particular, does not use templates.  It seems difficult to make template
classes work with OpenGL functions that incorporate the type into the function name.  RTTI
might help but making all decision at compile-time requires template meta-programming.

\section ack Acknowledgements  

Parts of CUGL are based on publications by Mark Kilgard, Ken Shoemake, F.S. Hill, and others.
Thanks also to Liping Ye and Qingzhe Huang for significant contributions.

*/

/** \defgroup functions Miscellaneous Functions  */

/** \defgroup errors Error Messages           */

/** \defgroup materials Materials and Models     */

/** \defgroup constants Constants                */

/** \defgroup examples Example Programs         */

// The following definitions avoid the need to read windows.h.
// Adapted from glut.h

#if defined(_WIN32)
# if 0
   /* This would put tons of macros and crap in our clean name space. */
#  define  WIN32_LEAN_AND_MEAN
#  include <windows.h>
# else
   /* XXX This is from Win32's <windef.h> */
#  ifndef APIENTRY
#   define GLUT_APIENTRY_DEFINED
#   if (_MSC_VER >= 800) || defined(_STDCALL_SUPPORTED) || defined(__BORLANDC__) || defined(__LCC__)
#    define APIENTRY    __stdcall
#   else
#    define APIENTRY
#   endif
#  endif
   /* XXX This is from Win32's <winnt.h> */
#  ifndef CALLBACK
#   if (defined(_M_MRX000) || defined(_M_IX86) || defined(_M_ALPHA) || defined(_M_PPC)) && !defined(MIDL_PASS) || defined(__LCC__)
#    define CALLBACK __stdcall
#   else
#    define CALLBACK
#   endif
#  endif
   /* XXX Hack for lcc compiler.  It doesn't support __declspec(dllimport), just __stdcall. */
#  if defined( __LCC__ )
#   undef WINGDIAPI
#   define WINGDIAPI __stdcall
#  else
   /* XXX This is from Win32's <wingdi.h> and <winnt.h> */
#   ifndef WINGDIAPI
#    define GLUT_WINGDIAPI_DEFINED
#    define WINGDIAPI __declspec(dllimport)
#   endif
#  endif
   /* XXX This is from Win32's <ctype.h> */
#  ifndef _WCHAR_T_DEFINED
typedef unsigned short wchar_t;
#   define _WCHAR_T_DEFINED
#  endif
# endif
#endif

#include <GL/glu.h>
#include <iostream>
#include <string>
#include <cmath>

namespace cugl
{

/**
 * \addtogroup constants
 * \{
 */

/**
 * Indicates version and last compile date.
 */
const char version[] = "CUGL V2 2005.11.18";

/**
 * A well-known mathematical constant.
 */
const double HALF_PI = 2 * atan(1.0);

/**
 * A well-known mathematical constant.
 */
const double PI = 4 * atan(1.0);

/**
 * A well-known mathematical constant.
 */
const double TWO_PI = 8 * atan(1.0);

/**
 * This constant is used by the predicate \c eq when comparing
 * two numbers for equality or inequality.  See \ref eq.
 */
const double EQUALITY_TOLERANCE = 1e-6;

// \}

/**
 * \defgroup typedefs Type definitions
 * \{
 */

/**
 * The type of OpenGL projection and model view matrices.
 */
typedef GLfloat GL_Matrix[4][4];

/** \} */

/**
 * \addtogroup errors
 * \{
 */

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

/**
 * Return a message describing an error.
 */
const char *getErrorString(CUGLErrorType err);

/**
 * Display a message if a 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.
 *
 * \note 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.
 *
 * \note This function is intended mainly for development and
 * should not normally be included in a production program.
 */
void checkOpenGLStatus();

/** \} */

/**
 * \addtogroup functions
 * \{
 */

/**
 * Compare two numbers for approximate equality.  See \ref eq.
 * The precision of the comparison is determined by the value of \c EQUALITY_TOLERANCE.
 * The numbers are consider equal if either one is zero and the other is less than
 * \c EQUALITY_TOLERANCE or if x/y differs from 1 by less than \c EQUALITY_TOLERANCE.
 */
template <typename T>
bool testEqual(T x, T y)
{
    if (x == 0)
        return fabs(y) < EQUALITY_TOLERANCE;
    else if (y == 0)
        return fabs(x) < EQUALITY_TOLERANCE;
    else
        return fabs(x/y - 1) < EQUALITY_TOLERANCE;
}

/** \} */

// Forward declarations.
class Vector;
class Quaternion;

/** 
 * An instance is a point represented by four homogeneous coordinates.
 * A point is represented by \c (x,y,z,w) in which each component is a \c GLfloat.
 * If \c 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 \c w=1.  A Point may be normalized if \c 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 the array \c 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.
    */
   explicit Point (const Quaternion & q);

   /**
    * Access coordinates: p[0] = x, p[1] = y, p[2] = z, p[3] = w.
    * Error BAD_INDEX reported if \c i is not one of {0,1,2, 3}.
    */
   GLfloat & operator[](int i);

   /**
    * Access coordinates: p[0] = x, p[1] = y, p[2] = z, p[3] = w.
    * This form is used with \c const instances (\c p[i] cannot appear on the left of \c =).
    * Error BAD_INDEX reported if \c i is not one of {0,1,2, 3}.
    */
   const GLfloat & operator[](int i) const;

   /**
    * Normalize this point.
    * A Point is in normal form if \c w=1.
    * Error ZERO_DIVISOR reported if \c w=0.
    * \note The value of this Point is changed by this function.
    */
   void normalize();

   /**
    * Return this Point in normalized form.
    * A Point is in normal form if \c w=1.
    * Error ZERO_DIVISOR reported if \c 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:
    * it must be one of \c GL_LIGHT0, \c GL_LIGHT1, ...
    */
   void light(GLenum lightNum) const;

   /**
    * Apply a translation to the current OpenGL matrix.
    * The translation moves the origin to this point.
    * \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 \c s.
    * This is implemented by multiplying the \c w component of the Point
    * by \c s.  If \c s=0, the result is a Point at infinity.
    * \return the Point at \c (x,y,z,s*w).
    */
   Point operator/(GLfloat s) const;

   /**
    * Use a Vector to move a Point.
    * The new coordinates of this Point are (x+w*v.x, y+w*v.y, z+w*v.z, w).
    */
   Point & operator+=(const Vector & v);

   /**
    * Use a Vector to move a Point.
    * The new coordinates of this Point are \c (x-w*v.x, y-w*v.y, z-w*v.z, w).
    */
   Point & operator-=(const Vector & v);

   /**
    * Scale a point.
    */
   friend Point operator*(const Point & p, GLfloat s);

   /**
    * Scale a 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(Line & k, Plane & p);

   /**
    * Compare two points.
    * This function returns \c true only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator==(const Point & p, const Point & q);

   /**
    * Compare two points.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Point & p, const Point & q);

  /**
   * Compare two points approximately.
   * See \ref eq.
   * \return \c true if the points are close and \c false otherwise.
   */
  friend bool eq(const Point & p, const Point & q);

  /**
   * Compare two points approximately.
   * See \ref eq.
   * \return \c false if the points are close and \c true otherwise.
   */
  friend bool neq(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.
    * \param os is the stream to write to.
    * \param p is the Point to be written.
    */
   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);

private:
   GLfloat x, y, z, 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(Point & p, Point & q);

   /**
    * Construct the Line from Point \c p to Point \c p+v.
    */
   Line(Point & p, Vector & v);

   /**
    * Find the point where this line meets the plane p.
    */
   friend Point meet(Line & k, 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.
   * See \ref eq.
    */
   friend bool operator==(const Line & x, const Line & y);

   /**
    * Compare two lines.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Line & x, const Line & y);

  /**
   * Compare two lines approximately.
   * See \ref eq.
   * \return \c true if the lines are close and \c false otherwise.
   */
  friend bool eq(const Line & p, const Line & q);

  /**
   * Compare two lines approximately.
   * See \ref eq.
   * \return \c false if the lines are close and \c true otherwise.
   */
  friend bool neq(const Line & p, const Line & q);

   /**
    * Write a representation of the line to the output stream.
    * The format is \c pt->pt where "\c pt" is the format used for points.
    * If \c setw is used to set the width, it is passed to the inserter
    * for Point.
    * \param os is the stream to write to.
    * \param k is the Line to be written.
    */
   friend std::ostream & operator<<(std::ostream & os, const Line & k);

private:
   Point s; // start point
   Point f; // finish point
};


/**
 * An instance is a plane defined by its equation \c Ax+By+Cx+D=0.
 *
 * Planes are used for clipping, shadows, and reflections
 * (see class Matrix).
 *
 * A homogeneous point (x,y,z,w) lies on the plane if Ax+By+Cz+Dw=0.
 * The notation for a plane is \c <A,B,C,D>.  The plane \c <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, \c y=0).  
 *
 * A Plane is in normal form if the Vector \c (A,B,C) is a unit vector. 
 * A Plane can be normalized by scaling \c A,B,C,D: see Plane::normalize().
 */

class Plane
{
   friend class Matrix;

public:
   /**
    * Construct a plane given the coefficients of its equation \c Ax+By+Cx+D=0.
    * If no arguments are provided, construct the plane \c <0,1,00>, corresponding to \c y=0.
    */
   Plane(GLfloat a = 0, GLfloat b = 1, GLfloat c = 0, GLfloat d = 0);

  /**
   * Construct a plane containing the three given points.
   */
   Plane(Point & p, Point & q, Point & r);

   /**
    * Construct a plane containing the line \c s and the point \c p.
    */
   Plane(Line & s, Point & p);

   /**
    * Construct the plane orthogonal to \c v and passing through \c p.
    */
   Plane(Vector & v, Point & p);

   /**
    * Normalize this plane.
    * The Plane \c (A,B,C,D) is in normal form when \c (A,B,C) is a unit vector.
    * Error ZERO_DIVISOR reported if \c |(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 \c (A,B,C,D) is in normal form when \c (A,B,C) is a unit vector.
    * Error ZERO_DIVISOR reported if \c |(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(Line & k, Plane & p);

   /**
    * Compare two planes.
    * This function returns \c true only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator==(const Plane & x, const Plane & y);

   /**
    * Compare two planes.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Plane & x, const Plane & y);

  /**
   * Compare two pplanes approximately.
   * See \ref eq.
   * \return \c true if the planes are close and \c false otherwise.
   */
  friend bool eq(const Plane & p, const Plane & q);

  /**
   * Compare two planes approximately.
   * See \ref eq.
   * \return \c false if the planes are close and \c true otherwise.
   */
  friend bool neq(const Plane & p, const Plane & q);

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

   /**
    * Write a description of the plane to the output stream as \c <a,b,c,d>.
    * Inserts the components of the plane into the output 
    * stream \c os in the format \c <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.
    * \param os is the stream to write to.
    * \param p is the Plane to be written.
    */
   friend std::ostream & operator<<(std::ostream & os, const Plane & p);

private:
   GLfloat a, b, c, d;
};


/**
 * An instance is a vector with 3 real components.
 * A vector v is a 3-vector represented by three orthogonal components 
 * \c vx, \c vy, and \c vz.  The norm of the vector \c 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:
   /**
    * The default constructor constructs the zero Vector: \c (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 \c points.
    * \note The vector returned is \b not a unit vector.
    */
   Vector(Point points[], int numPoints);

   /**
    * Construct a vector from two points.
    * Vector(p, q) is equivalent to p - q.
    */
   Vector(Point & p, 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 \c 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 \c 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 \c 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 apply a translation to the current OpenGL matrix.
    */
   void translate() const;

   /**
    * Use this vector to apply a rotation to the current OpenGL matrix.
    * The direction of the vector provides the axis of rotation.
    * \param angle determines the amount of rotation in degrees.
    */
   void rotate(GLfloat angle) 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;

   /**
    * Draw this vector as a line in the graphics window.
    * The line joins \c P and \c P+V, where \c P is the point provided.
    * \param p is the start point for the vector.
    */
   void draw(const Point & p) const;

   /**
    * Get or set a reference to a component of this vector.
    * \c v[0] is the \c x component; \c v[1] is the \c y component; \c v[2] is the \c 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 \c x component; \c v[1] is the \c y component; \c v[2] is the \c 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.
    * See also Vector::operator*().
    */
   friend Vector cross(const Vector & u, const Vector & v);

   /**
    * Return cross product of vectors \c u and \c v.
    * See also Vector::cross().
    */
   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 \c s*v.
    */
   friend Vector operator*(const Vector & v, GLfloat s);

   /**
    * Return Vector \c s*v.
    */
   friend Vector operator*(GLfloat s, const Vector & v);

   /**
    * Return Vector \c u+v.
    */
   friend Vector operator+(const Vector & u, const Vector & v);

   /**
    * Return Vector \c 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);

   /**
    * 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.
   * See \ref eq.
    */
   friend bool operator==(const Vector & x, const Vector & y);

   /**
    * Compare two vectors.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Vector & x, const Vector & y);

  /**
   * Compare two vectors approximately.
   * See \ref eq.
   * \return \c true if the vectors are close and \c false otherwise.
   */
  friend bool eq(const Vector & p, const Vector & q);

  /**
   * Compare two vectors approximately.
   * See \ref eq.
   * \return \c false if the vectors are close and \c true otherwise.
   */
  friend bool neq(const Vector & p, const Vector & q);

   /**
    * Write vector to output stream as \c (x,y,z).
    * Inserts the components of the vector into the output 
    * stream \c os in the format \c (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.
    * \param os is the stream to write to.
    * \param v is the Vector to be written.
    */
   friend std::ostream & operator<<(std::ostream & os, const Vector & v);

private:
   GLfloat x, y, z;
};

/**
 * \addtogroup constants
 * \{
 */

/**
 * Unit vector parallel to \c X axis.
 */
const Vector I(1, 0, 0);

/**
 * Unit vector parallel to \c Y axis.
 */
const Vector J(0, 1, 0);


/**
 * Unit vector parallel to \c 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(const 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 \c u, will yield \c v.
    * \pre The vectors \c u and \c v must be unit vectors.
    */
   Matrix(const Vector & u, const Vector & v); 

   /**
    * 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 determinant of this matrix.
    */
   GLfloat det() const;

   /**
    * Return the inverse of this Matrix.
    * The function uses Cramer's method: compute the matrix of cofactors
    * and divide each element by the determinant.
    * If the Matrix is singular, report \c SINGULAR_MATRIX and return the zero Matrix.
    * \return the inverse of this Matrix.
    */
   Matrix inv() const;

   /**
    * Multiply the current OpenGL matrix by this Matrix.
    */
   void apply() const;

   /**
    * Apply this Matrix to a Point.
    * Applying a Matrix to a Point is not really a meaningful operation.
    * The result is a Point \c q such that the Matrix transforms \c op to \c oq,
    * where \c o is the origin of the coordinate system.
    */
   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 two matrices.
    * \return \c lhs+rhs.
    */
   friend Matrix operator+(const Matrix & lhs, const Matrix & rhs);

   /**
    * Subtract two matrices.
    * \return \c lhs-rhs.
    */
   friend Matrix operator-(const Matrix & lhs, const Matrix & rhs);

   /**
    * Multiply two matrices.
    * \return \c lhs*rhs.
    */
   friend Matrix operator*(const Matrix & lhs, const Matrix & rhs);

   /**
    * Divide two matrices.
    * Matrixc "division" is performed by multiplying the first amtrix by the inverse of the second Matrix.
    * \return \c lhs*rhs.inv().
    */
   friend Matrix operator/(const Matrix & lhs, const Matrix & rhs);

   /**
    * Add a Matrix to this Matrix.
    */
   Matrix & operator+=(const Matrix & rhs);

   /**
    * Subtract a Matrix from this Matrix.
    */
   Matrix & operator-=(const Matrix & rhs);

   /**
    * Multiply this Matrix by another Matrix.
    */
   Matrix & operator*=(const Matrix & rhs);

   /**
    * Divide this Matrix by another Matrix.
    */
   Matrix & operator/=(const Matrix & rhs);


   /**
    * Compare two matrices.
    * This function returns \c true only if corresponding components are exactly equal.
    * See \ref eq.
    */
   friend bool operator==(const Matrix & x, const Matrix & y);

   /**
    * Compare two matrices.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Matrix & x, const Matrix & y);

  /**
   * Compare two matrices approximately.
   * See \ref eq.
   * \return \c true if the matrices are close and \c false otherwise.
   */
  friend bool eq(const Matrix & p, const Matrix & q);

  /**
   * Compare two matrices approximately.
   * See \ref eq.
   * \return \c false if the matrices are close and \c true otherwise.
   */
  friend bool neq(const Matrix & p, const Matrix & q);

   /**
    * 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.
    * \param os is the stream to write to.
    * \param m is the Matrix to be written.
    */
   friend std::ostream & operator<<(std::ostream & os, const Matrix & m);

private:
   GL_Matrix m;
};



/**
 * An instance is a quaternion represented as a (scalar, vector) pair.
 * We use the notation \c (s,v) for a quaternion with scalar component
 *   \c s and vector component \c 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 \c (0,v).
    * The result will be a unit quaternion if \c v is a unit vector,
    * in which case the quaternion represents a rotation through
    * 90 degrees about the axis \c v.
    */
   explicit Quaternion(const Vector & v) : s(0), v(v) {}

   /**
    * Construct the quaternion \c (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.
    */
   explicit Quaternion(Matrix m);

   /**
    * Construct a quaternion from Euler angles.
    * \param xr is the Euler angle for the X axis.
    * \param yr is the Euler angle for the Y axis.
    * \param zr is the Euler angle for the Z axis.
    */
   Quaternion(double xr, double yr, double zr);

   /**
    * Construct a Quaternion from a Point.
    * The quaternion consists of the vector \c (p[0],p[1],p[2]) and the scalar \c p[3].
    * \note This is a unusual operation and should not normally be used.
    * It is intended for experiments with non-linear transformations.
    */
   explicit 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 \c u, will yield \c v.
    * \pre The vectors \c u and \c v must be unit vectors.
    */
   Quaternion(const Vector & u, const Vector & v); 

   /**
    * Add another Quaternion to this Quaternion.
    */
   Quaternion & operator+=(const Quaternion & q);

   /**
    * Subtract another Quaternion from this Quaternion.
    */
   Quaternion & operator-=(const Quaternion & q);


   /**
    * Promote the vector \c v to a quaternion \c qv and 
    * return the quaternion product \c qv*q.
    */
   friend Quaternion operator*(const Vector & v, const Quaternion & q);

   /**
    * Promote the vector \c v to a quaternion \c qv and 
    * return the quaternion product \c q*qv.
    */
   friend Quaternion operator*(const Quaternion & q, const Vector & v);

   /**
    * Return the quaternion ratio \c q/r.
    * If q and r represent 
    * rotations, then \c q/r represents the rotation \c q followed by
    * the rotation that is the inverse of \c r.
    * \return the Quaternion ratio \c 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 \c a*q, where \c a is a scalar.
    * \c a is a scalar and a*(s,v) = (a*s, a*v).
    * \return the Quaternion \c a*q. 
    */
   friend Quaternion operator*(const Quaternion & q, GLfloat a);

   /**
    * Return the quaternion \c a*q, where \c a is a scalar.
    * \c a is a scalar and a*(s,v) = (a*s, a*v).
    * \return the Quaternion \c a*q. 
    */
   friend Quaternion operator*(GLfloat a, const Quaternion & q);

   /**
    * Return the quaternion \c q/a, where \c a is a scalar.
    * \c a is a scalar and (s,v)/a = (s/a, v/a).
    * Report error ZERO_DIVISOR if \c a = 0.
    * \return the Quaternion \c 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.
    * \return the inverse of this quaternion.
    */
   Quaternion inv() 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 \c v of the quaternion \c q = \c (s,v).
    * The same effect can be achieved with the constructor \c Vector::Vector(q).
    */
   Vector vector() const;

   /**
    * Return Scalar component \c s of the quaternion \c q = \c (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(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 \c x, \c y, \c z order.
    * \param xr is the Euler angle for the X axis.
    * \param yr is the Euler angle for the Y axis.
    * \param zr is the Euler angle for the Z axis.
    * \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 \c 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 \c (x1,y1) to \c (x2,y2).  The coordinates should be normalized 
    * so that both \c x and \c 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);

   /**
    * 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, \c ln(q1)+ln(q2) makes sense only if \c q1 and \c 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, \c exp(v1+v2)=exp(v1)*exp(v2) only if \c v1 and \c 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 \c (s,u) and \c (t,v) is \c st+dot(u,v)
    * where \c dot(u,v) denotes the vector dot product of \c u and \c 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.
   * See \ref eq.
    */
   friend bool operator==(const Quaternion & x, const Quaternion & y);

   /**
    * Compare two quaternions.
    * This function returns \c false only if corresponding components are exactly equal.
   * See \ref eq.
    */
   friend bool operator!=(const Quaternion & x, const Quaternion & y);

  /**
   * Compare two quaternions approximately.
   * See \ref eq.
   * \return \c true if the quaternions are close and \c false otherwise.
   */
  friend bool eq(const Quaternion & p, const Quaternion & q);

  /**
   * Compare two quaternions approximately.
   * See \ref eq.
   * \return \c false if the quaternions are close and \c true otherwise.
   */
  friend bool neq(const Quaternion & p, const Quaternion & q);

   /**
    * Write the quaternion to the output stream using the format "s (x,y,z)".
    * \param os is the stream to write to.
    * \param q is the Quaternion to be written.
    */
   friend std::ostream & operator<<(std::ostream & os, const Quaternion & q);

private:
    /**
     * The scalar component of the Quaternion.
     */
    GLfloat s;

    /**
     * The Vector component of the Quaternion.
     */
    Vector v;
};

/**
* An abstract base class for interpolation.
* Interpolators are intended to capture the general idea of a smooth
* transition between two configurations.  Examples include a translation
* (linear movement betwen two points), a rotation (smooth turning about 
* an axis), and camera motions (pan, tilt, track, etc.).
*
* The base class provides the essential functions for initializing and
* continuing a transition.  Derived classes add details for particular
* kinds of interpolation.
*
* The transition is defined by a sequence of steps, 0,1,2,...,\c maxSteps.
* This sequence is converted to a parameter, \c t, that varies from 0 to 1
* during the transition.  \c begin() initializes the transition (sets to 0),
* \c idle() increments the step counter, and \c atEnd() returns \c true when
* the step counter attains its maximum value.  As its name implies, \c idle()
* is usually called from within the GLUT idle function.
*
* The pure virtual function \c apply() is defined in derived classes and
* determines the effect of the interpolator: that is, a translation, rotation,
* camera position, etc.  It is usually called from the GLUT display function.

*/
class Interpolator
{
public:

    /**
    * Construct an interpolator.
    * \param maxSteps determines the number of interpolation steps.
    * \param repeat determines whether the \a idle() function will repeat the motion.
    */
    Interpolator(unsigned int maxSteps = 100, bool repeat = false) : maxSteps(maxSteps), step(0), repeat(repeat) {}

    /**
    * A pure virtual function that must be defined in derived classes to provide the effect of the motion.
    */
    virtual void apply() = 0;

    /**
    * Set the interpolator to an arbitrary position in its cycle.
    * By default, the position is the initial position.
    * \param s gives the step number for the position, defaults to zero.
    * \pre 0 <= s <= maxSteps.
    */
    void begin(unsigned int s = 0)
    {
        step = s;
    }

    /**
    * Set the total number of steps in a cycle.
    */
    void setSteps(int max) 
    { 
        maxSteps = max; 
    }

    /**
    * Test for the end of a cycle.
    * \return \a true if the cycle has ended, \a false otherwise.
    */
    bool atEnd() const
    {
        return step == maxSteps;
    }

    /**
    * Move to the next step.
    * If \a repeat is \a true, reinitialize at the end of each cycle, giving repetitive motion.
    * This function is usually called from within the GLUT \a idle() function.
    */
    void idle() 
    {
        if (step < maxSteps)
            ++step;
        else if (repeat)
            begin();
    }

protected:

    /**
    * Enable derived classes to obtain the current position in the cycle.
    * \return a \a double value in the range [0,1].
    */
    double getParam() const
    { 
        return double (step) / double(maxSteps); 
    }

private:

    /**
    * The total number of steps in a cycle.
    */
    unsigned int maxSteps;

    /**
    * The current step in the cycle.
    */
    unsigned int step;

    /**
    * Determines whether cycles are repeated.
    */
    bool repeat;
};

/**
 * A class for rotation interpolation.
 * This class uses quaternions to provide smooth interpolation of rotations.
 * See class Interpolator for an explanation of how to use it.
 */
class Rotation : public Interpolator
{
public:

    /**
     * Construct a rotation with default endpoints \c I and \c J.
     * \param maxSteps is the number of steps for the rotation.
     * \param repeat determines \c idle() will reset parameters at the end of a transition.
     */
    Rotation(unsigned int maxSteps = 100, bool repeat = false);

    /**
     * Construct a rotation with given end points.
     * \param first is a unit vector giving the initial direction.
     * \param last is a unit vector giving the final direction.
     * \param maxSteps is the number of steps for the rotation.
     * \param repeat determines \c idle() will reset parameters at the end of a transition.
     * \pre length(first) = length(last) = 1.
     */
    Rotation(const Vector & first, const Vector & last, unsigned int maxSteps = 100, bool repeat = false);

    /**
     * Set the starting and finishing directions for the rotation.
     * \param first is a unit vector giving the initial direction.
     * \param last is a unit vector giving the final direction.
     * \pre length(first) = length(last) = 1.
     */
    void set(const Vector & first, const Vector & last);

    /**
     * Set the final direction for the rotation.
     * The current direction is used as the initial direction of the new rotation.
     * This implies that, if this function is called during a rotation, the new
     * rotation will start immediately from the current position, without a "jump"
     * in the motion.
     * \param last is a unit vector giving the final direction.
     * \pre length(last) = 1.
     */
    void set(const Vector & last);

    /**
     * Apply the current rotation to the OpenGL model view matrix.
     */
    void apply();

    /**
     * \return the quaternion corresponding to the current rotation.
     */
    Quaternion get() const;

private:

    /**
     * Compute values for the constants required for the rotation calculations.
     * \param first is a unit vector giving the initial direction.
     * \param last is a unit vector giving the final direction.
     * \pre length(first) = length(last) = 1.
     */
    void initialize(const Vector & first, const Vector & last);

    /**
     * A unit vector whose direction is the initial direction of the rotation.
     */
    Vector vFirst;

    /**
     * A quaternion corresponding to the initial direction of the rotation.
     */
    Quaternion first;

    /**
     * A quaternion corresponding to the final direction of the rotation.
     */
    Quaternion last;

    /**
     * The angle between the initial and final directions.
     */
    double omega;

    /**
     * The sine of the angle between the initial and final directions.
     */
    double sinOmega;
};

/**
 * A class for translation interpolation.
 * A translation interpolation is represented as a Point that moves.
 * See class Interpolator for an explanation of how to use it.
 */
class Translation : public Interpolator
{
public:

    /**
     * Construct a default translation.
     * The initial and final points are both at the origin.
     * \param maxSteps is the number of steps for the rotation.
     * \param repeat determines \c idle() will reset parameters at the end of a transition.
     */
    Translation(int maxSteps = 100, bool repeat = false);

    /**
     * Construct a translation with given starting and finishing points.
     * \param first is the initial Point of the translation.
     * \param last is the final Point of the translation.
     * \param maxSteps is the number of steps for the rotation.
     * \param repeat determines \c idle() will reset parameters at the end of a transition.
     */
    Translation(const Point & first, const Point & last, int maxSteps = 100, bool repeat = false);

    /**
     * Set the starting and finishing points of a translation.
     * \param first is the initial Point of the translation.
     * \param last is the final Point of the translation.
     */
    void set(const Point & first, const Point & last);

    /**
     * Set the finishing point of a translation.
     * The current position is used as the initial position of the new rotation.
     * This implies that, if this function is called during a translation, the new
     * translation will start immediately from the current position, without a "jump"
     * in the motion.
     * \param pNew is a point that defines the final position.
     */
    void set(const Point & pNew);

    /**
     * The current position is passed to \c glTranslate().
     */
    void apply();

    /**
     * Return the Point corresponding to the current position of the translation.
     */
    Point get() const;

private:

    /**
     * The initial point of the translation.
     */
    Point first;

    /**
     * The final point of the translation.
     */
    Point last;
};

/**
 * A class for scalar interpolation.
 * A translation interpolation is represented as a double value that changes.
 * See class Interpolator for an explanation of how to use it.
 */
class Scaling : public Interpolator
{
public:

    /**
     * Construct a default scaling.
     * The initial value is 0 and the final value is 1.
     * \param maxSteps is the number of steps for the scaling.
     * \param repeat determines \c idle() will reset parameters at the end of a scaling.
     */
    Scaling(int maxSteps = 100, bool repeat = false);

    /**
     * Construct a translation with given starting and finishing values.
     * \param first is the initial value of the scaling.
     * \param last is the final value of the scaling.
     * \param maxSteps is the number of steps for the rotation.
     * \param repeat determines \c idle() will reset parameters at the end of a scaling.
     */
    Scaling(double first, double last, int maxSteps = 100, bool repeat = false);

    /**
     * Set the starting and finishing points of a scaling.
     * \param first is the initial value of the scaling.
     * \param last is the final value of the scaling.
     */
    void set(double first, double last);

    /**
     * Set the finishing point of a scaling.
     * The current value is used as the initial position of the new scaling.
     * This implies that, if this function is called during a scaling, the new
     * scaling will start immediately from the current position, without a "jump"
     * in the motion.
     * \param vNew is the nwe final value.
     */
    void set(double vNew);

    /**
     * Return the current value.
     */
    double get() const;

private:

    /**
     * The initial value of the scaling.
     */
    double first;

    /**
     * The final value of the scaling.
     */
    double last;
};

/**
 * A class for camera simulation.
 * This class provides standard "classical" camera movements: tracking
 * forwards/backwards, left/right, and up/down; panning left and right,
 * and tilting up and down.  The vertical axis remains vertical throughout
 * all movements.  The base class Interpolator provides smooth motion.
 *
 * The camera position is defined by an eye position (represented by a Point)
 * and a direction (represented by a Quaternion).
 */
class Camera : public Interpolator
{
public:

    /**
     * Construct a camera.
     * Since all parameters have default values, this is a default constructor.
     * \param eye is the initial position of the camera.
     * \param dir is the initial direction of the camera.
     * \param up is the initial up vector for the camera.
     * \param distance is the distance between the camera and the model.
     * \param maxSteps is the number of steps used for a transition.
     * \param repeat is \c true if the motion is to be repeated.
     */
    Camera(  Point eye = Point(), 
             Vector dir = K, 
             Vector up = J, 
             GLfloat distance = 1, 
             int maxSteps = 100, 
             bool repeat = false ) : 
       eyeFirst(eye), 
       eyeLast(eye), 
       dirFirst(dir), 
       dirLast(dir), 
       up(J),
       distance(distance),
       Interpolator(maxSteps, repeat)  
     {
     }

    /**
     * Apply a transformation that represents the camera position and orientation  to the current OpenGL model view matrix.
     * The transformation is implemented by calling \c gluLookAt.
     */
    void apply();

    /**
     * Set the distance from the camera position to the point of interest in the model.
     * By default, this distance is 1.
     */
    void set(GLfloat d) { distance = d; }

    /**
     * Set the position of the camera.
     * \param p is the new position of the camera.
     */
    void set(Point p);

    /**
     * Set the direction of the camera.
     * \param q is a quaternion specfying the rotation required to move the camera
     * from its current orientation to the desired orientation.  For example, to
     * pan righthrough an angle alpha, call \c camera.set(Quaternion(up,alpha)).
     */
    void set(Quaternion q);

    /**
     * Move the camera forwards, that is, towards the model.
     * \param distance is the distance moved.
     */
    void moveForward(GLfloat distance);

    /**
     * Move the camera backwards, that is, away from the model.
     * \param distance is the distance moved.
     */
    void moveBackward(GLfloat distance);

    /**
     * Move the camera vertically upwards, in the direction of the up vector.
     * \param distance is the distance moved.
     */
    void moveUp(GLfloat distance);

    /**
     * Move the camera vertically downwards, in the direction opposite to the up vector.
     * \param distance is the distance moved.
     */
    void moveDown(GLfloat distance);

    /**
     * Move the camera to the left.
     * The direction of motion makes a right angle with the line from
     * the camera to the model and with the up vector.
     * \param distance is the distance moved.
     */
    void moveLeft(GLfloat distance);

    /**
     * Move the camera to the right.
     * The direction of motion makes a right angle with the line from
     * the camera to the model and with the up vector.
     * \param distance is the distance moved.
     */
    void moveRight(GLfloat distance);

    /**
     * Pan the camera left.
     * The camera rotates using the up vector as an axis.
     * The rotation is to the left, the scenery moves to the right.
     * \param angle gives the amount of rotation in radians.
     */
    void panLeft(double angle);

    /**
     * Pan the camera right.
     * The camera rotates using the up vector as an axis.
     * The rotation is to the left, the scenery moves to the right.
     * \param angle gives the amount of rotation in radians.
     */
    void panRight(double angle);

    /**
     * Tilt the camera up.
     * The axis of rotation makes a right angle with the line joining
     * the camera to the model and with the up vector.
     * \param angle gives the amount of rotation in radians.
     */
    void tiltUp(double angle);

    /**
     * Tilt the camera down.
     * The axis of rotation makes a right angle with the line joining
     * the camera to the model and with the up vector.
     * \param angle gives the amount of rotation in radians.
     */
    void tiltDown(double angle);

    /**
     * Write a camera description.
     * This is mainly intended for debugging.  The information written
     * to the stream corresponds to a call to \c gluLookAt, with
     * the eye and model points, and the up vector.
     * \param os is the stream to write to.
     * \param cam is the camera.
     */
    friend std::ostream & operator<<(std::ostream & os, const Camera & cam);

private:

    /**
     * The kind of camera motion currently in effect.
     */
    enum { TRANSLATE, ROTATE } mode;

    /**
     * Return the current position of the camera.
     */
    Point eyeGet() const;

    /**
     * Return the direction in which the camera is currently pointing.
     */
    Vector dirGet() const;

    /**
     * The initial position of the camera.
     */
    Point eyeFirst;

    /**
     * The final direction of the camera.
     */
    Point eyeLast;

    /**
     * The initial direction of the camera.
     */
    Vector dirFirst;

    /**
     * The final direction of the camera.
     */
    Vector dirLast;

    /**
     * The 'up' vector for the scene.
     * The default value is \c J (the positive Y axis).
     */
    Vector up;

    /**
     * The distance from the camera to the model.
     */
    GLfloat distance;

    /**
     * The quaternion corresponding to the initial direction of the camera.
     */
    Quaternion qFirst;

    /**
     * The quaternion corresponding to the final direction of the camera.
     */
    Quaternion qLast;

    /**
     * sin(omega), where omega is the change of direction, in radians.
     */
    double sinOmega;

    /**
     * The change of direction, in radians.
     */
    double omega;
};


/**
  * 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 \c BMP image as a texture, the width and height
  * must be powers of two.  Thus a 256 x 512 \c 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.
    * 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 a path to a file with extension ".bmp".
    */
   PixelMap(const std::string & bmpFileName);

   /**
    * Construct a pixel map from a region of the OpenGL framebuffer.
    * \param x is the \c X coordinate of the lower left corner of the region.
    * \param y is the \c 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 a path to a file with extension ".bmp".
    */
   void read(const std::string & 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 \c X coordinate of the lower left corner of the region.
    * \param y is the \c 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 the constructor PixelMap(x, y, w, h).
    */
   void read(GLint x, GLint y, GLsizei w, GLsizei h);

   /**
    * Write a pixel map to an output stream as a \c BMP file.
    * \param bmpFileName is the name of the file to be written, with the extension '.bmp'.
    */
   void write(const std::string & bmpFileName) const;

   /**
    * 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() const;

   /**
    * 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() const;
      
   /**
    * 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);

   /**
    * 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.
    */
   std::string getName() const { return fileName; }

   /**
    * Return \c true if a pixel map has been read successfully.
    */
   bool ready() const { return pixels != NULL; }

   /**
    * Write a description of the pixel map to the output stream.
    * \param os is the stream to write to.
    * \param pm is the PixelMap to be written.
    */
   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);

   unsigned long numRows;
   unsigned long numCols;
   unsigned long size;
   std::string fileName;
   unsigned char *pixels;
};

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 \c (h,r) in which \c h is measured along
 * the axis of revolution and \c r is the distance from the
 * axis.  (If the axis is assumed to be vertical, then \c h
 * suggests 'height' and \c r suggests 'radius'.)
 * The surface is generated by rotating the profile
 * about the \c 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::Revolute(const Revolute&);

   /**
    * The assignment operator is private and not implemented: 
    * instances cannot be assigned.
    */
   const Revolute& Revolute::operator=(const Revolute&);

   int numSteps;
   int numSlices;
   double eccentricity;
   bool ready;
   GLfloatArray *coor;
   GLfloat *texCoor;
   Point *points;
   Vector *normals;
   Vector *faceNormals;
};


/**
 * \addtogroup 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;
}

/**
 * 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.
 * \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);

/**
 * Constructs a surface of revolution.
 * Constructs and renders a surface of revolution obtained by rotating
 * a curve about the \c 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]
   );
}

/** \} */

/**
 * \addtogroup materials
 * \{
 */

/**
 * \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: \c X = red, \c Y = green, \c 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 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);
}

/**
 * Use a Vector to move a Point.
 * \return the point \c p+v.
 */
inline Point operator+(const Point & p, const Vector & v)
{
    Point q(p);
    q += v;
    return q;
    // return Point(p.x + p.w * v.x, p.y + p.w * v.y, p.z + p.w * v.z, p.w);
}

/**
 * Use a Vector to move a Point.
 * \return the point \c p-v.
 */
inline Point operator-(const Point & p, const Vector & v)
{
    Point q(p);
    q -= v;
    return q;
    // return Point(p.x - p.w * v.x, p.y - p.w * v.y, p.z - p.w * v.z, p.w);
}

/**
 * Use a Vector to move a Point.
 * The corresponding function \c operator- is not provided because
 * there is no obvious interpretation of \c v-p.
 * However, there is an \c operator- for \c p-v.
 * \return Point at \c p+v.
 */
inline Point operator+(const Vector & v, const Point & p)
{
    Point q(p);
    q += v;
    return q;
    // 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 bool eq(const Point & p, const Point & q)
{ 
   return testEqual(p.x, q.x) && testEqual(p.y, q.y) && testEqual(p.z, q.z) && testEqual(p.w, q.w);
}

inline bool neq(const Point & p, const Point & q)
{ 
   return ! eq(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(Point & p, Point & q)
   : s(p), f(q)
{ }

inline Line::Line(Point & p, 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);
}

inline bool eq(const Line & m, const Line & n)
{ 
   return eq(m.s, n.s) && eq(n.s, n.f);
}

inline bool neq(const Line & m, const Line & n)
{ 
   return ! eq(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);
}

inline bool eq(const Plane & p, const Plane & q)
{ 
   return testEqual(p.a, q.a) && testEqual(p.b, q.b) && testEqual(p.c, q.c) && testEqual(p.d, q.d);
}

inline bool neq(const Plane & p, const Plane & q)
{ 
   return ! eq(p, q);
}



// class Vector: inlined constructors

inline Vector::Vector (GLfloat coordinates[])
{
   x = coordinates[0];
   y = coordinates[1];
   z = coordinates[2];
}

inline Vector::Vector(Point & p, 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 void Vector::rotate(GLfloat angle) const
{
    glRotatef(angle, x, y, z);
}

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

inline bool eq(const Vector & u, const Vector & v)
{ 
   return testEqual(u.x, v.x) && testEqual(u.y, v.y) && testEqual(u.z, v.z);
}

inline bool neq(const Vector & u, const Vector & v)
{ 
   return ! eq(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);
}

// 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 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 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 bool eq(const Matrix & m, const Matrix & n)
{ 
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         if (! testEqual(m(i,j), n(i,j)))
            return false;
   return true;
}

inline bool neq(const Matrix & m, const Matrix & n)
{ 
   for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
         if (! testEqual(m(i,j), n(i,j)))
            return true;
   return false;
}



// Class Quaternion: inlined member functions.

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 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 & 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;
}

/**
 * Return the quaternion \c 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 \c q+r. 
 */
inline Quaternion operator+(const Quaternion & q, const Quaternion & r)
{
    Quaternion res(q);
    res += r;
    return res;
}

/**
 * Return the quaternion \c q-r.
 * \note The difference of two unit quaternions is not in general a unit quaternion.  
 * \return the Quaternion \c q-r. 
 */
inline Quaternion operator-(const Quaternion & q, const Quaternion & r)
{
    Quaternion res(q);
    res -= r;
    return res;
}

/**
 * Return the quaternion product \c q*r.
 * If \c q and \c r represent 
 * rotations, then \c q*r represents the rotation \c q followed by the rotation \c r.
 * \note Quaternion multiplication is not commutative: \c q*r is not equal to \c r*q.
 * \return the Quaternion product \c q*r.  
 */
inline Quaternion operator*(const Quaternion & q, const Quaternion & r)
{
    Quaternion res(q);
    res *= r;
    return res;
}

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 operator/(const Quaternion & q, const Quaternion & r)
{
   GLfloat den = r.norm();
   return Quaternion
   (
      (q.s * r.s + q.v[0] * r.v[0] + q.v[1] * r.v[1] + q.v[2] * r.v[2]) / den, 
      (r.s * q.v[0] - q.s * r.v[0] - q.v[1] * r.v[2] + q.v[2] * r.v[1]) / den,
      (r.s * q.v[1] - q.s * r.v[1] - q.v[2] * r.v[0] + q.v[0] * r.v[2]) / den, 
      (r.s * q.v[2] - q.s * r.v[2] - q.v[0] * r.v[1] + q.v[1] * r.v[0]) / den 
   );
}

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

inline bool eq(const Quaternion & q, const Quaternion & r)
{ 
   return testEqual(q.s, r.s) && eq(q.v, r.v);
}

inline bool neq(const Quaternion & q, const Quaternion & r)
{ 
   return ! eq(q, r);
}

}; // end of namespace cugl

#endif

