# HG changeset patch
# User Eris Caffee <discordia@eldalin.com>
# Date 1318012428 18000
# Node ID 398894c52b697ee927fe6be88f33f4b5150da1e4
# Parent  139d189290e0ddd2ab124a324cf590b67672f17a
Changed mind.  Sticking with current design.  Matrix22<T> is done.

diff -r 139d189290e0 -r 398894c52b69 include/Math.h
--- a/include/Math.h	Fri Oct 07 11:09:15 2011 -0500
+++ b/include/Math.h	Fri Oct 07 13:33:48 2011 -0500
@@ -22,41 +22,42 @@
    ////////////////////////////////////////////////////////////////////////////////
    // Angles
 
-   namespace Angles
+   namespace Math
       {
       // This many decimals is enough to cover even IEEE quad precision (128 bit)
       // reals.  Run on PA-RISC lately?  :)
       const double PI = 3.14159265358979323846264338327950288;
-      }
 
-   ////////////////////////////////////////////////////////////////////////////////
-   // Mixed Vector / Matrix operations.
-   // * *=
-   //     Defined for Vector * Matrix and Matrix * Vector, but these
-   //     operators are not part of the class.  They are defined independently.
+      //////////////////////////////////////////////////////////////////////////
+      // Mixed Vector / Matrix operations.
+      // * *=
+      //     Defined for Vector * Matrix and Matrix * Vector, but these
+      //     operators are not part of the class.  They are defined independently.
+      //
 
-   // Vector * Matrix
-   template <typename T> inline Vector2<T>& operator*=(Vector2<T>& v, const Matrix22<T>& m)
-      {
-      Vector2<T> vres(0);
-      vres[0] = v[0]*m[0] + v[1]*m[1];
-      vres[1] = v[0]*m[2] + v[1]*m[3];
-      v = vres;
-      return v;
-      }
-   template <typename T> inline Vector2<T> operator*(const Vector2<T>& v, const Matrix22<T>& m)
-      { return Vector2<T>(v) *= m; }
+      // Vector * Matrix
+      template <typename T> inline Vector2<T>& operator*=(Vector2<T>& v, const Matrix22<T>& m)
+	 {
+	 Vector2<T> vres(0);
+	 vres[0] = v[0]*m[0] + v[1]*m[1];
+	 vres[1] = v[0]*m[2] + v[1]*m[3];
+	 v = vres;
+	 return v;
+	 }
+      template <typename T> inline Vector2<T> operator*(const Vector2<T>& v, const Matrix22<T>& m)
+	 { return Vector2<T>(v) *= m; }
 
-   // Matrix * Vector
-   // Note: can't define Matrix *= Vector since the result is a Vector.
-   template <typename T> inline Vector2<T> operator*(const Matrix22<T>& m, const Vector2<T>& v)
-      { 
-      Vector2<T> vres(0);
-      vres[0] = m[0]*v[0] + m[2]*v[1];
-      vres[1] = m[1]*v[0] + m[3]*v[1];
-      return vres;
-      }
+      // Matrix * Vector
+      // Note: can't define Matrix *= Vector since the result is a Vector.
+      template <typename T> inline Vector2<T> operator*(const Matrix22<T>& m, const Vector2<T>& v)
+	 { 
+	 Vector2<T> vres(0);
+	 vres[0] = m[0]*v[0] + m[2]*v[1];
+	 vres[1] = m[1]*v[0] + m[3]*v[1];
+	 return vres;
+	 }
 
+      } // namespace Math
 
    } // namespace arda
 
