Index: src/Helpers/fast_functions.hpp
===================================================================
--- src/Helpers/fast_functions.hpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ src/Helpers/fast_functions.hpp	(revision 753f026a17e1d05846e903705efe2ac582c0d091)
@@ -8,4 +8,6 @@
 #ifndef FAST_FUNCTIONS_HPP_
 #define FAST_FUNCTIONS_HPP_
+
+#include "defs.hpp"
 
 /**
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ src/Makefile.am	(revision 753f026a17e1d05846e903705efe2ac582c0d091)
@@ -9,5 +9,4 @@
 			   gslvector.cpp \
 			   linearsystemofequations.cpp \
-			   SingleVector.cpp \
 			   vector.cpp
 			   
@@ -15,5 +14,4 @@
 			   gslvector.hpp \
 			   linearsystemofequations.hpp \
-			   SingleVector.hpp \
 			   vector.hpp
 			   
Index: c/SingleVector.cpp
===================================================================
--- src/SingleVector.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ 	(revision )
@@ -1,537 +1,0 @@
-/** \file vector.cpp
- *
- * Function implementations for the class vector.
- *
- */
-
-
-#include "defs.hpp"
-#include "gslmatrix.hpp"
-#include "leastsquaremin.hpp"
-#include "memoryallocator.hpp"
-#include "SingleVector.hpp"
-#include "Helpers/fast_functions.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 ************************************/
-
-/** Constructor of class vector.
- */
-SingleVector::SingleVector() :
-  Vector(Baseconstructor())
-{
-  x[0] = x[1] = x[2] = 0.;
-};
-
-/** Constructor of class vector.
- */
-SingleVector::SingleVector(const double x1, const double x2, const double x3) :
-  Vector(Baseconstructor())
-{
-  x[0] = x1;
-  x[1] = x2;
-  x[2] = x3;
-};
-
-/**
- * Copy constructor
- */
-SingleVector::SingleVector(const Vector& src) :
-  Vector(Baseconstructor())
-{
-  x[0] = src[0];
-  x[1] = src[1];
-  x[2] = src[2];
-}
-
-SingleVector* SingleVector::clone() const{
-  return new SingleVector(x[0],x[1],x[2]);
-}
-
-/** Desctructor of class vector.
- */
-SingleVector::~SingleVector() {};
-
-/** Calculates square of distance between this and another vector.
- * \param *y array to second vector
- * \return \f$| x - y |^2\f$
- */
-double SingleVector::DistanceSquared(const Vector &y) const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += (x[i]-y[i])*(x[i]-y[i]);
-  return (res);
-};
-
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |\f$
- */
-double SingleVector::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[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);
-};
-
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |^2\f$
- */
-double SingleVector::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[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);
-};
-
-/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
- * \param *out ofstream for debugging messages
- * Tries to translate a vector into each adjacent neighbouring cell.
- */
-void SingleVector::KeepPeriodic(const double * const matrix)
-{
-//  int N[NDIM];
-//  bool flag = false;
-  //vector Shifted, TranslationVector;
-//  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
-//  Log() << Verbose(2) << "Vector is: ";
-//  Output(out);
-//  Log() << Verbose(0) << endl;
-  SingleVector TestVector(*this);
-  TestVector.InverseMatrixMultiplication(matrix);
-  for(int i=NDIM;i--;) { // correct periodically
-    if (TestVector[i] < 0) {  // get every coefficient into the interval [0,1)
-      TestVector[i] += ceil(TestVector[i]);
-    } else {
-      TestVector[i] -= floor(TestVector[i]);
-    }
-  }
-  TestVector.MatrixMultiplication(matrix);
-  CopyVector(TestVector);
-//  Log() << Verbose(2) << "New corrected vector is: ";
-//  Output(out);
-//  Log() << Verbose(0) << endl;
-//  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
-};
-
-/** Calculates scalar product between this and another vector.
- * \param *y array to second vector
- * \return \f$\langle x, y \rangle\f$
- */
-double SingleVector::ScalarProduct(const Vector &y) const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += x[i]*y[i];
-  return (res);
-};
-
-
-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
- *  -# ATTENTION: Only three dim.
- *  \param *y array to vector with which to calculate crossproduct
- *  \return \f$ x \times y \f&
- */
-void SingleVector::VectorProduct(const Vector &y)
-{
-  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]);
-  CopyVector(tmp);
-};
-
-
-/** projects this vector onto plane defined by \a *y.
- * \param *y normal vector of plane
- * \return \f$\langle x, y \rangle\f$
- */
-void SingleVector::ProjectOntoPlane(const Vector &y)
-{
-  Vector tmp;
-  tmp = y;
-  tmp.Normalize();
-  tmp.Scale(ScalarProduct(tmp));
-  *this -= tmp;
-};
-
-/** Calculates the minimum distance of this vector to the plane.
- * \param *out output stream for debugging
- * \param *PlaneNormal normal of plane
- * \param *PlaneOffset offset of plane
- * \return distance to plane
- */
-double SingleVector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
-{
-  // first create part that is orthonormal to PlaneNormal with withdraw
-  Vector temp = VecFromRep(this)- PlaneOffset;
-  temp.MakeNormalTo(PlaneNormal);
-  temp.Scale(-1.);
-  // then add connecting vector from plane to point
-  temp += VecFromRep(this)-PlaneOffset;
-  double sign = temp.ScalarProduct(PlaneNormal);
-  if (fabs(sign) > MYEPSILON)
-    sign /= fabs(sign);
-  else
-    sign = 0.;
-
-  return (temp.Norm()*sign);
-};
-
-/** Calculates the projection of a vector onto another \a *y.
- * \param *y array to second vector
- */
-void SingleVector::ProjectIt(const Vector &y)
-{
-  Vector helper = y;
-  helper.Scale(-(ScalarProduct(y)));
-  AddVector(helper);
-};
-
-/** Calculates the projection of a vector onto another \a *y.
- * \param *y array to second vector
- * \return Vector
- */
-Vector SingleVector::Projection(const Vector &y) const
-{
-  Vector helper = y;
-  helper.Scale((ScalarProduct(y)/y.NormSquared()));
-
-  return helper;
-};
-
-/** Checks whether vector has all components zero.
- * @return true - vector is zero, false - vector is not
- */
-bool SingleVector::IsZero() const
-{
-  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
-};
-
-/** Checks whether vector has length of 1.
- * @return true - vector is normalized, false - vector is not
- */
-bool SingleVector::IsOne() const
-{
-  return (fabs(Norm() - 1.) < MYEPSILON);
-};
-
-/** Checks whether vector is normal to \a *normal.
- * @return true - vector is normalized, false - vector is not
- */
-bool SingleVector::IsNormalTo(const Vector &normal) const
-{
-  if (ScalarProduct(normal) < MYEPSILON)
-    return true;
-  else
-    return false;
-};
-
-/** Checks whether vector is normal to \a *normal.
- * @return true - vector is normalized, false - vector is not
- */
-bool SingleVector::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;
-};
-
-/** Calculates the angle between this and another vector.
- * \param *y array to second vector
- * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
- */
-double SingleVector::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);
-};
-
-
-double& SingleVector::operator[](size_t i){
-  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
-}
-
-const double& SingleVector::operator[](size_t i) const{
-  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
-  return x[i];
-}
-
-double* SingleVector::get(){
-  return x;
-}
-/** Scales each atom coordinate by an individual \a factor.
- * \param *factor pointer to scaling factor
- */
-void SingleVector::ScaleAll(const double *factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= factor[i];
-};
-
-void SingleVector::Scale(const double factor)
-{
-  for (int i=NDIM;i--;)
-    x[i] *= factor;
-};
-
-/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
- * \param *M matrix of box
- * \param *Minv inverse matrix
- */
-void SingleVector::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);
-};
-
-/** Do a matrix multiplication.
- * \param *matrix NDIM_NDIM array
- */
-void SingleVector::MatrixMultiplication(const double * const M)
-{
-  SingleVector C;
-  // do the matrix multiplication
-  C[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
-  C[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
-  C[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
-
-  CopyVector(C);
-};
-
-/** Do a matrix multiplication with the \a *A' inverse.
- * \param *matrix NDIM_NDIM array
- */
-bool SingleVector::InverseMatrixMultiplication(const double * const A)
-{
-  SingleVector 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[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
-    C[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
-    C[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
-
-    CopyVector(C);
-    return true;
-  } else {
-    return false;
-  }
-};
-
-
-/** Mirrors atom against a given plane.
- * \param n[] normal vector of mirror plane.
- */
-void SingleVector::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];
-};
-
-
-/** Calculates orthonormal vector to one given vector.
- * Just subtracts the projection onto the given vector from this vector.
- * The removed part of the vector is Vector::Projection()
- * \param *x1 vector
- * \return true - success, false - vector is zero
- */
-bool SingleVector::MakeNormalTo(const Vector &y1)
-{
-  bool result = false;
-  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
-  Vector x1;
-  x1 = factor * y1;
-  SubtractVector(x1);
-  for (int i=NDIM;i--;)
-    result = result || (fabs(x[i]) > MYEPSILON);
-
-  return result;
-};
-
-/** Creates this vector as one of the possible orthonormal ones to the given one.
- * Just scan how many components of given *vector are unequal to zero and
- * try to get the skp of both to be zero accordingly.
- * \param *vector given vector
- * \return true - success, false - failure (null vector given)
- */
-bool SingleVector::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;
-  }
-};
-
-/**
- * Checks whether this vector is within the parallelepiped defined by the given three vectors and
- * their offset.
- *
- * @param offest for the origin of the parallelepiped
- * @param three vectors forming the matrix that defines the shape of the parallelpiped
- */
-bool SingleVector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
-{
-  Vector a = VecFromRep(this);
-  a -= offset;
-  a.InverseMatrixMultiplication(parallelepiped);
-  bool isInside = true;
-
-  for (int i=NDIM;i--;)
-    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
-
-  return isInside;
-}
-
-
-void SingleVector::CopyVector(SingleVector& vec) {
-  for(int i =NDIM; i--;)
-    x[i]=vec.x[i];
-}
-
-bool SingleVector::isBaseClass() const{
-  return false;
-}
Index: c/SingleVector.hpp
===================================================================
--- src/SingleVector.hpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ 	(revision )
@@ -1,79 +1,0 @@
-#ifndef SINGLEVECTOR_HPP_
-#define SINGLEVECTOR_HPP_
-
-using namespace std;
-
-/*********************************************** includes ***********************************/
-
-// include config.h
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <iostream>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_multimin.h>
-#include "vector.hpp"
-
-#include "defs.hpp"
-
-/********************************************** declarations *******************************/
-
-/** Single vector.
- * basically, just a x[3] but with helpful functions
- */
-class SingleVector : public Vector {
-public:
-
-  SingleVector();
-  SingleVector(const double x1, const double x2, const double x3);
-  SingleVector(const Vector& src);
-  virtual ~SingleVector();
-
-  virtual double DistanceSquared(const Vector &y) const;
-  virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
-  virtual double PeriodicDistance(const Vector &y, const double * const cell_size) const;
-  virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
-  virtual double ScalarProduct(const Vector &y) const;
-  virtual double Angle(const Vector &y) const;
-  virtual bool IsZero() const;
-  virtual bool IsOne() const;
-  virtual bool IsNormalTo(const Vector &normal) const;
-  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 Mirror(const Vector &x);
-  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 bool GetOneNormalVector(const Vector &x1);
-  virtual bool MakeNormalTo(const Vector &y1);
-  virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
-  virtual void WrapPeriodically(const double * const M, const double * const Minv);
-
-  // Accessors ussually come in pairs... and sometimes even more than that
-  virtual double& operator[](size_t i);
-  virtual const double& operator[](size_t i) const;
-
-  // 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];
-
-};
-
-#endif /*SINGLEVECTOR_HPP_*/
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ src/vector.cpp	(revision 753f026a17e1d05846e903705efe2ac582c0d091)
@@ -6,6 +6,7 @@
 
 
-#include "SingleVector.hpp"
+#include "vector.hpp"
 #include "Helpers/Assert.hpp"
+#include "Helpers/fast_functions.hpp"
 
 #include <iostream>
@@ -18,31 +19,28 @@
 /** Constructor of class vector.
  */
-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) :
-  rep(new SingleVector(x1,x2,x3))
-{};
+Vector::Vector()
+{
+  x[0] = x[1] = x[2] = 0.;
+};
 
 /**
  * Copy constructor
  */
