--- trunk/src/math/Vector.hpp	2005/01/13 19:40:37	253
+++ trunk/src/math/Vector.hpp	2006/05/17 21:51:42	963
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -53,440 +53,442 @@
 #include <cmath>
 #include <iostream>
 #include <math.h>
+#include "config.h"
 namespace oopse {
 
-    static const double epsilon = 0.000001;
+  static const RealType epsilon = 0.000001;
 
-    template<typename T>
-    inline bool equal(T e1, T e2) {
-        return e1 == e2;
-    }
+  template<typename T>
+  inline bool equal(T e1, T e2) {
+    return e1 == e2;
+  }
 
-    template<>
-    inline bool equal(float e1, float e2) {
-        return fabs(e1 - e2) < epsilon;
-    }
+  //template<>
+  //inline bool equal(float e1, float e2) {
+  //  return fabs(e1 - e2) < epsilon;
+  //}
 
-    template<>
-    inline bool equal(double e1, double e2) {
-        return fabs(e1 - e2) < epsilon;
-    }
+  template<>
+  inline bool equal(RealType e1, RealType e2) {
+    return fabs(e1 - e2) < epsilon;
+  }
 
     
-    /**
-     * @class Vector Vector.hpp "math/Vector.hpp"
-     * @brief Fix length vector class
-     */
-    template<typename Real, unsigned int Dim>
-    class Vector{
-        public:
+  /**
+   * @class Vector Vector.hpp "math/Vector.hpp"
+   * @brief Fix length vector class
+   */
+  template<typename Real, unsigned int Dim>
+  class Vector{
+  public:
 
-            typedef Real ElemType;
-            typedef Real* ElemPoinerType;
+    typedef Real ElemType;
+    typedef Real* ElemPoinerType;
 
-            /** default constructor */
-            inline Vector(){
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = 0;
-            }
+    /** default constructor */
+    inline Vector(){
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = 0;
+    }
 
-            /** Constructs and initializes a Vector from a vector */
-            inline Vector(const Vector<Real, Dim>& v) {
-                *this  = v;
-            }
+    /** Constructs and initializes a Vector from a vector */
+    inline Vector(const Vector<Real, Dim>& v) {
+      *this  = v;
+    }
 
-            /** copy assignment operator */
-            inline Vector<Real, Dim>& operator=(const Vector<Real, Dim>& v) {
-                if (this == &v)
-                    return *this;
+    /** copy assignment operator */
+    inline Vector<Real, Dim>& operator=(const Vector<Real, Dim>& v) {
+      if (this == &v)
+	return *this;
                 
-            	for (unsigned int i = 0; i < Dim; i++)            
-                    data_[i] = v[i];
+      for (unsigned int i = 0; i < Dim; i++)            
+	this->data_[i] = v[i];
                 
-                return *this;
-            }
+      return *this;
+    }
 
-            template<typename T>
-            inline Vector(const T& s){
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = s;
-            }
+    template<typename T>
+    inline Vector(const T& s){
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = s;
+    }
             
-            /** Constructs and initializes a Vector from an array */            
-            inline Vector( Real* v) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = v[i];
-            }
+    /** Constructs and initializes a Vector from an array */            
+    inline Vector( Real* v) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = v[i];
+    }
 
-            /** 
-             * Returns reference of ith element.
-             * @return reference of ith element
-             * @param i index
-             */
-            inline Real& operator[](unsigned int  i) {
-                assert( i < Dim);
-                return data_[i];
-            }
+    /** 
+     * Returns reference of ith element.
+     * @return reference of ith element
+     * @param i index
+     */
+    inline Real& operator[](unsigned int  i) {
+      assert( i < Dim);
+      return this->data_[i];
+    }
 
-            /** 
-             * Returns reference of ith element.
-             * @return reference of ith element
-             * @param i index
-             */
-            inline Real& operator()(unsigned int  i) {
-                assert( i < Dim);
-                return data_[i];
-            }
+    /** 
+     * Returns reference of ith element.
+     * @return reference of ith element
+     * @param i index
+     */
+    inline Real& operator()(unsigned int  i) {
+      assert( i < Dim);
+      return this->data_[i];
+    }
 
-            /** 
-             * Returns constant reference of ith element.
-             * @return reference of ith element
-             * @param i index
-             */
-            inline  const Real& operator[](unsigned int i) const {
-                assert( i < Dim);
-                return data_[i];
-            }
+    /** 
+     * Returns constant reference of ith element.
+     * @return reference of ith element
+     * @param i index
+     */
+    inline  const Real& operator[](unsigned int i) const {
+      assert( i < Dim);
+      return this->data_[i];
+    }
 
-            /** 
-             * Returns constant reference of ith element.
-             * @return reference of ith element
-             * @param i index
-             */
-            inline  const Real& operator()(unsigned int i) const {
-                assert( i < Dim);
-                return data_[i];
-            }
+    /** 
+     * Returns constant reference of ith element.
+     * @return reference of ith element
+     * @param i index
+     */
+    inline  const Real& operator()(unsigned int i) const {
+      assert( i < Dim);
+      return this->data_[i];
+    }
 
-            /** Copy the internal data to an array*/
-            void getArray(Real* array) {
-                for (unsigned int i = 0; i < Dim; i ++) {
-                    array[i] = data_[i];
-                }                
-            }
+    /** Copy the internal data to an array*/
+    void getArray(Real* array) {
+      for (unsigned int i = 0; i < Dim; i ++) {
+	array[i] = this->data_[i];
+      }                
+    }
 
-            /** Returns the pointer of internal array */
-            Real* getArrayPointer() {
-                return data_;
-            }
+    /** Returns the pointer of internal array */
+    Real* getArrayPointer() {
+      return this->data_;
+    }
             
-            /**
-             * Tests if this vetor is equal to other vector
-             * @return true if equal, otherwise return false
-             * @param v vector to be compared
-             */
-             inline bool operator ==(const Vector<Real, Dim>& v) {
-
-                for (unsigned int i = 0; i < Dim; i ++) {
-                    if (!equal(data_[i], v[i])) {
-                        return false;
-                    }
-                }
+    /**
+     * Tests if this vetor is equal to other vector
+     * @return true if equal, otherwise return false
+     * @param v vector to be compared
+     */
+    inline bool operator ==(const Vector<Real, Dim>& v) {
+
+      for (unsigned int i = 0; i < Dim; i ++) {
+	if (!equal(this->data_[i], v[i])) {
+	  return false;
+	}
+      }
                 
-                return true;
-            }
+      return true;
+    }
 
-            /**
-             * Tests if this vetor is not equal to other vector
-             * @return true if equal, otherwise return false
-             * @param v vector to be compared
-             */
-            inline bool operator !=(const Vector<Real, Dim>& v) {
-                return !(*this == v);
-            }
+    /**
+     * Tests if this vetor is not equal to other vector
+     * @return true if equal, otherwise return false
+     * @param v vector to be compared
+     */
+    inline bool operator !=(const Vector<Real, Dim>& v) {
+      return !(*this == v);
+    }
              
-            /** Negates the value of this vector in place. */           
-            inline void negate() {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = -data_[i];
-            }
+    /** Negates the value of this vector in place. */           
+    inline void negate() {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = -this->data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the negation of vector v1.
-            * @param v1 the source vector
-            */
-            inline void negate(const Vector<Real, Dim>& v1) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = -v1.data_[i];
+    /**
+     * Sets the value of this vector to the negation of vector v1.
+     * @param v1 the source vector
+     */
+    inline void negate(const Vector<Real, Dim>& v1) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = -v1.data_[i];
 
-            }
+    }
             
