Index: src/Exceptions/MultipleSolutionsException.hpp
===================================================================
--- src/Exceptions/MultipleSolutionsException.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
+++ src/Exceptions/MultipleSolutionsException.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -0,0 +1,38 @@
+/*
+ * InfiniteSolutionsException.hpp
+ *
+ *  Created on: May 28, 2010
+ *      Author: crueger
+ */
+
+#ifndef INFINITESOLUTIONSEXCEPTION_HPP_
+#define INFINITESOLUTIONSEXCEPTION_HPP_
+
+template<class T>
+class MultipleSolutionsException : public MathException
+{
+public:
+  MultipleSolutionsException(std::string _file, int _line, const T& _solution) :
+    MathException(_file,_line),
+    solution(_solution)
+  {}
+
+  virtual ~MultipleSolutionsException() throw(){}
+
+  virtual const char* what() const throw(){
+    return "Multiple solutions possible....";
+  }
+
+  /**
+   * use this method in case you don't care that there is more than one
+   * solution and just any solution is fine
+   */
+  virtual T getAnySolution(){
+    return solution;
+  }
+
+private:
+  T solution;
+};
+
+#endif /* MULTIPLESOLUTIONSEXCEPTION_HPP_ */
Index: src/Exceptions/SkewException.cpp
===================================================================
--- src/Exceptions/SkewException.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
+++ src/Exceptions/SkewException.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -0,0 +1,22 @@
+/*
+ * SkewException.cpp
+ *
+ *  Created on: May 27, 2010
+ *      Author: crueger
+ */
+
+#include "SkewException.hpp"
+
+using namespace std;
+
+SkewException::SkewException(string _file, int _line) :
+  MathException(_file,_line)
+{}
+
+SkewException::~SkewException() throw()
+{}
+
+const char* SkewException::what() const throw()
+{
+  return "Given lines were askew.";
+}
Index: src/Exceptions/SkewException.hpp
===================================================================
--- src/Exceptions/SkewException.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
+++ src/Exceptions/SkewException.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -0,0 +1,22 @@
+/*
+ * SkewException.hpp
+ *
+ *  Created on: May 27, 2010
+ *      Author: crueger
+ */
+
+#ifndef SKEWEXCEPTION_HPP_
+#define SKEWEXCEPTION_HPP_
+
+#include "Exceptions/MathException.hpp"
+
+class SkewException : public MathException
+{
+public:
+  SkewException(std::string _file, int _line);
+  virtual ~SkewException() throw ();
+
+  virtual const char* what() const throw();
+};
+
+#endif /* SKEWEXCEPTION_HPP_ */
Index: src/Legacy/oldmenu.cpp
===================================================================
--- src/Legacy/oldmenu.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Legacy/oldmenu.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -26,4 +26,5 @@
 #include "vector_ops.hpp"
 #include "Plane.hpp"
+#include "Line.hpp"
 
 #include "UIElements/UIFactory.hpp"
@@ -185,5 +186,5 @@
           // rotate vector around first angle
           first->x = x;
-          first->x = RotateVector(first->x,z,b - M_PI);
+          first->x = Line(zeroVec,z).rotateVector(first->x,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
           // remove the projection onto the rotation plane of the second angle
@@ -201,5 +202,5 @@
           // rotate another vector around second angle
           n = y;
-          n = RotateVector(n,x,c - M_PI);
+          n = Line(zeroVec,x).rotateVector(n,c - M_PI);
           Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 
Index: src/Line.cpp
===================================================================
--- src/Line.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Line.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -11,11 +11,25 @@
 
 #include "vector.hpp"
-
-Line::Line(Vector &_origin, Vector &_direction) :
-  origin(new Vector(_origin)),
+#include "log.hpp"
+#include "verbose.hpp"
+#include "gslmatrix.hpp"
+#include "info.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
+#include "Plane.hpp"
+
+using namespace std;
+
+Line::Line(const Vector &_origin, const Vector &_direction) :
   direction(new Vector(_direction))
 {
   direction->Normalize();
-}
+  origin.reset(new Vector(_origin.partition(*direction).second));
+}
+
+Line::Line(const Line &src) :
+  origin(new Vector(*src.origin)),
+  direction(new Vector(*src.direction))
+{}
 
 Line::~Line()