-Vector::Vector(const Vector& src) :
-  rep(src.rep->clone())
-{}
+
+Vector::Vector(const Vector& src)
+{
+  x[0] = src[0];
+  x[1] = src[1];
+  x[2] = src[2];
+}
+
+/** Constructor of class vector.
+ */
+Vector::Vector(const double x1, const double x2, const double x3)
+{
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
 
 /**
@@ -50,8 +48,9 @@
  */
 Vector& Vector::operator=(const Vector& src){
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   // check for self assignment
   if(&src!=this){
-    rep.reset(src.rep->clone());
+    x[0] = src[0];
+    x[1] = src[1];
+    x[2] = src[2];
   }
   return *this;
@@ -68,6 +67,8 @@
 double Vector::DistanceSquared(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->DistanceSquared(y);
+  double res = 0.;
+  for (int i=NDIM;i--;)
+    res += (x[i]-y[i])*(x[i]-y[i]);
+  return (res);
 };
 
@@ -88,6 +89,32 @@
 double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->PeriodicDistance(y,cell_size);
+  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[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);
 };
 
@@ -99,6 +126,32 @@
 double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->PeriodicDistanceSquared(y,cell_size);
+  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[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);
 };
 
@@ -109,6 +162,24 @@
 void Vector::KeepPeriodic(const double * const matrix)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->KeepPeriodic(matrix);