-            /**
-            * Sets the value of this vector to the sum of itself and v1 (*this += v1).
-            * @param v1 the other vector
-            */
-            inline void add( const Vector<Real, Dim>& v1 ) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] += v1.data_[i];
-            }
+    /**
+     * Sets the value of this vector to the sum of itself and v1 (*this += v1).
+     * @param v1 the other vector
+     */
+    inline void add( const Vector<Real, Dim>& v1 ) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] += v1.data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the sum of v1 and v2 (*this = v1 + v2).
-            * @param v1 the first vector
-            * @param v2 the second vector
-            */
-            inline void add( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = v1.data_[i] + v2.data_[i];
-            }
+    /**
+     * Sets the value of this vector to the sum of v1 and v2 (*this = v1 + v2).
+     * @param v1 the first vector
+     * @param v2 the second vector
+     */
+    inline void add( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = v1.data_[i] + v2.data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the difference  of itself and v1 (*this -= v1).
-            * @param v1 the other vector
-            */
-            inline void sub( const Vector<Real, Dim>& v1 ) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] -= v1.data_[i];
-            }
+    /**
+     * Sets the value of this vector to the difference  of itself and v1 (*this -= v1).
+     * @param v1 the other vector
+     */
+    inline void sub( const Vector<Real, Dim>& v1 ) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] -= v1.data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the difference of vector v1 and v2 (*this = v1 - v2).
-            * @param v1 the first vector
-            * @param v2 the second vector
-            */
-            inline void sub( const Vector<Real, Dim>& v1, const Vector  &v2 ){
-            	for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = v1.data_[i] - v2.data_[i];
-            }
+    /**
+     * Sets the value of this vector to the difference of vector v1 and v2 (*this = v1 - v2).
+     * @param v1 the first vector
+     * @param v2 the second vector
+     */
+    inline void sub( const Vector<Real, Dim>& v1, const Vector  &v2 ){
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = v1.data_[i] - v2.data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the scalar multiplication of itself (*this *= s).
-            * @param s the scalar value
-            */
-            inline void mul( Real s ) {
-                for (unsigned int i = 0; i < Dim; i++)
-                   data_[i] *= s;
-            }
+    /**
+     * Sets the value of this vector to the scalar multiplication of itself (*this *= s).
+     * @param s the scalar value
+     */
+    inline void mul( Real s ) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] *= s;
+    }
 
