Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/Makefile.am	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -9,4 +9,5 @@
 			   gslvector.cpp \
 			   linearsystemofequations.cpp \
+			   SingleVector.cpp \
 			   vector.cpp
 			   
@@ -14,4 +15,5 @@
 			   gslvector.hpp \
 			   linearsystemofequations.hpp \
+			   SingleVector.hpp \
 			   vector.hpp
 			   
@@ -122,5 +124,4 @@
 		 periodentafel.cpp \
 		 Plane.cpp \
-		 SingleVector.cpp \
 		 tesselation.cpp \
 		 tesselationhelpers.cpp \
@@ -158,5 +159,4 @@
          periodentafel.hpp \
          Plane.hpp \
-         SingleVector.hpp \
          stackclass.hpp \
          tesselation.hpp \
Index: src/SingleVector.cpp
===================================================================
--- src/SingleVector.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/SingleVector.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -21,9 +21,14 @@
 #include <gsl/gsl_vector.h>
 
+#include <iostream>
+
+using namespace std;
+
 /************************************ Functions for class vector ************************************/
 
 /** Constructor of class vector.
  */
-SingleVector::SingleVector()
+SingleVector::SingleVector() :
+  Vector(Baseconstructor())
 {
   x[0] = x[1] = x[2] = 0.;
@@ -32,5 +37,6 @@
 /** Constructor of class vector.
  */
-SingleVector::SingleVector(const double x1, const double x2, const double x3)
+SingleVector::SingleVector(const double x1, const double x2, const double x3) :
+  Vector(Baseconstructor())
 {
   x[0] = x1;
@@ -42,5 +48,6 @@
  * Copy constructor
  */
-SingleVector::SingleVector(const Vector& src)
+SingleVector::SingleVector(const Vector& src) :
+  Vector(Baseconstructor())
 {
   x[0] = src[0];
@@ -49,15 +56,6 @@
 }
 
-/**
- * Assignment operator
- */
-Vector& SingleVector::operator=(const Vector& src){
-  // check for self assignment
-  if(&src!=this){
-    x[0] = src[0];
-    x[1] = src[1];
-    x[2] = src[2];
-  }
-  return *this;
+SingleVector* SingleVector::clone() const{
+  return new SingleVector(x[0],x[1],x[2]);
 }
 
@@ -76,16 +74,4 @@
     res += (x[i]-y[i])*(x[i]-y[i]);
   return (res);
-};
-
-/** Calculates distance between this and another vector.
- * \param *y array to second vector
- * \return \f$| x - y |\f$
- */
-double SingleVector::Distance(const Vector &y) const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += (x[i]-y[i])*(x[i]-y[i]);
-  return (sqrt(res));
 };
 
@@ -173,10 +159,9 @@
 //  bool flag = false;
   //vector Shifted, TranslationVector;
-  Vector TestVector;
 //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
 //  Log() << Verbose(2) << "Vector is: ";
 //  Output(out);
 //  Log() << Verbose(0) << endl;
-  TestVector = *this;
+  SingleVector TestVector(*this);
   TestVector.InverseMatrixMultiplication(matrix);
   for(int i=NDIM;i--;) { // correct periodically
@@ -208,4 +193,14 @@
 
 
+void SingleVector::AddVector(const Vector &y) {
+  for(int i=NDIM;i--;)
+    x[i] += y[i];
+}
+
+void SingleVector::SubtractVector(const Vector &y){
+  for(int i=NDIM;i--;)
+    x[i] -= y[i];
+}
+
 /** Calculates VectorProduct between this and another vector.
  *  -# returns the Product in place of vector from which it was initiated
@@ -216,9 +211,9 @@
 void SingleVector::VectorProduct(const Vector &y)
 {
-  Vector tmp;
+  SingleVector tmp;
   tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
   tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
   tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
-  (*this) = tmp;
+  CopyVector(tmp);
 };
 
@@ -245,12 +240,10 @@
 double SingleVector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
 {
-  Vector temp;
-
   // first create part that is orthonormal to PlaneNormal with withdraw
-  temp = (*this )- PlaneOffset;
+  Vector temp = VecFromRep(this)- PlaneOffset;
   temp.MakeNormalTo(PlaneNormal);
   temp.Scale(-1.);
   // then add connecting vector from plane to point
-  temp += (*this)-PlaneOffset;
+  temp += VecFromRep(this)-PlaneOffset;
   double sign = temp.ScalarProduct(PlaneNormal);
   if (fabs(sign) > MYEPSILON)
@@ -282,60 +275,4 @@
 
   return helper;
-};
-
-/** Calculates norm of this vector.
- * \return \f$|x|\f$
- */
-double SingleVector::Norm() const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += this->x[i]*this->x[i];
-  return (sqrt(res));
-};
-
-/** Calculates squared norm of this vector.
- * \return \f$|x|^2\f$
- */
-double SingleVector::NormSquared() const
-{
-  return (ScalarProduct(*this));
-};
-
-/** Normalizes this vector.
- */
-void SingleVector::Normalize()
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += this->x[i]*this->x[i];
-  if (fabs(res) > MYEPSILON)
-    res = 1./sqrt(res);
-  Scale(&res);
-};
-
-/** Zeros all components of this vector.
- */
-void SingleVector::Zero()
-{
-  for (int i=NDIM;i--;)
-    this->x[i] = 0.;
-};
-
-/** Zeros all components of this vector.
- */
-void SingleVector::One(const double one)
-{
-  for (int i=NDIM;i--;)
-    this->x[i] = one;
-};
-
-/** Initialises all components of this vector.
- */
-void SingleVector::Init(const double x1, const double x2, const double x3)
-{
-  x[0] = x1;
-  x[1] = x2;
-  x[2] = x3;
 };
 
@@ -410,12 +347,4 @@
 }
 