+  //  int N[NDIM];
+  //  bool flag = false;
+    //vector Shifted, TranslationVector;
+  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
+  //  Log() << Verbose(2) << "Vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+    InverseMatrixMultiplication(matrix);
+    for(int i=NDIM;i--;) { // correct periodically
+      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
+        at(i) += ceil(at(i));
+      } else {
+        at(i) -= floor(at(i));
+      }
+    }
+    MatrixMultiplication(matrix);
+  //  Log() << Verbose(2) << "New corrected vector is: ";
+  //  Output(out);
+  //  Log() << Verbose(0) << endl;
+  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
 };
 
@@ -119,6 +190,8 @@
 double Vector::ScalarProduct(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->ScalarProduct(y);
+  double res = 0.;
+  for (int i=NDIM;i--;)
+    res += x[i]*y[i];
+  return (res);
 };
 
@@ -132,6 +205,9 @@
 void Vector::VectorProduct(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->VectorProduct(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;
 };
 
@@ -143,6 +219,9 @@
 void Vector::ProjectOntoPlane(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ProjectOntoPlane(y);
+  Vector tmp;
+  tmp = y;
+  tmp.Normalize();
+  tmp.Scale(ScalarProduct(tmp));
+  *this -= tmp;
 };
 
@@ -155,6 +234,17 @@
 double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->DistanceToPlane(PlaneNormal,PlaneOffset);
+  // first create part that is orthonormal to PlaneNormal with withdraw
+  Vector temp = (*this) - PlaneOffset;
+  temp.MakeNormalTo(PlaneNormal);
+  temp.Scale(-1.);
+  // then add connecting vector from plane to point
+  temp += (*this)-PlaneOffset;
+  double sign = temp.ScalarProduct(PlaneNormal);
+  if (fabs(sign) > MYEPSILON)
+    sign /= fabs(sign);
+  else
+    sign = 0.;
+
+  return (temp.Norm()*sign);
 };
 