@@ -24,9 +38,183 @@
 
 double Line::distance(const Vector &point) const{
-  return fabs(point.ScalarProduct(*direction) - origin->ScalarProduct(*direction));
+  // get any vector from line to point
+  Vector helper = point - *origin;
+  // partition this vector along direction
+  // the residue points from the line to the point
+  return helper.partition(*direction).second.Norm();
 }
 
 Vector Line::getClosestPoint(const Vector &point) const{
-  double factor = point.ScalarProduct(*direction) - origin->ScalarProduct(*direction);
-  return (point - (factor * (*direction)));
-}
+  // get any vector from line to point
+  Vector helper = point - *origin;
+  // partition this vector along direction
+  // add only the part along the direction
+  return *origin + helper.partition(*direction).first;
+}
+
+Vector Line::getDirection() const{
+  return *direction;
+}
+
+Vector Line::getOrigin() const{
+  return *origin;
+}
+
+vector<Vector> Line::getPointsOnLine() const{
+  vector<Vector> res;
+  res.reserve(2);
+  res.push_back(*origin);
+  res.push_back(*origin+*direction);
+  return res;
+}
+
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
+Vector Line::getIntersection(const Line& otherLine) const{
+  Info FunctionInfo(__func__);
+
+  pointset line1Points = getPointsOnLine();
+
+  Vector Line1a = line1Points[0];
+  Vector Line1b = line1Points[1];
+
+  pointset line2Points = otherLine.getPointsOnLine();
+
+  Vector Line2a = line2Points[0];
+  Vector Line2b = line2Points[1];
+
+  Vector res;
+
+  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
+
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a[i]);
+    M->Set(1, i, Line1b[i]);
+    M->Set(2, i, Line2a[i]);
+    M->Set(3, i, Line2b[i]);
+  }
+
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
+    throw SkewException(__FILE__,__LINE__);
+  }
+
+  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
+
+
+  // constuct a,b,c
+  Vector a = Line1b - Line1a;
+  Vector b = Line2b - Line2a;
+  Vector c = Line2a - Line1a;
+  Vector d = Line2b - Line1b;
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   res.Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel = Line1a - Line2a;
+    factor = parallel.ScalarProduct(a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      res = Line2a;
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return res;
+    } else {
+      parallel = Line1a - Line2b;
+      factor = parallel.ScalarProduct(a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        res = Line2b;
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return res;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    res.Zero();
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1 = c;
+  temp1.VectorProduct(b);
+  temp2 = a;
+  temp2.VectorProduct(b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+
+  // construct intersection
+  res = a;
+  res.Scale(s);
+  res += Line1a;
+  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
+
+  return res;
+}
+
+/** Rotates the vector by an angle of \a alpha around this line.
+ * \param rhs Vector to rotate
+ * \param alpha rotation angle in radian
+ */
+Vector Line::rotateVector(const Vector &rhs, double alpha) const{
+  Vector helper = rhs;
+
+  // translate the coordinate system so that the line goes through (0,0,0)
+  helper -= *origin;
+
+  // partition the vector into a part that gets rotated and a part that lies along the line
+  pair<Vector,Vector> parts = helper.partition(*direction);
+
+  // we just keep anything that is along the axis
+  Vector res = parts.first;
+
+  // the rest has to be rotated
+  Vector a = parts.second;
+  // we only have to do the rest, if we actually could partition the vector
+  if(!a.IsZero()){
+    // construct a vector that is orthogonal to a and direction and has length |a|
+    Vector y = a;
+    // direction is normalized, so the result has length |a|
+    y.VectorProduct(*direction);
+
+    res += cos(alpha) * a + sin(alpha) * y;
+  }
+
+  // translate the coordinate system back
+  res += *origin;
+  return res;
+}
+
+Plane Line::getOrthogonalPlane(const Vector &origin) const{
+  return Plane(getDirection(),origin);
+}
+
+Line makeLineThrough(const Vector &x1, const Vector &x2){
+  if(x1==x2){
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+  return Line(x1,x1-x2);
+}
Index: src/Line.hpp
===================================================================
--- src/Line.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Line.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -12,15 +12,29 @@
 
 #include <memory>
+#include <vector>
 
 class Vector;
+class Plane;
 
 class Line : public Space
 {
 public:
-  Line(Vector &_origin, Vector &_direction);
+  Line(const Vector &_origin, const Vector &_direction);
+  Line(const Line& _src);
   virtual ~Line();
 
-  virtual double distance(const Vector &point) const=0;
-  virtual Vector getClosestPoint(const Vector &point) const=0;
+  virtual double distance(const Vector &point) const;
+  virtual Vector getClosestPoint(const Vector &point) const;
+
+  Vector getDirection() const;
+  Vector getOrigin() const;
+
+  std::vector<Vector> getPointsOnLine() const;
+
+  Vector getIntersection(const Line& otherLine) const;
+
+  Vector rotateVector(const Vector &rhs, double alpha) const;
+
+  Plane getOrthogonalPlane(const Vector &origin) const;
 
 private:
@@ -29,3 +43,8 @@
 };
 
+/**
+ * Named constructor to make a line through two points
+ */
+Line makeLineThrough(const Vector &x1, const Vector &x2);
+
 #endif /* LINE_HPP_ */
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Makefile.am	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -125,12 +125,14 @@
 				   
 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
-  Exceptions/LinearDependenceException.cpp \
-  Exceptions/MathException.cpp \
-  Exceptions/ZeroVectorException.cpp
+				  Exceptions/LinearDependenceException.cpp \
+				  Exceptions/MathException.cpp \
+				  Exceptions/SkewException.cpp \
+				  Exceptions/ZeroVectorException.cpp
 				  
 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
-  Exceptions/LinearDependenceException.hpp \
-  Exceptions/MathException.hpp \
-  Exceptions/ZeroVectorException.hpp
+				  Exceptions/LinearDependenceException.hpp \
+				  Exceptions/MathException.hpp \
+				  Exceptions/SkewException.hpp \
+				  Exceptions/ZeroVectorException.hpp
 
 SOURCE = \
Index: src/Plane.cpp
===================================================================
--- src/Plane.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Plane.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -14,4 +14,6 @@
 #include "Helpers/Assert.hpp"
 #include <cmath>
+#include "Line.hpp"
+#include "Exceptions/MultipleSolutionsException.hpp"
 
 /**
@@ -42,5 +44,4 @@
 /**
  * Constructs a plane from two direction vectors and a offset.
- * If no offset is given a plane through origin is assumed
  */
 Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException) :
@@ -113,9 +114,9 @@
 }
 
-Vector Plane::getOffsetVector() {
+Vector Plane::getOffsetVector() const {
   return getOffset()*getNormal();
 }
 
-vector<Vector> Plane::getPointsOnPlane(){
+vector<Vector> Plane::getPointsOnPlane() const{
   std::vector<Vector> res;
   res.reserve(3);
@@ -147,39 +148,41 @@
  * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
  */
-Vector Plane::GetIntersection(const Vector &Origin, const Vector &LineVector)
+Vector Plane::GetIntersection(const Line& line) const
 {
   Info FunctionInfo(__func__);
   Vector res;
 
-  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
-  Vector Direction = LineVector - Origin;
-  Direction.Normalize();
-  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
-  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
-  double factor1 = Direction.ScalarProduct(*normalVector.get());
-  if (fabs(factor1) < MYEPSILON) { // Uniqueness: line parallel to plane?
-    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  double factor2 = Origin.ScalarProduct(*normalVector.get());
-  if (fabs(factor2-offset) < MYEPSILON) { // Origin is in-plane
-    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
-    res = Origin;
-    return res;
-  }
-
+  double factor1 = getNormal().ScalarProduct(line.getDirection());
+  if(fabs(factor1)<MYEPSILON){
+    // the plane is parallel... under all circumstances this is bad luck
+    // we no have either no or infinite solutions
+    if(isContained(line.getOrigin())){
+      throw MultipleSolutionsException<Vector>(__FILE__,__LINE__,line.getOrigin());
+    }
+    else{
+      throw LinearDependenceException(__FILE__,__LINE__);
+    }
+  }
+
+  double factor2 = getNormal().ScalarProduct(line.getOrigin());
   double scaleFactor = (offset-factor2)/factor1;
 
-  //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
-  Direction.Scale(scaleFactor);
-  res = Origin + Direction;
-  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
-
-  // test whether resulting vector really is on plane
-  ASSERT(fabs(res.ScalarProduct(*normalVector) - offset) < MYEPSILON,
-         "Calculated line-Plane intersection does not lie on plane.");
+  res = line.getOrigin() + scaleFactor * line.getDirection();
+
+  // tests to make sure the resulting vector really is on plane and line
+  ASSERT(isContained(res),"Calculated line-Plane intersection does not lie on plane.");
+  ASSERT(line.isContained(res),"Calculated line-Plane intersection does not lie on line.");
   return res;
 };
+
+Vector Plane::mirrorVector(const Vector &rhs) const {
+  Vector helper = getVectorToPoint(rhs);
+  // substract twice the Vector to the plane
+  return rhs+2*helper;
+}
+
+Line Plane::getOrthogonalLine(const Vector &origin) const{
+  return Line(origin,getNormal());
+}
 
 /************ Methods inherited from Space ****************/
Index: src/Plane.hpp
===================================================================
--- src/Plane.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Plane.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -17,4 +17,5 @@
 
 class Vector;
+class Line;
 
 class Plane : public Space
@@ -42,13 +43,23 @@
    * Same as getOffset()*getNormal();
    */
-  Vector getOffsetVector();
+  Vector getOffsetVector() const;
 
   /**
    * returns three seperate points on this plane
    */
-  std::vector<Vector> getPointsOnPlane();
+  std::vector<Vector> getPointsOnPlane() const;
 
   // some calculations
-  Vector GetIntersection(const Vector &Origin, const Vector &LineVector);
+  Vector GetIntersection(const Line &Line) const;
+
+  Vector mirrorVector(const Vector &rhs) const;
+
+  /**
+   * get a Line that is orthogonal to this plane, going through a chosen
+   * point.
+   *
+   * The point does not have to lie on the plane itself.
+   */
+  Line getOrthogonalLine(const Vector &origin) const;
 
   /****** Methods inherited from Space ***********/
Index: src/Space.cpp
===================================================================
--- src/Space.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Space.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -19,4 +19,9 @@
 }
 
+Vector Space::getVectorToPoint(const Vector &origin) const{
+  Vector support = getClosestPoint(origin);
+  return support-origin;
+}
+
 bool Space::isContained(const Vector &point) const{
   return (distance(point)) < MYEPSILON;
Index: src/Space.hpp
===================================================================
--- src/Space.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/Space.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -17,7 +17,33 @@
   virtual ~Space();
 
+  /**
+   * Calculates the distance between a Space and a Vector.
+   */
   virtual double distance(const Vector &point) const=0;
+
+  /**
+   * get the closest point inside the space to another point
+   */
   virtual Vector getClosestPoint(const Vector &point) const=0;
+
+  /**
+   * get the shortest Vector from a point to a Space.
+   *
+   * The Vector always points from the given Vector to the point in space
+   * returned by Plane::getClosestPoint().
+   */
+  virtual Vector getVectorToPoint(const Vector &point) const;
+
+  /**
+   * Test wether a point is contained in the space.
+   *
+   * returns true, when the point lies inside and false
+   * otherwise.
+   */
   virtual bool isContained(const Vector &point) const;
+
+  /**
+   * Tests if this space contains the center of the coordinate system.
+   */
   virtual bool hasZero() const;
 
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/atom.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -67,11 +67,13 @@
 atom *atom::GetTrueFather()
 {
-  atom *walker = this;
-  do {
-    if (walker == walker->father) // top most father is the one that points on itself
-      break;
-    walker = walker->father;
-  } while (walker != NULL);
-  return walker;
+  if(father == this){ // top most father is the one that points on itself
+    return this;
+  }
+  else if(!father) {
+    return 0;
+  }
+  else {
+    return father->GetTrueFather();
+  }
 };
 
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/molecule_geometry.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -16,4 +16,7 @@
 #include "molecule.hpp"
 #include "World.hpp"
+#include "Plane.hpp"
+#include <boost/foreach.hpp>
+
 
 /************************************* Functions for class molecule *********************************/
@@ -256,5 +259,9 @@
 void molecule::Mirror(const Vector *n)
 {
-  ActOnAllVectors( &Vector::Mirror, *n );
+  OBSERVE;
+  Plane p(*n,0);
+  BOOST_FOREACH( atom* iter, atoms ){
+    (*iter->node) = p.mirrorVector(*iter->node);
+  }
 };
 
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/tesselation.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -17,4 +17,5 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -441,34 +442,27 @@
 
   try {
-    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
-  }
-  catch (LinearDependenceException &excp) {
-    Log() << Verbose(1) << excp;
-    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
-    return false;
-  }
-
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
-
-  if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
-    return true;
-  }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
-    return true;
-  }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
-    return true;
-  }
-  // Calculate cross point between one baseline and the line from the third endpoint to intersection
-  int i = 0;
-  do {
-    try {
-      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
-                                                    *(endpoints[(i+1)%3]->node->node),
-                                                    *(endpoints[(i+2)%3]->node->node),
-                                                    *Intersection);
+    Line centerLine = makeLineThrough(*MolCenter, *x);
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
+
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
+
+    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
+      return true;
+    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
+      return true;
+    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
+      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
+      return true;
+    }
+    // Calculate cross point between one baseline and the line from the third endpoint to intersection
+    int i = 0;
+    do {
+      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
+      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
+      CrossPoint = line1.getIntersection(line2);
       helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
       CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
@@ -477,17 +471,14 @@
       if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
         DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        i=4;
-        break;
+        return false;
       }
       i++;
-    } catch (LinearDependenceException &excp){
-      break;
-    }
-  } while (i < 3);
-  if (i == 3) {
+    } while (i < 3);
     DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
     return true;
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);
+  }
+  catch (MathException &excp) {
+    Log() << Verbose(1) << excp;
+    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
     return false;
   }