-double& SingleVector::at(size_t i){
-  return (*this)[i];
-}
-
-const double& SingleVector::at(size_t i) const{
-  return (*this)[i];
-}
-
 double* SingleVector::get(){
   return x;
@@ -424,14 +353,8 @@
  * \param *factor pointer to scaling factor
  */
-void SingleVector::Scale(const double ** const factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= (*factor)[i];
-};
-
-void SingleVector::Scale(const double * const factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= *factor;
+void SingleVector::ScaleAll(const double *factor)
+{
+  for (int i=NDIM;i--;)
+    x[i] *= factor[i];
 };
 
@@ -440,13 +363,4 @@
   for (int i=NDIM;i--;)
     x[i] *= factor;
-};
-
-/** Translate atom by given vector.
- * \param trans[] translation vector.
- */
-void SingleVector::Translate(const Vector &trans)
-{
-  for (int i=NDIM;i--;)
-    x[i] += trans[i];
 };
 
@@ -474,5 +388,5 @@
 void SingleVector::MatrixMultiplication(const double * const M)
 {
-  Vector C;
+  SingleVector C;
   // do the matrix multiplication
   C[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
@@ -480,5 +394,5 @@
   C[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
 
-  *this = C;
+  CopyVector(C);
 };
 
@@ -488,5 +402,5 @@
 bool SingleVector::InverseMatrixMultiplication(const double * const A)
 {
-  Vector C;
+  SingleVector C;
   double B[NDIM*NDIM];
   double detA = RDET3(A);
@@ -511,5 +425,5 @@
     C[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
 
-    *this = C;
+    CopyVector(C);
     return true;
   } else {
@@ -518,17 +432,4 @@
 };
 
-
-/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
- * this vector = x1*factors[0] + x2* factors[1] + x3*factors[2]
- * \param *x1 first vector
- * \param *x2 second vector
- * \param *x3 third vector
- * \param *factors three-component vector with the factor for each given vector
- */
-void SingleVector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
-{
-  for(int i=NDIM;i--;)
-    x[i] = factors[0]*x1[i] + factors[1]*x2[i] + factors[2]*x3[i];
-};
 
 /** Mirrors atom against a given plane.
@@ -615,5 +516,6 @@
 bool SingleVector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  Vector a = (*this) - offset;
+  Vector a = VecFromRep(this);
+  a -= offset;
   a.InverseMatrixMultiplication(parallelepiped);
   bool isInside = true;
@@ -624,2 +526,12 @@
   return isInside;
 }
+
+
+void SingleVector::CopyVector(SingleVector& vec) {
+  for(int i =NDIM; i--;)
+    x[i]=vec.x[i];
+}
+
+bool SingleVector::isBaseClass() const{
+  return false;
+}
Index: src/SingleVector.hpp
===================================================================
--- src/SingleVector.hpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/SingleVector.hpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -31,5 +31,4 @@
   virtual ~SingleVector();
 
-  virtual double Distance(const Vector &y) const;
   virtual double DistanceSquared(const Vector &y) const;
   virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
@@ -37,6 +36,4 @@
   virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
   virtual double ScalarProduct(const Vector &y) const;
-  virtual double Norm() const;
-  virtual double NormSquared() const;
   virtual double Angle(const Vector &y) const;
   virtual bool IsZero() const;
@@ -45,24 +42,18 @@
   virtual bool IsEqualTo(const Vector &a) const;
 
+  virtual void AddVector(const Vector &y);
+  virtual void SubtractVector(const Vector &y);
   virtual void VectorProduct(const Vector &y);
   virtual void ProjectOntoPlane(const Vector &y);
   virtual void ProjectIt(const Vector &y);
   virtual Vector Projection(const Vector &y) const;
-  virtual void Zero();
-  virtual void One(const double one);
-  virtual void Init(const double x1, const double x2, const double x3);
-  virtual void Normalize();
-  virtual void Translate(const Vector &x);
   virtual void Mirror(const Vector &x);
-  virtual void Scale(const double ** const factor);
-  virtual void Scale(const double * const factor);
+  virtual void ScaleAll(const double *factor);
   virtual void Scale(const double factor);
   virtual void MatrixMultiplication(const double * const M);
   virtual bool InverseMatrixMultiplication(const double * const M);
   virtual void KeepPeriodic(const double * const matrix);
-  virtual void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
   virtual bool GetOneNormalVector(const Vector &x1);
   virtual bool MakeNormalTo(const Vector &y1);
-  //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
   virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   virtual void WrapPeriodically(const double * const M, const double * const Minv);
@@ -71,14 +62,15 @@
   virtual double& operator[](size_t i);
   virtual const double& operator[](size_t i) const;
-  virtual double& at(size_t i);
-  virtual const double& at(size_t i) const;
-
-  // Assignment operator
-  virtual Vector &operator=(const Vector& src);
 
   // Access to internal structure
   virtual double* get();
+protected:
 
 private:
+  // method used for protection, i.e. to avoid infinite recursion
+  // when our internal rep becomes messed up
+  virtual bool isBaseClass() const;
+  virtual SingleVector *clone() const;
+  void CopyVector(SingleVector &rhs);
   double x[NDIM];
 
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/boundary.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -827,5 +827,5 @@
 
   // calculate filler grid in [0,1]^3
-  FillerDistance.Init(distance[0], distance[1], distance[2]);
+  FillerDistance = Vector(distance[0], distance[1], distance[2]);
   FillerDistance.InverseMatrixMultiplication(M);
   for(int i=0;i<NDIM;i++)
@@ -841,5 +841,5 @@
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
-        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         CurrentPosition.MatrixMultiplication(M);
         Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/builder.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -1644,13 +1644,13 @@
                 first = World::getInstance().createAtom();
                 first->type = periode->FindElement(1);
-                first->x.Init(0.441, -0.143, 0.);
+                first->x = Vector(0.441, -0.143, 0.);
                 filler->AddAtom(first);
                 second = World::getInstance().createAtom();
                 second->type = periode->FindElement(1);
-                second->x.Init(-0.464, 1.137, 0.0);
+                second->x = Vector(-0.464, 1.137, 0.0);
                 filler->AddAtom(second);
                 third = World::getInstance().createAtom();
                 third->type = periode->FindElement(8);
-                third->x.Init(-0.464, 0.177, 0.);
+                third->x = Vector(-0.464, 0.177, 0.);
                 filler->AddAtom(third);
                 filler->AddBond(first, third, 1);
Index: src/ellipsoid.cpp
===================================================================
--- src/ellipsoid.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/ellipsoid.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -59,5 +59,5 @@
   Matrix[8] = cos(theta);
   helper.MatrixMultiplication(Matrix);
-  helper.Scale(InverseLength);
+  helper.ScaleAll(InverseLength);
   //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
 
@@ -70,5 +70,5 @@
   theta = -EllipsoidAngle[1];
   phi = -EllipsoidAngle[2];
-  helper.Scale(EllipsoidLength);
+  helper.ScaleAll(EllipsoidLength);
   Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
   Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/molecule.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -239,5 +239,5 @@
     matrix = ReturnFullMatrixforSymmetric(cell_size);
     Orthovector1.MatrixMultiplication(matrix);
-    InBondvector.SubtractVector(Orthovector1); // subtract just the additional translation
+    InBondvector -= Orthovector1; // subtract just the additional translation
     Free(&matrix);
     bondlength = InBondvector.Norm();
@@ -272,5 +272,5 @@
         FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
       }
-      InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
+      InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
       FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
       FirstOtherAtom->x = InBondvector;  // ... and add distance vector to replacement atom
@@ -356,6 +356,6 @@
         SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
       }
-      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
-      SecondOtherAtom->x.Scale(&BondRescale);
+      FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
+      SecondOtherAtom->x *= BondRescale;
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
       for(int i=NDIM;i--;) { // and make relative to origin atom
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/molecule_geometry.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -191,4 +191,9 @@
 /** Scales all atoms by \a *factor.
  * \param *factor pointer to scaling factor
+ *
+ * TODO: Is this realy what is meant, i.e.
+ * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
+ * or rather
+ * x=(**factor) * x (as suggested by comment)
  */
 void molecule::Scale(const double ** const factor)
@@ -199,6 +204,6 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Scale(factor);
-    ptr->x.Scale(factor);
+      ptr->Trajectory.R.at(j).ScaleAll(*factor);
+    ptr->x.ScaleAll(*factor);
   }
 };
@@ -214,6 +219,6 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Translate(*trans);
-    ptr->x.Translate(*trans);
+      ptr->Trajectory.R.at(j) += (*trans);
+    ptr->x += (*trans);
   }
 };
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/tesselationhelpers.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -427,10 +427,10 @@
   double volume;
 
