diff --git a/src/matrix.cpp b/src/matrix.cpp
new file mode 100644
--- /dev/null
+++ b/src/matrix.cpp
@@ -0,0 +1,178 @@
+/*
+ *
+ *  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 <ostream>
+#include <cmath>
+#include "vector.h"
+#include "matrix.h"
+#include "tiny.h"
+
+static const std::string _module_id("$Id$");
+
+Matrix::Matrix(const Vector &c1, const Vector &c2, const Vector &c3) {
+  
+  Alloc();
+
+  mat[0][0]=c1.x; mat[0][1]=c2.x; mat[0][2]=c3.x;
+  mat[1][0]=c1.y; mat[1][1]=c2.y; mat[1][2]=c3.y;
+  mat[2][0]=c1.z; mat[2][1]=c2.z; mat[2][2]=c3.z;
+
+    
+}
+
+void Matrix::Alloc(void) {
+
+  // constructor
+//  mat=(double **)malloc(3*sizeof(double *));
+//  mat[0]=(double *)malloc(9*sizeof(double));
+  mat = new double*[3];
+  mat[0] = new double[9];
+  for (int i=1;i<3;i++)
+    mat[i]=mat[i-1]+3;
+  
+}
+
+Matrix::~Matrix() {
+ 
+  // destructor
+  //free(mat[0]);
+  //free(mat);
+  delete[] mat[0];
+  delete[] mat;
+}
+
+Matrix::Matrix(void) {
+  
+  // constructor
+  Alloc();
+
+  // clear matrix
+  for (int i=0;i<9;i++) {
+    mat[0][i]=0.;
+  }
+ 
+}
+
+Matrix::Matrix(const Matrix &source) {
+ 
+  // copy constructor
+  Alloc();
+  
+  for (int i=0;i<9;i++) {
+    mat[0][i]=source.mat[0][i];
+  }
+  
+}
+
+
+void Matrix::operator=(const Matrix &source) {
+  
+  // assignment
+  
+  // don't assign to self
+  if (this==&source) return;
+  
+  // copy 
+  for (int i=0;i<9;i++)
+    mat[0][i]=source.mat[0][i];
+  
+}
+
+
+void Matrix::print(ostream *os) {
+  
+  *os << "{ { " << mat[0][0] << "," << mat[0][1] << "," << mat[0][2] << "},{" << mat[1][0] << "," << mat[1][1] << "," << mat[1][2] << "},{" << mat[2][0] << "," << mat[2][1] << "," << mat[2][2] << "} }";
+
+}
+
+ostream &operator<<(ostream &os, Matrix &v) {
+  v.print(&os);
+  return os;
+}
+
+
+Vector Matrix::operator*(const Vector &v) const {
+  
+  // matrix * vector
+  Vector result;
+  
+  result.x = mat[0][0]*v.x+mat[0][1]*v.y+mat[0][2]*v.z;
+  result.y = mat[1][0]*v.x+mat[1][1]*v.y+mat[1][2]*v.z;
+  result.z = mat[2][0]*v.x+mat[2][1]*v.y+mat[2][2]*v.z;
+  
+  return result;
+  
+  
+}
+
+
+bool Matrix::operator==(Matrix &m) const {
+   
+  for (int i=0;i<9;i++) {
+    if ((mat[0][i]-m.mat[0][i])>TINY)
+      return false;
+  }
+  return true;
+  
+}
+
+double Matrix::Det(void) const {
+  
+  return 
+    - mat[0][2]*mat[0][4]*mat[0][6]
+    + mat[0][1]*mat[0][5]*mat[0][6] 
+    + mat[0][2]*mat[0][3]*mat[0][7]
+    - mat[0][0]*mat[0][5]*mat[0][7]
+    - mat[0][1]*mat[0][3]*mat[0][8]
+    + mat[0][0]*mat[0][4]*mat[0][8];
+
+}
+
+Matrix Matrix::Inverse(void) const {
+  
+  // return the Inverse of this matrix
+  double rd=1./Det(); // Reciproce Det;
+  Matrix inverse;
+  inverse.mat[0][0]=rd*(-mat[0][5]*mat[0][7]+mat[0][4]*mat[0][8]);
+  inverse.mat[0][1]=rd*( mat[0][2]*mat[0][7]-mat[0][1]*mat[0][8]);
+  inverse.mat[0][2]=rd*(-mat[0][2]*mat[0][4]+mat[0][1]*mat[0][5]);
+  inverse.mat[0][3]=rd*( mat[0][5]*mat[0][6]-mat[0][3]*mat[0][8]);
+  inverse.mat[0][4]=rd*(-mat[0][2]*mat[0][6]+mat[0][0]*mat[0][8]);
+  inverse.mat[0][5]=rd*( mat[0][2]*mat[0][3]-mat[0][0]*mat[0][5]);
+  inverse.mat[0][6]=rd*(-mat[0][4]*mat[0][6]+mat[0][3]*mat[0][7]);
+  inverse.mat[0][7]=rd*( mat[0][1]*mat[0][6]-mat[0][0]*mat[0][7]);
+  inverse.mat[0][8]=rd*(-mat[0][1]*mat[0][3]+mat[0][0]*mat[0][4]);
+  
+
+  return inverse;
+}
+
+void Matrix::Rot2D(double theta) { 
+  
+  // make this matrix a rotation matrix over theta
+  // see http://mathworld.wolfram.com/RotationMatrix.html
+ 
+  mat[0][0] = cos(theta); mat[0][1]=sin(theta);
+  mat[1][0] = -sin(theta); mat[1][1]=cos(theta);
+  mat[0][2] = mat[1][2] = mat[2][0] = mat[2][1] = mat[2][2] = 0.;
+
+}