diff -r 139d189290e0 -r 398894c52b69 include/Matrix.h
--- a/include/Matrix.h	Fri Oct 07 11:09:15 2011 -0500
+++ b/include/Matrix.h	Fri Oct 07 13:33:48 2011 -0500
@@ -10,182 +10,283 @@
 
 namespace arda 
    {
-   //////////////////////////////////////////////////////////////////////////
-   // Maxtrices
-   //
-   // Supported operations on matrices
-   //
-   // []
-   //     If you have a Matrix named m you can access it's member elements as
-   //     m[0], m[1], m[2], etc.  The indexing is COLUMN MAJOR.  So for a 2x2
-   //     matrix, m[0] and m[1] are the elements of the first column, not the
-   //     first row.
-   // == != + += - -=
-   //     Defined for operations on two matrices of the same type.
-   // * *= / /=
-   //     Defined for scalar multiplication/division with int, float, and double.
-   //     * is defined as a standalone template operator for the scalar * Matrix form.
-   // * *=
-   //     Also defined for Matrix * Matrix
-   //
-   // assign()          Assign the value of one matrix to another.
-   // det()             Calculates the determinant of a matrix.
-   // getstring()       Returns a printable string representation of the matrix.
-   // transpose()       Returns the transpose of a matrix.
-   // setidentity()     Sets a matrix to the identity matrix.
+   namespace Math
+      {
+      //////////////////////////////////////////////////////////////////////////
+      // Maxtrices
+      //
+      // Supported operations on matrices
+      //
+      // []
+      //     If you have a Matrix named m you can access it's member elements as
+      //     m[0], m[1], m[2], etc.  The indexing is COLUMN MAJOR.  So for a 2x2
+      //     matrix, m[0] and m[1] are the elements of the first column, not the
+      //     first row.
+      // == != + += - -=
+      //     Defined for operations on two matrices of the same type.
+      // * *= / /=
+      //     Defined for scalar multiplication/division with int, float, and double.
+      //     * is defined as a standalone template operator for the scalar * Matrix
+      //     form.
+      // * *=
+      //     Also defined for Matrix * Matrix
+      //
+      // assign()          Assign the value of one matrix to another.
+      // getcol()          Returns the vector that is column n of the matrix.
+      // getstring()       Returns a printable string representation of the matrix.
+      // setcol()          Set the vector that is column n of the matrix.
+      //                   Overloaded so that you can set from a vector or from
+      //                   specific values.
+      // setidentity()     Sets a matrix to the identity matrix.
+      //
+      //
+      //
+      // The following are functions, not class methods.  This is because I think
+      // that the notation det(m) and transpose(m), for example, is more natural
+      // to write and read than m.det() and m.transpose().  Using functional
+      // notation just more naturally mimics the standard mathematical notation.
+      //
+      // det()             Calculates the determinant of a matrix.
+      // transpose()       Returns the transpose of a matrix.
+      
+      //////////////////////////////////////////////////////////////////////////
+      template <typename T> class Matrix22
+	 {
+	 public:
+	 T m[4];
 
-   // inverse()         Returns the inverse of a non-singular matrix.
-   // getcol()          Returns the vector that is column n of the matrix.
-   // setcol()          Set the vector that is column n of the matrix.
+	 // Constructors
+	 Matrix22() {}
+	 explicit Matrix22(T a) { m[0] = m[1] = m[2] = m[3] = a; }
+	 Matrix22(T a0, T a1, T a2, T a3) 
+	    { 
+	    m[0] = a0; m[2] = a2;
+	    m[1] = a1; m[3] = a3; 
+	    } 
+
+	 // Array indexing
+	 // Remember: column major order is used.
+	 inline T& operator[](unsigned int i) 
+	    { assert (i<4); return m[i]; }
+	 inline T operator[](unsigned int i) const
+	    { assert (i<4); return m[i]; }
+
+	 // Assignment
+	 inline Matrix22<T>& operator=(const Matrix22<T>& m2)
+	    { 
+	    m[0] = m2[0]; m[2] = m2[2]; 
+	    m[1] = m2[1]; m[3] = m2[3]; 
+	    return *this;
+	    }
+	 inline Matrix22<T>& assign(T a0, T a1, T a2, T a3)
+	    { 
+	    m[0] = a0; m[2] = a2;
+	    m[1] = a1; m[3] = a3; 
+	    return *this;
+	    }
+
+	 // Comparison
+	 inline bool operator==(Matrix22<T>& m2) const
+	    {
+	    return  
+		  m[0] == m2[0] && m[2] == m2[2] &&
+		  m[1] == m2[1] && m[3] == m2[3]; 
+	    }
+	 inline bool operator!=(Matrix22<T>& m2) const
+	    { return ! (*this == m2); } 
+
+	 // Matrix addition
+	 inline Matrix22<T>& operator+=(const Matrix22<T>& m2)
+	    {
+	    m[0] += m2[0]; m[2] += m2[2]; 
+	    m[1] += m2[1]; m[3] += m2[3];
+	    return *this;
+	    }
+	 inline Matrix22<T> operator+(const Matrix22<T>& m2) const
+	    { return Matrix22<T>(*this) += m2; }
+
+	 // Matrix subtraction
+	 inline Matrix22<T>& operator-=(const Matrix22<T>& m2)
+	    {
+	    m[0] -= m2[0]; m[2] -= m2[2]; 
+	    m[1] -= m2[1]; m[3] -= m2[3];
+	    return *this;
+	    }
+	 inline Matrix22<T> operator-(const Matrix22<T>& m2) const
+	    { return Matrix22<T>(*this) -= m2; }
       
-   //////////////////////////////////////////////////////////////////////////
-   template <typename T> class Matrix22
-      {
-      public:
-      T m[4];
+	 // Scalar multiplication
+	 inline Matrix22<T>& operator*=(const int a)
+	    {
+	    m[0] *= a; m[2] *= a;
+	    m[1] *= a; m[3] *= a;
+	    return *this;
+	    }
+	 inline Matrix22<T>& operator*=(const float a)
+	    {
+	    m[0] *= a; m[2] *= a;
+	    m[1] *= a; m[3] *= a;
+	    return *this;
+	    }
+	 inline Matrix22<T>& operator*=(const double a)
+	    {
+	    m[0] *= a; m[2] *= a;
+	    m[1] *= a; m[3] *= a;
+	    return *this;
+	    }
 
-      // Constructors
-      Matrix22() {}
-      Matrix22(T a) { m[0] = m[1] = m[2] = m[3] = a; }
-      Matrix22(T a0, T a1, T a2, T a3) 
-	 { 
-	 m[0] = a0; m[2] = a2;
-	 m[1] = a1; m[3] = a3; 
-	 } 
+	 // Scalar division
+	 inline Matrix22<T>& operator/=(const int a)
+	    {
+	    assert(a!=0);
+	    m[0] /= a; m[2] /= a;
+	    m[1] /= a; m[3] /= a;
+	    return *this;
+	    }
+	 inline Matrix22<T>& operator/=(const float a)
+	    {
+	    assert(a!=0);
+	    m[0] /= a; m[2] /= a;
+	    m[1] /= a; m[3] /= a;
+	    return *this;
+	    }
+	 inline Matrix22<T>& operator/=(const double a)
+	    {
+	    assert(a!=0);
+	    m[0] /= a; m[2] /= a;
+	    m[1] /= a; m[3] /= a;
+	    return *this;
+	    }
+	 inline Matrix22<T> operator/(const int a)
+	    { return Matrix22<T>(*this) /= a; }
+	 inline Matrix22<T> operator/(const float a)
+	    { return Matrix22<T>(*this) /= a; }
+	 inline Matrix22<T> operator/(const double a)
+	    { return Matrix22<T>(*this) /= a; }
 
-      // Array indexing
-      // Remember: column major order is used.
-      inline T& operator[](unsigned int i) 
-	 { assert (i<4); return m[i]; }
-      inline const T& operator[](unsigned int i) const
-	 { assert (i<4); return m[i]; }
+	 // Matrix multiplication
+	 // Not sure if this should be inlined at all.  Sure the compiler will
+	 // probably ignore the inline request here, but maybe it won't and maybe
+	 // that would be bad.  Needs real testing.
+	 inline Matrix22<T>& operator*=(Matrix22<T> m2)
+	    {
+	    const int size=2;
+	    Matrix22<T> mres(0);
+	    int i, j, k;
+	    for (i=0; i<size; ++i)
+	       for (j=0; j<size; ++j)
+		  for (k=0; k<size; ++k)
+		     mres[size*i+j] += (*this)[size*k+j] * m2[size*i+k];
+	    *this = mres;
+	    return *this;
+	    }
+	 inline Matrix22<T> operator*(Matrix22<T> m2) const
+	    { return Matrix22<T>(*this) *= m2; }
 
-      // Assignment
-      inline Matrix22<T>& operator=(const Matrix22<T>& m2)
-	 { 
-	 m[0] = m2[0]; m[2] = m2[2]; 
-	 m[1] = m2[1]; m[3] = m2[3]; 
-	 return *this;
-	 }
-      inline Matrix22<T>& assign(T a0, T a1, T a2, T a3)
-	 { 
-	 m[0] = a0; m[2] = a2;
-	 m[1] = a1; m[3] = a3; 
-	 return *this;
-	 }
 
-      // Comparison
-      inline bool operator==(Matrix22<T>& m2) const
-	 {
-	 return  
-	       m[0] == m2[0] && m[2] == m2[2] &&
-	       m[1] == m2[1] && m[3] == m2[3]; 
-	 }
-      inline bool operator!=(Matrix22<T>& m2) const
-	 { return ! (*this == m2); } 
+	 inline Vector2<T> getcol(const unsigned int i) const
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    return Vector2<T>(m[i*size], m[i*size+1]);
+	    }
 
-      // Matrix addition
-      inline Matrix22<T>& operator+=(const Matrix22<T>& m2)
-	 {
-	 m[0] += m2[0]; m[2] += m2[2]; 
-	 m[1] += m2[1]; m[3] += m2[3];
-	 return *this;
-	 }
-      inline Matrix22<T> operator+(const Matrix22<T>& m2) const
-	 { return Matrix22<T>(*this) += m2; }
+	 inline Vector2<T> getrow(const unsigned int i) const
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    return Vector2<T>(m[i], m[i+size]);
+	    }
+	 
+	 inline Matrix22<T>& setcol(const unsigned int i, const Vector2<T>& v)
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    m[i*size] = v[0];
+	    m[i*size+1] = v[1];
+	    return *this;
+	    }
 
-      // Matrix subtraction
-      inline Matrix22<T>& operator-=(const Matrix22<T>& m2)
-	 {
-	 m[0] -= m2[0]; m[2] -= m2[2]; 
-	 m[1] -= m2[1]; m[3] -= m2[3];
-	 return *this;
-	 }
-      inline Matrix22<T> operator-(const Matrix22<T>& m2) const
-	 { return Matrix22<T>(*this) -= m2; }
+	 inline Matrix22<T>& setcol(const unsigned int i, const T a, const T b)
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    m[i*size] = a;
+	    m[i*size+1] = b;
+	    return *this;
+	    }
+
+	 inline Matrix22<T>& setrow(const unsigned int i, const Vector2<T>& v)
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    m[i] = v[0];
+	    m[i+size] = v[1];
+	    return *this;
+	    }
+
+	 inline Matrix22<T>& setrow(const unsigned int i, const T a, const T b)
+	    {
+	    assert(i<2);
+	    const int size=2;
+	    m[i] = a;
+	    m[i+size] = b;
+	    return *this;
+	    }
+
+	 std::string & getstring(std::string & s) const;
+
+	 inline Matrix22<T>& setidentity()
+	    {
+	    memset(m, 0, sizeof(m));
+	    m[0] = m[3] = 1;
+	    return *this;
+	    }
+
+	 };
+
+      // Scalar multiplication continued
+      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const int a)
+	 { return Matrix22<T>(m) *= a; }
+      template <typename T> inline Matrix22<T> operator*(const int a, const Matrix22<T>& m)
+	 { return Matrix22<T>(m) *= a; }
+      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const float a)
+	 { return Matrix22<T>(m) *= a; }
+      template <typename T> inline Matrix22<T> operator*(const float a, const Matrix22<T>& m)
+	 { return Matrix22<T>(m) *= a; }
+      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const double a)
+	 { return Matrix22<T>(m) *= a; }
+      template <typename T> inline Matrix22<T> operator*(const double a, const Matrix22<T>& m)
+	 { return Matrix22<T>(m) *= a; }
+
+
+
+
+      //////////////////////////////////////////////////////////////////////////
+      typedef Matrix22<int> Matrix22i;
+      typedef Matrix22<float> Matrix22f;
+      typedef Matrix22<double> Matrix22d;
+
+#ifdef NOTHING   
+      typedef Matrix33<int> Matrix33i;
+      typedef Matrix33<float> Matrix33f;
+      typedef Matrix33<double> Matrix33d;
       