-  TetraederVector[0].CopyVector(a);
-  TetraederVector[1].CopyVector(b);
-  TetraederVector[2].CopyVector(c);
+  TetraederVector[0] = a;
+  TetraederVector[1] = b;
+  TetraederVector[2] = c;
   for (int j=0;j<3;j++)
     TetraederVector[j].SubtractVector(d);
-  Point.CopyVector(TetraederVector[0]);
+  Point = TetraederVector[0];
   Point.VectorProduct(TetraederVector[1]);
   volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
Index: src/unittests/ActOnAllUnitTest.cpp
===================================================================
--- src/unittests/ActOnAllUnitTest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/ActOnAllUnitTest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -70,15 +70,8 @@
 
   // scaling by value
-  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 2. );
+  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, factor );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
 
-  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, 0.5 );
-  CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
-
-  // scaling by ref
-  VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&factor );
-  CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
-
-  VL.ActOnAll( (void (Vector::*)(const double * const)) &Vector::Scale, (const double * const)&inverse );
+  VL.ActOnAll( (void (Vector::*)(const double)) &Vector::Scale, inverse );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
 
@@ -90,8 +83,8 @@
     inverses[i] = 1./factors[i];
   }
-  VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&factors );
+  VL.ActOnAll<Vector,void,const double*>(&Vector::ScaleAll, factors );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
 
-  VL.ActOnAll( (void (Vector::*)(const double ** const)) &Vector::Scale, (const double ** const)&inverses );
+  VL.ActOnAll<Vector,void,const double*>(&Vector::ScaleAll, inverses );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
   Free(factors);