-            /**
-            * Sets the value of this vector to the scalar multiplication of vector v1  
-            * (*this = s * v1).
-            * @param v1 the vector            
-            * @param s the scalar value
-            */
-            inline void mul( const Vector<Real, Dim>& v1, Real s) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = s * v1.data_[i];
-            }
+    /**
+     * Sets the value of this vector to the scalar multiplication of vector v1  
+     * (*this = s * v1).
+     * @param v1 the vector            
+     * @param s the scalar value
+     */
+    inline void mul( const Vector<Real, Dim>& v1, Real s) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = s * v1.data_[i];
+    }
 
-            /**
-            * Sets the value of this vector to the scalar division of itself  (*this /= s ).
-            * @param s the scalar value
-            */             
-            inline void div( Real s) {
-                for (unsigned int i = 0; i < Dim; i++)            
-                    data_[i] /= s;
-            }
-
-            /**
-            * Sets the value of this vector to the scalar division of vector v1  (*this = v1 / s ).
-            * @param v1 the source vector
-            * @param s the scalar value
-            */                         
-            inline void div( const Vector<Real, Dim>& v1, Real s ) {
-                for (unsigned int i = 0; i < Dim; i++)
-                    data_[i] = v1.data_[i] / s;
-            }
+    /**
+     * Sets the value of this vector to the scalar division of itself  (*this /= s ).
+     * @param s the scalar value
+     */             
+    inline void div( Real s) {
+      for (unsigned int i = 0; i < Dim; i++)            
+	this->data_[i] /= s;
+    }
 
-            /** @see #add */
-            inline Vector<Real, Dim>& operator +=( const Vector<Real, Dim>& v1 ) {
-                add(v1);
-                return *this;
-            }
+    /**
+     * Sets the value of this vector to the scalar division of vector v1  (*this = v1 / s ).
+     * @param v1 the source vector
+     * @param s the scalar value
+     */                         
+    inline void div( const Vector<Real, Dim>& v1, Real s ) {
+      for (unsigned int i = 0; i < Dim; i++)
+	this->data_[i] = v1.data_[i] / s;
+    }
 
-            /** @see #sub */
-            inline Vector<Real, Dim>& operator -=( const Vector<Real, Dim>& v1 ) {
-                sub(v1);
-                return *this;
-            }
+    /** @see #add */
+    inline Vector<Real, Dim>& operator +=( const Vector<Real, Dim>& v1 ) {
+      add(v1);
+      return *this;
+    }
 
-            /** @see #mul */
-            inline Vector<Real, Dim>& operator *=( Real s) {
-                mul(s);
-                return *this;
-            }
+    /** @see #sub */
+    inline Vector<Real, Dim>& operator -=( const Vector<Real, Dim>& v1 ) {
+      sub(v1);
+      return *this;
+    }
 
-            /** @see #div */
-            inline Vector<Real, Dim>& operator /=( Real s ) {
-                div(s);
-                return *this;
-            }
+    /** @see #mul */
+    inline Vector<Real, Dim>& operator *=( Real s) {
+      mul(s);
+      return *this;
+    }
 
-            /**
-             * Returns the length of this vector.
-             * @return the length of this vector
-             */
-             inline Real length() {
-                return sqrt(lengthSquare());  
-            }
+    /** @see #div */
+    inline Vector<Real, Dim>& operator /=( Real s ) {
+      div(s);
+      return *this;
+    }
+
+    /**
+     * Returns the length of this vector.
+     * @return the length of this vector
+     */
+    inline Real length() {
+      return sqrt(lengthSquare());  
+    }
             
-            /**
-             * Returns the squared length of this vector.
-             * @return the squared length of this vector
-             */
-             inline Real lengthSquare() {
-                return dot(*this, *this);
-            }
+    /**
+     * Returns the squared length of this vector.
+     * @return the squared length of this vector
+     */
+    inline Real lengthSquare() {
+      return dot(*this, *this);
+    }
             