-      // Scalar multiplication
-      inline Matrix22<T>& operator*=(const int a)
-	 {
-	 m[0] *= a; m[2] *= a;
-	 m[1] *= a; m[3] *= a;
-	 return *this;
-	 }
-      inline Matrix22<T>& operator*=(const float a)
-	 {
-	 m[0] *= a; m[2] *= a;
-	 m[1] *= a; m[3] *= a;
-	 return *this;
-	 }
-      inline Matrix22<T>& operator*=(const double a)
-	 {
-	 m[0] *= a; m[2] *= a;
-	 m[1] *= a; m[3] *= a;
-	 return *this;
-	 }
+      typedef Matrix44<int> Matrix44i;
+      typedef Matrix44<float> Matrix44f;
+      typedef Matrix44<double> Matrix44d;
+      
+#endif
 
-      // Scalar division
-      inline Matrix22<T>& operator/=(const int a)
-	 {
-	 assert(a!=0);
-	 m[0] /= a; m[2] /= a;
-	 m[1] /= a; m[3] /= a;
-	 return *this;
-	 }
-      inline Matrix22<T>& operator/=(const float a)
-	 {
-	 assert(a!=0);
-	 m[0] /= a; m[2] /= a;
-	 m[1] /= a; m[3] /= a;
-	 return *this;
-	 }
-      inline Matrix22<T>& operator/=(const double a)
-	 {
-	 assert(a!=0);
-	 m[0] /= a; m[2] /= a;
-	 m[1] /= a; m[3] /= a;
-	 return *this;
-	 }
-      inline Matrix22<T> operator/(const int a)
-	 { return Matrix22<T>(*this) /= a; }
-      inline Matrix22<T> operator/(const float a)
-	 { return Matrix22<T>(*this) /= a; }
-      inline Matrix22<T> operator/(const double a)
-	 { return Matrix22<T>(*this) /= a; }
 
-      // Matrix multiplication
-      // Not sure if this should be inlined at all.  Sure the compiler will
-      // probably ignoe the inline request here, but maybe it won't and maybe
-      // that would be bad.  Needs real testing.
-      inline Matrix22<T>& operator*=(Matrix22<T> m2)
-	 {
-	 const int size=2;
-	 Matrix22<T> mres(0);
-	 int i, j, k;
-	 for (i=0; i<size; ++i)
-	    for (j=0; j<size; ++j)
-	       for (k=0; k<size; ++k)
-		  mres[size*i+j] += (*this)[size*k+j] * m2[size*i+k];
-	 *this = mres;
-	 return *this;
-	 }
-      inline Matrix22<T> operator*(Matrix22<T> m2) const
-	 { return Matrix22<T>(*this) *= m2; }
 
-      // methods
-      inline double det() const
+      //////////////////////////////////////////////////////////////////////////
+      // Matrix functions
+
+      template <typename T> inline double det(const Matrix22<T> m)
 	 { return (double) (m[0] * m[3]) - (double) (m[1] * m[2]); }
 
-      std::string & getstring(std::string & s) const;
-
-      inline Matrix22<T>& setidentity()
-	 {
-	 memset(m, 0, sizeof(m));
-	 m[0] = m[3] = 1;
-	 return *this;
-	 }
-
-      inline Matrix22<T> transpose() const
+      template <typename T> inline Matrix22<T> transpose(const Matrix22<T>& m)
 	 {
 	 const int size=2;
 	 Matrix22<T> mres(0);
@@ -195,44 +296,15 @@
 	       mres[size*i+j] = m[size*j+i];
 	 return mres;
 	 }
-      };
 
-   // Scalar multiplication continued
-   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const int a)
-      { return Matrix22<T>(m) *= a; }
-   template <typename T> inline Matrix22<T> operator*(const int a, const Matrix22<T>& m)
-      { return Matrix22<T>(m) *= a; }
-   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const float a)
-      { return Matrix22<T>(m) *= a; }
-   template <typename T> inline Matrix22<T> operator*(const float a, const Matrix22<T>& m)
-      { return Matrix22<T>(m) *= a; }
-   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const double a)
-      { return Matrix22<T>(m) *= a; }
-   template <typename T> inline Matrix22<T> operator*(const double a, const Matrix22<T>& m)
-      { return Matrix22<T>(m) *= a; }
-
-
-
-
-   //////////////////////////////////////////////////////////////////////////
-   typedef Matrix22<int> Matrix22i;
-   typedef Matrix22<float> Matrix22f;
-   typedef Matrix22<double> Matrix22d;
-
-#ifdef NOTHING   
-   typedef Matrix33<int> Matrix33i;
-   typedef Matrix33<float> Matrix33f;
-   typedef Matrix33<double> Matrix33d;
-      
-   typedef Matrix44<int> Matrix44i;
-   typedef Matrix44<float> Matrix44f;
-   typedef Matrix44<double> Matrix44d;
-      
-#endif
+      } // namespace Math
    } // namespace arda
 
+
 ////////////////////////////////////////////////////////////////////////////////