Index: src/unittests/AnalysisCorrelationToPointUnitTest.cpp
===================================================================
--- src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -64,17 +64,17 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
+  *Walker->node = Vector(1., 0., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
+  *Walker->node = Vector(0., 1., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
+  *Walker->node = Vector(1., 1., 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
===================================================================
--- src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -40,5 +40,5 @@
 void AnalysisCorrelationToSurfaceUnitTest::setUp()
 {
-  ASSERT_DO(Assert::Throw);
+  //ASSERT_DO(Assert::Throw);
 
   atom *Walker = NULL;
@@ -73,17 +73,17 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(1., 0., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(0., 1., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(1., 1., 0. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
 
@@ -106,18 +106,18 @@
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
-  Walker->node->Init(4., 0., 4. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = carbon;
-  Walker->node->Init(0., 4., 4. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = carbon;
-  Walker->node->Init(4., 4., 0. );
+  *Walker->node = Vector(4., 0., 4. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = carbon;
+  *Walker->node = Vector(0., 4., 4. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = carbon;
+  *Walker->node = Vector(4., 4., 0. );
   TestMolecule->AddAtom(Walker);
   // add inner atoms
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
-  Walker->node->Init(0.5, 0.5, 0.5 );
+  *Walker->node = Vector(0.5, 0.5, 0.5 );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/AnalysisPairCorrelationUnitTest.cpp
===================================================================
--- src/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -63,17 +63,17 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
+  *Walker->node = Vector(1., 0., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
+  *Walker->node = Vector(0., 1., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
+  *Walker->node = Vector(1., 1., 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/analysisbondsunittest.cpp
===================================================================
--- src/unittests/analysisbondsunittest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/analysisbondsunittest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -69,21 +69,21 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1.5, 0., 1.5 );
+  *Walker->node = Vector(1.5, 0., 1.5 );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 1.5, 1.5 );
+  *Walker->node = Vector(0., 1.5, 1.5 );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1.5, 1.5, 0. );
+  *Walker->node = Vector(1.5, 1.5, 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = carbon;
-  Walker->node->Init(0.5, 0.5, 0.5 );
+  *Walker->node = Vector(0.5, 0.5, 0.5 );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/bondgraphunittest.cpp
===================================================================
--- src/unittests/bondgraphunittest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/bondgraphunittest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -65,17 +65,17 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
+  *Walker->node = Vector(1., 0., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
+  *Walker->node = Vector(0., 1., 1. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
+  *Walker->node = Vector(1., 1., 0. );
   TestMolecule->AddAtom(Walker);
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/listofbondsunittest.cpp
===================================================================
--- src/unittests/listofbondsunittest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/listofbondsunittest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -58,17 +58,17 @@
   Walker = World::getInstance().createAtom();
   Walker->type = hydrogen;
-  Walker->node->Init(1., 0., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 1., 1. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(1., 1., 0. );
-  TestMolecule->AddAtom(Walker);
-  Walker = World::getInstance().createAtom();
-  Walker->type = hydrogen;
-  Walker->node->Init(0., 0., 0. );
+  *Walker->node = Vector(1., 0., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(0., 1., 1. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(1., 1., 0. );
+  TestMolecule->AddAtom(Walker);
+  Walker = World::getInstance().createAtom();
+  Walker->type = hydrogen;
+  *Walker->node = Vector(0., 0., 0. );
   TestMolecule->AddAtom(Walker);
 
Index: src/unittests/tesselation_boundarytriangleunittest.cpp
===================================================================
--- src/unittests/tesselation_boundarytriangleunittest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/tesselation_boundarytriangleunittest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -86,49 +86,49 @@
 
   // simple test on y line
-  Point.Init(-1.,0.5,0.);
-  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.5,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(-4.,0.5,0.);
+  Point = Vector(-1.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point = Vector(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(-4.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.5,0.);
+  Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on x line
-  Point.Init(0.5,-1.,0.);
-  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(0.5,-6.,0.);
+  Point = Vector(0.5,-1.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point = Vector(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(0.5,-6.,0.);
   CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.,0.);
+  Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on slanted line
-  Point.Init(1.,1.,0.);
+  Point = Vector(1.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.5,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(5.,5.,0.);
+  Point = Vector(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(5.,5.,0.);
   CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.5,0.);
+  Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on first node
-  Point.Init(-1.,-1.,0.);
+  Point = Vector(-1.,-1.,0.);
   CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.,0.);
+  Point = Vector(0.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on second node
-  Point.Init(0.,2.,0.);
-  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,1.,0.);
+  Point = Vector(0.,2.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point = Vector(0.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on third node
-  Point.Init(2.,0.,0.);
-  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(1.,0.,0.);
+  Point = Vector(2.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point = Vector(1.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 };
@@ -142,59 +142,59 @@
 
   // straight down/up
-  Point.Init(1./3.,1./3.,+5.);
+  Point = Vector(1./3.,1./3.,+5.);
   CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(1./3.,1./3.,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(1./3.,1./3.,-5.);
+  Point = Vector(1./3.,1./3.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(1./3.,1./3.,-5.);
   CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(1./3.,1./3.,0.);
+  Point = Vector(1./3.,1./3.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on y line
-  Point.Init(-1.,0.5,+2.);
+  Point = Vector(-1.,0.5,+2.);
   CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.5,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(-1.,0.5,-3.);
+  Point = Vector(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(-1.,0.5,-3.);
   CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.5,0.);
+  Point = Vector(0.,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on x line
-  Point.Init(0.5,-1.,+1.);
+  Point = Vector(0.5,-1.,+1.);
   CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(0.5,-1.,-2.);
+  Point = Vector(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(0.5,-1.,-2.);
   CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.,0.);
+  Point = Vector(0.5,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on slanted line
-  Point.Init(1.,1.,+3.);
+  Point = Vector(1.,1.,+3.);
   CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.5,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-  Point.Init(1.,1.,-4.);
+  Point = Vector(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point = Vector(1.,1.,-4.);
   CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.5,0.5,0.);
+  Point = Vector(0.5,0.5,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on first node
-  Point.Init(-1.,-1.,5.);
+  Point = Vector(-1.,-1.,5.);
   CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,0.,0.);
+  Point = Vector(0.,0.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on second node
-  Point.Init(0.,2.,5.);
+  Point = Vector(0.,2.,5.);
   CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(0.,1.,0.);
+  Point = Vector(0.,1.,0.);
   CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
 
   // simple test on third node
-  Point.Init(2.,0.,5.);
+  Point = Vector(2.,0.,5.);
   CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
-  Point.Init(1.,0.,0.);
-  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
-};
+  Point = Vector(1.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+};
Index: src/unittests/vectorunittest.cpp
===================================================================
--- src/unittests/vectorunittest.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/unittests/vectorunittest.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -26,4 +26,8 @@
 #endif /*HAVE_TESTRUNNER*/
 
+#include <iostream>
+
+using namespace std;
+
 /********************************************** Test classes **************************************/
 
@@ -34,9 +38,9 @@
 void VectorTest::setUp()
 {
-  zero.Init(0.,0.,0.);
-  unit.Init(1.,0.,0.);
-  otherunit.Init(0.,1.,0.);
-  notunit.Init(0.,1.,1.);
-  two.Init(2.,1.,0.);
+  zero  = Vector(0.,0.,0.);
+  unit = Vector(1.,0.,0.);
+  otherunit = Vector(0.,1.,0.);
+  notunit = Vector(0.,1.,1.);
+  two = Vector(2.,1.,0.);
 };
 
@@ -69,5 +73,5 @@
   double factor;
   // copy vector
-  fixture.Init(2.,3.,4.);
+  fixture = Vector(2.,3.,4.);
   CPPUNIT_ASSERT_EQUAL( Vector(2.,3.,4.), fixture );
   // summation and scaling
@@ -223,5 +227,5 @@
 void VectorTest::VectorRotationTest()
 {
-  fixture.Init(-1.,0.,0.);
+  fixture = Vector(-1.,0.,0.);
 
   // zero vector does not change
@@ -266,18 +270,18 @@
   parallelepiped[8] = 1;
 
-  fixture.CopyVector(zero);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(2.5,2.5,2.5);
+  fixture = zero;
+  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
+  fixture = Vector(2.5,2.5,2.5);
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(1.,1.,1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(3.5,3.5,3.5);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(2.,2.,2.);
+  fixture = Vector(1.,1.,1.);
+  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
+  fixture = Vector(3.5,3.5,3.5);
+  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
+  fixture = Vector(2.,2.,2.);
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(2.,3.,2.);
+  fixture = Vector(2.,3.,2.);
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture.Init(-2.,2.,-1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-}
-
+  fixture = Vector(-2.,2.,-1.);
+  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
+}
+
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/vector.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -6,18 +6,11 @@
 
 
-#include "defs.hpp"
-#include "gslmatrix.hpp"
-#include "leastsquaremin.hpp"
-#include "memoryallocator.hpp"
-#include "vector.hpp"
-#include "Helpers/fast_functions.hpp"
+#include "SingleVector.hpp"
 #include "Helpers/Assert.hpp"
-#include "Plane.hpp"
-#include "Exceptions/LinearDependenceException.hpp"
-
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_permutation.h>
-#include <gsl/gsl_vector.h>
+
+#include <iostream>
+
+using namespace std;
+
 
 /************************************ Functions for class vector ************************************/
@@ -25,27 +18,31 @@
 /** Constructor of class vector.
  */
-Vector::Vector()
-{
-  x[0] = x[1] = x[2] = 0.;
-};
+Vector::Vector() :
+  rep(new SingleVector())
+{};
+
+Vector::Vector(Baseconstructor)  // used by derived objects to construct their bases
+{}
+
+Vector::Vector(Baseconstructor,const Vector* v) :
+  rep(v->clone())
+{}
+
+Vector Vector::VecFromRep(const Vector* v){
+  return Vector(Baseconstructor(),v);
+}
 
 /** Constructor of class vector.
  */
-Vector::Vector(const double x1, const double x2, const double x3)
-{
-  x[0] = x1;
-  x[1] = x2;
-  x[2] = x3;
-};
+Vector::Vector(const double x1, const double x2, const double x3) :
+  rep(new SingleVector(x1,x2,x3))
+{};
 
 /**
  * Copy constructor
  */
-Vector::Vector(const Vector& src)
-{
-  x[0] = src[0];
-  x[1] = src[1];
-  x[2] = src[2];
-}
+Vector::Vector(const Vector& src) :
+  rep(src.rep->clone())
+{}
 
 /**
@@ -53,9 +50,8 @@
  */
 Vector& Vector::operator=(const Vector& src){
+  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   // check for self assignment
   if(&src!=this){
-    x[0] = src[0];
-    x[1] = src[1];
-    x[2] = src[2];
+    rep.reset(src.rep->clone());
   }
   return *this;
@@ -72,8 +68,6 @@
 double Vector::DistanceSquared(const Vector &y) const
 {
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += (x[i]-y[i])*(x[i]-y[i]);
-  return (res);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->DistanceSquared(y);
 };
 
@@ -94,32 +88,6 @@
 double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
 {
-  double res = Distance(y), tmp, matrix[NDIM*NDIM];
-  Vector Shiftedy, TranslationVector;
-  int N[NDIM];
-  matrix[0] = cell_size[0];
-  matrix[1] = cell_size[1];
-  matrix[2] = cell_size[3];
-  matrix[3] = cell_size[1];
-  matrix[4] = cell_size[2];
-  matrix[5] = cell_size[4];
-  matrix[6] = cell_size[3];
-  matrix[7] = cell_size[4];
-  matrix[8] = cell_size[5];
-  // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-  for (N[0]=-1;N[0]<=1;N[0]++)
-    for (N[1]=-1;N[1]<=1;N[1]++)
-      for (N[2]=-1;N[2]<=1;N[2]++) {
-        // create the translation vector
-        TranslationVector.Zero();
-        for (int i=NDIM;i--;)
-          TranslationVector.x[i] = (double)N[i];
-        TranslationVector.MatrixMultiplication(matrix);
-        // add onto the original vector to compare with
-        Shiftedy = y + TranslationVector;
-        // get distance and compare with minimum so far
-        tmp = Distance(Shiftedy);
-        if (tmp < res) res = tmp;
-      }
-  return (res);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->PeriodicDistance(y,cell_size);
 };
 
@@ -131,32 +99,6 @@
 double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
 {
-  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
-  Vector Shiftedy, TranslationVector;
-  int N[NDIM];
-  matrix[0] = cell_size[0];
-  matrix[1] = cell_size[1];
-  matrix[2] = cell_size[3];
-  matrix[3] = cell_size[1];
-  matrix[4] = cell_size[2];
-  matrix[5] = cell_size[4];
-  matrix[6] = cell_size[3];
-  matrix[7] = cell_size[4];
-  matrix[8] = cell_size[5];
-  // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-  for (N[0]=-1;N[0]<=1;N[0]++)
-    for (N[1]=-1;N[1]<=1;N[1]++)
-      for (N[2]=-1;N[2]<=1;N[2]++) {
-        // create the translation vector
-        TranslationVector.Zero();
-        for (int i=NDIM;i--;)
-          TranslationVector.x[i] = (double)N[i];
-        TranslationVector.MatrixMultiplication(matrix);
-        // add onto the original vector to compare with
-        Shiftedy = y + TranslationVector;
-        // get distance and compare with minimum so far
-        tmp = DistanceSquared(Shiftedy);
-        if (tmp < res) res = tmp;
-      }
-  return (res);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->PeriodicDistanceSquared(y,cell_size);
 };
 
@@ -167,27 +109,6 @@
 void Vector::KeepPeriodic(const double * const matrix)
 {
-//  int N[NDIM];
-//  bool flag = false;
-  //vector Shifted, TranslationVector;
-  Vector TestVector;
-//  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
-//  Log() << Verbose(2) << "Vector is: ";
-//  Output(out);
-//  Log() << Verbose(0) << endl;
-  TestVector = (*this);
-  TestVector.InverseMatrixMultiplication(matrix);
-  for(int i=NDIM;i--;) { // correct periodically
-    if (TestVector.x[i] < 0) {  // get every coefficient into the interval [0,1)
-      TestVector.x[i] += ceil(TestVector.x[i]);
-    } else {
-      TestVector.x[i] -= floor(TestVector.x[i]);
-    }
-  }
-  TestVector.MatrixMultiplication(matrix);
-  (*this) = TestVector;
-//  Log() << Verbose(2) << "New corrected vector is: ";
-//  Output(out);
-//  Log() << Verbose(0) << endl;
-//  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->KeepPeriodic(matrix);
 };
 
@@ -198,8 +119,6 @@
 double Vector::ScalarProduct(const Vector &y) const
 {
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += x[i]*y[i];
-  return (res);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->ScalarProduct(y);
 };
 
@@ -213,9 +132,6 @@
 void Vector::VectorProduct(const Vector &y)
 {
-  Vector tmp;
-  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
-  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
-  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
-  (*this) = tmp;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->VectorProduct(y);
 };
 
@@ -227,8 +143,6 @@
 void Vector::ProjectOntoPlane(const Vector &y)
 {
-  Vector tmp = y;
-  tmp.Normalize();
-  tmp *= ScalarProduct(tmp);
-  (*this) -= tmp;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->ProjectOntoPlane(y);
 };
 
@@ -241,18 +155,6 @@
 double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
 {
-  // first create part that is orthonormal to PlaneNormal with withdraw
-  Vector temp = (*this) - PlaneOffset;
-  temp.MakeNormalTo(PlaneNormal);
-  temp *= -1.;
-  // then add connecting vector from plane to point
-  temp += (*this);
-  temp -= PlaneOffset;
-  double sign = temp.ScalarProduct(PlaneNormal);
-  if (fabs(sign) > MYEPSILON)
-    sign /= fabs(sign);
-  else
-    sign = 0.;
-
-  return (temp.Norm()*sign);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->DistanceToPlane(PlaneNormal,PlaneOffset);
 };
 
@@ -262,7 +164,6 @@
 void Vector::ProjectIt(const Vector &y)
 {
-  Vector helper = y;
-  helper.Scale(-(ScalarProduct(y)));
-  AddVector(helper);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->ProjectIt(y);
 };
 
@@ -273,8 +174,6 @@
 Vector Vector::Projection(const Vector &y) const
 {
-  Vector helper = y;
-  helper.Scale((ScalarProduct(y)/y.NormSquared()));
-
-  return helper;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->Projection(y);
 };
 
@@ -299,10 +198,6 @@
 void Vector::Normalize()
 {
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += this->x[i]*this->x[i];
-  if (fabs(res) > MYEPSILON)
-    res = 1./sqrt(res);
-  Scale(&res);
+  double factor = Norm();
+  (*this) *= 1/factor;
 };
 
@@ -311,6 +206,5 @@
 void Vector::Zero()
 {
-  for (int i=NDIM;i--;)
-    this->x[i] = 0.;
+  rep.reset(new SingleVector());
 };
 
@@ -319,15 +213,5 @@
 void Vector::One(const double one)
 {
-  for (int i=NDIM;i--;)
-    this->x[i] = one;
-};
-
-/** Initialises all components of this vector.
- */
-void Vector::Init(const double x1, const double x2, const double x3)
-{
-  x[0] = x1;
-  x[1] = x2;
-  x[2] = x3;
+  rep.reset(new SingleVector(one,one,one));
 };
 
@@ -337,5 +221,6 @@
 bool Vector::IsZero() const
 {
-  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->IsZero();
 };
 
@@ -345,5 +230,6 @@
 bool Vector::IsOne() const
 {
-  return (fabs(Norm() - 1.) < MYEPSILON);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->IsOne();
 };
 
@@ -353,8 +239,6 @@
 bool Vector::IsNormalTo(const Vector &normal) const
 {
-  if (ScalarProduct(normal) < MYEPSILON)
-    return true;
-  else
-    return false;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->IsNormalTo(normal);
 };
 
@@ -364,10 +248,6 @@
 bool Vector::IsEqualTo(const Vector &a) const
 {
-  bool status = true;
-  for (int i=0;i<NDIM;i++) {
-    if (fabs(x[i] - a[i]) > MYEPSILON)
-      status = false;
-  }
-  return status;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->IsEqualTo(a);
 };
 
@@ -378,26 +258,17 @@
 double Vector::Angle(const Vector &y) const
 {
-  double norm1 = Norm(), norm2 = y.Norm();
-  double angle = -1;
-  if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
-    angle = this->ScalarProduct(y)/norm1/norm2;
-  // -1-MYEPSILON occured due to numerical imprecision, catch ...
-  //Log() << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
-  if (angle < -1)
-    angle = -1;
-  if (angle > 1)
-    angle = 1;
-  return acos(angle);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->Angle(y);
 };
 
 
 double& Vector::operator[](size_t i){
-  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return (*rep)[i];
 }
 
 const double& Vector::operator[](size_t i) const{
-  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return (*rep)[i];
 }
 
@@ -411,5 +282,6 @@
 
 double* Vector::get(){
-  return x;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->get();
 }
 
@@ -421,8 +293,6 @@
 bool Vector::operator==(const Vector& b) const
 {
-  bool status = true;
-  for (int i=0;i<NDIM;i++)
-    status = status && (fabs((*this)[i] - b[i]) < MYEPSILON);
-  return status;
+  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
+  return IsEqualTo(b);
 };
 
@@ -467,4 +337,5 @@
 Vector const Vector::operator+(const Vector& b) const
 {
+  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.AddVector(b);
@@ -479,4 +350,5 @@
 Vector const Vector::operator-(const Vector& b) const
 {
+  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.SubtractVector(b);
@@ -520,32 +392,17 @@
 };
 
-/** Scales each atom coordinate by an individual \a factor.
- * \param *factor pointer to scaling factor
- */
-void Vector::Scale(const double ** const factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= (*factor)[i];
-};
-
-void Vector::Scale(const double * const factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= *factor;
-};
+
+void Vector::ScaleAll(const double *factor)
+{
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->ScaleAll(factor);
+};
+
+
 
 void Vector::Scale(const double factor)
 {
-  for (int i=NDIM;i--;)
-    x[i] *= factor;
-};
-
-/** Translate atom by given vector.
- * \param trans[] translation vector.
- */
-void Vector::Translate(const Vector &trans)
-{
-  for (int i=NDIM;i--;)
-    x[i] += trans[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->Scale(factor);
 };
 
@@ -556,14 +413,6 @@
 void Vector::WrapPeriodically(const double * const M, const double * const Minv)
 {
-  MatrixMultiplication(Minv);
-  // truncate to [0,1] for each axis
-  for (int i=0;i<NDIM;i++) {
-    x[i] += 0.5;  // set to center of box
-    while (x[i] >= 1.)
-      x[i] -= 1.;
-    while (x[i] < 0.)
-      x[i] += 1.;
-  }
-  MatrixMultiplication(M);
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->WrapPeriodically(M,Minv);
 };
 
@@ -573,12 +422,6 @@
 void Vector::MatrixMultiplication(const double * const M)
 {
-  Vector C;
-  // do the matrix multiplication
-  C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
-  C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
-  C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
-  // transfer the result into this
-  for (int i=NDIM;i--;)
-    x[i] = C.x[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->MatrixMultiplication(M);
 };
 
@@ -588,33 +431,6 @@
 bool Vector::InverseMatrixMultiplication(const double * const A)
 {
-  Vector C;
-  double B[NDIM*NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-
-    // do the matrix multiplication
-    C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
-    C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
-    C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
-    // transfer the result into this
-    for (int i=NDIM;i--;)
-      x[i] = C.x[i];
-    return true;
-  } else {
-    return false;
-  }
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->InverseMatrixMultiplication(A);
 };
 
@@ -639,9 +455,6 @@
 void Vector::Mirror(const Vector &n)
 {
-  double projection;
-  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
-  // withdraw projected vector twice from original one
-  for (int i=NDIM;i--;)
-    x[i] -= 2.*projection*n[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->Mirror(n);
 };
 
@@ -655,12 +468,6 @@
 bool Vector::MakeNormalTo(const Vector &y1)
 {
-  bool result = false;
-  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
-  Vector x1 = factor * y1 ;
-  SubtractVector(x1);
-  for (int i=NDIM;i--;)
-    result = result || (fabs(x[i]) > MYEPSILON);
-
-  return result;
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->MakeNormalTo(y1);
 };
 
@@ -673,36 +480,6 @@
 bool Vector::GetOneNormalVector(const Vector &GivenVector)
 {
-  int Components[NDIM]; // contains indices of non-zero components
-  int Last = 0;   // count the number of non-zero entries in vector
-  int j;  // loop variables
-  double norm;
-
-  for (j=NDIM;j--;)
-    Components[j] = -1;
-  // find two components != 0
-  for (j=0;j<NDIM;j++)
-    if (fabs(GivenVector[j]) > MYEPSILON)
-      Components[Last++] = j;
-
-  switch(Last) {
-    case 3:  // threecomponent system
-    case 2:  // two component system
-      norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
-      x[Components[2]] = 0.;
-      // in skp both remaining parts shall become zero but with opposite sign and third is zero
-      x[Components[1]] = -1./GivenVector[Components[1]] / norm;
-      x[Components[0]] = 1./GivenVector[Components[0]] / norm;
-      return true;
-      break;
-    case 1: // one component system
-      // set sole non-zero component to 0, and one of the other zero component pendants to 1
-      x[(Components[0]+2)%NDIM] = 0.;
-      x[(Components[0]+1)%NDIM] = 1.;
-      x[Components[0]] = 0.;
-      return true;
-      break;
-    default:
-      return false;
-  }
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->GetOneNormalVector(GivenVector);
 };
 
@@ -712,6 +489,6 @@
 void Vector::AddVector(const Vector &y)
 {
-  for (int i=NDIM;i--;)
-    this->x[i] += y[i];
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->AddVector(y);
 }
 
@@ -721,195 +498,7 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  for (int i=NDIM;i--;)
-    this->x[i] -= y[i];
-}
-
-/** Copy vector \a y componentwise.
- * \param y vector
- */
-void Vector::CopyVector(const Vector &y)
-{
-  // check for self assignment
-  if(&y!=this) {
-    for (int i=NDIM;i--;)
-      this->x[i] = y.x[i];
-  }
-}
-
-// this function is completely unused so it is deactivated until new uses arrive and a new
-// place can be found
-#if 0
-/** Solves a vectorial system consisting of two orthogonal statements and a norm statement.
- * This is linear system of equations to be solved, however of the three given (skp of this vector\
- * with either of the three hast to be zero) only two are linear independent. The third equation
- * is that the vector should be of magnitude 1 (orthonormal). This all leads to a case-based solution
- * where very often it has to be checked whether a certain value is zero or not and thus forked into
- * another case.
- * \param *x1 first vector
- * \param *x2 second vector
- * \param *y third vector
- * \param alpha first angle
- * \param beta second angle
- * \param c norm of final vector
- * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
- * \bug this is not yet working properly
- */
-bool Vector::SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c)
-{
-  double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
-  double ang; // angle on testing
-  double sign[3];
-  int i,j,k;
-  A = cos(alpha) * x1->Norm() * c;
-  B1 = cos(beta + M_PI/2.) * y->Norm() * c;
-  B2 = cos(beta) * x2->Norm() * c;
-  C = c * c;
-  Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
-  int flag = 0;
-  if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
-    if (fabs(x1->x[1]) > MYEPSILON) {
-      flag = 1;
-    } else if (fabs(x1->x[2]) > MYEPSILON) {
-       flag = 2;
-    } else {
-      return false;
-    }
-  }
-  switch (flag) {
-    default:
-    case 0:
-      break;
-    case 2:
-      flip(x1->x[0],x1->x[1]);
-      flip(x2->x[0],x2->x[1]);
-      flip(y->x[0],y->x[1]);
-      //flip(x[0],x[1]);
-      flip(x1->x[1],x1->x[2]);
-      flip(x2->x[1],x2->x[2]);
-      flip(y->x[1],y->x[2]);
-      //flip(x[1],x[2]);
-    case 1:
-      flip(x1->x[0],x1->x[1]);
-      flip(x2->x[0],x2->x[1]);
-      flip(y->x[0],y->x[1]);
-      //flip(x[0],x[1]);
-      flip(x1->x[1],x1->x[2]);
-      flip(x2->x[1],x2->x[2]);
-      flip(y->x[1],y->x[2]);
-      //flip(x[1],x[2]);
-      break;
-  }
-  // now comes the case system
-  D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
-  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
-  D3 = y->x[0]/x1->x[0]*A-B1;
-  Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
-  if (fabs(D1) < MYEPSILON) {
-    Log() << Verbose(2) << "D1 == 0!\n";
-    if (fabs(D2) > MYEPSILON) {
-      Log() << Verbose(3) << "D2 != 0!\n";
-      x[2] = -D3/D2;
-      E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
-      E2 = -x1->x[1]/x1->x[0];
-      Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
-      F1 = E1*E1 + 1.;
-      F2 = -E1*E2;
-      F3 = E1*E1 + D3*D3/(D2*D2) - C;
-      Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
-      if (fabs(F1) < MYEPSILON) {
-        Log() << Verbose(4) << "F1 == 0!\n";
-        Log() << Verbose(4) << "Gleichungssystem linear\n";
-        x[1] = F3/(2.*F2);
-      } else {
-        p = F2/F1;
-        q = p*p - F3/F1;
-        Log() << Verbose(4) << "p " << p << "\tq " << q << endl;
-        if (q < 0) {
-          Log() << Verbose(4) << "q < 0" << endl;
-          return false;
-        }
-        x[1] = p + sqrt(q);
-      }
-      x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
-    } else {
-      Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n";
-      return false;
-    }
-  } else {
-    E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
-    E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
-    Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
-    F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
-    F2 = -(E1*E2 + D2*D3/(D1*D1));
-    F3 = E1*E1 + D3*D3/(D1*D1) - C;
-    Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
-    if (fabs(F1) < MYEPSILON) {
-      Log() << Verbose(3) << "F1 == 0!\n";
-      Log() << Verbose(3) << "Gleichungssystem linear\n";
-      x[2] = F3/(2.*F2);
-    } else {
-      p = F2/F1;
-      q = p*p - F3/F1;
-      Log() << Verbose(3) << "p " << p << "\tq " << q << endl;
-      if (q < 0) {
-        Log() << Verbose(3) << "q < 0" << endl;
-        return false;
-      }
-      x[2] = p + sqrt(q);
-    }
-    x[1] = (-D2 * x[2] - D3)/D1;
-    x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
-  }
-  switch (flag) { // back-flipping
-    default:
-    case 0:
-      break;
-    case 2:
-      flip(x1->x[0],x1->x[1]);
-      flip(x2->x[0],x2->x[1]);
-      flip(y->x[0],y->x[1]);
-      flip(x[0],x[1]);
-      flip(x1->x[1],x1->x[2]);
-      flip(x2->x[1],x2->x[2]);
-      flip(y->x[1],y->x[2]);
-      flip(x[1],x[2]);
-    case 1:
-      flip(x1->x[0],x1->x[1]);
-      flip(x2->x[0],x2->x[1]);
-      flip(y->x[0],y->x[1]);
-      //flip(x[0],x[1]);
-      flip(x1->x[1],x1->x[2]);
-      flip(x2->x[1],x2->x[2]);
-      flip(y->x[1],y->x[2]);
-      flip(x[1],x[2]);
-      break;
-  }
-  // one z component is only determined by its radius (without sign)
-  // thus check eight possible sign flips and determine by checking angle with second vector
-  for (i=0;i<8;i++) {
-    // set sign vector accordingly
-    for (j=2;j>=0;j--) {
-      k = (i & pot(2,j)) << j;
-      Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
-      sign[j] = (k == 0) ? 1. : -1.;
-    }
-    Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
-    // apply sign matrix
-    for (j=NDIM;j--;)
-      x[j] *= sign[j];
-    // calculate angle and check
-    ang = x2->Angle (this);
-    Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
-    if (fabs(ang - cos(beta)) < MYEPSILON) {
-      break;
-    }
-    // unapply sign matrix (is its own inverse)
-    for (j=NDIM;j--;)
-      x[j] *= sign[j];
-  }
-  return true;
-};
-
-#endif
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  rep->SubtractVector(y);
+}
 
 /**
@@ -922,11 +511,14 @@
 bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  Vector a = (*this) - offset;
-  a.InverseMatrixMultiplication(parallelepiped);
-  bool isInside = true;
-
-  for (int i=NDIM;i--;)
-    isInside = isInside && ((a.x[i] <= 1) && (a.x[i] >= 0));
-
-  return isInside;
-}
+  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
+  return rep->IsInParallelepiped(offset, parallelepiped);
+}
+
+bool Vector::isBaseClass() const{
+  return true;
+}
+
+Vector* Vector::clone() const{
+  ASSERT(false, "Cannot clone a base Vector object");
+  return 0;
+}
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/vector.hpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -15,4 +15,6 @@
 #include <gsl/gsl_multimin.h>
 
+#include <memory>
+
 #include "defs.hpp"
 
@@ -26,5 +28,5 @@
   // this struct is used to indicate calls to the Baseconstructor from inside vectors.
   struct Baseconstructor{};
-  public:
+public:
 
   Vector();
@@ -33,5 +35,6 @@
   virtual ~Vector();
 
-  virtual double Distance(const Vector &y) const;
+  // Method implemented by forwarding to the Representation
+
   virtual double DistanceSquared(const Vector &y) const;
   virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
@@ -39,6 +42,4 @@
   virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
   virtual double ScalarProduct(const Vector &y) const;
-  virtual double Norm() const;
-  virtual double NormSquared() const;
   virtual double Angle(const Vector &y) const;
   virtual bool IsZero() const;
@@ -49,25 +50,16 @@
   virtual void AddVector(const Vector &y);
   virtual void SubtractVector(const Vector &y);
-  virtual void CopyVector(const Vector &y);
   virtual void VectorProduct(const Vector &y);
   virtual void ProjectOntoPlane(const Vector &y);
   virtual void ProjectIt(const Vector &y);
   virtual Vector Projection(const Vector &y) const;
-  virtual void Zero();
-  virtual void One(const double one);
-  virtual void Init(const double x1, const double x2, const double x3);
-  virtual void Normalize();
-  virtual void Translate(const Vector &x);
   virtual void Mirror(const Vector &x);
-  virtual void Scale(const double ** const factor);
-  virtual void Scale(const double * const factor);
+  virtual void ScaleAll(const double *factor);
   virtual void Scale(const double factor);
   virtual void MatrixMultiplication(const double * const M);
   virtual bool InverseMatrixMultiplication(const double * const M);
   virtual void KeepPeriodic(const double * const matrix);
-  virtual void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
   virtual bool GetOneNormalVector(const Vector &x1);
   virtual bool MakeNormalTo(const Vector &y1);
-  //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
   virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   virtual void WrapPeriodically(const double * const M, const double * const Minv);
@@ -76,6 +68,6 @@
   virtual double& operator[](size_t i);
   virtual const double& operator[](size_t i) const;
-  virtual double& at(size_t i);
-  virtual const double& at(size_t i) const;
+  double& at(size_t i);
+  const double& at(size_t i) const;
 
   // Assignment operator
@@ -84,4 +76,13 @@
   // Access to internal structure
   virtual double* get();
+
+  // Methods that are derived directly from other methods
+  double Distance(const Vector &y) const;
+  double Norm() const;
+  double NormSquared() const;
+  void Normalize();
+  void Zero();
+  void One(const double one);
+  void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
 
   // operators for mathematical operations
@@ -92,6 +93,17 @@
   Vector const operator-(const Vector& b) const;
 
+protected:
+  typedef std::auto_ptr<Vector> rep_ptr;
+  Vector(Baseconstructor);
+  Vector(Baseconstructor,const Vector*);
+  static Vector VecFromRep(const Vector*);
+
 private:
-  double x[NDIM];
+  // method used for protection, i.e. to avoid infinite recursion
+  // when our internal rep becomes messed up
+  virtual bool isBaseClass() const;
+  virtual Vector* clone() const;
+  // this is used to represent the vector internally
+  rep_ptr rep;
 
 };
Index: src/vector_ops.cpp
===================================================================
--- src/vector_ops.cpp	(revision 2733825a2c3f421fa61f61e58f954fa295e6e23d)
+++ src/vector_ops.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
@@ -120,5 +120,5 @@
   Vector res;
   // normalise this vector with respect to axis
-  a.CopyVector(vec);
+  a = vec;
   a.ProjectOntoPlane(axis);
   // construct normal vector
