diff --git a/src/vector.cpp b/src/vector.cpp
new file mode 100644
--- /dev/null
+++ b/src/vector.cpp
@@ -0,0 +1,341 @@
+/*
+ *
+ *  This file is part of the Virtual Leaf.
+ *
+ *  The Virtual Leaf is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  The Virtual Leaf is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  Copyright 2010 Roeland Merks.
+ *
+ */
+
+#include <string>
+#include <stddef.h>
+#include <ostream>
+#include <limits.h>
+#include <stdio.h>
+//#include <math.h>
+#include <stdarg.h>
+#include "sqr.h"
+#include "pi.h"
+#include "vector.h"
+#include "tiny.h"
+
+static const std::string _module_id("$Id$");
+
+void Vector::operator=(const Vector &source) {
+  
+  // assignment
+  
+  // don't assign to self
+  if (this==&source) return;
+  
+  x=source.x;
+  y=source.y;
+  z=source.z;
+  
+}
+
+
+
+
+
+ostream &Vector::print(ostream &os) const {
+  os << "(" << x << ", " << y << ", " << z << ")";
+  return os;
+}
+
+ostream &operator<<(ostream &os, const Vector &v) {
+  v.print(os);
+  return os;
+}
+
+
+Vector Vector::operator+(const Vector &v) const {
+  
+  Vector result;
+  result.x=x+v.x;
+  result.y=y+v.y;
+  result.z=z+v.z;
+  
+  return result;
+
+}
+
+
+Vector& Vector::operator-=(const Vector &v) {
+  
+  x-=v.x;
+  y-=v.y;
+  z-=v.z;
+  
+  return *this;
+
+}
+
+Vector Vector::operator/(const double divisor) const {
+  
+  
+  Vector result;
+  
+  result.x=x/divisor;
+  result.y=y/divisor;
+  result.z=z/divisor;
+
+  return result;
+  
+  
+}
+
+
+Vector Vector::operator*(const double multiplier) const {
+
+  Vector result;
+  
+  result.x=x*multiplier;
+  result.y=y*multiplier;
+  result.z=z*multiplier;
+
+  return result;
+
+  
+}
+
+
+Vector operator*(const double multiplier, const Vector &v) {
+
+  Vector result;
+  
+  result.x=v.x*multiplier;
+  result.y=v.y*multiplier;
+  result.z=v.z*multiplier;
+
+  return result;
+  
+}
+
+Vector &Vector::operator/=(const double divisor) {
+  x/=divisor;
+  y/=divisor;
+  z/=divisor;
+  
+  return *this;
+
+}
+
+Vector &Vector::operator*=(const double multiplier) {
+
+  x*=multiplier;
+  y*=multiplier;
+  z*=multiplier;
+
+  return *this;
+
+}
+
+Vector Vector::operator*(const Vector &v) const {
+  
+  // cross product ("uitproduct")
+  Vector result;
+  
+  result.x=y*v.z-z*v.y;
+  result.y=z*v.x-x*v.z;
+  result.z=x*v.y-y*v.x;
+  
+  return result;
+  
+  
+}
+
+
+double InnerProduct(const Vector &v1, const Vector &v2) {
+
+  // Inner product ("inproduct")
+  double result;
+  result=v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;
+  return result;
+
+}
+
+double Vector::Angle(const Vector &v) const {
+  
+  // angle between this vector and another vector
+  
+  // angle is within range of [0,pi] radians
+  
+  // angle is arccosine of the inner product over the product of the norms of the vectors
+  
+  double cos_angle=InnerProduct(*this,v)/(Norm()*v.Norm());
+  
+  // check for computational inaccuracies
+  if (cos_angle<=-1)
+    return Pi;
+  
+  if (cos_angle>=1)
+    return 0.;
+  
+  double angle=acos(cos_angle);
+  
+  return angle;
+
+}
+
+double Vector::SignedAngle(const Vector &v) const {
+  
+  // angle between this vector and another vector
+  
+  // angle is within range of [-pi,pi] radians
+  
+  // angle is arccosine of the inner product over the product of the norms of the vectors
+  
+  double cos_angle=InnerProduct(*this,v)/(Norm()*v.Norm());
+  
+  // check for computational inaccuracies
+  if (cos_angle<=-1)
+    return Pi;
+  
+  if (cos_angle>=1)
+    return 0.;
+  
+  double angle=acos(cos_angle);
+  
+  double sign = (InnerProduct ( Perp2D(), v ) )>0.?1.:-1.;
+  return angle * sign;
+
+}
+
+
+bool Vector::operator==(const Vector &v) const {
+   
+  // "sloppy equal"
+  if ((fabs(x-v.x)<TINY) && (fabs(y-v.y)<TINY) && (fabs(z-v.z)<TINY))
+    return true;
+  else
+    return false;
+  
+}
+
+bool Vector::operator< (const Vector &v) const {
+  
+  // Compare x coordinate
+  if (x<v.x) 
+    return true;
+  else 
+    return false;
+  
+}
+
+
+double Vector::SqrNorm(void) const {
+  
+  // return the square of the norm
+  // added this function to avoid taking the square root and 
+  // the square again if this value is needed
+  return DSQR(x)+DSQR(y)+DSQR(z);
+
+}
+
+void Vector::Normalise(void) {
+  
+  double norm;
+  norm=Norm(); // Absolute value;
+  
+  if (norm>0.) { // if the norm is 0, don't normalise 
+    // (otherwise division by zero occurs)
+    
+    (*this)/=norm;
+  }
+  
+}
+
+Vector Vector::Normalised(void) const {
+  
+  double norm;
+  norm=Norm(); // Absolute value;
+  
+  if (norm>0.) { // if the norm is 0, don't normalise 
+    // (otherwise division by zero occurs)
+    return (*this)/norm;
+  } else {
+    return *this;
+  }
+    
+}
+
+bool Vector::SameDirP(const Vector &v) {
+  
+  // return true if the two (parallel) vectors point in the same direction
+  //  if (x*v.x>=0 && y*v.y>=0 && z*v.z>=0)
+  double angle=Angle(v);
+  
+  if (angle<(Pi/2))
+    return true;
+  else
+    return false;
+ 
+}
+
+double Vector::Max(void) {
+  
+  // Find maximum value of vector
+  double max;
+  
+  max=x > y ? x : y;
+  max=max > z ? max : z;
+  
+  return max;
+}
+
+double Vector::Min(void) {
+  
+  // Find minimum value of vector
+  double min;
+  
+  min=x < y ? x : y;
+  min=min < z ? min : z;
+  
+  return min;
+}
+
+// test
+
+#ifdef TEST
+
+void main() {
+  
+  Vector v(1.,2.,3.);
+  
+  Vector d;
+  
+  d=5*v;
+  
+//  cerr << d << "\n";
+  
+  d=v*5;
+  
+//  cerr << d << "\n";
+  
+  d=v/5;
+  
+//  cerr << d << "\n";
+
+//  cerr << d.Norm() << "\n";
+  
+  d.Normalise();
+  
+//  cerr << d << "  " << d.Norm() << "\n";
+  
+  
+}
+
+#endif