-            /** Normalizes this vector in place */
-            inline void normalize() {
-                Real len;
+    /** Normalizes this vector in place */
+    inline void normalize() {
+      Real len;
 
-                len = length();
+      len = length();
                 
-                //if (len < oopse:epsilon)
-                //  throw();
+      //if (len < oopse:epsilon)
+      //  throw();
                 
-                *this /= len;
-            }
+      *this /= len;
+    }
 
-            /**
-             * Tests if this vector is normalized
-             * @return true if this vector is normalized, otherwise return false
-             */
-            inline bool isNormalized() {
-                return equal(lengthSquare(), 1.0);
-            }           
-            
-        protected:
-            Real data_[Dim];
+    /**
+     * Tests if this vector is normalized
+     * @return true if this vector is normalized, otherwise return false
+     */
+    inline bool isNormalized() {
+      return equal(lengthSquare(), (RealType)1);
+    }           
+
+    unsigned int size() {return Dim;}
+  protected:
+    Real data_[Dim];
         
-    };
+  };
 
-    /** unary minus*/
-    template<typename Real, unsigned int Dim>    
-    inline Vector<Real, Dim> operator -(const Vector<Real, Dim>& v1){
-        Vector<Real, Dim> tmp(v1);
-        tmp.negate();
-        return tmp;
-    }
+  /** unary minus*/
+  template<typename Real, unsigned int Dim>    
+  inline Vector<Real, Dim> operator -(const Vector<Real, Dim>& v1){
+    Vector<Real, Dim> tmp(v1);
+    tmp.negate();
+    return tmp;
+  }
 
-    /**
-     * Return the sum of two vectors  (v1 - v2). 
-     * @return the sum of two vectors
-     * @param v1 the first vector
-     * @param v2 the second vector
-     */   
-    template<typename Real, unsigned int Dim>    
-    inline Vector<Real, Dim> operator +(const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2) {
-        Vector<Real, Dim> result;
+  /**
+   * Return the sum of two vectors  (v1 - v2). 
+   * @return the sum of two vectors
+   * @param v1 the first vector
+   * @param v2 the second vector
+   */   
+  template<typename Real, unsigned int Dim>    
+  inline Vector<Real, Dim> operator +(const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2) {
+    Vector<Real, Dim> result;
         
-        result.add(v1, v2);
-        return result;        
-    }
+    result.add(v1, v2);
+    return result;        
+  }
 
-    /**
-     * Return the difference of two vectors  (v1 - v2). 
-     * @return the difference of two vectors
-     * @param v1 the first vector
-     * @param v2 the second vector
-     */  
-    template<typename Real, unsigned int Dim>    
-    Vector<Real, Dim> operator -(const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2) {
-        Vector<Real, Dim> result;
-        result.sub(v1, v2);
-        return result;        
-    }
+  /**
+   * Return the difference of two vectors  (v1 - v2). 
+   * @return the difference of two vectors
+   * @param v1 the first vector
+   * @param v2 the second vector
+   */  
+  template<typename Real, unsigned int Dim>    
+  Vector<Real, Dim> operator -(const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2) {
+    Vector<Real, Dim> result;
+    result.sub(v1, v2);
+    return result;        
+  }
     
-    /**
-     * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). 
-     * @return  the vaule of scalar multiplication of this vector
-     * @param v1 the source vector
-     * @param s the scalar value
-     */ 
-    template<typename Real, unsigned int Dim>                 
-    Vector<Real, Dim> operator * ( const Vector<Real, Dim>& v1, Real s) {       
-        Vector<Real, Dim> result;
-        result.mul(v1,s);
-        return result;           
-    }
+  /**
+   * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). 
+   * @return  the vaule of scalar multiplication of this vector
+   * @param v1 the source vector
+   * @param s the scalar value
+   */ 
+  template<typename Real, unsigned int Dim>                 
+  Vector<Real, Dim> operator * ( const Vector<Real, Dim>& v1, Real s) {       
+    Vector<Real, Dim> result;
+    result.mul(v1,s);
+    return result;           
+  }
     