@@ -516,7 +507,8 @@
   GetCenter(&Direction);
   try {
-    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
-  }
-  catch (LinearDependenceException &excp) {
+    Line l = makeLineThrough(*x, Direction);
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
+  }
+  catch (MathException &excp) {
     (*ClosestPoint) = (*x);
   }
@@ -541,5 +533,6 @@
     Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
+    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
     CrossDirection[i] = CrossPoint[i] - InPlane;
     CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/tesselationhelpers.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -14,4 +14,5 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -666,7 +667,7 @@
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
-                                                   *(Base->endpoints[1]->node->node),
-                                                     NewOffset, NewDirection);
+  Line line1 = makeLineThrough(*(Base->endpoints[0]->node->node),*(Base->endpoints[1]->node->node));
+  Line line2 = makeLineThrough(NewOffset, NewDirection);
+  *Intersection = line1.getIntersection(line2);
   Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);
Index: src/unittests/LineUnittest.cpp
===================================================================
--- src/unittests/LineUnittest.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
+++ src/unittests/LineUnittest.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -0,0 +1,353 @@
+/*
+ * LineUnittest.cpp
+ *
+ *  Created on: May 27, 2010
+ *      Author: crueger
+ */
+
+#include "LineUnittest.hpp"
+
+#include "vector.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <iostream>
+
+using namespace std;
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+CPPUNIT_TEST_SUITE_REGISTRATION( LineUnittest );
+
+void LineUnittest::setUp(){
+  // three lines along the axes
+  la1 = new Line(zeroVec,e1);
+  la2 = new Line(zeroVec,e2);
+  la3 = new Line(zeroVec,e3);
+
+  // the lines along the planes defined by two coordinate axes
+  lp1 = new Line(e1,e1-e2);
+  lp2 = new Line(e2,e2-e3);
+  lp3 = new Line(e3,e3-e1);
+}
+void LineUnittest::tearDown(){
+  delete la1;
+  delete la2;
+  delete la3;
+
+  delete lp1;
+  delete lp2;
+  delete lp3;
+}
+
+void LineUnittest::constructionErrorTest(){
+  // test some constructions
+
+  // direction+origin should never fail
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e1));
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e2));
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,e3));
+
+  // two points fails if both points are the same
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e1,e2));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e2,e3));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e3,e1));
+  // for zerovectors
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e1,zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e2,zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e3,zeroVec));
+  // now we pass two times the same point
+  CPPUNIT_ASSERT_THROW(makeLineThrough(zeroVec,zeroVec),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(e1,e1),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(e2,e2),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(e3,e3),LinearDependenceException);
+
+}
+
+bool testDirection(const Vector &dir1,const Vector &dir2){
+  return (dir1==dir2) || (dir1==-1*dir2);
+}
+
+void LineUnittest::constructionResultTest(){
+  // test all directions
+  CPPUNIT_ASSERT(testDirection(la1->getDirection(),e1));
+  CPPUNIT_ASSERT(testDirection(la2->getDirection(),e2));
+  CPPUNIT_ASSERT(testDirection(la3->getDirection(),e3));
+
+  // test origins
+  CPPUNIT_ASSERT_EQUAL(la1->getOrigin(),zeroVec);
+  CPPUNIT_ASSERT_EQUAL(la2->getOrigin(),zeroVec);
+  CPPUNIT_ASSERT_EQUAL(la2->getOrigin(),zeroVec);
+
+  // test if desired points are on the lines
+  CPPUNIT_ASSERT(la1->isContained(zeroVec));
+  CPPUNIT_ASSERT(la2->isContained(zeroVec));
+  CPPUNIT_ASSERT(la3->isContained(zeroVec));
+
+  CPPUNIT_ASSERT(la1->isContained(e1));
+  CPPUNIT_ASSERT(la2->isContained(e2));
+  CPPUNIT_ASSERT(la3->isContained(e3));
+
+  CPPUNIT_ASSERT(lp1->isContained(e1));
+  CPPUNIT_ASSERT(lp2->isContained(e2));
+  CPPUNIT_ASSERT(lp3->isContained(e3));
+
+  CPPUNIT_ASSERT(lp1->isContained(e2));
+  CPPUNIT_ASSERT(lp2->isContained(e3));
+  CPPUNIT_ASSERT(lp3->isContained(e1));
+}
+
+void LineUnittest::isContainedTest(){
+  // Zerovector on the axes lines
+  CPPUNIT_ASSERT(la1->isContained(zeroVec));
+  CPPUNIT_ASSERT(la2->isContained(zeroVec));
+  CPPUNIT_ASSERT(la3->isContained(zeroVec));
+
+  // multiples of the second support vector
+  CPPUNIT_ASSERT(la1->isContained(e1));
+  CPPUNIT_ASSERT(la2->isContained(e2));
+  CPPUNIT_ASSERT(la3->isContained(e3));
+
+  CPPUNIT_ASSERT(la1->isContained(2*e1));
+  CPPUNIT_ASSERT(la2->isContained(2*e2));
+  CPPUNIT_ASSERT(la3->isContained(2*e3));
+
+  CPPUNIT_ASSERT(la1->isContained(3*e1));
+  CPPUNIT_ASSERT(la2->isContained(3*e2));
+  CPPUNIT_ASSERT(la3->isContained(3*e3));
+
+  // negative multiples
+  CPPUNIT_ASSERT(la1->isContained(-1*e1));
+  CPPUNIT_ASSERT(la2->isContained(-1*e2));
+  CPPUNIT_ASSERT(la3->isContained(-1*e3));
+
+  CPPUNIT_ASSERT(la1->isContained(-2*e1));
+  CPPUNIT_ASSERT(la2->isContained(-2*e2));
+  CPPUNIT_ASSERT(la3->isContained(-2*e3));
+
+  // points that should not be on the lines
+  CPPUNIT_ASSERT(!la1->isContained(e2));
+  CPPUNIT_ASSERT(!la2->isContained(e3));
+  CPPUNIT_ASSERT(!la3->isContained(e1));
+
+  CPPUNIT_ASSERT(!la1->isContained(2*e2));
+  CPPUNIT_ASSERT(!la2->isContained(2*e3));
+  CPPUNIT_ASSERT(!la3->isContained(2*e1));
+
+  CPPUNIT_ASSERT(!la1->isContained(-1*e2));
+  CPPUNIT_ASSERT(!la2->isContained(-1*e3));
+  CPPUNIT_ASSERT(!la3->isContained(-1*e1));
+
+  // For the plane lines
+  CPPUNIT_ASSERT(lp1->isContained(e1));
+  CPPUNIT_ASSERT(lp2->isContained(e2));
+  CPPUNIT_ASSERT(lp3->isContained(e3));
+
+  CPPUNIT_ASSERT(lp1->isContained(e2));
+  CPPUNIT_ASSERT(lp2->isContained(e3));
+  CPPUNIT_ASSERT(lp3->isContained(e1));
+
+  CPPUNIT_ASSERT(lp1->isContained(e1+2*(e1-e2)));
+  CPPUNIT_ASSERT(lp2->isContained(e2+2*(e2-e3)));
+  CPPUNIT_ASSERT(lp3->isContained(e3+2*(e3-e1)));
+
+  CPPUNIT_ASSERT(lp1->isContained(e1-2*(e1-e2)));
+  CPPUNIT_ASSERT(lp2->isContained(e2-2*(e2-e3)));
+  CPPUNIT_ASSERT(lp3->isContained(e3-2*(e3-e1)));
+}
+
+void LineUnittest::intersectionTest(){
+  Vector fixture;
+
+  // intersection of the axis lines
+  fixture = la1->getIntersection(*la2);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la2->getIntersection(*la3);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la3->getIntersection(*la1);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+
+  // axes and plane lines
+  fixture = la1->getIntersection(*lp1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->getIntersection(*lp2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->getIntersection(*lp3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = la1->getIntersection(*lp3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->getIntersection(*lp1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->getIntersection(*lp2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // two plane lines
+  fixture = lp1->getIntersection(*lp2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = lp2->getIntersection(*lp3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = lp3->getIntersection(*lp1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+
+  // When we have two times the same line, we check if the point is on the line
+  fixture = la1->getIntersection(*la1);
+  CPPUNIT_ASSERT(la1->isContained(fixture));
+  fixture = la2->getIntersection(*la2);
+  CPPUNIT_ASSERT(la2->isContained(fixture));
+  fixture = la3->getIntersection(*la3);
+  CPPUNIT_ASSERT(la3->isContained(fixture));
+
+  fixture = lp1->getIntersection(*lp1);
+  CPPUNIT_ASSERT(lp1->isContained(fixture));
+  fixture = lp2->getIntersection(*lp2);
+  CPPUNIT_ASSERT(lp2->isContained(fixture));
+  fixture = lp3->getIntersection(*lp3);
+  CPPUNIT_ASSERT(lp3->isContained(fixture));
+
+  // lines that are askew should produce an Error
+  CPPUNIT_ASSERT_THROW(lp1->getIntersection(*la3),SkewException);
+  CPPUNIT_ASSERT_THROW(lp2->getIntersection(*la1),SkewException);
+  CPPUNIT_ASSERT_THROW(lp3->getIntersection(*la2),SkewException);
+
+  CPPUNIT_ASSERT_THROW(la1->getIntersection(*lp2),SkewException);
+  CPPUNIT_ASSERT_THROW(la2->getIntersection(*lp3),SkewException);
+  CPPUNIT_ASSERT_THROW(la3->getIntersection(*lp1),SkewException);
+}
+
+void LineUnittest::rotationTest(){
+  Vector fixture;
+
+  // rotate zero Vector along the axes lines by various degrees
+  fixture = la1->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la2->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la3->rotateVector(zeroVec,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+
+  fixture = la1->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la2->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = la3->rotateVector(zeroVec,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+
+  // rotate vectors on the axis around their lines
+  fixture = la1->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = la1->rotateVector(e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la2->rotateVector(e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la3->rotateVector(e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // more vectors on the axis
+  fixture = la1->rotateVector(2*e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
+  fixture = la2->rotateVector(2*e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
+  fixture = la3->rotateVector(2*e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
+
+  fixture = la1->rotateVector(2*e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e1);
+  fixture = la2->rotateVector(2*e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e2);
+  fixture = la3->rotateVector(2*e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*e3);
+
+  // negative factors
+  fixture = la1->rotateVector(-1*e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(-1*e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(-1*e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+
+  fixture = la1->rotateVector(-1*e1,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(-1*e2,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(-1*e3,2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+
+
+
+  // now the real rotations
+  // e2 around e1
+  fixture = la1->rotateVector(e2,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = la1->rotateVector(e2,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  fixture = la1->rotateVector(e2,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la1->rotateVector(e2,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+
+  // e3 around e2
+  fixture = la2->rotateVector(e3,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = la2->rotateVector(e3,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la2->rotateVector(e3,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  fixture = la2->rotateVector(e3,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // e1 around e3
+  fixture = la3->rotateVector(e1,0);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = la3->rotateVector(e1,1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = la3->rotateVector(e1,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+  fixture = la3->rotateVector(e1,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+
+
+  // and some rotation around the plane lines
+
+  // Vectors on the line
+  fixture = lp1->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = lp1->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+
+  fixture = lp2->rotateVector(e2,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = lp2->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = lp3->rotateVector(e3,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  fixture = lp3->rotateVector(e1,1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+
+  // the real stuff
+  fixture = lp1->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,1,0));
+  fixture = lp2->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(0,1,1));
+  fixture = lp3->rotateVector(zeroVec,M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,Vector(1,0,1));
+
+  fixture = lp1->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = lp2->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = lp3->rotateVector(zeroVec,2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+}
Index: src/unittests/LineUnittest.hpp
===================================================================
--- src/unittests/LineUnittest.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
+++ src/unittests/LineUnittest.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -0,0 +1,46 @@
+/*
+ * LineUnittest.hpp
+ *
+ *  Created on: May 27, 2010
+ *      Author: crueger
+ */
+
+#ifndef LINEUNITTEST_HPP_
+#define LINEUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Line.hpp"
+
+class LineUnittest : public CppUnit::TestFixture
+{
+
+  CPPUNIT_TEST_SUITE( LineUnittest) ;
+  CPPUNIT_TEST ( constructionErrorTest );
+  CPPUNIT_TEST ( constructionResultTest );
+  CPPUNIT_TEST ( isContainedTest );
+  CPPUNIT_TEST ( intersectionTest );
+  CPPUNIT_TEST ( rotationTest );
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+  void setUp();
+  void tearDown();
+
+  void constructionErrorTest();
+  void constructionResultTest();
+  void isContainedTest();
+  void intersectionTest();
+  void rotationTest();
+
+private:
+  Line *la1;
+  Line *la2;
+  Line *la3;
+  Line *lp1;
+  Line *lp2;
+  Line *lp3;
+
+};
+
+#endif /* LINEUNITTEST_HPP_ */
Index: src/unittests/Makefile.am
===================================================================
--- src/unittests/Makefile.am	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/unittests/Makefile.am	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -26,4 +26,5 @@
   InfoUnitTest \
   LinearSystemOfEquationsUnitTest \
+  LineUnittest \
   LinkedCellUnitTest \
   ListOfBondsUnitTest \
@@ -69,4 +70,5 @@
   infounittest.cpp \
   linearsystemofequationsunittest.cpp \
+  LineUnittest.cpp \
   LinkedCellUnitTest.cpp \
   listofbondsunittest.cpp \
@@ -104,4 +106,5 @@
   infounittest.hpp \
   linearsystemofequationsunittest.hpp \
+  LineUnittest.hpp \
   LinkedCellUnitTest.hpp \
   listofbondsunittest.hpp \
@@ -170,4 +173,7 @@
 LinearSystemOfEquationsUnitTest_LDADD = ${ALLLIBS}
 
+LineUnittest_SOURCES = UnitTestMain.cpp LineUnittest.cpp LineUnittest.hpp
+LineUnittest_LDADD = ${ALLLIBS}
+
 LinkedCellUnitTest_SOURCES = UnitTestMain.cpp LinkedCellUnitTest.cpp LinkedCellUnitTest.hpp 
 LinkedCellUnitTest_LDADD = ${ALLLIBS}
Index: src/unittests/PlaneUnittest.cpp
===================================================================
--- src/unittests/PlaneUnittest.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/unittests/PlaneUnittest.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -17,4 +17,5 @@
 
 #include "vector.hpp"
+#include "Line.hpp"
 
 CPPUNIT_TEST_SUITE_REGISTRATION( PlaneUnittest );
@@ -153,5 +154,63 @@
   CPPUNIT_ASSERT(fabs(p4->distance(e1)-1) < MYEPSILON);
   CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(e1));
-
-
-}
+}
+
+void PlaneUnittest::mirrorTest(){
+  Vector fixture;
+
+  // some Vectors that lie on the planes
+  fixture = p1->mirrorVector(e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = p1->mirrorVector(e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = p1->mirrorVector(e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = p2->mirrorVector(zeroVec);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = p2->mirrorVector(e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = p2->mirrorVector(e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+
+  fixture = p3->mirrorVector(zeroVec);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = p3->mirrorVector(e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  fixture = p3->mirrorVector(e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  fixture = p4->mirrorVector(zeroVec);
+  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+  fixture = p4->mirrorVector(e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  fixture = p4->mirrorVector(e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+
+  // some Vectors outside of the planes
+  {
+    Vector t = (2./3.)*(e1+e2+e3);
+    fixture = p1->mirrorVector(zeroVec);
+    CPPUNIT_ASSERT_EQUAL(fixture,t);
+  }
+
+  fixture = p2->mirrorVector(e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
+  fixture = p3->mirrorVector(e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
+  fixture = p4->mirrorVector(e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
+}
+
+void PlaneUnittest::LineIntersectionTest(){
+  Vector fixture;
+  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
+  Line l1 = makeLineThrough(zeroVec,Vector(2,1,0));
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e1, zeroVec).GetIntersection(l1) );
+  CPPUNIT_ASSERT_EQUAL( zeroVec, fixture );
+
+  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
+  Line l2 = makeLineThrough(e1,Vector(0,1,1));
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e2, Vector(2,1,0)).GetIntersection(l2) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
+}
Index: src/unittests/PlaneUnittest.hpp
===================================================================
--- src/unittests/PlaneUnittest.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/unittests/PlaneUnittest.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -20,4 +20,6 @@
   CPPUNIT_TEST ( pointsTest );
   CPPUNIT_TEST ( operationsTest );
+  CPPUNIT_TEST ( mirrorTest );
+  CPPUNIT_TEST ( LineIntersectionTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -30,4 +32,6 @@
   void pointsTest();
   void operationsTest();
+  void mirrorTest();
+  void LineIntersectionTest();
 
 private:
Index: src/unittests/vectorunittest.cpp
===================================================================
--- src/unittests/vectorunittest.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/unittests/vectorunittest.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -215,71 +215,4 @@
 }
 
-/** UnitTest for line intersections.
- */
-void VectorTest::LineIntersectionTest()
-{
-  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unit, zero).GetIntersection(zero, two) );
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(otherunit, two).GetIntersection( unit, notunit) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
-
-  // four vectors equal to zero
-  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(zero, zero, zero, zero), LinearDependenceException);
-  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // four vectors equal to unit
-  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, unit, unit, unit), LinearDependenceException);
-  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // two equal lines
-  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, two));
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // line from (1,0,0) to (2,1,0) cuts line from (1,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, otherunit) );
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // line from (1,0,0) to (0,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, zero, zero, two) );
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // line from (1,0,0) to (2,1,0) cuts line from (0,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, zero, otherunit) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., -1., 0.), fixture );
-};
-
-/** UnitTest for vector rotations.
- */
-void VectorTest::VectorRotationTest()
-{
-  fixture = Vector(-1.,0.,0.);
-
-  // zero vector does not change
-  fixture = RotateVector(zero,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  fixture = RotateVector(zero, two, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
-
-  // vector on axis does not change
-  fixture = RotateVector(unit,unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // rotations
-  fixture = RotateVector(otherunit, unit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( Vector(0.,-1.,0.), fixture );
-
-  fixture = RotateVector(otherunit, unit, 2. * M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(otherunit,unit, 0);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-
-  fixture = RotateVector(Vector(0.,0.,1.), notunit, M_PI);
-  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
-}
 
 /**
Index: src/unittests/vectorunittest.hpp
===================================================================
--- src/unittests/vectorunittest.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/unittests/vectorunittest.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -27,6 +27,4 @@
     CPPUNIT_TEST ( ProjectionTest );
     CPPUNIT_TEST ( NormalsTest );
-    CPPUNIT_TEST ( LineIntersectionTest );
-    CPPUNIT_TEST ( VectorRotationTest );
     CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/vector.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -213,7 +213,7 @@
 {
   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]);
+  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;
 };
@@ -232,30 +232,4 @@
   *this -= tmp;
 };
-
-/** Calculates the minimum distance vector 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
- * \return distance vector onto to plane
- */
-Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
-{
-  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.;
-
-  temp.Normalize();
-  temp.Scale(sign);
-  return temp;
-};
-
 
 /** Calculates the minimum distance of this vector to the plane.
@@ -551,4 +525,21 @@
   MatrixMultiplication(M);
 };
+
+std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
+  double factor = ScalarProduct(rhs)/rhs.NormSquared();
+  Vector res= factor * rhs;
+  return make_pair(res,(*this)-res);
+}
+
+std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
+  Vector helper = *this;
+  pointset res;
+  for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
+    pair<Vector,Vector> currPart = helper.partition(*iter);
+    res.push_back(currPart.first);
+    helper = currPart.second;
+  }
+  return make_pair(res,helper);
+}
 
 /** Do a matrix multiplication.
@@ -611,16 +602,4 @@
 };
 
-/** Mirrors atom against a given plane.
- * \param n[] normal vector of mirror plane.
- */
-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--;)
-    at(i) -= 2.*projection*n[i];
-};
-
 /** Calculates orthonormal vector to one given vectors.
  * Just subtracts the projection onto the given vector from this vector.
@@ -633,6 +612,5 @@
   bool result = false;
   double factor = y1.ScalarProduct(*this)/y1.NormSquared();
-  Vector x1;
-  x1 = factor * y1;
+  Vector x1 = factor * y1;
   SubtractVector(x1);
   for (int i=NDIM;i--;)
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/vector.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -16,4 +16,5 @@
 
 #include <memory>
+#include <vector>
 
 #include "defs.hpp"
@@ -21,4 +22,8 @@
 
 /********************************************** declarations *******************************/
+
+class Vector;
+
+typedef std::vector<Vector> pointset;
 
 /** Single vector.
@@ -39,5 +44,4 @@
 
   double DistanceSquared(const Vector &y) const;
-  Vector GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
   double DistanceToSpace(const Space& space) const;
   double PeriodicDistance(const Vector &y, const double * const cell_size) const;
@@ -56,5 +60,4 @@
   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);
@@ -66,4 +69,6 @@
   bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   void WrapPeriodically(const double * const M, const double * const Minv);
+  std::pair<Vector,Vector> partition(const Vector&) const;
+  std::pair<pointset,Vector> partition(const pointset&) const;
 
   // Accessors ussually come in pairs... and sometimes even more than that
Index: src/vector_ops.cpp
===================================================================
--- src/vector_ops.cpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/vector_ops.cpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -15,4 +15,5 @@
 #include "Helpers/fast_functions.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
+#include "Exceptions/SkewException.hpp"
 
 #include <gsl/gsl_linalg.h>
@@ -110,132 +111,2 @@
   return true;
 };
-
-/** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
- * \param *axis rotation axis
- * \param alpha rotation angle in radian
- */
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha)
-{
-  Vector a,y;
-  Vector res;
-  // normalise this vector with respect to axis
-  a = vec;
-  a.ProjectOntoPlane(axis);
-  // construct normal vector
-  try {
-    y = Plane(axis,a,0).getNormal();
-  }
-  catch (MathException &excp) {
-    // The normal vector cannot be created if there is linar dependency.
-    // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
-    return vec;
-  }
-  y.Scale(vec.Norm());
-  // scale normal vector by sine and this vector by cosine
-  y.Scale(sin(alpha));
-  a.Scale(cos(alpha));
-  res = vec.Projection(axis);
-  // add scaled normal vector onto this vector
-  res += y;
-  // add part in axis direction
-  res += a;
-  return res;
-};
-
-/** Calculates the intersection of the two lines that are both on the same plane.
- * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
- * \param *out output stream for debugging
- * \param *Line1a first vector of first line
- * \param *Line1b second vector of first line
- * \param *Line2a first vector of second line
- * \param *Line2b second vector of second line
- * \return true - \a this will contain the intersection on return, false - lines are parallel
- */
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b)
-{
-  Info FunctionInfo(__func__);
-
-  Vector res;
-
-  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
-
-  M->SetAll(1.);
-  for (int i=0;i<3;i++) {
-    M->Set(0, i, Line1a[i]);
-    M->Set(1, i, Line1b[i]);
-    M->Set(2, i, Line2a[i]);
-    M->Set(3, i, Line2b[i]);
-  }
-
-  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    cout << "\t" << M->Get(i,j);
-  //  cout << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
-
-
-  // constuct a,b,c
-  Vector a = Line1b - Line1a;
-  Vector b = Line2b - Line2a;
-  Vector c = Line2a - Line1a;
-  Vector d = Line2b - Line1b;
-  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   res.Zero();
-   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
-   throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // check for parallelity
-  Vector parallel;
-  double factor = 0.;
-  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
-    parallel = Line1a - Line2a;
-    factor = parallel.ScalarProduct(a)/a.Norm();
-    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-      res = Line2a;
-      Log() << Verbose(1) << "Lines conincide." << endl;
-      return res;
-    } else {
-      parallel = Line1a - Line2b;
-      factor = parallel.ScalarProduct(a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        res = Line2b;
-        Log() << Verbose(1) << "Lines conincide." << endl;
-        return res;
-      }
-    }
-    Log() << Verbose(1) << "Lines are parallel." << endl;
-    res.Zero();
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // obtain s
-  double s;
-  Vector temp1, temp2;
-  temp1 = c;
-  temp1.VectorProduct(b);
-  temp2 = a;
-  temp2.VectorProduct(b);
-  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
-
-  // construct intersection
-  res = a;
-  res.Scale(s);
-  res += Line1a;
-  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
-
-  return res;
-};
Index: src/vector_ops.hpp
===================================================================
--- src/vector_ops.hpp	(revision e0582670720906c59f098684b31f4fd9f892fc21)
+++ src/vector_ops.hpp	(revision 76c0d685b7570bb1f520bbd307e7486df6502bb4)
@@ -10,6 +10,4 @@
 
 bool LSQdistance(Vector &res,const Vector **vectors, int num);
-Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha);
-Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b);
 
 #endif /* VECTOR_OPS_HPP_ */
