Index: src/Plane.cpp
===================================================================
--- src/Plane.cpp	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/Plane.cpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -44,5 +44,5 @@
  * If no offset is given a plane through origin is assumed
  */
-Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(LinearDependenceException) :
+Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException) :
     normalVector(new Vector()),
     offset(_offset)
@@ -50,5 +50,9 @@
   Vector x1 = y1;
   Vector x2 = y2;
-  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(x2)) < MYEPSILON)) {
+  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON)) {
+    throw ZeroVectorException(__FILE__,__LINE__);
+  }
+
+  if((fabs(x1.Angle(x2)) < MYEPSILON)) {
     throw LinearDependenceException(__FILE__,__LINE__);
   }
@@ -92,9 +96,9 @@
 
 
-Vector Plane::getNormal(){
+Vector Plane::getNormal() const{
   return *normalVector;
 }
 
-double Plane::getOffset(){
+double Plane::getOffset() const{
   return offset;
 }
@@ -106,15 +110,15 @@
 vector<Vector> Plane::getPointsOnPlane(){
   std::vector<Vector> res;
+  res.reserve(3);
   // first point on the plane
-  res[0] = getOffsetVector();
-  // first is orthogonal to the plane...
-  // an orthogonal vector to this one lies on the plane
+  res.push_back(getOffsetVector());
+  // get a vector that has direction of plane
   Vector direction;
-  direction.GetOneNormalVector(res[0]);
-  res[1] = res[0]+direction;
-  // get an orthogonal vector to direction and offset (lies on the plane)
-  direction.VectorProduct(res[0]);
+  direction.GetOneNormalVector(getNormal());
+  res.push_back(res[0]+direction);
+  // get an orthogonal vector to direction and normal (has direction of plane)
+  direction.VectorProduct(getNormal());
   direction.Normalize();
-  res[2] = res[0] +direction;
+  res.push_back(res[0] +direction);
   return res;
 }
@@ -178,13 +182,17 @@
 
 Vector Plane::getClosestPoint(const Vector &point) const{
-  Vector difference = distance(point) * (*normalVector);
-  if(difference.IsZero()){
+  double factor = point.ScalarProduct(*normalVector)-offset;
+  if(fabs(factor) < MYEPSILON){
     // the point itself lies on the plane
     return point;
   }
-  // get the direction this vector is pointing
-  double sign = difference.ScalarProduct(*normalVector);
-  // sign cannot be zero, since normalVector and difference are both != zero
-  sign = sign/fabs(sign);
-  return (point - (sign * difference));
+  Vector difference = factor * (*normalVector);
+  return (point - difference);
 }
+
+// Operators
+
+ostream &operator << (ostream &ost,const Plane &p){
+  ost << "<" << p.getNormal() << ";x> - " << p.getOffset() << "=0";
+  return ost;
+}
Index: src/Plane.hpp
===================================================================
--- src/Plane.hpp	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/Plane.hpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -11,4 +11,5 @@
 #include <memory>
 #include <vector>
+#include <iostream>
 #include "Space.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
@@ -22,6 +23,6 @@
 public:
   Plane(const Vector &y1, const Vector &y2, const Vector &y3) throw(LinearDependenceException);
-  Plane(const Vector &y1, const Vector &y2, double _offset) throw(LinearDependenceException);
-  Plane(const Vector &_normalVector, double _offset=0)  throw(ZeroVectorException);
+  Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException);
+  Plane(const Vector &_normalVector, double _offset)  throw(ZeroVectorException);
   Plane(const Vector &_normalVector, const Vector &_offsetVector) throw(ZeroVectorException);
   virtual ~Plane();
@@ -31,9 +32,9 @@
    * returns normal Vector for a plane
    */
-  Vector getNormal();
+  Vector getNormal() const;
   /**
    * returns the distance of the plane from the origin
    */
-  double getOffset();
+  double getOffset() const;
   /**
    * returns a vector that points on the plane.
@@ -60,3 +61,5 @@
 };
 
+std::ostream &operator<< (std::ostream &ost,const Plane& p);
+
 #endif /* PLANE_HPP_ */
Index: src/unittests/Makefile.am
===================================================================
--- src/unittests/Makefile.am	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/unittests/Makefile.am	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -30,4 +30,5 @@
   MemoryAllocatorUnitTest \
   MoleculeDescriptorTest \
+  PlaneUnittest \
   ObserverTest \
   SingletonTest \
@@ -70,4 +71,5 @@
   memoryusageobserverunittest.cpp \
   MoleculeDescriptorTest.cpp \
+  PlaneUnittest.cpp \
   ObserverTest.cpp \
   SingletonTest.cpp \
@@ -102,4 +104,5 @@
   memoryusageobserverunittest.hpp \
   MoleculeDescriptorTest.hpp \
+  PlaneUnittest.hpp \
   ObserverTest.hpp \
   SingletonTest.hpp \
@@ -168,4 +171,7 @@
 MoleculeDescriptorTest_LDADD = ${ALLLIBS}
 
+PlaneUnittest_SOURCES = UnitTestMain.cpp PlaneUnittest.cpp PlaneUnittest.hpp
+PlaneUnittest_LDADD = ${ALLLIBS}
+
 SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
 SingletonTest_LDADD = $(BOOST_LIB) ${BOOST_THREAD_LIB}
Index: src/unittests/PlaneUnittest.cpp
===================================================================
--- src/unittests/PlaneUnittest.cpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
+++ src/unittests/PlaneUnittest.cpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -0,0 +1,157 @@
+/*
+ * PlaneUnittest.cpp
+ *
+ *  Created on: Apr 30, 2010
+ *      Author: crueger
+ */
+
+#include "PlaneUnittest.hpp"
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#ifdef HAVE_TESTRUNNER
+#include "UnitTestMain.hpp"
+#endif /*HAVE_TESTRUNNER*/
+
+#include "vector.hpp"
+
+CPPUNIT_TEST_SUITE_REGISTRATION( PlaneUnittest );
+
+void PlaneUnittest::setUp(){
+  p1 = new Plane(e1,e2,e3);
+  p2 = new Plane(e1,e2,zeroVec);
+  p3 = new Plane(e1,zeroVec,e3);
+  p4 = new Plane(zeroVec,e2,e3);
+}
+
+void PlaneUnittest::tearDown(){
+  delete p1;
+  delete p2;
+  delete p3;
+  delete p4;
+}
+
+void PlaneUnittest::constructionErrorTest(){
+  // try several method of construction..
+  // see if error checking works
+
+  // three points
+  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,e3));
+  // when only two points are differnt this gives an error
+  CPPUNIT_ASSERT_THROW(Plane(e1,e2,e2),LinearDependenceException);
+  // same with only one point
+  CPPUNIT_ASSERT_THROW(Plane(e1,e1,e1),LinearDependenceException);
+
+  // use two vector giving two directions
+  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,0));
+  // and again this is actually only one vector
+  CPPUNIT_ASSERT_THROW(Plane(e1,e1,0),LinearDependenceException);
+  // Zero vector does not give a good direction
+  CPPUNIT_ASSERT_THROW(Plane(e1,zeroVec,0),ZeroVectorException);
+
+  // use a normalvector and an scalar offset
+  CPPUNIT_ASSERT_NO_THROW(Plane(e1,0));
+  // The zero vector is no good as a normalvector
+  CPPUNIT_ASSERT_THROW(Plane(zeroVec,0),ZeroVectorException);
+
+  // use a normalvector and an offset vector
+  CPPUNIT_ASSERT_NO_THROW(Plane(e1,zeroVec));
+  // and the bad zeroVector again
+  CPPUNIT_ASSERT_THROW(Plane(zeroVec,zeroVec),ZeroVectorException);
+}
+
+
+// we need to test normals independent of the direction
+bool testNormal(const Vector &normal1, const Vector &normal2){
+  return (normal1==normal2) || (normal1==-1*normal2);
+}
+
+void PlaneUnittest::constructionResultTest(){
+  {
+    // construct with three points on plane
+    Plane p1(e1,e2,zeroVec);
+    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
+
+    Plane p2(e1,e3,zeroVec);
+    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
+
+    Plane p3(e2,e3,zeroVec);
+    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
+  }
+  {
+    // construct with two directions + offset
+    Plane p1(e1,e2,0);
+    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
+
+    Plane p2(e1,e3,0);
+    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
+
+    Plane p3(e2,e3,0);
+    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
+  }
+}
+
+void PlaneUnittest::pointsTest(){
+  std::vector<Vector> points1 = p1->getPointsOnPlane();
+  CPPUNIT_ASSERT(p1->isContained(points1[0]));
+  CPPUNIT_ASSERT(p1->isContained(points1[1]));
+  CPPUNIT_ASSERT(p1->isContained(points1[2]));
+  // check that the three points differ
+  CPPUNIT_ASSERT(points1[0]!=points1[1]);
+  CPPUNIT_ASSERT(points1[0]!=points1[2]);
+  CPPUNIT_ASSERT(points1[1]!=points1[2]);
+
+
+  std::vector<Vector> points2 = p2->getPointsOnPlane();
+  CPPUNIT_ASSERT(p2->isContained(points2[0]));
+  CPPUNIT_ASSERT(p2->isContained(points2[1]));
+  CPPUNIT_ASSERT(p2->isContained(points2[2]));
+  // check that the three points differ
+  CPPUNIT_ASSERT(points2[0]!=points2[1]);
+  CPPUNIT_ASSERT(points2[0]!=points2[2]);
+  CPPUNIT_ASSERT(points2[1]!=points2[2]);
+
+  std::vector<Vector> points3 = p3->getPointsOnPlane();
+  CPPUNIT_ASSERT(p3->isContained(points3[0]));
+  CPPUNIT_ASSERT(p3->isContained(points3[1]));
+  CPPUNIT_ASSERT(p3->isContained(points3[2]));
+  // check that the three points differ
+  CPPUNIT_ASSERT(points3[0]!=points3[1]);
+  CPPUNIT_ASSERT(points3[0]!=points3[2]);
+  CPPUNIT_ASSERT(points3[1]!=points3[2]);
+
+  std::vector<Vector> points4 = p4->getPointsOnPlane();
+  CPPUNIT_ASSERT(p4->isContained(points4[0]));
+  CPPUNIT_ASSERT(p4->isContained(points4[1]));
+  CPPUNIT_ASSERT(p4->isContained(points4[2]));
+  // check that the three points differ
+  CPPUNIT_ASSERT(points4[0]!=points4[1]);
+  CPPUNIT_ASSERT(points4[0]!=points4[2]);
+  CPPUNIT_ASSERT(points4[1]!=points4[2]);
+}
+
+
+void PlaneUnittest::operationsTest(){
+  {
+    Vector t = (1./3.)*(e1+e2+e3);
+    CPPUNIT_ASSERT(fabs(p1->distance(zeroVec)-t.Norm()) < MYEPSILON);
+    CPPUNIT_ASSERT_EQUAL(t,p1->getClosestPoint(zeroVec));
+  }
+
+  CPPUNIT_ASSERT(fabs(p2->distance(e3)-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p2->getClosestPoint(e3));
+  CPPUNIT_ASSERT(fabs(p3->distance(e2)-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p3->getClosestPoint(e2));
+  CPPUNIT_ASSERT(fabs(p4->distance(e1)-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(e1));
+
+
+}
Index: src/unittests/PlaneUnittest.hpp
===================================================================
--- src/unittests/PlaneUnittest.hpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
+++ src/unittests/PlaneUnittest.hpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -0,0 +1,40 @@
+/*
+ * PlaneUnittest.hpp
+ *
+ *  Created on: Apr 30, 2010
+ *      Author: crueger
+ */
+
+#ifndef PLANEUNITTEST_HPP_
+#define PLANEUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "Plane.hpp"
+
+class PlaneUnittest : public CppUnit::TestFixture
+{
+  CPPUNIT_TEST_SUITE( PlaneUnittest) ;
+  CPPUNIT_TEST ( constructionErrorTest );
+  CPPUNIT_TEST ( constructionResultTest );
+  CPPUNIT_TEST ( pointsTest );
+  CPPUNIT_TEST ( operationsTest );
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+  void setUp();
+  void tearDown();
+
+  void constructionErrorTest();
+  void constructionResultTest();
+  void pointsTest();
+  void operationsTest();
+
+private:
+  Plane *p1;
+  Plane *p2;
+  Plane *p3;
+  Plane *p4;
+};
+
+#endif /* PLANEUNITTEST_HPP_ */
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/vector.cpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -421,4 +421,9 @@
 };
 
+bool Vector::operator!=(const Vector& b) const
+{
+  return !IsEqualTo(b);
+}
+
 /** Sums vector \a to this lhs component-wise.
  * \param a base vector
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/vector.hpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -89,4 +89,5 @@
   // operators for mathematical operations
   bool operator==(const Vector& b) const;
+  bool operator!=(const Vector& b) const;
   const Vector& operator+=(const Vector& b);
   const Vector& operator-=(const Vector& b);
Index: src/vector_ops.cpp
===================================================================
--- src/vector_ops.cpp	(revision 1513a7402fe19aa16ecbdf33eddd7a4e2e403ff8)
+++ src/vector_ops.cpp	(revision fa5a6afe5abec6b36502dfd0ee72e520831287b7)
@@ -126,5 +126,5 @@
     y = Plane(axis,a,0).getNormal();
   }
-  catch (LinearDependenceException &excp) {
+  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.