@@ -164,6 +254,5 @@
 void Vector::ProjectIt(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ProjectIt(y);
+  (*this) += (-ScalarProduct(y))*y;
 };
 
@@ -174,6 +263,8 @@
 Vector Vector::Projection(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->Projection(y);
+  Vector helper = y;
+  helper.Scale((ScalarProduct(y)/y.NormSquared()));
+
+  return helper;
 };
 
@@ -206,5 +297,5 @@
 void Vector::Zero()
 {
-  rep.reset(new SingleVector());
+  at(0)=at(1)=at(2)=0;
 };
 
@@ -213,5 +304,5 @@
 void Vector::One(const double one)
 {
-  rep.reset(new SingleVector(one,one,one));
+  at(0)=at(1)=at(2)=one;
 };
 
@@ -221,6 +312,5 @@
 bool Vector::IsZero() const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsZero();
+  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
 };
 
@@ -230,6 +320,5 @@
 bool Vector::IsOne() const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsOne();
+  return (fabs(Norm() - 1.) < MYEPSILON);
 };
 
@@ -239,6 +328,8 @@
 bool Vector::IsNormalTo(const Vector &normal) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsNormalTo(normal);
+  if (ScalarProduct(normal) < MYEPSILON)
+    return true;
+  else
+    return false;
 };
 