-template <typename T> std::string & arda::Matrix22<T>::getstring(std::string & s) const
+// Matrix methods
+
+template <typename T> std::string & arda::Math::Matrix22<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
diff -r 139d189290e0 -r 398894c52b69 include/Vector.h
--- a/include/Vector.h	Fri Oct 07 11:09:15 2011 -0500
+++ b/include/Vector.h	Fri Oct 07 13:33:48 2011 -0500
@@ -10,449 +10,453 @@
 
 namespace arda 
    {
-   ////////////////////////////////////////////////////////////////////////////////
-   // Vectors
-   //
-   // Templated classes for 2, 3, and 4 dimensional vectors of any type,
-   //     although these types are really intended to be int, float, and double 
-   //     only.  (i.e. if you roll your own base type, it had better act like
-   //     a numeric scalar.)
-   //
-   // Supported operations on vectors
-   //
-   // Construction
-   //     Vector v;             // Default zero vector.
-   //     Vector v2(v);         // Construct from other vector.
-   //     T a, b;
-   //     Vector v(a, b)        // Construct from values.
-   // []
-   //     If you have a Vector named v you can access it's member elements as
-   //     v[0], v[1], v[2], and v[3]. in addtion to v.x, v.y, v.z, and v.w
-   // == != + += - -=
-   //     Defined for operations on two vectors of the same type.
-   // * *= / /=
-   //     Defined for scalar multiplication/division with int, float, and double.
-   //     * is defined as a standalone template operator for the scalar * Vector form.
-   // assign(T a, T b [, T b [, T d]])
-   //     Directly assign values to the vector.
-   // cross(Vector3& v2)
-   // cross(Vector3& v2, Vector3& vres)
-   //     Cross product of the vector with v2.  Only defined for Vector3 types.
-   //     The version that takes two arguments stores the result in the second
-   //     argument.  This may be faster since no temporary object is created
-   //     during the operation.
-   // dot(Vector& v2)
-   //     The dot product of the vector with v2.
-   // get_angle(Vector& v2)
-   //     Returns the angle (in radians) between the vector and v2.
-   //     Not meaningful for int vectors.
-   // get_anglen(Vector& v2)
-   //     Like get_angle(), but the vectors must already be normalized.
-   //     Not meaningful for int vectors.
-   // getstring(std::string & s)
-   //     Returns a C string representation of the vector.
-   //     This is one of the few non-inline Vector functions.
-   // length()
-   //     Get the length of a vector.
-   //     Returns a double for all vector types.
-   // normalize()
-   //     Normalizes the vecor.
-   //     Not meaningful for int vectors.
-   // proj(Vector& v2, Vector& vres)
-   //     Calculates the projection of the vector onto v2 and returns the result 
-   //     in vres. 
-   //     Not meaningful for int vectors.
-   //     Returns vres.
-   //
-   //
-   // Note that the following typedefs are defined for convenience:
-   //
-   // typedef Vector2<int> Vector2i;
-   // typedef Vector2<float> Vector2f;
-   // typedef Vector2<double> Vector2d;
-   // typedef Vector3<int> Vector3i;
-   // typedef Vector3<float> Vector3f;
-   // typedef Vector3<double> Vector3d;
-   // typedef Vector4<int> Vector4i;
-   // typedef Vector4<float> Vector4f;
-   // typedef Vector4<double> Vector4d;
-   //
-   //
-   //
-   // Scalar multiplication and division are overloaded.  Why?
-   // I'm not sure I need this.  The idea is to support all likely scalar
-   // types without having to rely on implicit or explicit type conversion,
-   // especially not a conversion that might result in loss of precision.
-   // Multiplying a float Vector by a double scalar, for example.
-   //
-   // Is there a better way to do this?  Am I worrying about nothing?
-   //
-   //
-   //
-   // I suppose I could make this into a single Vector template on two parameters
-   // template <typename T, unsigned int N> class Vector
-   //    {
-   //    T v[N];
-   //
-   //    etc...
-   //
-   // This might buy perfect generality, but it would come at the expense of
-   // loops in all of the code, so I don't think it would really be worth it.
+   namespace Math
+      {
+      ////////////////////////////////////////////////////////////////////////////////
+      // Vectors
+      //
+      // Templated classes for 2, 3, and 4 dimensional vectors of any type,
+      //     although these types are really intended to be int, float, and double 
+      //     only.  (i.e. if you roll your own base type, it had better act like
+      //     a numeric scalar.)
+      //
+      // Supported operations on vectors
+      //
+      // Construction
+      //     Vector v;             // Default zero vector.
+      //     Vector v2(v);         // Construct from other vector.
+      //     T a, b;
+      //     Vector v(a, b)        // Construct from values.
+      // []
+      //     If you have a Vector named v you can access it's member elements as
+      //     v[0], v[1], v[2], and v[3]. in addtion to v.x, v.y, v.z, and v.w
+      // == != + += - -=
+      //     Defined for operations on two vectors of the same type.
+      // * *= / /=
+      //     Defined for scalar multiplication/division with int, float, and double.
+      //     * is defined as a standalone template operator for the scalar * Vector form.
+      // assign(T a, T b [, T b [, T d]])
+      //     Directly assign values to the vector.
+      // cross(Vector3& v2)
+      // cross(Vector3& v2, Vector3& vres)
+      //     Cross product of the vector with v2.  Only defined for Vector3 types.
+      //     The version that takes two arguments stores the result in the second
+      //     argument.  This may be faster since no temporary object is created
+      //     during the operation.
+      // dot(Vector& v2)
+      //     The dot product of the vector with v2.
+      // get_angle(Vector& v2)
+      //     Returns the angle (in radians) between the vector and v2.
+      //     Not meaningful for int vectors.
+      // get_anglen(Vector& v2)
+      //     Like get_angle(), but the vectors must already be normalized.
+      //     Not meaningful for int vectors.
+      // getstring(std::string & s)
+      //     Returns a C string representation of the vector.
+      //     This is one of the few non-inline Vector functions.
+      // length()
+      //     Get the length of a vector.
+      //     Returns a double for all vector types.
+      // normalize()
+      //     Normalizes the vecor.
+      //     Not meaningful for int vectors.
+      // proj(Vector& v2, Vector& vres)
+      //     Calculates the projection of the vector onto v2 and returns the result 
+      //     in vres. 
+      //     Not meaningful for int vectors.
+      //     Returns vres.
+      //
+      //
+      // Note that the following typedefs are defined for convenience:
+      //
+      // typedef Vector2<int> Vector2i;
+      // typedef Vector2<float> Vector2f;
+      // typedef Vector2<double> Vector2d;
+      // typedef Vector3<int> Vector3i;
+      // typedef Vector3<float> Vector3f;
+      // typedef Vector3<double> Vector3d;
+      // typedef Vector4<int> Vector4i;
+      // typedef Vector4<float> Vector4f;
+      // typedef Vector4<double> Vector4d;
+      //
+      //
+      //
+      // Scalar multiplication and division are overloaded.  Why?
+      // I'm not sure I need this.  The idea is to support all likely scalar
+      // types without having to rely on implicit or explicit type conversion,
+      // especially not a conversion that might result in loss of precision.
+      // Multiplying a float Vector by a double scalar, for example.
+      //
+      // Is there a better way to do this?  Am I worrying about nothing?
+      //
+      //
+      //
+      // I suppose I could make this into a single Vector template on two parameters
+      // template <typename T, unsigned int N> class Vector
+      //    {
+      //    T v[N];
+      //
+      //    etc...
+      //
+      // This might buy perfect generality, but it would come at the expense of
+      // loops in all of the code, so I don't think it would really be worth it.
+      
+      
+      /////////////////////////////////////////////////////////////////////////////
+      template <typename T> class Vector2
+	 {
+	 public:
+	 T x, y;
+	 
+	 // Constructors
+	 Vector2() {}
+	 explicit Vector2(const T a) : x (a), y(a) {}
+	 Vector2(const T a, const T b) : x (a), y(b) {}
 
+	 // Array indexing
+	 inline T& operator[](unsigned int i)
+	    { assert (i<2); if (i==0) return x; return y; }
+	 inline T operator[](unsigned int i) const
+	    { assert (i<2); if (i==0) return x; return y; }
 
-   /////////////////////////////////////////////////////////////////////////////
-   template <typename T> class Vector2
-      {
-      public:
-      T x, y;
+	 // Assignment
+	 inline Vector2<T>& operator=(const Vector2<T>& v2)
+	    { x = v2.x; y = v2.y; return *this; }
+	 inline Vector2<T>& assign(const T a, const T b)
+	    { x = a; y = b; return *this; }
 
-      // Constructors
-      Vector2() {}
-      Vector2(T a) : x (a), y(a) {}
-      Vector2(T a, T b) : x (a), y(b) {}
+	 // Comparison
+	 inline bool operator==(Vector2<T>& v2) const
+	    { return ((x == v2.x) && (y == v2.y)); }
+	 inline bool operator!=(Vector2<T>& v2) const
+	    { return ! (*this == v2); }
 
-      // Array indexing
-      inline T& operator[](unsigned int i)
-	 { assert (i<2); if (i==0) return x; return y; }
-      inline const T& operator[](unsigned int i) const
-	 { assert (i<2); if (i==0) return x; return y; }
+	 // Vector addition
+	 inline Vector2<T>& operator+=(const Vector2<T>& v2)
+	    { x += v2.x; y += v2.y; return *this; }
+	 inline Vector2<T> operator+(const Vector2<T>& v2) const
+	    { return Vector2<T>(*this) += v2; }
 
-      // Assignment
-      inline Vector2<T>& operator=(const Vector2<T>& v2)
-	 { x = v2.x; y = v2.y; return *this; }
-      inline Vector2<T>& assign(const T a, const T b)
-	 { x = a; y = b; return *this; }
+	 // Vector subtraction
+	 inline Vector2<T>& operator-=(const Vector2<T>& v2)
+	    { x -= v2.x; y -= v2.y; return *this; }
+	 inline Vector2<T> operator-(const Vector2<T>& v2) const
+	    { return Vector2<T>(*this) -= v2; }
 
-      // Comparison
-      inline bool operator==(Vector2<T>& v2) const
-	 { return ((x == v2.x) && (y == v2.y)); }
-      inline bool operator!=(Vector2<T>& v2) const
-	 { return ! (*this == v2); }
+	 // Scalar multiplication
+	 inline Vector2<T>& operator*=(const int a)
+	    { x *= a; y *= a; return *this; }
+	 inline Vector2<T>& operator*=(const float a)
+	    { x *= a; y *= a; return *this; }
+	 inline Vector2<T>& operator*=(const double a)
+	    { x *= a; y *= a; return *this; }
+	 inline Vector2<T> operator*(const int a)
+	    { return Vector2<T>(*this) *= a;}
+	 inline Vector2<T> operator*(const float a)
+	    { return Vector2<T>(*this) *= a;}
+	 inline Vector2<T> operator*(const double a)
+	    { return Vector2<T>(*this) *= a;}
 
-      // Vector addition
-      inline Vector2<T>& operator+=(const Vector2<T>& v2)
-	 { x += v2.x; y += v2.y; return *this; }
-      inline Vector2<T> operator+(const Vector2<T>& v2) const
-	 { return Vector2<T>(*this) += v2; }
+	 // Scalar division
+	 inline Vector2<T>& operator/=(const int a)
+	    { assert(a!=0); x /= a; y /= a; return *this; }
+	 inline Vector2<T>& operator/=(const float a)
+	    { assert(a!=0); x /= a; y /= a; return *this; }
+	 inline Vector2<T>& operator/=(const double a)
+	    { assert(a!=0); x /= a; y /= a; return *this; }
+	 inline Vector2<T> operator/(const int a)
+	    { return Vector2<T>(*this) /= a;}
+	 inline Vector2<T> operator/(const float a)
+	    { return Vector2<T>(*this) /= a;}
+	 inline Vector2<T> operator/(const double a)
+	    { return Vector2<T>(*this) /= a;}
 
-      // Vector subtraction
-      inline Vector2<T>& operator-=(const Vector2<T>& v2)
-	 { x -= v2.x; y -= v2.y; return *this; }
-      inline Vector2<T> operator-(const Vector2<T>& v2) const
-	 { return Vector2<T>(*this) -= v2; }
 
-      // Scalar multiplication
-      inline Vector2<T>& operator*=(const int a)
-	 { x *= a; y *= a; return *this; }
-      inline Vector2<T>& operator*=(const float a)
-	 { x *= a; y *= a; return *this; }
-      inline Vector2<T>& operator*=(const double a)
-	 { x *= a; y *= a; return *this; }
-      inline Vector2<T> operator*(const int a)
-	 { return Vector2<T>(*this) *= a;}
-      inline Vector2<T> operator*(const float a)
-	 { return Vector2<T>(*this) *= a;}
-      inline Vector2<T> operator*(const double a)
-	 { return Vector2<T>(*this) *= a;}
+	 // methods
 
-      // Scalar division
-      inline Vector2<T>& operator/=(const int a)
-	 { assert(a!=0); x /= a; y /= a; return *this; }
-      inline Vector2<T>& operator/=(const float a)
-	 { assert(a!=0); x /= a; y /= a; return *this; }
-      inline Vector2<T>& operator/=(const double a)
-	 { assert(a!=0); x /= a; y /= a; return *this; }
-      inline Vector2<T> operator/(const int a)
-	 { return Vector2<T>(*this) /= a;}
-      inline Vector2<T> operator/(const float a)
-	 { return Vector2<T>(*this) /= a;}
-      inline Vector2<T> operator/(const double a)
-	 { return Vector2<T>(*this) /= a;}
+	 inline T dot(const Vector2<T>& v2) const
+	    { return x*v2.x + y*v2.y; }
 
+	 inline double get_angle(Vector2<T>& v2)
+	    { 
+	    double tmp = dot(v2); 
+	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
+	    }
 
-      // methods
+	 inline double get_anglen(Vector2<T>& v2)
+	    { return acos((double) dot(v2)); }
 
-      inline T dot(const Vector2<T>& v2) const
-	 { return x*v2.x + y*v2.y; }
+	 std::string & getstring(std::string & s) const;
 
-      inline double get_angle(Vector2<T>& v2)
-	 { 
-	 double tmp = dot(v2); 
-	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
-	 }
+	 inline double length(void) const
+	    { return sqrt((dot(*this))); }
 
-      inline double get_anglen(Vector2<T>& v2)
-	 { return acos((double) dot(v2)); }
+	 inline Vector2<T>& normalize()
+	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; return *this; }
 
-      std::string & getstring(std::string & s) const;
+	 inline Vector2<T> proj(Vector2<T>& v2)
+	    { Vector2<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 inline Vector2<T>& proj(Vector2<T>& v2, Vector2<T>& vres)
+	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 };
 
-      inline double length(void) const
-	 { return sqrt((dot(*this))); }
+      // Scalar multiplication, continued
+      template <typename T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
+	 { return Vector2<T>(v) *= a;}
+      template <typename T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
+	 { return Vector2<T>(v) *= a;}
+      template <typename T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
+	 { return Vector2<T>(v) *= a;}
 
-      inline Vector2<T>& normalize()
-	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; return *this; }
 
-      inline Vector2<T> proj(Vector2<T>& v2)
-	 { Vector2<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      inline Vector2<T>& proj(Vector2<T>& v2, Vector2<T>& vres)
-	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      };
+      /////////////////////////////////////////////////////////////////////////////
+      template <typename T> class Vector3
+	 {
+	 public:
+	 T x, y, z;
 
-   // Scalar multiplication, continued
-   template <typename T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
-      { return Vector2<T>(v) *= a;}
-   template <typename T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
-      { return Vector2<T>(v) *= a;}
-   template <typename T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
-      { return Vector2<T>(v) *= a;}
+	 // Constructors
+	 Vector3() {}
+	 explicit Vector3(const T a) : x (a), y(a), z(a) {}
+	 Vector3(const T a, const T b, const T c) : x (a), y(b), z(c) {}
 
+	 // Array indexing
+	 inline T& operator[](unsigned int i)
+	    { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
+	 inline T operator[](unsigned int i) const
+	    { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
 
-   /////////////////////////////////////////////////////////////////////////////
-   template <typename T> class Vector3
-      {
-      public:
-      T x, y, z;
+	 // Assignment
+	 inline Vector3<T>& operator=(const Vector3<T>& v2)
+	    { x = v2.x; y = v2.y; z = v2.z; return *this; }
+	 inline Vector3<T>& assign(const T a, const T b, const T c)
+	    { x = a; y = b; z = c; return *this; }
 
-      // Constructors
-      Vector3() {}
-      Vector3(T a) : x (a), y(a), z(a) {}
-      Vector3(T a, T b, T c) : x (a), y(b), z(c) {}
+	 // Comparison
+	 inline bool operator==(Vector3<T>& v2) const
+	    { return ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
+	 inline bool operator!=(Vector3<T>& v2) const
+	    { return ! (*this == v2); }
 
-      // Array indexing
-      inline T& operator[](unsigned int i)
-	 { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
-      inline const T& operator[](unsigned int i) const
-	 { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
+	 // Vector addition
+	 inline Vector3<T>& operator+=(const Vector3<T>& v2)
+	    { x += v2.x; y += v2.y; z += v2.z; return *this; }
+	 inline Vector3<T> operator+(const Vector3<T>& v2) const
+	    { return Vector3<T>(*this) += v2; }
 
-      // Assignment
-      inline Vector3<T>& operator=(const Vector3<T>& v2)
-	 { x = v2.x; y = v2.y; z = v2.z; return *this; }
-      inline Vector3<T>& assign(const T a, const T b, const T c)
-	 { x = a; y = b; z = c; return *this; }
+	 // Vector subtraction
+	 inline Vector3<T>& operator-=(const Vector3<T>& v2)
+	    { x -= v2.x; y -= v2.y; z -= v2.z; return *this; }
+	 inline Vector3<T> operator-(const Vector3<T>& v2) const
+	    { return Vector3<T>(*this) -= v2; }
 
-      // Comparison
-      inline bool operator==(Vector3<T>& v2) const
-	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
-      inline bool operator!=(Vector3<T>& v2) const
-	 { return ! (*this == v2); }
+	 // Scalar multiplication
+	 inline Vector3<T>& operator*=(const int a)
+	    { x *= a; y *= a; z *= a; return *this; }
+	 inline Vector3<T>& operator*=(const float a)
+	    { x *= a; y *= a; z *= a; return *this; }
+	 inline Vector3<T>& operator*=(const double a)
+	    { x *= a; y *= a; z *= a; return *this; }
+	 inline Vector3<T> operator*(const int a)
+	    { return Vector3<T>(*this) *= a;}
+	 inline Vector3<T> operator*(const float a)
+	    { return Vector3<T>(*this) *= a;}
+	 inline Vector3<T> operator*(const double a)
+	    { return Vector3<T>(*this) *= a;}
 
-      // Vector addition
-      inline Vector3<T>& operator+=(const Vector3<T>& v2)
-	 { x += v2.x; y += v2.y; z += v2.z; return *this; }
-      inline Vector3<T> operator+(const Vector3<T>& v2) const
-	 { return Vector3<T>(*this) += v2; }
+	 // Scalar division
+	 inline Vector3<T>& operator/=(const int a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
+	 inline Vector3<T>& operator/=(const float a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
+	 inline Vector3<T>& operator/=(const double a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
+	 inline Vector3<T> operator/(const int a)
+	    { return Vector3<T>(*this) /= a;}
+	 inline Vector3<T> operator/(const float a)
+	    { return Vector3<T>(*this) /= a;}
+	 inline Vector3<T> operator/(const double a)
+	    { return Vector3<T>(*this) /= a;}
 
-      // Vector subtraction
-      inline Vector3<T>& operator-=(const Vector3<T>& v2)
-	 { x -= v2.x; y -= v2.y; z -= v2.z; return *this; }
-      inline Vector3<T> operator-(const Vector3<T>& v2) const
-	 { return Vector3<T>(*this) -= v2; }
 
-      // Scalar multiplication
-      inline Vector3<T>& operator*=(const int a)
-	 { x *= a; y *= a; z *= a; return *this; }
-      inline Vector3<T>& operator*=(const float a)
-	 { x *= a; y *= a; z *= a; return *this; }
-      inline Vector3<T>& operator*=(const double a)
-	 { x *= a; y *= a; z *= a; return *this; }
-      inline Vector3<T> operator*(const int a)
-	 { return Vector3<T>(*this) *= a;}
-      inline Vector3<T> operator*(const float a)
-	 { return Vector3<T>(*this) *= a;}
-      inline Vector3<T> operator*(const double a)
-	 { return Vector3<T>(*this) *= a;}
+	 // methods
 
-      // Scalar division
-      inline Vector3<T>& operator/=(const int a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
-      inline Vector3<T>& operator/=(const float a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
-      inline Vector3<T>& operator/=(const double a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
-      inline Vector3<T> operator/(const int a)
-	 { return Vector3<T>(*this) /= a;}
-      inline Vector3<T> operator/(const float a)
-	 { return Vector3<T>(*this) /= a;}
-      inline Vector3<T> operator/(const double a)
-	 { return Vector3<T>(*this) /= a;}
+	 inline Vector3<T>& cross(Vector3<T>& v2, Vector3<T>& vres)
+	    { 
+	    vres.x =  y*v2.z - v2.y*z;
+	    vres.y = -x*v2.z + v2.x*z;
+	    vres.z =  x*v2.y - v2.x*y;
+	    return vres;
+	    }
+	 inline Vector3<T> cross(Vector3<T>& v2)
+	    { 
+	    Vector3<T> vres;
+	    vres.x =  y*v2.z - v2.y*z;
+	    vres.y = -x*v2.z + v2.x*z;
+	    vres.z =  x*v2.y - v2.x*y;
+	    return vres;
+	    }
 
+	 inline T dot(const Vector3<T>& v2) const
+	    { return x*v2.x + y*v2.y + z*v2.z; }
 
-      // methods
+	 inline double get_angle(Vector3<T>& v2)
+	    { 
+	    double tmp = dot(v2); 
+	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
+	    }
 
-      inline Vector3<T>& cross(Vector3<T>& v2, Vector3<T>& vres)
-	 { 
-	 vres.x =  y*v2.z - v2.y*z;
-	 vres.y = -x*v2.z + v2.x*z;
-	 vres.z =  x*v2.y - v2.x*y;
-	 return vres;
-	 }
-      inline Vector3<T> cross(Vector3<T>& v2)
-	 { 
-	 Vector3<T> vres;
-	 vres.x =  y*v2.z - v2.y*z;
-	 vres.y = -x*v2.z + v2.x*z;
-	 vres.z =  x*v2.y - v2.x*y;
-	 return vres;
-	 }
+	 inline double get_anglen(Vector3<T>& v2)
+	    { return acos((double) dot(v2)); }
 
-      inline T dot(const Vector3<T>& v2) const
-	 { return x*v2.x + y*v2.y + z*v2.z; }
+	 std::string & getstring(std::string & s) const;
 
-      inline double get_angle(Vector3<T>& v2)
-	 { 
-	 double tmp = dot(v2); 
-	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
-	 }
+	 inline double length(void) const
+	    { return sqrt((dot(*this))); }
 
-      inline double get_anglen(Vector3<T>& v2)
-	 { return acos((double) dot(v2)); }
+	 inline Vector3<T>& normalize()
+	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; return *this; }
 
-      std::string & getstring(std::string & s) const;
+	 inline Vector3<T> proj(Vector3<T>& v2)
+	    { Vector3<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 inline Vector3<T>& proj(Vector3<T>& v2, Vector3<T>& vres)
+	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 };
 
-      inline double length(void) const
-	 { return sqrt((dot(*this))); }
+      // Scalar multiplication, continued
+      template <typename T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
+	 { return Vector3<T>(v) *= a;}
+      template <typename T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
+	 { return Vector3<T>(v) *= a;}
+      template <typename T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
+	 { return Vector3<T>(v) *= a;}
 
-      inline Vector3<T>& normalize()
-	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; return *this; }
+      /////////////////////////////////////////////////////////////////////////////
+      template <typename T> class Vector4
+	 {
+	 public:
+	 T x, y, z, w;
 
-      inline Vector3<T> proj(Vector3<T>& v2)
-	 { Vector3<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      inline Vector3<T>& proj(Vector3<T>& v2, Vector3<T>& vres)
-	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      };
+	 // Constructors
+	 Vector4() {}
+	 explicit Vector4(const T a) : x (a), y(a), z(a), w(a) {}
+	 Vector4(const T a, const T b, const T c, const T d) : x (a), y(b), z(c), w(d) {}
 
-   // Scalar multiplication, continued
-   template <typename T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
-      { return Vector3<T>(v) *= a;}
-   template <typename T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
-      { return Vector3<T>(v) *= a;}
-   template <typename T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
-      { return Vector3<T>(v) *= a;}
+	 // Array indexing
+	 inline T& operator[](unsigned int i)
+	    { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
+	 inline T operator[](unsigned int i) const
+	    { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
 
-   /////////////////////////////////////////////////////////////////////////////
-   template <typename T> class Vector4
-      {
-      public:
-      T x, y, z, w;
+	 // Assignment
+	 inline Vector4<T>& operator=(const Vector4<T>& v2)
+	    { x = v2.x; y = v2.y; z = v2.z; w = v2.w; return *this; }
+	 inline Vector4<T>& assign(const T a, const T b, const T c, const T d)
+	    { x = a; y = b; z = c; w = d; return *this; }
 
-      // Constructors
-      Vector4() {}
-      Vector4(T a) : x (a), y(a), z(a), w(a) {}
-      Vector4(T a, T b, T c, T d) : x (a), y(b), z(c), w(d) {}
+	 // Comparison
+	 inline bool operator==(Vector4<T>& v2) const
+	    { return ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
+	 inline bool operator!=(Vector4<T>& v2) const
+	    { return ! (*this == v2); }
 
-      // Array indexing
-      inline T& operator[](unsigned int i)
-	 { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
-      inline const T& operator[](unsigned int i) const
-	 { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
+	 // Vector addition
+	 inline Vector4<T>& operator+=(const Vector4<T>& v2)
+	    { x += v2.x; y += v2.y; z += v2.z; w += v2.w; return *this; }
+	 inline Vector4<T> operator+(const Vector4<T>& v2) const
+	    { return Vector4<T>(*this) += v2; }
 
-      // Assignment
-      inline Vector4<T>& operator=(const Vector4<T>& v2)
-	 { x = v2.x; y = v2.y; z = v2.z; w = v2.w; return *this; }
-      inline Vector4<T>& assign(const T a, const T b, const T c, const T d)
-	 { x = a; y = b; z = c; w = d; return *this; }
+	 // Vector subtraction
+	 inline Vector4<T>& operator-=(const Vector4<T>& v2)
+	    { x -= v2.x; y -= v2.y; z -= v2.z; w -= v2.w; return *this; }
+	 inline Vector4<T> operator-(const Vector4<T>& v2) const
+	    { return Vector4<T>(*this) -= v2; }
 
-      // Comparison
-      inline bool operator==(Vector4<T>& v2) const
-	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
-      inline bool operator!=(Vector4<T>& v2) const
-	 { return ! (*this == v2); }
+	 // Scalar multiplication
+	 inline Vector4<T>& operator*=(const int a)
+	    { x *= a; y *= a; z *= a; w *= a; return *this; }
+	 inline Vector4<T>& operator*=(const float a)
+	    { x *= a; y *= a; z *= a; w *= a; return *this; }
+	 inline Vector4<T>& operator*=(const double a)
+	    { x *= a; y *= a; z *= a; w *= a; return *this; }
+	 inline Vector4<T> operator*(const int a)
+	    { return Vector4<T>(*this) *= a;}
+	 inline Vector4<T> operator*(const float a)
+	    { return Vector4<T>(*this) *= a;}
+	 inline Vector4<T> operator*(const double a)
+	    { return Vector4<T>(*this) *= a;}
 
-      // Vector addition
-      inline Vector4<T>& operator+=(const Vector4<T>& v2)
-	 { x += v2.x; y += v2.y; z += v2.z; w += v2.w; return *this; }
-      inline Vector4<T> operator+(const Vector4<T>& v2) const
-	 { return Vector4<T>(*this) += v2; }
+	 // Scalar division
+	 inline Vector4<T>& operator/=(const int a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
+	 inline Vector4<T>& operator/=(const float a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
+	 inline Vector4<T>& operator/=(const double a)
+	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
+	 inline Vector4<T> operator/(const int a)
+	    { return Vector4<T>(*this) /= a;}
+	 inline Vector4<T> operator/(const float a)
+	    { return Vector4<T>(*this) /= a;}
+	 inline Vector4<T> operator/(const double a)
+	    { return Vector4<T>(*this) /= a;}
 
-      // Vector subtraction
-      inline Vector4<T>& operator-=(const Vector4<T>& v2)
-	 { x -= v2.x; y -= v2.y; z -= v2.z; w -= v2.w; return *this; }
-      inline Vector4<T> operator-(const Vector4<T>& v2) const
-	 { return Vector4<T>(*this) -= v2; }
 
-      // Scalar multiplication
-      inline Vector4<T>& operator*=(const int a)
-	 { x *= a; y *= a; z *= a; w *= a; return *this; }
-      inline Vector4<T>& operator*=(const float a)
-	 { x *= a; y *= a; z *= a; w *= a; return *this; }
-      inline Vector4<T>& operator*=(const double a)
-	 { x *= a; y *= a; z *= a; w *= a; return *this; }
-      inline Vector4<T> operator*(const int a)
-	 { return Vector4<T>(*this) *= a;}
-      inline Vector4<T> operator*(const float a)
-	 { return Vector4<T>(*this) *= a;}
-      inline Vector4<T> operator*(const double a)
-	 { return Vector4<T>(*this) *= a;}
+	 // methods
 
-      // Scalar division
-      inline Vector4<T>& operator/=(const int a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
-      inline Vector4<T>& operator/=(const float a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
-      inline Vector4<T>& operator/=(const double a)
-	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
-      inline Vector4<T> operator/(const int a)
-	 { return Vector4<T>(*this) /= a;}
-      inline Vector4<T> operator/(const float a)
-	 { return Vector4<T>(*this) /= a;}
-      inline Vector4<T> operator/(const double a)
-	 { return Vector4<T>(*this) /= a;}
+	 inline T dot(const Vector4<T>& v2) const
+	    { return x*v2.x + y*v2.y + z*v2.z + w*v2.w; }
 
+	 inline double get_angle(Vector4<T>& v2)
+	    { 
+	    double tmp = dot(v2); 
+	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
+	    }
 
-      // methods
+	 inline double get_anglen(Vector4<T>& v2)
+	    { return acos((double) dot(v2)); }
 
-      inline T dot(const Vector4<T>& v2) const
-	 { return x*v2.x + y*v2.y + z*v2.z + w*v2.w; }
+	 std::string & getstring(std::string & s) const;
 
-      inline double get_angle(Vector4<T>& v2)
-	 { 
-	 double tmp = dot(v2); 
-	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
-	 }
+	 inline double length(void) const
+	    { return sqrt((dot(*this))); }
 
-      inline double get_anglen(Vector4<T>& v2)
-	 { return acos((double) dot(v2)); }
+	 inline Vector4<T>& normalize()
+	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; w /= l; return *this; }
 
-      std::string & getstring(std::string & s) const;
+	 inline Vector4<T> proj(Vector4<T>& v2)
+	    { Vector4<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 inline Vector4<T>& proj(Vector4<T>& v2, Vector4<T>& vres)
+	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
+	 };
 
-      inline double length(void) const
-	 { return sqrt((dot(*this))); }
+      // Scalar multiplication, continued
+      template <typename T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
+	 { return Vector4<T>(v) *= a;}
+      template <typename T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
+	 { return Vector4<T>(v) *= a;}
+      template <typename T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
+	 { return Vector4<T>(v) *= a;}
 
-      inline Vector4<T>& normalize()
-	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; w /= l; return *this; }
+      /////////////////////////////////////////////////////////////////////////////
 
-      inline Vector4<T> proj(Vector4<T>& v2)
-	 { Vector4<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      inline Vector4<T>& proj(Vector4<T>& v2, Vector4<T>& vres)
-	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
-      };
+      typedef Vector2<int> Vector2i;
+      typedef Vector2<float> Vector2f;
+      typedef Vector2<double> Vector2d;
 
-   // Scalar multiplication, continued
-   template <typename T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
-      { return Vector4<T>(v) *= a;}
-   template <typename T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
-      { return Vector4<T>(v) *= a;}
-   template <typename T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
-      { return Vector4<T>(v) *= a;}
+      typedef Vector3<int> Vector3i;
+      typedef Vector3<float> Vector3f;
+      typedef Vector3<double> Vector3d;
 
-   /////////////////////////////////////////////////////////////////////////////
+      typedef Vector4<int> Vector4i;
+      typedef Vector4<float> Vector4f;
+      typedef Vector4<double> Vector4d;
 
-   typedef Vector2<int> Vector2i;
-   typedef Vector2<float> Vector2f;
-   typedef Vector2<double> Vector2d;
-
-   typedef Vector3<int> Vector3i;
-   typedef Vector3<float> Vector3f;
-   typedef Vector3<double> Vector3d;
-
-   typedef Vector4<int> Vector4i;
-   typedef Vector4<float> Vector4f;
-   typedef Vector4<double> Vector4d;
+      } // namespace Math
 
    } // namespace arda
 
 ////////////////////////////////////////////////////////////////////////////////
-template <typename T> std::string & arda::Vector2<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Math::Vector2<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
@@ -462,7 +466,7 @@
    }
 
 ////////////////////////////////////////////////////////////////////////////////
-template <typename T> std::string & arda::Vector3<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Math::Vector3<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
@@ -472,7 +476,7 @@
    }
 
 ////////////////////////////////////////////////////////////////////////////////
-template <typename T> std::string & arda::Vector4<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Math::Vector4<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
diff -r 139d189290e0 -r 398894c52b69 src/Matrix.cpp
--- a/src/Matrix.cpp	Fri Oct 07 11:09:15 2011 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,515 +0,0 @@
-#include "Matrix.h"
-
-#include <iostream>
-using namespace std;
-#include <string>
-#include <sstream>
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring2m(Matrix2i m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << " ], "
-	 "[ " << m[2] << ", " << m[3] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring2m(Matrix2f m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << " ], "
-	 "[ " << m[2] << ", " << m[3] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring2m(Matrix2d m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << " ], "
-	 "[ " << m[2] << ", " << m[3] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring3m(Matrix3i m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
-	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
-	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring3m(Matrix3f m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
-	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
-	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring3m(Matrix3d m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
-	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
-	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring4m(Matrix4i m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
-	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
-	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
-	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring4m(Matrix4f m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
-	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
-	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
-	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Matrix::getstring4m(Matrix4d m, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "[ "
-	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
-	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
-	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
-	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec2mv(Matrix2i m, Vector2i v, Vector2i vres)
-   {
-   vres[0] = m[0]*v[0] + m[2]*v[1];
-   vres[1] = m[1]*v[0] + m[3]*v[1];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec2mv(Matrix2f m, Vector2f v, Vector2f vres)
-   {
-   vres[0] = m[0]*v[0] + m[2]*v[1];
-   vres[1] = m[1]*v[0] + m[3]*v[1];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec2mv(Matrix2d m, Vector2d v, Vector2d vres)
-   {
-   vres[0] = m[0]*v[0] + m[2]*v[1];
-   vres[1] = m[1]*v[0] + m[3]*v[1];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec3mv(Matrix3i m, Vector3i v, Vector3i vres)
-   {
-   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
-   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
-   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec3mv(Matrix3f m, Vector3f v, Vector3f vres)
-   {
-   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
-   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
-   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec3mv(Matrix3d m, Vector3d v, Vector3d vres)
-   {
-   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
-   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
-   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec4mv(Matrix4i m, Vector4i v, Vector4i vres)
-   {
-   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
-   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
-   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
-   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec4mv(Matrix4f m, Vector4f v, Vector4f vres)
-   {
-   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
-   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
-   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
-   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
-   return vres;
-  
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec4mv(Matrix4d m, Vector4d v, Vector4d vres)
-   {
-   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
-   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
-   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
-   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
-   return vres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec2vm(Vector2i v, Matrix2i m, Vector2i vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1];
-   vres[1] = v[0]*m[2] + v[1]*m[3];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec2vm(Vector2f v, Matrix2f m, Vector2f vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1];
-   vres[1] = v[0]*m[2] + v[1]*m[3];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec2vm(Vector2d v, Matrix2d m, Vector2d vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1];
-   vres[1] = v[0]*m[2] + v[1]*m[3];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec3vm(Vector3i v, Matrix3i m, Vector3i vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
-   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
-   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec3vm(Vector3f v, Matrix3f m, Vector3f vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
-   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
-   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec3vm(Vector3d v, Matrix3d m, Vector3d vres)
-   {
-   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
-   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
-   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multvec4vm(Vector4i v, Matrix4i m, Vector4i vres)
-   {
-   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
-   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
-   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
-   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multvec4vm(Vector4f v, Matrix4f m, Vector4f vres)
-   {
-   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
-   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
-   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
-   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
-   return vres;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multvec4vm(Vector4d v, Matrix4d m, Vector4d vres)
-   {
-   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
-   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
-   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
-   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
-   return vres;   
-   }
-
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multmat2(Matrix2i m1, Matrix2i m2, Matrix2i mres)
-   {
-   const int size=2;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multmat2(Matrix2f m1, Matrix2f m2, Matrix2f mres)
-   {
-   const int size=2;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multmat2(Matrix2d m1, Matrix2d m2, Matrix2d mres)
-   {
-   const int size=2;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multmat3(Matrix3i m1, Matrix3i m2, Matrix3i mres)
-   {
-   const int size=3;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multmat3(Matrix3f m1, Matrix3f m2, Matrix3f mres)
-   {
-   const int size=3;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multmat3(Matrix3d m1, Matrix3d m2, Matrix3d mres)
-   {
-   const int size=3;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::multmat4(Matrix4i m1, Matrix4i m2, Matrix4i mres)
-   {
-   const int size=4;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::multmat4(Matrix4f m1, Matrix4f m2, Matrix4f mres)
-   {
-   const int size=4;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::multmat4(Matrix4d m1, Matrix4d m2, Matrix4d mres)
-   {
-   cout << "here" << endl;
-   const int size=4;
-   int i, j, k;
-   for (i=0; i<size*size; ++i)
-      mres[i] = 0;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; ++j)
-	 for (k=0; k<size; ++k)
-	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::transpose2(Matrix2i m, Matrix2i mres)
-   {
-   const int size=2;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::transpose2(Matrix2f m, Matrix2f mres)
-   {
-   const int size=2;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::transpose2(Matrix2d m, Matrix2d mres)
-   {
-   const int size=2;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::transpose3(Matrix3i m, Matrix3i mres)
-   {
-   const int size=3;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::transpose3(Matrix3f m, Matrix3f mres)
-   {
-   const int size=3;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::transpose3(Matrix3d m, Matrix3d mres)
-   {
-   const int size=3;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-int * arda::Matrix::transpose4(Matrix4i m, Matrix4i mres)
-   {
-   const int size=4;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-float * arda::Matrix::transpose4(Matrix4f m, Matrix4f mres)
-   {
-   const int size=4;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-double * arda::Matrix::transpose4(Matrix4d m, Matrix4d mres)
-   {
-   const int size=4;
-   int i, j;
-   for (i=0; i<size; ++i)
-      for (j=0; j<size; j++)
-	 mres[size*i+j] = m[size*j+i];
-   return mres;
-   }
-
-
diff -r 139d189290e0 -r 398894c52b69 src/main.cpp
--- a/src/main.cpp	Fri Oct 07 11:09:15 2011 -0500
+++ b/src/main.cpp	Fri Oct 07 13:33:48 2011 -0500
@@ -6,7 +6,7 @@
 using namespace std;
 
 #include "Math.h"
-using namespace arda;
+using namespace arda::Math;
 
 ////////////////////////////////////////////////////////////////////////////////
 void test_vector(void)
@@ -588,8 +588,8 @@
    double d;
 
    cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
-   d = m22i_1.det();
-   cout << setw(40) << "d = m22i_1.det(): " << d << endl;
+   d = det(m22i_1);
+   cout << setw(40) << "d = det(m22i_1): " << d << endl;
    }
 
    {
@@ -599,8 +599,8 @@
    Matrix22i m22i_2;
 
    cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
-   m22i_2 = m22i_1.transpose();
-   cout << setw(40) << "m22i_2 = m22i_1.transpose(): " << m22i_2.getstring(s1) << endl;
+   m22i_2 = transpose(m22i_1);
+   cout << setw(40) << "m22i_2 = transpose(m22i_1): " << m22i_2.getstring(s1) << endl;
    }
 
    {
@@ -612,6 +612,55 @@
    m22i_1.setidentity();
    cout << setw(40) << "m22i_1.setidentity(): " << m22i_1.getstring(s1) << endl;
    }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix getrow() and getcol()" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Vector2i v2i_1(0);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   v2i_1 = m22i_1.getrow(0);
+   cout << setw(40) << "v2i_1 = m22i_1.getrow(0): " << v2i_1.getstring(s1) << endl;
+   v2i_1 = m22i_1.getrow(1);
+   cout << setw(40) << "v2i_1 = m22i_1.getrow(1): " << v2i_1.getstring(s1) << endl;
+   v2i_1 = m22i_1.getcol(0);
+   cout << setw(40) << "v2i_1 = m22i_1.getcol(0): " << v2i_1.getstring(s1) << endl;
+   v2i_1 = m22i_1.getcol(1);
+   cout << setw(40) << "v2i_1 = m22i_1.getcol(1): " << v2i_1.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix setrow() and setcol()" << endl;
+   Matrix22i m22i_1;
+   Vector2i v2i_1(2, 3);
+
+   m22i_1.setidentity();
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   m22i_1.setrow(0, v2i_1);
+   cout << setw(40) << "m22i_1.setrow(0, v2i_1): " << m22i_1.getstring(s1) << endl;
+
+   m22i_1.setidentity();
+   m22i_1.setrow(1, v2i_1);
+   cout << setw(40) << "m22i_1.setrow(1, v2i_1): " << m22i_1.getstring(s1) << endl;
+
+   m22i_1.setidentity();
+   m22i_1.setrow(1, 4, 5);
+   cout << setw(40) << "m22i_1.setrow(1, 4, 5): " << m22i_1.getstring(s1) << endl;
+
+   m22i_1.setidentity();
+   m22i_1.setcol(0, v2i_1);
+   cout << setw(40) << "m22i_1.setcol(0, v2i_1): " << m22i_1.getstring(s1) << endl;
+
+   m22i_1.setidentity();
+   m22i_1.setcol(1, v2i_1);
+   cout << setw(40) << "m22i_1.setcol(1, v2i_1): " << m22i_1.getstring(s1) << endl;
+
+   m22i_1.setidentity();
+   m22i_1.setcol(1, 4, 5);
+   cout << setw(40) << "m22i_1.setcol(1, 4, 5): " << m22i_1.getstring(s1) << endl;
+   }
    }
 
 ////////////////////////////////////////////////////////////////////////////////