-    /**
-     * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). 
-     * @return  the vaule of scalar multiplication of this vector
-     * @param s the scalar value
-     * @param v1 the source vector
-     */  
-    template<typename Real, unsigned int Dim>
-    Vector<Real, Dim> operator * ( Real s, const Vector<Real, Dim>& v1 ) {
-        Vector<Real, Dim> result;
-        result.mul(v1, s);
-        return result;           
-    }
+  /**
+   * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). 
+   * @return  the vaule of scalar multiplication of this vector
+   * @param s the scalar value
+   * @param v1 the source vector
+   */  
+  template<typename Real, unsigned int Dim>
+  Vector<Real, Dim> operator * ( Real s, const Vector<Real, Dim>& v1 ) {
+    Vector<Real, Dim> result;
+    result.mul(v1, s);
+    return result;           
+  }
 
-    /**
-     * Returns the  value of division of a vector by a scalar. 
-     * @return  the vaule of scalar division of this vector
-     * @param v1 the source vector
-     * @param s the scalar value
-     */
-    template<typename Real, unsigned int Dim>    
-    Vector<Real, Dim> operator / ( const Vector<Real, Dim>& v1, Real s) {       
-        Vector<Real, Dim> result;
-        result.div( v1,s);
-        return result;           
-    }
+  /**
+   * Returns the  value of division of a vector by a scalar. 
+   * @return  the vaule of scalar division of this vector
+   * @param v1 the source vector
+   * @param s the scalar value
+   */
+  template<typename Real, unsigned int Dim>    
+  Vector<Real, Dim> operator / ( const Vector<Real, Dim>& v1, Real s) {       
+    Vector<Real, Dim> result;
+    result.div( v1,s);
+    return result;           
+  }
     
-    /**
-     * Returns the dot product of two Vectors
-     * @param v1 first vector
-     * @param v2 second vector
-     * @return the dot product of v1 and v2
-     */
-    template<typename Real, unsigned int Dim>    
-    inline Real dot( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
-        Real tmp;
-        tmp = 0;
+  /**
+   * Returns the dot product of two Vectors
+   * @param v1 first vector
+   * @param v2 second vector
+   * @return the dot product of v1 and v2
+   */
+  template<typename Real, unsigned int Dim>    
+  inline Real dot( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
+    Real tmp;
+    tmp = 0;
 
-        for (unsigned int i = 0; i < Dim; i++)
-            tmp += v1[i] * v2[i];
+    for (unsigned int i = 0; i < Dim; i++)
+      tmp += v1[i] * v2[i];
 
-        return tmp;
-    }
+    return tmp;
+  }
 
-    /**
-     * Returns the distance between  two Vectors
-     * @param v1 first vector
-     * @param v2 second vector
-     * @return the distance between v1 and v2
-     */	
-    template<typename Real, unsigned int Dim>    
-    inline Real distance( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
-        Vector<Real, Dim> tempVector = v1 - v2;
-        return tempVector.length();
-    }
+  /**
+   * Returns the distance between  two Vectors
+   * @param v1 first vector
+   * @param v2 second vector
+   * @return the distance between v1 and v2
+   */	
+  template<typename Real, unsigned int Dim>    
+  inline Real distance( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
+    Vector<Real, Dim> tempVector = v1 - v2;
+    return tempVector.length();
+  }
 
-    /**
-     * Returns the squared distance between  two Vectors
-     * @param v1 first vector
-     * @param v2 second vector
-     * @return the squared distance between v1 and v2
-     */
-    template<typename Real, unsigned int Dim>
-    inline Real distanceSquare( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
-        Vector<Real, Dim> tempVector = v1 - v2;
-        return tempVector.lengthSquare();
-    }
+  /**
+   * Returns the squared distance between  two Vectors
+   * @param v1 first vector
+   * @param v2 second vector
+   * @return the squared distance between v1 and v2
+   */
+  template<typename Real, unsigned int Dim>
+  inline Real distanceSquare( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) {
+    Vector<Real, Dim> tempVector = v1 - v2;
+    return tempVector.lengthSquare();
+  }
 
-    /**
-     * Write to an output stream
-     */
-    template<typename Real, unsigned int Dim>
-    std::ostream &operator<< ( std::ostream& o, const Vector<Real, Dim>& v) {
+  /**
+   * Write to an output stream
+   */
+  template<typename Real, unsigned int Dim>
+  std::ostream &operator<< ( std::ostream& o, const Vector<Real, Dim>& v) {
 
-        o << "[ ";
+    o << "[ ";
         
-        for (unsigned int i = 0 ; i< Dim; i++) {
-            o << v[i];
+    for (unsigned int i = 0 ; i< Dim; i++) {
+      o << v[i];
 
-            if (i  != Dim -1) {
-                o<< ", ";
-            }
-        }
-
-        o << " ]";
-        return o;        
+      if (i  != Dim -1) {
+	o<< ", ";
+      }
     }
+
+    o << " ]";
+    return o;        
+  }
     
 }
 #endif