@@ -248,6 +339,10 @@
 bool Vector::IsEqualTo(const Vector &a) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->IsEqualTo(a);
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };
 
@@ -258,17 +353,26 @@
 double Vector::Angle(const Vector &y) const
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->Angle(y);
+  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);
 };
 
 
 double& Vector::operator[](size_t i){
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return (*rep)[i];
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
 }
 
 const double& Vector::operator[](size_t i) const{
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return (*rep)[i];
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
 }
 
@@ -282,6 +386,5 @@
 
 double* Vector::get(){
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->get();
+  return x;
 }
 
@@ -293,5 +396,4 @@
 bool Vector::operator==(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   return IsEqualTo(b);
 };
@@ -337,5 +439,4 @@
 Vector const Vector::operator+(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.AddVector(b);
@@ -350,5 +451,4 @@
 Vector const Vector::operator-(const Vector& b) const
 {
-  ASSERT(isBaseClass(),"Operator used on Derived Vector object");
   Vector x = *this;
   x.SubtractVector(b);
@@ -395,6 +495,6 @@
 void Vector::ScaleAll(const double *factor)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->ScaleAll(factor);
+  for (int i=NDIM;i--;)
+    x[i] *= factor[i];
 };
 
@@ -403,6 +503,6 @@
 void Vector::Scale(const double factor)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->Scale(factor);
+  for (int i=NDIM;i--;)
+    x[i] *= factor;
 };
 
@@ -413,6 +513,14 @@
 void Vector::WrapPeriodically(const double * const M, const double * const Minv)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->WrapPeriodically(M,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);
 };
 
@@ -422,6 +530,8 @@
 void Vector::MatrixMultiplication(const double * const M)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->MatrixMultiplication(M);
+  // do the matrix multiplication
+  at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
+  at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
+  at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
 };
 
@@ -431,6 +541,30 @@
 bool Vector::InverseMatrixMultiplication(const double * const A)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->InverseMatrixMultiplication(A);
+  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
+    at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
+    at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
+    at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
+
+    return true;
+  } else {
+    return false;
+  }
 };
 
@@ -455,6 +589,9 @@
 void Vector::Mirror(const Vector &n)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->Mirror(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];
 };
 
@@ -468,6 +605,13 @@
 bool Vector::MakeNormalTo(const Vector &y1)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->MakeNormalTo(y1);
+  bool result = false;
+  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
+  Vector x1;
+  x1 = factor * y1;
+  SubtractVector(x1);
+  for (int i=NDIM;i--;)
+    result = result || (fabs(x[i]) > MYEPSILON);
+
+  return result;
 };
 
@@ -480,6 +624,36 @@
 bool Vector::GetOneNormalVector(const Vector &GivenVector)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  return rep->GetOneNormalVector(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;
+  }
 };
 
@@ -489,6 +663,6 @@
 void Vector::AddVector(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->AddVector(y);
+  for(int i=NDIM;i--;)
+    x[i] += y[i];
 }
 
@@ -498,6 +672,6 @@
 void Vector::SubtractVector(const Vector &y)
 {
-  ASSERT((rep.get()) && (!rep->isBaseClass()),"Representation stored in vector Object was not of derived type");
-  rep->SubtractVector(y);
+  for(int i=NDIM;i--;)
+    x[i] -= y[i];
 }
 
@@ -511,14 +685,11 @@
 bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  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;
-}
+  Vector a = (*this)-offset;
+  a.InverseMatrixMultiplication(parallelepiped);
+  bool isInside = true;
+
+  for (int i=NDIM;i--;)
+    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
+
+  return isInside;
+}
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 1bd79efa209ce2203a329de83169d52ea290fa17)
+++ src/vector.hpp	(revision 753f026a17e1d05846e903705efe2ac582c0d091)
@@ -37,43 +37,43 @@
   // Method implemented by forwarding to the Representation
 
-  virtual double DistanceSquared(const Vector &y) const;
-  virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
-  virtual double PeriodicDistance(const Vector &y, const double * const cell_size) const;
-  virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
-  virtual double ScalarProduct(const Vector &y) const;
-  virtual double Angle(const Vector &y) const;
-  virtual bool IsZero() const;
-  virtual bool IsOne() const;
-  virtual bool IsNormalTo(const Vector &normal) const;
-  virtual bool IsEqualTo(const Vector &a) const;
+  double DistanceSquared(const Vector &y) const;
+  double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
+  double PeriodicDistance(const Vector &y, const double * const cell_size) const;
+  double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
+  double ScalarProduct(const Vector &y) const;
+  double Angle(const Vector &y) const;
+  bool IsZero() const;
+  bool IsOne() const;
+  bool IsNormalTo(const Vector &normal) const;
+  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 Mirror(const Vector &x);
-  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 bool GetOneNormalVector(const Vector &x1);
-  virtual bool MakeNormalTo(const Vector &y1);
-  virtual bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
-  virtual void WrapPeriodically(const double * const M, const double * const Minv);
+  void AddVector(const Vector &y);
+  void SubtractVector(const Vector &y);
+  void VectorProduct(const Vector &y);
+  void ProjectOntoPlane(const Vector &y);
+  void ProjectIt(const Vector &y);
+  Vector Projection(const Vector &y) const;
+  void Mirror(const Vector &x);
+  void ScaleAll(const double *factor);
+  void Scale(const double factor);
+  void MatrixMultiplication(const double * const M);
+  bool InverseMatrixMultiplication(const double * const M);
+  void KeepPeriodic(const double * const matrix);
+  bool GetOneNormalVector(const Vector &x1);
+  bool MakeNormalTo(const Vector &y1);
+  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
+  void WrapPeriodically(const double * const M, const double * const Minv);
 
   // Accessors ussually come in pairs... and sometimes even more than that
-  virtual double& operator[](size_t i);
-  virtual const double& operator[](size_t i) const;
+  double& operator[](size_t i);
+  const double& operator[](size_t i) const;
   double& at(size_t i);
   const double& at(size_t i) const;
 
   // Assignment operator
-  virtual Vector &operator=(const Vector& src);
+  Vector &operator=(const Vector& src);
 
   // Access to internal structure
-  virtual double* get();
+  double* get();
 
   // Methods that are derived directly from other methods
@@ -94,16 +94,7 @@
 
 protected:
-  typedef std::auto_ptr<Vector> rep_ptr;
-  Vector(Baseconstructor);
-  Vector(Baseconstructor,const Vector*);
-  static Vector VecFromRep(const Vector*);
 
 private:
-  // 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;
+  double x[NDIM];
 
 };
