Index: src/Box.cpp
===================================================================
--- src/Box.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Box.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -14,4 +14,6 @@
 #include "Matrix.hpp"
 #include "vector.hpp"
+#include "Shapes/BaseShapes.hpp"
+#include "Shapes/ShapeOps.hpp"
 
 #include "Helpers/Assert.hpp"
@@ -139,4 +141,9 @@
 }
 
+Shape Box::getShape() const{
+  return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
+
+}
+
 Box &Box::operator=(const Box &src){
   if(&src!=this){
Index: src/Box.hpp
===================================================================
--- src/Box.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Box.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -11,4 +11,5 @@
 class Matrix;
 class Vector;
+class Shape;
 
 #include <list>
@@ -89,4 +90,6 @@
   double periodicDistance(const Vector &point1,const Vector &point2) const;
 
+  Shape getShape() const;
+
 private:
   Matrix *M;    //!< Defines the layout of the box
Index: src/Descriptors/AtomShapeDescriptor.cpp
===================================================================
--- src/Descriptors/AtomShapeDescriptor.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Descriptors/AtomShapeDescriptor.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,26 @@
+/*
+ * AtomShapeDescriptor.cpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#include "AtomShapeDescriptor.hpp"
+#include "AtomShapeDescriptor_impl.hpp"
+
+#include "atom.hpp"
+
+AtomShapeDescriptor_impl::AtomShapeDescriptor_impl(const Shape &_shape) :
+  shape(_shape)
+{}
+
+AtomShapeDescriptor_impl::~AtomShapeDescriptor_impl(){}
+
+bool AtomShapeDescriptor_impl::predicate(std::pair<atomId_t,atom*> atom){
+  return shape.isInside(*atom.second->node);
+}
+
+
+AtomDescriptor AtomByShape(const Shape &shape){
+  return AtomDescriptor(AtomDescriptor::impl_ptr(new AtomShapeDescriptor_impl(shape)));
+}
Index: src/Descriptors/AtomShapeDescriptor.hpp
===================================================================
--- src/Descriptors/AtomShapeDescriptor.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Descriptors/AtomShapeDescriptor.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,17 @@
+/*
+ * AtomShapeDescriptor.hpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#ifndef ATOMSHAPEDESCRIPTOR_HPP_
+#define ATOMSHAPEDESCRIPTOR_HPP_
+
+#include "Descriptors/AtomDescriptor.hpp"
+
+class Shape;
+
+AtomDescriptor AtomByShape(const Shape &shape);
+
+#endif /* ATOMSHAPEDESCRIPTOR_HPP_ */
Index: src/Descriptors/AtomShapeDescriptor_impl.hpp
===================================================================
--- src/Descriptors/AtomShapeDescriptor_impl.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Descriptors/AtomShapeDescriptor_impl.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,26 @@
+/*
+ * AtomShapeDescriptor_impl.hpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#ifndef ATOMSHAPEDESCRIPTOR_IMPL_HPP_
+#define ATOMSHAPEDESCRIPTOR_IMPL_HPP_
+
+#include "Descriptors/AtomDescriptor_impl.hpp"
+
+#include "Shapes/Shape.hpp"
+
+class AtomShapeDescriptor_impl : public AtomDescriptor_impl
+{
+public:
+  AtomShapeDescriptor_impl(const Shape&);
+  virtual  ~AtomShapeDescriptor_impl();
+
+  bool predicate(std::pair<atomId_t,atom*> atom);
+private:
+  Shape shape;
+};
+
+#endif /* ATOMSHAPEDESCRIPTOR_IMPL_HPP_ */
Index: src/Exceptions/NotOnSurfaceException.cpp
===================================================================
--- src/Exceptions/NotOnSurfaceException.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Exceptions/NotOnSurfaceException.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,21 @@
+/*
+ * NotOnSurfaceException.cpp
+ *
+ *  Created on: Jul 14, 2010
+ *      Author: crueger
+ */
+
+#include "NotOnSurfaceException.hpp"
+
+using namespace std;
+
+NotOnSurfaceException::NotOnSurfaceException(string _file, int _line) :
+  ShapeException(_file,_line)
+{}
+
+NotOnSurfaceException::~NotOnSurfaceException() throw ()
+{}
+
+const char* NotOnSurfaceException::what() const throw(){
+  return "Given point did not lie on surface of Shape";
+}
Index: src/Exceptions/NotOnSurfaceException.hpp
===================================================================
--- src/Exceptions/NotOnSurfaceException.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Exceptions/NotOnSurfaceException.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,22 @@
+/*
+ * NotOnSurfaceException.hpp
+ *
+ *  Created on: Jul 14, 2010
+ *      Author: crueger
+ */
+
+#ifndef NOTONSURFACEEXCEPTION_HPP_
+#define NOTONSURFACEEXCEPTION_HPP_
+
+#include "Exceptions/ShapeException.hpp"
+
+class NotOnSurfaceException : public ShapeException
+{
+public:
+  NotOnSurfaceException(std::string _file, int _line);
+  virtual ~NotOnSurfaceException() throw();
+
+  virtual const char* what() const throw();
+};
+
+#endif /* NOTONSURFACEEXCEPTION_HPP_ */
Index: src/Exceptions/ShapeException.cpp
===================================================================
--- src/Exceptions/ShapeException.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Exceptions/ShapeException.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,15 @@
+/*
+ * ShapeException.cpp
+ *
+ *  Created on: Jul 14, 2010
+ *      Author: crueger
+ */
+
+#include "ShapeException.hpp"
+
+ShapeException::ShapeException(std::string file, int line) :
+  CustomException(file,line)
+{}
+
+ShapeException::~ShapeException() throw ()
+{}
Index: src/Exceptions/ShapeException.hpp
===================================================================
--- src/Exceptions/ShapeException.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/Exceptions/ShapeException.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,20 @@
+/*
+ * ShapeException.hpp
+ *
+ *  Created on: Jul 14, 2010
+ *      Author: crueger
+ */
+
+#ifndef SHAPEEXCEPTION_HPP_
+#define SHAPEEXCEPTION_HPP_
+
+#include "Exceptions/CustomException.hpp"
+
+class ShapeException : public CustomException
+{
+public:
+  ShapeException(std::string file, int line);
+  virtual ~ShapeException() throw();
+};
+
+#endif /* SHAPEEXCEPTION_HPP_ */
Index: src/Line.cpp
===================================================================
--- src/Line.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Line.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -37,4 +37,12 @@
 Line::~Line()
 {}
+
+Line &Line::operator=(const Line& rhs){
+  if(this!=&rhs){
+    origin.reset(new Vector(*rhs.origin));
+    direction.reset(new Vector(*rhs.direction));
+  }
+  return *this;
+}
 
 
@@ -211,4 +219,31 @@
 }
 
+Line Line::rotateLine(const Line &rhs, double alpha) const{
+  Vector lineOrigin = rotateVector(rhs.getOrigin(),alpha);
+  Vector helper = rhs.getDirection();
+  // rotate the direction without considering the ofset
+  pair<Vector,Vector> parts = helper.partition(*direction);
+  Vector lineDirection = parts.first;
+  Vector a = parts.second;
+  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);
+
+    lineDirection += cos(alpha) * a + sin(alpha) * y;
+  }
+  return Line(lineOrigin,lineDirection);
+}
+
+Plane Line::rotatePlane(const Plane &rhs, double alpha) const{
+  vector<Vector> points = rhs.getPointsOnPlane();
+  transform(points.begin(),
+            points.end(),
+            points.begin(),
+            boost::bind(&Line::rotateVector,this,_1,alpha));
+  return Plane(points[0],points[1],points[2]);
+}
+
 Plane Line::getOrthogonalPlane(const Vector &origin) const{
   return Plane(getDirection(),origin);
@@ -234,4 +269,22 @@
 }
 
+LinePoint Line::getLinePoint(const Vector &point) const{
+  ASSERT(isContained(point),"Line point queried for point not on line");
+  Vector helper = point - (*origin);
+  double param = helper.ScalarProduct(*direction);
+  return LinePoint(*this,param);
+}
+
+LinePoint Line::posEndpoint() const{
+  return LinePoint(*this, numeric_limits<double>::infinity());
+}
+LinePoint Line::negEndpoint() const{
+  return LinePoint(*this,-numeric_limits<double>::infinity());
+}
+
+bool operator==(const Line &x,const Line &y){
+  return *x.origin == *y.origin && *x.direction == *y.direction;
+}
+
 Line makeLineThrough(const Vector &x1, const Vector &x2){
   if(x1==x2){
@@ -240,2 +293,48 @@
   return Line(x1,x1-x2);
 }
+
+
+/******************************** Points on the line ********************/
+
+LinePoint::LinePoint(const LinePoint &src) :
+  line(src.line),param(src.param)
+{}
+
+LinePoint::LinePoint(const Line &_line, double _param) :
+  line(_line),param(_param)
+{}
+
+LinePoint& LinePoint::operator=(const LinePoint &src){
+  line=src.line;
+  param=src.param;
+  return *this;
+}
+
+Vector LinePoint::getPoint() const{
+  ASSERT(!isInfinite(),"getPoint() on infinite LinePoint called");
+  return (*line.origin)+param*(*line.direction);
+}
+
+Line LinePoint::getLine() const{
+  return line;
+}
+
+bool LinePoint::isInfinite() const{
+  return isPosInfinity() || isNegInfinity();
+}
+bool LinePoint::isPosInfinity() const{
+  return param == numeric_limits<double>::infinity();
+}
+bool LinePoint::isNegInfinity() const{
+  return param ==-numeric_limits<double>::infinity();
+}
+
+bool operator==(const LinePoint &x, const LinePoint &y){
+  ASSERT(x.line==y.line,"Operation on two points of different lines");
+  return x.param == y.param;
+
+}
+bool operator<(const LinePoint &x, const LinePoint &y){
+  ASSERT(x.line==y.line,"Operation on two points of different lines");
+  return x.param<y.param;
+}
Index: src/Line.hpp
===================================================================
--- src/Line.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Line.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -16,11 +16,16 @@
 class Vector;
 class Plane;
+class LinePoint;
 
 class Line : public Space
 {
+  friend bool operator==(const Line&,const Line&);
+  friend class LinePoint;
 public:
   Line(const Vector &_origin, const Vector &_direction);
   Line(const Line& _src);
   virtual ~Line();
+
+  Line &operator=(const Line& rhs);
 
   virtual double distance(const Vector &point) const;
@@ -35,8 +40,16 @@
 
   Vector rotateVector(const Vector &rhs, double alpha) const;
+  Line rotateLine(const Line &rhs, double alpha) const;
+  Plane rotatePlane(const Plane &rhs, double alpha) const;
 
   Plane getOrthogonalPlane(const Vector &origin) const;
 
   std::vector<Vector> getSphereIntersections() const;
+
+  LinePoint getLinePoint(const Vector&) const;
+  LinePoint posEndpoint() const;
+  LinePoint negEndpoint() const;
+
+
 
 private:
@@ -45,4 +58,6 @@
 };
 
+bool operator==(const Line&,const Line&);
+
 /**
  * Named constructor to make a line through two points
@@ -50,3 +65,36 @@
 Line makeLineThrough(const Vector &x1, const Vector &x2);
 
+/**
+ * Class for representing points on a line
+ * These objects allow comparison of points on the same line as well as specifying the
+ * infinite "endpoints" of a line.
+ */
+class LinePoint{
+  friend class Line;
+  friend bool operator==(const LinePoint&, const LinePoint&);
+  friend bool operator<(const LinePoint&, const LinePoint&);
+public:
+  LinePoint(const LinePoint&);
+  LinePoint& operator=(const LinePoint&);
+  Vector getPoint() const;
+  Line getLine() const;
+  bool isInfinite() const;
+  bool isPosInfinity() const;
+  bool isNegInfinity() const;
+
+private:
+  LinePoint(const Line&,double);
+  Line line;
+  double param;
+};
+
+bool operator==(const LinePoint&, const LinePoint&);
+bool operator<(const LinePoint&, const LinePoint&);
+
+inline bool operator!= (const LinePoint& x, const LinePoint& y) { return !(x==y); }
+inline bool operator>  (const LinePoint& x, const LinePoint& y) { return y<x; }
+inline bool operator<= (const LinePoint& x, const LinePoint& y) { return !(y<x); }
+inline bool operator>= (const LinePoint& x, const LinePoint& y) { return !(x<y); }
+
+
 #endif /* LINE_HPP_ */
Index: src/LineSegment.cpp
===================================================================
--- src/LineSegment.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/LineSegment.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,90 @@
+/*
+ * LineSegment.cpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#include "LineSegment.hpp"
+#include "Helpers/Assert.hpp"
+
+#include "Line.hpp"
+#include "vector.hpp"
+
+using namespace std;
+
+LineSegment::LineSegment(const Vector &point1,const Vector &point2){
+  line.reset(new Line(point1,point1-point2));
+  begin.reset(new LinePoint(line->getLinePoint(point1)));
+  end.reset(new LinePoint(line->getLinePoint(point2)));
+}
+
+LineSegment::LineSegment(const LinePoint &lp1,const LinePoint &lp2){
+  ASSERT(lp1.getLine()==lp2.getLine(),"Cannot make a lineSegment with two different lines");
+  line.reset(new Line(lp1.getLine()));
+  begin.reset(new LinePoint(lp1));
+  end.reset(new LinePoint(lp2));
+}
+LineSegment::LineSegment(const LineSegment &src)
+{
+  line.reset(new Line(*src.line));
+  begin.reset(new LinePoint(*src.begin));
+  end.reset(new LinePoint(*src.end));
+}
+LineSegment::~LineSegment(){
+}
+LineSegment &LineSegment::operator=(const LineSegment &src){
+  if(this!=&src){
+    lineptr _line(new Line(*src.line));
+    pointptr _begin(new LinePoint(*src.begin));
+    pointptr _end(new LinePoint(*src.end));
+    line = _line;
+    begin = _begin;
+    end = _end;
+  }
+  return *this;
+}
+
+bool LineSegment::isContained(const Vector &point) const{
+  if(!line->isContained(point)){
+    return false;
+  }
+  LinePoint lp = line->getLinePoint(point);
+  return *begin <= lp && lp <= *end;
+}
+
+bool LineSegment::isContained(const LinePoint &lp) const{
+  ASSERT(lp.getLine()==*line,"Cannot compare Linepoints from different Lines");
+  return *begin <= lp && lp <= *end;
+}
+
+bool LineSegment::overlaps(const LineSegment &lineSeg) const{
+  ASSERT(*lineSeg.line == *line,"Cannot compare Linesegments from different Lines");
+  // mutual comparison of endpoints
+  return lineSeg.isContained(*end) || isContained(*lineSeg.end);
+}
+
+bool LineSegment::hasZeroLength() const{
+  return *begin==*end;
+}
+
+LinePoint LineSegment::getBegin() const{
+  return *begin;
+}
+
+LinePoint LineSegment::getEnd() const{
+  return *end;
+}
+
+Line LineSegment::getLine() const{
+  return *line;
+}
+
+bool operator==(const LineSegment &x,const LineSegment &y){
+  ASSERT(*x.line==*y.line,"Comparison of Linesegments from different lines");
+  return *x.begin == *y.begin && *x.end == *y.end;
+}
+bool operator<(const LineSegment &x,const LineSegment &y){
+  ASSERT(*x.line==*y.line,"Comparison of Linesegments from different lines");
+  return (*x.begin != *y.begin) ? *x.begin<*y.begin:*x.end<*y.end;
+}
Index: src/LineSegment.hpp
===================================================================
--- src/LineSegment.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/LineSegment.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,48 @@
+/*
+ * LineSegment.hpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#ifndef LINESEGMENT_HPP_
+#define LINESEGMENT_HPP_
+
+#include <memory>
+
+class Vector;
+class Line;
+class LinePoint;
+
+class LineSegment
+{
+  friend bool operator==(const LineSegment&,const LineSegment&);
+  friend bool operator<(const LineSegment&,const LineSegment&);
+
+  typedef std::auto_ptr<LinePoint> pointptr;
+  typedef std::auto_ptr<Line> lineptr;
+public:
+  LineSegment(const Vector&,const Vector&);
+  LineSegment(const LinePoint&,const LinePoint&);
+  LineSegment(const LineSegment&);
+  virtual ~LineSegment();
+  LineSegment &operator=(const LineSegment&);
+
+  bool isContained(const Vector &) const;
+  bool isContained(const LinePoint &) const;
+  bool overlaps(const LineSegment&) const;
+  bool hasZeroLength() const;
+  LinePoint getBegin() const;
+  LinePoint getEnd() const;
+  Line getLine() const;
+
+private:
+  pointptr begin;
+  pointptr end;
+  lineptr line;
+};
+
+bool operator==(const LineSegment&,const LineSegment&);
+bool operator<(const LineSegment&,const LineSegment&);
+
+#endif /* LINESEGMENT_HPP_ */
Index: src/LineSegmentSet.cpp
===================================================================
--- src/LineSegmentSet.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/LineSegmentSet.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,223 @@
+/*
+ * LineSegmentSet.cpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#include "LineSegmentSet.hpp"
+
+#include "Helpers/Assert.hpp"
+#include "Line.hpp"
+#include "LineSegment.hpp"
+
+#include "Helpers/Range.hpp"
+
+#include <algorithm>
+#include <boost/bind.hpp>
+
+using namespace std;
+
+LineSegmentSet::LineSegmentSet(const Line &_line) :
+  line(new Line(_line))
+{}
+
+LineSegmentSet::LineSegmentSet(const LineSegmentSet &src)
+{
+  line.reset(new Line(*src.line));
+  segments=src.segments;
+}
+
+LineSegmentSet& LineSegmentSet::operator=(const LineSegmentSet &src){
+  if(this!=&src){
+    lineptr _line(new Line(*src.line));
+    segments=src.segments;
+    line = _line;
+  }
+  return *this;
+}
+
+void LineSegmentSet::insert(const LineSegment &lineSeg){
+  ASSERT(*line==lineSeg.getLine(),"Inserted LineSegment was not from the line of the set");
+  // we might get messed up sets
+  if(lineSeg.hasZeroLength()){
+    // nothing to do here... there is practically nothing that is being inserted
+    return;
+  }
+  if(lineSeg.getBegin()>lineSeg.getEnd()){
+    // just insert the reversed set
+    insert(LineSegment(lineSeg.getEnd(),lineSeg.getBegin()));
+    return;
+  }
+  // get the first overlapping segment
+  set_t::iterator begin = find_if(segments.begin(),
+                                  segments.end(),
+                                  boost::bind(&LineSegment::overlaps,lineSeg,_1));
+  // see if we have overlapping segments
+  if(begin==segments.end()){
+    // no overlapping segments... we can just insert and get going
+    segments.insert(lineSeg);
+    return;
+  }
+  // get the last overlapping segment
+  set_t::iterator last = begin;
+  while(last!=segments.end() && lineSeg.overlaps(*last))++last;
+  last--;
+  // all segments between last and end now overlap and both are valid iterators
+  ASSERT(begin!=segments.end() && last!=segments.end(),"ERROR: last or begin are behind end of set");
+  // get the endpoints of the section
+  LinePoint beginPoint = min(begin->getBegin(),lineSeg.getBegin());
+  LinePoint endPoint = max(last->getEnd(),lineSeg.getEnd());
+  LineSegment insertSeg(beginPoint,endPoint);
+  // delete all overlapping segments
+  last++;
+  segments.erase(begin,last);
+  // insert the new combined segments
+  segments.insert(insertSeg);
+}
+
+void LineSegmentSet::erase(const LineSegment &lineSeg){
+  ASSERT(*line==lineSeg.getLine(),"Erased LineSegment was not from the line of the set");
+  // we might get messed up sets
+  if(lineSeg.hasZeroLength()){
+    // nothing to do here... there is practically nothing that is being erase
+    return;
+  }
+  if(lineSeg.getBegin()>lineSeg.getEnd()){
+    // just erase the reversed set
+    erase(LineSegment(lineSeg.getEnd(),lineSeg.getBegin()));
+    return;
+  }
+  // get the first overlapping segment
+  set_t::iterator begin = find_if(segments.begin(),
+                                  segments.end(),
+                                  boost::bind(&LineSegment::overlaps,lineSeg,_1));
+  if(begin==segments.end()){
+    // nothing to erase
+    return;
+  }
+  // get the last overlapping segment
+  set_t::iterator last = begin;
+  while(last!=segments.end() && lineSeg.overlaps(*last))++last;
+  last--;
+  // at the ends there might be two overlapping lineSegments of which parts remain
+  LineSegment beginSegment(begin->getBegin(),lineSeg.getBegin());
+  LineSegment endSegment(lineSeg.getEnd(),last->getEnd());
+  // erase all elements that have been overlapped
+  last++;
+  segments.erase(begin,last);
+  // insert the two segments if they still have space left
+  if(beginSegment.getBegin()<beginSegment.getEnd())
+    segments.insert(beginSegment);
+  if(endSegment.getBegin()<endSegment.getEnd())
+      segments.insert(endSegment);
+}
+
+LineSegmentSet::iterator LineSegmentSet::begin(){
+  return segments.begin();
+}
+LineSegmentSet::const_iterator LineSegmentSet::begin() const{
+  return segments.begin();
+}
+LineSegmentSet::iterator LineSegmentSet::end(){
+  return segments.end();
+}
+LineSegmentSet::const_iterator LineSegmentSet::end() const{
+  return segments.end();
+}
+
+bool LineSegmentSet::isContained(const Vector &point) const{
+  if(!line->isContained(point)){
+    return false;
+  }
+  LinePoint lp = line->getLinePoint(point);
+  for(set_t::iterator iter=segments.begin();iter!=segments.end();++iter){
+    if(iter->isContained(point)){
+      return true;
+    }
+    if(iter->getBegin()>lp){
+      break;
+    }
+  }
+  return false;
+}
+
+Line LineSegmentSet::getLine(){
+  return *line;
+}
+
+LineSegmentSet merge(const LineSegmentSet &x,const LineSegmentSet &y){
+  typedef range<LineSegmentSet::set_t::iterator> itRange;
+  ASSERT(*x.line==*y.line,"Lines of merges LineSegmentSets don't match");
+  LineSegmentSet res(*x.line);
+  // standard merge algorithm with a few modifications to deal with overlapping lines
+  itRange small(x.segments.begin(),x.segments.end());
+  itRange large(y.segments.begin(),y.segments.end());
+  while(small.first!=small.last && large.first != large.last){
+    if(large.first->getBegin()<small.first->getBegin()){
+      itRange tmp = small;
+      small = large;
+      large = tmp;
+    }
+
+    if(small.first->overlaps(*large.first)){
+      // get all overlapping segments from larger part
+      LineSegmentSet::set_t::iterator last = large.first;
+      while(last!=large.last && small.first->overlaps(*last))++last;
+      last--;
+      // get the endpoint of the new segment to be inserted
+      LinePoint lpEnd = max(small.first->getEnd(),last->getEnd());
+      res.segments.insert(LineSegment(small.first->getBegin(),lpEnd));
+      large.first=last; ++large.first; // all overlapping segments of the larger part have been dealt with
+    }
+    else{
+      // no overlap, we can just insert
+      res.segments.insert(*small.first);
+    }
+    ++small.first;
+  }
+  // deal with the rest of the iterator
+  // we don't know which one is done, so we just do both
+  res.segments.insert(small.first,small.last);
+  res.segments.insert(large.first,large.last);
+  return res;
+}
+
+LineSegmentSet intersect(const LineSegmentSet &x, const LineSegmentSet &y){
+  typedef range<LineSegmentSet::set_t::iterator> itRange;
+  ASSERT(*x.line==*y.line,"Lines of merges LineSegmentSets don't match");
+  LineSegmentSet res(*x.line);
+  itRange small(x.segments.begin(),x.segments.end());
+  itRange large(y.segments.begin(),y.segments.end());
+
+  while(small.first!=small.last && large.first != large.last){
+    if(large.first->getBegin()<small.first->getBegin()){
+      itRange tmp = small;
+      small = large;
+      large = tmp;
+    }
+    // add overlapping parts of segments from large
+    while(large.first!=large.last && small.first->overlaps(*large.first)){
+      LinePoint lpEnd = min(small.first->getEnd(),large.first->getEnd());
+      res.segments.insert(LineSegment(large.first->getBegin(),lpEnd));
+      ++large.first;
+    }
+    ++small.first;
+  }
+  return res;
+}
+
+LineSegmentSet invert(const LineSegmentSet &x){
+  LineSegmentSet res(*x.line);
+  LinePoint lpBegin = x.line->negEndpoint();
+  for(LineSegmentSet::set_t::iterator iter = x.segments.begin();iter!=x.segments.end();++iter){
+    if(lpBegin!=iter->getBegin())
+      res.segments.insert(LineSegment(lpBegin,iter->getBegin()));
+    // the end of the current scanned segment is the beginning of the next one to insert
+    lpBegin = iter->getEnd();
+  }
+  // we might need to extend to infinity
+  if(!lpBegin.isPosInfinity())
+    res.segments.insert(LineSegment(lpBegin,x.line->posEndpoint()));
+  return res;
+}
Index: src/LineSegmentSet.hpp
===================================================================
--- src/LineSegmentSet.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
+++ src/LineSegmentSet.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -0,0 +1,52 @@
+/*
+ * LineSegmentSet.hpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#ifndef LINESEGMENTSET_HPP_
+#define LINESEGMENTSET_HPP_
+
+#include <memory>
+#include <set>
+
+class Vector;
+class Line;
+class LineSegment;
+
+class LineSegmentSet {
+  friend LineSegmentSet merge(const LineSegmentSet&,const LineSegmentSet&);
+  friend LineSegmentSet intersect(const LineSegmentSet&, const LineSegmentSet&);
+  friend LineSegmentSet invert(const LineSegmentSet &);
+
+  typedef std::set<LineSegment> set_t;
+  typedef std::auto_ptr<Line> lineptr;
+public:
+  typedef set_t::iterator iterator;
+  typedef set_t::const_iterator const_iterator;
+
+  LineSegmentSet(const Line&);
+  LineSegmentSet(const LineSegmentSet&);
+  LineSegmentSet& operator=(const LineSegmentSet&);
+
+  void insert(const LineSegment&);
+  void erase(const LineSegment&);
+
+  iterator       begin();
+  const_iterator begin() const;
+  iterator       end();
+  const_iterator end()   const;
+
+  bool isContained(const Vector &) const;
+  Line getLine();
+private:
+  lineptr line;
+  set_t segments;
+};
+
+LineSegmentSet merge(const LineSegmentSet&,const LineSegmentSet&);
+LineSegmentSet intersect(const LineSegmentSet&, const LineSegmentSet&);
+LineSegmentSet invert(const LineSegmentSet &);
+
+#endif /* LINESEGMENTSET_HPP_ */
Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Makefile.am	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -84,4 +84,6 @@
   Exceptions/MissingValueException.cpp \
   Exceptions/NotInvertibleException.cpp \
+  Exceptions/NotOnSurfaceException.cpp \
+  Exceptions/ShapeException.cpp \
   Exceptions/ParseError.cpp \
   Exceptions/SkewException.cpp \
@@ -95,4 +97,6 @@
   Exceptions/MissingValueException.hpp \
   Exceptions/NotInvertibleException.hpp \
+  Exceptions/NotOnSurfaceException.hpp \
+  Exceptions/ShapeException.hpp \
   Exceptions/ParseError.hpp \
   Exceptions/SkewException.hpp \
@@ -153,4 +157,5 @@
   Descriptors/AtomIdDescriptor.cpp \
   Descriptors/AtomSelectionDescriptor.cpp \
+  Descriptors/AtomShapeDescriptor.cpp \
   Descriptors/AtomTypeDescriptor.cpp \
   Descriptors/MoleculeDescriptor.cpp \
@@ -165,4 +170,5 @@
   Descriptors/AtomIdDescriptor.hpp \
   Descriptors/AtomSelectionDescriptor.hpp \
+  Descriptors/AtomShapeDescriptor.hpp \
   Descriptors/AtomTypeDescriptor.hpp \
   Descriptors/MoleculeDescriptor.hpp \
@@ -214,4 +220,6 @@
   leastsquaremin.cpp \
   Line.cpp \
+  LineSegment.cpp \
+  LineSegmentSet.cpp \
   linkedcell.cpp \
   log.cpp \
@@ -265,4 +273,6 @@
   leastsquaremin.hpp \
   Line.hpp \
+  LineSegment.hpp \
+  LineSegmentSet.hpp \
   linkedcell.hpp \
   lists.hpp \
Index: src/Matrix.cpp
===================================================================
--- src/Matrix.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Matrix.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -213,4 +213,11 @@
 }
 
+Matrix Matrix::transpose() const{
+  MatrixContent *newContent = new MatrixContent();
+  gsl_matrix_transpose_memcpy(newContent->content, content->content);
+  Matrix res = Matrix(newContent);
+  return res;
+}
+
 Matrix &Matrix::operator*=(const double factor){
   gsl_matrix_scale(content->content, factor);
Index: src/Matrix.hpp
===================================================================
--- src/Matrix.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Matrix.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -96,4 +96,6 @@
   Matrix invert() const;
 
+  Matrix transpose() const;
+
   // operators
   Matrix &operator=(const Matrix&);
Index: src/Plane.cpp
===================================================================
--- src/Plane.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Plane.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -15,4 +15,5 @@
 #include "verbose.hpp"
 #include "Helpers/Assert.hpp"
+#include "helpers.hpp"
 #include <cmath>
 #include "Line.hpp"
@@ -106,4 +107,12 @@
 Plane::~Plane()
 {}
+
+Plane &Plane::operator=(const Plane &rhs){
+  if(&rhs!=this){
+    normalVector.reset(new Vector(*rhs.normalVector));
+    offset = rhs.offset;
+  }
+  return *this;
+}
 
 
@@ -188,4 +197,9 @@
 }
 
+bool Plane::onSameSide(const Vector &point1,const Vector &point2) const{
+  return sign(point1.ScalarProduct(*normalVector)-offset) ==
+         sign(point2.ScalarProduct(*normalVector)-offset);
+}
+
 /************ Methods inherited from Space ****************/
 
@@ -207,4 +221,8 @@
 // Operators
 
+bool operator==(const Plane &x,const Plane &y){
+  return *x.normalVector == *y.normalVector && x.offset == y.offset;
+}
+
 ostream &operator << (ostream &ost,const Plane &p){
   ost << "<" << p.getNormal() << ";x> - " << p.getOffset() << "=0";
Index: src/Plane.hpp
===================================================================
--- src/Plane.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Plane.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -21,4 +21,5 @@
 class Plane : public Space
 {
+  friend bool operator==(const Plane&,const Plane&);
   typedef std::auto_ptr<Vector> vec_ptr;
 public:
@@ -29,4 +30,6 @@
   Plane(const Plane& plane);
   virtual ~Plane();
+
+  Plane &operator=(const Plane&);
 
   // Accessor Functions
@@ -63,4 +66,9 @@
   Line getOrthogonalLine(const Vector &origin) const;
 
+  /**
+   * Test if two points are on the same side of the plane
+   */
+  bool onSameSide(const Vector&,const Vector&) const;
+
   /****** Methods inherited from Space ***********/
 
@@ -73,4 +81,6 @@
 };
 
+bool operator==(const Plane&,const Plane&);
+
 std::ostream &operator<< (std::ostream &ost,const Plane& p);
 
Index: src/Shapes/BaseShapes.cpp
===================================================================
--- src/Shapes/BaseShapes.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/BaseShapes.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -8,9 +8,43 @@
 #include "Shapes/BaseShapes.hpp"
 #include "Shapes/BaseShapes_impl.hpp"
+#include "Shapes/ShapeOps.hpp"
 
 #include "vector.hpp"
+#include "Helpers/Assert.hpp"
+
+#include "Line.hpp"
+#include "Plane.hpp"
+#include "LineSegment.hpp"
+#include "LineSegmentSet.hpp"
+
+#include <cmath>
+#include <algorithm>
 
 bool Sphere_impl::isInside(const Vector &point){
   return point.NormSquared()<=1;
+}
+
+bool Sphere_impl::isOnSurface(const Vector &point){
+  return fabs(point.NormSquared()-1)<MYEPSILON;
+}
+
+Vector Sphere_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  return point;
+}
+
+LineSegmentSet Sphere_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  std::vector<Vector> intersections = line.getSphereIntersections();
+  if(intersections.size()==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Sphere_impl::toString(){
+  return "Sphere()";
 }
 
@@ -20,10 +54,96 @@
 }
 
+Shape Sphere(const Vector &center,double radius){
+  return translate(resize(Sphere(),radius),center);
+}
+
+Shape Ellipsoid(const Vector &center, const Vector &radius){
+  return translate(stretch(Sphere(),radius),center);
+}
+
 bool Cuboid_impl::isInside(const Vector &point){
-  return point[0]<=1 && point[1]<=1 && point[2]<=1;
+  return fabs(point[0])<=1 && fabs(point[1])<=1 && fabs(point[2])<=1;
+}
+
+bool Cuboid_impl::isOnSurface(const Vector &point){
+  bool retVal = isInside(point);
+  // test all borders of the cuboid
+  // double fabs
+  retVal = retVal &&
+           ((fabs(fabs(point[0])-1)  < MYEPSILON) ||
+            (fabs(fabs(point[1])-1)  < MYEPSILON) ||
+            (fabs(fabs(point[2])-1)  < MYEPSILON));
+  return retVal;
+}
+
+Vector Cuboid_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res;
+  // figure out on which sides the Vector lies (maximum 3, when it is in a corner)
+  for(int i=NDIM;i--;){
+    if(fabs(fabs(point[i])-1)<MYEPSILON){
+      // add the scaled (-1/+1) Vector to the set of surface vectors
+      res[i] = point[i];
+    }
+  }
+  ASSERT(res.NormSquared()>=1 && res.NormSquared()<=3,"To many or to few sides found for this Vector");
+
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet Cuboid_impl::getLineIntersections(const Line &line){
+  LineSegmentSet res(line);
+  // get the intersection on each of the six faces
+  vector<Vector> intersections;
+  intersections.resize(2);
+  int c=0;
+  int x[2]={-1,+1};
+  for(int i=NDIM;i--;){
+    for(int p=0;p<2;++p){
+      if(c==2) goto end; // I know this sucks, but breaking two loops is stupid
+      Vector base;
+      base[i]=x[p];
+      // base now points to the surface and is normal to it at the same time
+      Plane p(base,base);
+      Vector intersection = p.GetIntersection(line);
+      if(isInside(intersection)){
+        // if we have a point on the edge it might already be contained
+        if(c==1 && intersections[0]==intersection)
+          continue;
+        intersections[c++]=intersection;
+      }
+    }
+  }
+  end:
+  if(c==2){
+    res.insert(LineSegment(intersections[0],intersections[1]));
+  }
+  return res;
+}
+
+string Cuboid_impl::toString(){
+  return "Cuboid()";
 }
 
 Shape Cuboid(){
-  Shape::impl_ptr impl = Shape::impl_ptr(new Sphere_impl());
+  Shape::impl_ptr impl = Shape::impl_ptr(new Cuboid_impl());
   return Shape(impl);
 }
+
+Shape Cuboid(const Vector &corner1, const Vector &corner2){
+  // make sure the two edges are upper left front and lower right back
+  Vector sortedC1;
+  Vector sortedC2;
+  for(int i=NDIM;i--;){
+    sortedC1[i] = min(corner1[i],corner2[i]);
+    sortedC2[i] = max(corner1[i],corner2[i]);
+    ASSERT(corner1[i]!=corner2[i],"Given points for cuboid edges did not define a valid space");
+  }
+  // get the middle point
+  Vector middle = (1./2.)*(sortedC1+sortedC2);
+  Vector factors = sortedC2-middle;
+  return translate(stretch(Cuboid(),factors),middle);
+}
Index: src/Shapes/BaseShapes.hpp
===================================================================
--- src/Shapes/BaseShapes.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/BaseShapes.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -12,5 +12,8 @@
 
 Shape Sphere();
+Shape Sphere(const Vector &center,double radius);
+Shape Ellipsoid(const Vector &center, const Vector &radius);
 Shape Cuboid();
+Shape Cuboid(const Vector &corner1, const Vector &corner2);
 
 #endif /* BASESHAPES_HPP_ */
Index: src/Shapes/BaseShapes_impl.hpp
===================================================================
--- src/Shapes/BaseShapes_impl.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/BaseShapes_impl.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -13,8 +13,16 @@
 class Sphere_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
 };
 
 class Cuboid_impl : public Shape_impl {
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
 };
 
Index: src/Shapes/Shape.cpp
===================================================================
--- src/Shapes/Shape.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/Shape.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -9,4 +9,10 @@
 #include "Shape_impl.hpp"
 
+#include "Helpers/Assert.hpp"
+
+#include <string>
+
+using namespace std;
+
 Shape::Shape(const Shape& src) :
   impl(src.getImpl())
@@ -17,4 +23,16 @@
 bool Shape::isInside(const Vector &point) const{
   return impl->isInside(point);
+}
+
+bool Shape::isOnSurface(const Vector &point) const{
+  return impl->isOnSurface(point);
+}
+
+Vector Shape::getNormal(const Vector &point) const throw (NotOnSurfaceException){
+  return impl->getNormal(point);
+}
+
+LineSegmentSet Shape::getLineIntersections(const Line &line){
+  return impl->getLineIntersections(line);
 }
 
@@ -30,4 +48,8 @@
 }
 
+std::string Shape::toString() const{
+  return impl->toString();
+}
+
 Shape::impl_ptr Shape::getImpl() const{
   return impl;
@@ -65,4 +87,55 @@
 }
 
+bool AndShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector AndShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet AndShape_impl::getLineIntersections(const Line &line){
+  return intersect(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string AndShape_impl::toString(){
+  return string("(") + lhs->toString() + string("&&") + rhs->toString() + string(")");
+}
+
 Shape operator&&(const Shape &lhs,const Shape &rhs){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new AndShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
@@ -82,4 +155,55 @@
 }
 
+bool OrShape_impl::isOnSurface(const Vector &point){
+  // check the number of surfaces that this point is on
+  int surfaces =0;
+  surfaces += lhs->isOnSurface(point);
+  surfaces += rhs->isOnSurface(point);
+
+  switch(surfaces){
+    case 0:
+      return false;
+      // no break necessary
+    case 1:
+      // if it is inside for the object where it does not lie on
+      // the surface the whole point lies inside
+      return (lhs->isOnSurface(point) && !rhs->isInside(point)) ||
+             (rhs->isOnSurface(point) && !lhs->isInside(point));
+      // no break necessary
+    case 2:
+      {
+        // it lies on both Shapes... could be an edge or an inner point
+        // test the direction of the normals
+        Vector direction=lhs->getNormal(point)+rhs->getNormal(point);
+        // if the directions are opposite we lie on the inside
+        return !direction.IsZero();
+      }
+      // no break necessary
+    default:
+      // if this happens there is something wrong
+      ASSERT(0,"Default case should have never been used");
+  }
+  return false; // never reached
+}
+
+Vector OrShape_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector res;
+  if(!isOnSurface(point)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  res += lhs->isOnSurface(point)?lhs->getNormal(point):zeroVec;
+  res += rhs->isOnSurface(point)?rhs->getNormal(point):zeroVec;
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet OrShape_impl::getLineIntersections(const Line &line){
+  return merge(lhs->getLineIntersections(line),rhs->getLineIntersections(line));
+}
+
+string OrShape_impl::toString(){
+  return string("(") + lhs->toString() + string("||") + rhs->toString() + string(")");
+}
+
 Shape operator||(const Shape &lhs,const Shape &rhs){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new OrShape_impl(getShapeImpl(lhs),getShapeImpl(rhs)));
@@ -99,6 +223,28 @@
 }
 
+bool NotShape_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(point);
+}
+
+Vector NotShape_impl::getNormal(const Vector &point) throw(NotOnSurfaceException){
+  return -1*arg->getNormal(point);
+}
+
+LineSegmentSet NotShape_impl::getLineIntersections(const Line &line){
+  return invert(arg->getLineIntersections(line));
+}
+
+string NotShape_impl::toString(){
+  return string("!") + arg->toString();
+}
+
 Shape operator!(const Shape &arg){
   Shape::impl_ptr newImpl = Shape::impl_ptr(new NotShape_impl(getShapeImpl(arg)));
   return Shape(newImpl);
 }
+
+/**************** global operations *********************************/
+ostream &operator<<(ostream &ost,const Shape &shape){
+  ost << shape.toString();
+  return ost;
+}
Index: src/Shapes/Shape.hpp
===================================================================
--- src/Shapes/Shape.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/Shape.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -10,7 +10,12 @@
 
 #include <boost/shared_ptr.hpp>
+#include <iosfwd>
+
+#include "Exceptions/NotOnSurfaceException.hpp"
 
 class Vector;
 class Shape_impl;
+class LineSegmentSet;
+class Line;
 
 class Shape
@@ -25,7 +30,12 @@
 
   bool isInside(const Vector &point) const;
+  bool isOnSurface(const Vector &point) const;
+  Vector getNormal(const Vector &point) const throw(NotOnSurfaceException);
+
+  LineSegmentSet getLineIntersections(const Line&);
 
   Shape &operator=(const Shape& rhs);
 
+  std::string toString() const;
 protected:
   impl_ptr getImpl() const;
@@ -42,3 +52,5 @@
 Shape operator!(const Shape&);
 
+std::ostream &operator<<(std::ostream&,const Shape&);
+
 #endif /* SHAPE_HPP_ */
Index: src/Shapes/ShapeOps.cpp
===================================================================
--- src/Shapes/ShapeOps.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/ShapeOps.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -11,8 +11,56 @@
 #include "Helpers/Assert.hpp"
 
+/*************** Base case ***********************/
+
+ShapeOpsBase_impl::ShapeOpsBase_impl(const Shape::impl_ptr &_arg) :
+  arg(_arg){}
+
+ShapeOpsBase_impl::~ShapeOpsBase_impl(){}
+
+bool ShapeOpsBase_impl::isInside(const Vector &point){
+  return arg->isInside(translateIn(point));
+}
+
+bool ShapeOpsBase_impl::isOnSurface(const Vector &point){
+  return arg->isOnSurface(translateIn(point));
+}
+
+Vector ShapeOpsBase_impl::getNormal(const Vector &point) throw (NotOnSurfaceException){
+  Vector helper = translateIn(point);
+  if(!arg->isOnSurface(helper)){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  Vector res = translateOutNormal(arg->getNormal(helper));
+  res.Normalize();
+  return res;
+}
+
+LineSegmentSet ShapeOpsBase_impl::getLineIntersections(const Line &line){
+  Line newLine(translateIn(line.getOrigin()),translateIn(line.getDirection()));
+  LineSegmentSet res(line);
+  LineSegmentSet helper = getArg()->getLineIntersections(newLine);
+  for(LineSegmentSet::iterator iter = helper.begin();iter!=helper.end();++iter){
+    LinePoint lpBegin = iter->getBegin();
+    LinePoint lpEnd = iter->getBegin();
+    // translate both linepoints
+    lpBegin = lpBegin.isNegInfinity()?
+                line.negEndpoint():
+                line.getLinePoint(translateOutPos(lpBegin.getPoint()));
+    lpEnd = lpEnd.isPosInfinity()?
+              line.posEndpoint():
+              line.getLinePoint(translateOutPos(lpEnd.getPoint()));
+    res.insert(LineSegment(lpBegin,lpEnd));
+  }
+  return res;
+}
+
+Shape::impl_ptr ShapeOpsBase_impl::getArg(){
+  return arg;
+}
+
 /********************* Resize ********************/
 
 Resize_impl::Resize_impl(const Shape::impl_ptr &_arg,double _size) :
-  arg(_arg), size(_size)
+  ShapeOpsBase_impl(_arg), size(_size)
 {
   ASSERT(size>0,"Cannot resize a Shape to size zero or below");
@@ -21,6 +69,20 @@
 Resize_impl::~Resize_impl(){}
 
-bool Resize_impl::isInside(const Vector& point){
-  return arg->isInside((1/size) * point);
+Vector Resize_impl::translateIn(const Vector& point){
+  return (1/size) * point;
+}
+
+Vector Resize_impl::translateOutPos(const Vector& point){
+  return size * point;
+}
+
+Vector Resize_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Resize_impl::toString(){
+  stringstream sstr;
+  sstr << "resize(" << getArg()->toString() << "," << size << ")";
+  return sstr.str();
 }
 
@@ -33,11 +95,25 @@
 
 Translate_impl::Translate_impl(const Shape::impl_ptr &_arg, const Vector &_offset) :
-  arg(_arg),offset(_offset)
+  ShapeOpsBase_impl(_arg),offset(_offset)
 {}
 
 Translate_impl::~Translate_impl(){}
 
-bool Translate_impl::isInside(const Vector& point){
-  return arg->isInside(point-offset);
+Vector Translate_impl::translateIn(const Vector& point){
+  return point-offset;
+}
+
+Vector Translate_impl::translateOutPos(const Vector& point){
+  return point+offset;
+}
+
+Vector Translate_impl::translateOutNormal(const Vector& point){
+  return point;
+}
+
+string Translate_impl::toString(){
+  stringstream sstr;
+  sstr << "translate(" << getArg()->toString() << "," << offset << ")";
+  return sstr.str();
 }
 
@@ -50,5 +126,5 @@
 
 Stretch_impl::Stretch_impl(const Shape::impl_ptr &_arg, const Vector &_factors) :
-  arg(_arg),factors(_factors)
+  ShapeOpsBase_impl(_arg),factors(_factors)
 {
   ASSERT(factors[0]>0,"cannot stretch a shape by a negative amount");
@@ -62,8 +138,32 @@
 Stretch_impl::~Stretch_impl(){}
 
-bool Stretch_impl::isInside(const Vector& point){
+Vector Stretch_impl::translateIn(const Vector& point){
   Vector helper=point;
   helper.ScaleAll(reciFactors);
-  return arg->isInside(helper);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutPos(const Vector& point){
+  Vector helper=point;
+  helper.ScaleAll(factors);
+  return helper;
+}
+
+Vector Stretch_impl::translateOutNormal(const Vector& point){
+  Vector helper=point;
+  // the normalFactors are derived from appearances of the factors
+  // with in the vectorproduct
+  Vector normalFactors;
+  normalFactors[0]=factors[1]*factors[2];
+  normalFactors[1]=factors[0]*factors[2];
+  normalFactors[2]=factors[0]*factors[1];
+  helper.ScaleAll(normalFactors);
+  return helper;
+}
+
+string Stretch_impl::toString(){
+  stringstream sstr;
+  sstr << "stretch(" << getArg()->toString() << "," << factors << ")";
+  return sstr.str();
 }
 
@@ -76,5 +176,5 @@
 
 Transform_impl::Transform_impl(const Shape::impl_ptr &_arg, const Matrix &_transformation) :
-  arg(_arg),transformation(_transformation)
+  ShapeOpsBase_impl(_arg),transformation(_transformation)
 {
   transformationInv = transformation.invert();
@@ -83,6 +183,21 @@
 Transform_impl::~Transform_impl(){}
 
-bool Transform_impl::isInside(const Vector& point){
-  return arg->isInside(transformationInv * point);
+Vector Transform_impl::translateIn(const Vector& point){
+  return transformationInv * point;
+}
+
+Vector Transform_impl::translateOutPos(const Vector& point){
+  return transformation * point;
+}
+
+Vector Transform_impl::translateOutNormal(const Vector& point){
+  Matrix mat = transformation.invert().transpose();
+  return mat * point;
+}
+
+string Transform_impl::toString(){
+  stringstream sstr;
+  sstr << "transform(" << getArg()->toString() << "," << transformation << ")";
+  return sstr.str();
 }
 
Index: src/Shapes/ShapeOps_impl.hpp
===================================================================
--- src/Shapes/ShapeOps_impl.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/ShapeOps_impl.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -13,46 +13,75 @@
 #include "Matrix.hpp"
 
-class Resize_impl :  public Shape_impl
+class ShapeOpsBase_impl : public Shape_impl{
+public:
+  ShapeOpsBase_impl(const Shape::impl_ptr&);
+  virtual ~ShapeOpsBase_impl();
+  virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw (NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+protected:
+  virtual Vector translateIn(const Vector &point)=0;
+  virtual Vector translateOutPos(const Vector &point)=0;
+  virtual Vector translateOutNormal(const Vector &point)=0;
+  Shape::impl_ptr getArg();
+private:
+  Shape::impl_ptr arg;
+};
+
+class Resize_impl :  public ShapeOpsBase_impl
 {
 public:
   Resize_impl(const Shape::impl_ptr&,double);
   virtual ~Resize_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
 private:
-  Shape::impl_ptr arg;
   double size;
 };
 
-class Translate_impl :  public Shape_impl
+class Translate_impl :  public ShapeOpsBase_impl
 {
 public:
   Translate_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Translate_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
 private:
-  Shape::impl_ptr arg;
   Vector offset;
 };
 
-class Stretch_impl : public Shape_impl
+class Stretch_impl : public ShapeOpsBase_impl
 {
 public:
   Stretch_impl(const Shape::impl_ptr&, const Vector&);
   virtual ~Stretch_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
 private:
-  Shape::impl_ptr arg;
   Vector factors;
   Vector reciFactors;
 };
 
-class Transform_impl : public Shape_impl
+class Transform_impl : public ShapeOpsBase_impl
 {
 public:
   Transform_impl(const Shape::impl_ptr&, const Matrix&);
   virtual ~Transform_impl();
-  virtual bool isInside(const Vector& point);
+protected:
+  virtual Vector translateIn(const Vector &point);
+  virtual Vector translateOutPos(const Vector &point);
+  virtual Vector translateOutNormal(const Vector &point);
+  virtual std::string toString();
 private:
-  Shape::impl_ptr arg;
   Matrix transformation;
   Matrix transformationInv;
Index: src/Shapes/Shape_impl.hpp
===================================================================
--- src/Shapes/Shape_impl.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/Shapes/Shape_impl.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -10,4 +10,9 @@
 
 #include "Shapes/Shape.hpp"
+#include "vector.hpp"
+#include "Line.hpp"
+#include "LineSegment.hpp"
+#include "LineSegmentSet.hpp"
+
 
 class Shape_impl {
@@ -16,4 +21,8 @@
   virtual ~Shape_impl(){};
   virtual bool isInside(const Vector &point)=0;
+  virtual bool isOnSurface(const Vector &point)=0;
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException)=0;
+  virtual LineSegmentSet getLineIntersections(const Line&)=0;
+  virtual std::string toString()=0;
 };
 
@@ -23,4 +32,18 @@
     return true;
   }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    LineSegmentSet res(line);
+    res.insert(LineSegment(line.negEndpoint(),line.posEndpoint()));
+    return res;
+  }
+  virtual std::string toString(){
+    return "Everywhere()";
+  }
 };
 
@@ -28,4 +51,16 @@
   virtual bool isInside(const Vector &point){
     return false;
+  }
+  virtual bool isOnSurface(const Vector &point){
+    return false;
+  }
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException){
+    throw NotOnSurfaceException(__FILE__,__LINE__);
+  }
+  virtual LineSegmentSet getLineIntersections(const Line &line){
+    return LineSegmentSet(line);
+  }
+  virtual std::string toString(){
+    return "Nowhere()";
   }
 };
@@ -36,4 +71,8 @@
   virtual ~AndShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
 private:
   Shape::impl_ptr lhs;
@@ -46,4 +85,8 @@
   virtual ~OrShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
 private:
   Shape::impl_ptr lhs;
@@ -56,4 +99,8 @@
   virtual ~NotShape_impl();
   virtual bool isInside(const Vector &point);
+  virtual bool isOnSurface(const Vector &point);
+  virtual Vector getNormal(const Vector &point) throw(NotOnSurfaceException);
+  virtual LineSegmentSet getLineIntersections(const Line&);
+  virtual std::string toString();
 private:
   Shape::impl_ptr arg;
Index: src/helpers.cpp
===================================================================
--- src/helpers.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/helpers.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -141,2 +141,12 @@
   exit(255);
 }
+
+sign_t sign(double value){
+  if(fabs(value)<MYEPSILON){
+    return Zero;
+  }
+  if(value<0)
+    return Minus;
+  else
+    return Plus;
+}
Index: src/helpers.hpp
===================================================================
--- src/helpers.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/helpers.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -43,4 +43,10 @@
     const char *file, const int line);
   //__attribute__ ((__return__));
+
+typedef enum {
+  Minus = -1,
+  Zero = 0,
+  Plus = +1
+} sign_t;
 
 double ask_value(const char *text);
@@ -52,4 +58,5 @@
 int CompareDoubles (const void * a, const void * b);
 void performCriticalExit();
+sign_t sign(double value);
 
 /********************************************** helpful template functions *********************************/
Index: src/unittests/LineUnittest.cpp
===================================================================
--- src/unittests/LineUnittest.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/unittests/LineUnittest.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -29,12 +29,12 @@
 void LineUnittest::setUp(){
   // three lines along the axes
-  la1 = new Line(zeroVec,e1);
-  la2 = new Line(zeroVec,e2);
-  la3 = new Line(zeroVec,e3);
+  la1 = new Line(zeroVec,unitVec[0]);
+  la2 = new Line(zeroVec,unitVec[1]);
+  la3 = new Line(zeroVec,unitVec[2]);
 
   // 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);
+  lp1 = new Line(unitVec[0],unitVec[0]-unitVec[1]);
+  lp2 = new Line(unitVec[1],unitVec[1]-unitVec[2]);
+  lp3 = new Line(unitVec[2],unitVec[2]-unitVec[0]);
 }
 void LineUnittest::tearDown(){
@@ -52,21 +52,21 @@
 
   // 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));
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[0]));
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[1]));
+  CPPUNIT_ASSERT_NO_THROW(Line(zeroVec,unitVec[2]));
 
   // 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));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[0],unitVec[1]));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[1],unitVec[2]));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[2],unitVec[0]));
   // for zerovectors
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e1,zeroVec));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e2,zeroVec));
-  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(e3,zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[0],zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[1],zeroVec));
+  CPPUNIT_ASSERT_NO_THROW(makeLineThrough(unitVec[2],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);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[0],unitVec[0]),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[1],unitVec[1]),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(makeLineThrough(unitVec[2],unitVec[2]),LinearDependenceException);
 
 }
@@ -78,7 +78,7 @@
 void LineUnittest::constructionResultTest(){
   // test all directions
-  CPPUNIT_ASSERT(testDirection(la1->getDirection(),e1));
-  CPPUNIT_ASSERT(testDirection(la2->getDirection(),e2));
-  CPPUNIT_ASSERT(testDirection(la3->getDirection(),e3));
+  CPPUNIT_ASSERT(testDirection(la1->getDirection(),unitVec[0]));
+  CPPUNIT_ASSERT(testDirection(la2->getDirection(),unitVec[1]));
+  CPPUNIT_ASSERT(testDirection(la3->getDirection(),unitVec[2]));
 
   // test origins
@@ -92,15 +92,15 @@
   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));
+  CPPUNIT_ASSERT(la1->isContained(unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(unitVec[2]));
+
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]));
+
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[2]));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[0]));
 }
 
@@ -112,54 +112,54 @@
 
   // 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));
+  CPPUNIT_ASSERT(la1->isContained(unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(unitVec[2]));
+
+  CPPUNIT_ASSERT(la1->isContained(2*unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(2*unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(2*unitVec[2]));
+
+  CPPUNIT_ASSERT(la1->isContained(3*unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(3*unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(3*unitVec[2]));
 
   // 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));
+  CPPUNIT_ASSERT(la1->isContained(-1*unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(-1*unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(-1*unitVec[2]));
+
+  CPPUNIT_ASSERT(la1->isContained(-2*unitVec[0]));
+  CPPUNIT_ASSERT(la2->isContained(-2*unitVec[1]));
+  CPPUNIT_ASSERT(la3->isContained(-2*unitVec[2]));
 
   // 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));
+  CPPUNIT_ASSERT(!la1->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(!la2->isContained(unitVec[2]));
+  CPPUNIT_ASSERT(!la3->isContained(unitVec[0]));
+
+  CPPUNIT_ASSERT(!la1->isContained(2*unitVec[1]));
+  CPPUNIT_ASSERT(!la2->isContained(2*unitVec[2]));
+  CPPUNIT_ASSERT(!la3->isContained(2*unitVec[0]));
+
+  CPPUNIT_ASSERT(!la1->isContained(-1*unitVec[1]));
+  CPPUNIT_ASSERT(!la2->isContained(-1*unitVec[2]));
+  CPPUNIT_ASSERT(!la3->isContained(-1*unitVec[0]));
 
   // 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)));
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]));
+
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[1]));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[2]));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[0]));
+
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]+2*(unitVec[0]-unitVec[1])));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]+2*(unitVec[1]-unitVec[2])));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]+2*(unitVec[2]-unitVec[0])));
+
+  CPPUNIT_ASSERT(lp1->isContained(unitVec[0]-2*(unitVec[0]-unitVec[1])));
+  CPPUNIT_ASSERT(lp2->isContained(unitVec[1]-2*(unitVec[1]-unitVec[2])));
+  CPPUNIT_ASSERT(lp3->isContained(unitVec[2]-2*(unitVec[2]-unitVec[0])));
 }
 
@@ -177,24 +177,24 @@
   // axes and plane lines
   fixture = la1->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
   fixture = la2->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = la3->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   fixture = la1->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
   fixture = la2->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = la3->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   // two plane lines
   fixture = lp1->getIntersection(*lp2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
   fixture = lp2->getIntersection(*lp3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
   fixture = lp3->getIntersection(*lp1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
 
   // When we have two times the same line, we check if the point is on the line
@@ -242,80 +242,80 @@
 
   // 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);
+  fixture = la1->rotateVector(unitVec[0],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = la2->rotateVector(unitVec[1],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = la3->rotateVector(unitVec[2],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
+
+  fixture = la1->rotateVector(unitVec[0],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = la2->rotateVector(unitVec[1],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = la3->rotateVector(unitVec[2],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   // 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);
+  fixture = la1->rotateVector(2*unitVec[0],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[0]);
+  fixture = la2->rotateVector(2*unitVec[1],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[1]);
+  fixture = la3->rotateVector(2*unitVec[2],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[2]);
+
+  fixture = la1->rotateVector(2*unitVec[0],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[0]);
+  fixture = la2->rotateVector(2*unitVec[1],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[1]);
+  fixture = la3->rotateVector(2*unitVec[2],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,2*unitVec[2]);
 
   // 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);
+  fixture = la1->rotateVector(-1*unitVec[0],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
+  fixture = la2->rotateVector(-1*unitVec[1],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = la3->rotateVector(-1*unitVec[2],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
+
+  fixture = la1->rotateVector(-1*unitVec[0],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
+  fixture = la2->rotateVector(-1*unitVec[1],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = la3->rotateVector(-1*unitVec[2],2.);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
 
 
 
   // 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);
+  // unitVec[1] around unitVec[0]
+  fixture = la1->rotateVector(unitVec[1],0);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = la1->rotateVector(unitVec[1],1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
+  fixture = la1->rotateVector(unitVec[1],M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = la1->rotateVector(unitVec[1],2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+
+  // unitVec[2] around unitVec[1]
+  fixture = la2->rotateVector(unitVec[2],0);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
+  fixture = la2->rotateVector(unitVec[2],1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
+  fixture = la2->rotateVector(unitVec[2],M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
+  fixture = la2->rotateVector(unitVec[2],2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
+
+  // unitVec[0] around unitVec[2]
+  fixture = la3->rotateVector(unitVec[0],0);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = la3->rotateVector(unitVec[0],1./2.*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = la3->rotateVector(unitVec[0],M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
+  fixture = la3->rotateVector(unitVec[0],2*M_PI);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
 
 
@@ -323,18 +323,18 @@
 
   // 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);
+  fixture = lp1->rotateVector(unitVec[0],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = lp1->rotateVector(unitVec[1],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+
+  fixture = lp2->rotateVector(unitVec[1],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = lp2->rotateVector(unitVec[2],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
+
+  fixture = lp3->rotateVector(unitVec[2],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
+  fixture = lp3->rotateVector(unitVec[0],1.);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
 
   // the real stuff
@@ -358,6 +358,6 @@
     std::vector<Vector> res = la1->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e1));
-    CPPUNIT_ASSERT(testDirection(res[1],e1));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[0]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[0]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -366,6 +366,6 @@
     std::vector<Vector> res = la2->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e2));
-    CPPUNIT_ASSERT(testDirection(res[1],e2));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[1]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[1]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -374,6 +374,6 @@
     std::vector<Vector> res = la3->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT(testDirection(res[0],e3));
-    CPPUNIT_ASSERT(testDirection(res[1],e3));
+    CPPUNIT_ASSERT(testDirection(res[0],unitVec[2]));
+    CPPUNIT_ASSERT(testDirection(res[1],unitVec[2]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -382,6 +382,6 @@
     std::vector<Vector> res = lp1->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
-    CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
+    CPPUNIT_ASSERT((res[0]==unitVec[0]) || (res[0]==unitVec[1]));
+    CPPUNIT_ASSERT((res[1]==unitVec[0]) || (res[1]==unitVec[1]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -390,6 +390,6 @@
     std::vector<Vector> res = lp2->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
-    CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
+    CPPUNIT_ASSERT((res[0]==unitVec[1]) || (res[0]==unitVec[2]));
+    CPPUNIT_ASSERT((res[1]==unitVec[1]) || (res[1]==unitVec[2]));
     CPPUNIT_ASSERT(res[0]!=res[1]);
   }
@@ -398,7 +398,7 @@
     std::vector<Vector> res = lp3->getSphereIntersections();
     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
-    CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
-    CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
-    CPPUNIT_ASSERT(res[0]!=res[1]);
-  }
-}
+    CPPUNIT_ASSERT((res[0]==unitVec[2]) || (res[0]==unitVec[0]));
+    CPPUNIT_ASSERT((res[1]==unitVec[2]) || (res[1]==unitVec[0]));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+}
Index: src/unittests/MatrixUnittest.cpp
===================================================================
--- src/unittests/MatrixUnittest.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/unittests/MatrixUnittest.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -39,28 +39,28 @@
   }
   perm1 = new Matrix();
-  perm1->column(0) = e1;
-  perm1->column(1) = e3;
-  perm1->column(2) = e2;
+  perm1->column(0) = unitVec[0];
+  perm1->column(1) = unitVec[2];
+  perm1->column(2) = unitVec[1];
 
 
   perm2 = new Matrix();
-  perm2->column(0) = e2;
-  perm2->column(1) = e1;
-  perm2->column(2) = e3;
+  perm2->column(0) = unitVec[1];
+  perm2->column(1) = unitVec[0];
+  perm2->column(2) = unitVec[2];
 
   perm3 = new Matrix();
-  perm3->column(0) = e2;
-  perm3->column(1) = e3;
-  perm3->column(2) = e1;
+  perm3->column(0) = unitVec[1];
+  perm3->column(1) = unitVec[2];
+  perm3->column(2) = unitVec[0];
 
   perm4 = new Matrix();
-  perm4->column(0) = e3;
-  perm4->column(1) = e2;
-  perm4->column(2) = e1;
+  perm4->column(0) = unitVec[2];
+  perm4->column(1) = unitVec[1];
+  perm4->column(2) = unitVec[0];
 
   perm5 = new Matrix();
-  perm5->column(0) = e3;
-  perm5->column(1) = e1;
-  perm5->column(2) = e2;
+  perm5->column(0) = unitVec[2];
+  perm5->column(1) = unitVec[0];
+  perm5->column(2) = unitVec[1];
 
 }
@@ -108,10 +108,10 @@
 
   mat.one();
-  CPPUNIT_ASSERT_EQUAL(mat.row(0),e1);
-  CPPUNIT_ASSERT_EQUAL(mat.row(1),e2);
-  CPPUNIT_ASSERT_EQUAL(mat.row(2),e3);
-  CPPUNIT_ASSERT_EQUAL(mat.column(0),e1);
-  CPPUNIT_ASSERT_EQUAL(mat.column(1),e2);
-  CPPUNIT_ASSERT_EQUAL(mat.column(2),e3);
+  CPPUNIT_ASSERT_EQUAL(mat.row(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(mat.row(1),unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(mat.row(2),unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(1),unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(mat.column(2),unitVec[2]);
 
   Vector t1=Vector(1.,1.,1.);
@@ -204,23 +204,23 @@
 
   res = (*diagonal)*(*perm1);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),2*e2);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
   res = (*diagonal)*(*perm2);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),3*e3);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),3*unitVec[2]);
   res = (*diagonal)*(*perm3);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),e1);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
   res = (*diagonal)*(*perm4);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),2*e2);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),e1);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),unitVec[0]);
   res = (*diagonal)*(*perm5);
-  CPPUNIT_ASSERT_EQUAL(res.column(0),3*e3);
-  CPPUNIT_ASSERT_EQUAL(res.column(1),e1);
-  CPPUNIT_ASSERT_EQUAL(res.column(2),2*e2);
+  CPPUNIT_ASSERT_EQUAL(res.column(0),3*unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(res.column(1),unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(res.column(2),2*unitVec[1]);
 }
 
@@ -261,42 +261,42 @@
 
 void MatrixUnittest::VecMultTest(){
-  CPPUNIT_ASSERT_EQUAL((*zero)*e1,zeroVec);
-  CPPUNIT_ASSERT_EQUAL((*zero)*e2,zeroVec);
-  CPPUNIT_ASSERT_EQUAL((*zero)*e3,zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[0],zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[1],zeroVec);
+  CPPUNIT_ASSERT_EQUAL((*zero)*unitVec[2],zeroVec);
   CPPUNIT_ASSERT_EQUAL((*zero)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*one)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*one)*e2,e2);
-  CPPUNIT_ASSERT_EQUAL((*one)*e3,e3);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[1],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*one)*unitVec[2],unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*one)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e2,2*e2);
-  CPPUNIT_ASSERT_EQUAL((*diagonal)*e3,3*e3);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[1],2*unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*diagonal)*unitVec[2],3*unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*diagonal)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e1,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e2,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm1)*e3,e2);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[0],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[1],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm1)*unitVec[2],unitVec[1]);
   CPPUNIT_ASSERT_EQUAL((*perm1)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e1,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e2,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm2)*e3,e3);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[0],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[1],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm2)*unitVec[2],unitVec[2]);
   CPPUNIT_ASSERT_EQUAL((*perm2)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e1,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e2,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm3)*e3,e1);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[0],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[1],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm3)*unitVec[2],unitVec[0]);
   CPPUNIT_ASSERT_EQUAL((*perm3)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e1,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e2,e2);
-  CPPUNIT_ASSERT_EQUAL((*perm4)*e3,e1);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[0],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[1],unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL((*perm4)*unitVec[2],unitVec[0]);
   CPPUNIT_ASSERT_EQUAL((*perm4)*zeroVec,zeroVec);
 
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e1,e3);
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e2,e1);
-  CPPUNIT_ASSERT_EQUAL((*perm5)*e3,e2);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[0],unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[1],unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL((*perm5)*unitVec[2],unitVec[1]);
   CPPUNIT_ASSERT_EQUAL((*perm5)*zeroVec,zeroVec);
 
Index: src/unittests/PlaneUnittest.cpp
===================================================================
--- src/unittests/PlaneUnittest.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/unittests/PlaneUnittest.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -24,8 +24,8 @@
 
 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);
+  p1 = new Plane(unitVec[0],unitVec[1],unitVec[2]);
+  p2 = new Plane(unitVec[0],unitVec[1],zeroVec);
+  p3 = new Plane(unitVec[0],zeroVec,unitVec[2]);
+  p4 = new Plane(zeroVec,unitVec[1],unitVec[2]);
 }
 
@@ -42,24 +42,24 @@
 
   // three points
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,e3));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],unitVec[2]));
   // when only two points are differnt this gives an error
-  CPPUNIT_ASSERT_THROW(Plane(e1,e2,e2),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[1],unitVec[1]),LinearDependenceException);
   // same with only one point
-  CPPUNIT_ASSERT_THROW(Plane(e1,e1,e1),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],unitVec[0]),LinearDependenceException);
 
   // use two vector giving two directions
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,e2,0));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],unitVec[1],0));
   // and again this is actually only one vector
-  CPPUNIT_ASSERT_THROW(Plane(e1,e1,0),LinearDependenceException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],unitVec[0],0),LinearDependenceException);
   // Zero vector does not give a good direction
-  CPPUNIT_ASSERT_THROW(Plane(e1,zeroVec,0),ZeroVectorException);
+  CPPUNIT_ASSERT_THROW(Plane(unitVec[0],zeroVec,0),ZeroVectorException);
 
   // use a normalvector and an scalar offset
-  CPPUNIT_ASSERT_NO_THROW(Plane(e1,0));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],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));
+  CPPUNIT_ASSERT_NO_THROW(Plane(unitVec[0],zeroVec));
   // and the bad zeroVector again
   CPPUNIT_ASSERT_THROW(Plane(zeroVec,zeroVec),ZeroVectorException);
@@ -75,28 +75,28 @@
   {
     // construct with three points on plane
-    Plane p1(e1,e2,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    Plane p1(unitVec[0],unitVec[1],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
 
-    Plane p2(e1,e3,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    Plane p2(unitVec[0],unitVec[2],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
 
-    Plane p3(e2,e3,zeroVec);
-    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    Plane p3(unitVec[1],unitVec[2],zeroVec);
+    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
   }
   {
     // construct with two directions + offset
-    Plane p1(e1,e2,0);
-    CPPUNIT_ASSERT(testNormal(e3,p1.getNormal()));
+    Plane p1(unitVec[0],unitVec[1],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[2],p1.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p1.getOffset());
 
-    Plane p2(e1,e3,0);
-    CPPUNIT_ASSERT(testNormal(e2,p2.getNormal()));
+    Plane p2(unitVec[0],unitVec[2],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[1],p2.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p2.getOffset());
 
-    Plane p3(e2,e3,0);
-    CPPUNIT_ASSERT(testNormal(e1,p3.getNormal()));
+    Plane p3(unitVec[1],unitVec[2],0);
+    CPPUNIT_ASSERT(testNormal(unitVec[0],p3.getNormal()));
     CPPUNIT_ASSERT_EQUAL(0.,p3.getOffset());
   }
@@ -145,15 +145,15 @@
 void PlaneUnittest::operationsTest(){
   {
-    Vector t = (1./3.)*(e1+e2+e3);
+    Vector t = (1./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
     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));
+  CPPUNIT_ASSERT(fabs(p2->distance(unitVec[2])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p2->getClosestPoint(unitVec[2]));
+  CPPUNIT_ASSERT(fabs(p3->distance(unitVec[1])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p3->getClosestPoint(unitVec[1]));
+  CPPUNIT_ASSERT(fabs(p4->distance(unitVec[0])-1) < MYEPSILON);
+  CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(unitVec[0]));
 }
 
@@ -162,45 +162,45 @@
 
   // 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 = p1->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p1->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = p1->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   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 = p2->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p2->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
 
   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 = p3->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[0]);
+  fixture = p3->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   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);
+  fixture = p4->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[1]);
+  fixture = p4->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,unitVec[2]);
 
   // some Vectors outside of the planes
   {
-    Vector t = (2./3.)*(e1+e2+e3);
+    Vector t = (2./3.)*(unitVec[0]+unitVec[1]+unitVec[2]);
     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);
+  fixture = p2->mirrorVector(unitVec[2]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[2]);
+  fixture = p3->mirrorVector(unitVec[1]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[1]);
+  fixture = p4->mirrorVector(unitVec[0]);
+  CPPUNIT_ASSERT_EQUAL(fixture,-1*unitVec[0]);
 }
 
@@ -209,10 +209,10 @@
   // 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_NO_THROW(fixture = Plane(unitVec[0], 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) );
+  Line l2 = makeLineThrough(unitVec[0],Vector(0,1,1));
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unitVec[1], Vector(2,1,0)).GetIntersection(l2) );
   CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
 }
Index: src/unittests/ShapeUnittest.cpp
===================================================================
--- src/unittests/ShapeUnittest.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/unittests/ShapeUnittest.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -17,4 +17,5 @@
 
 #include "Shapes/Shape.hpp"
+#include "Shapes/BaseShapes.hpp"
 #include "vector.cpp"
 
@@ -24,31 +25,31 @@
 void ShapeUnittest::setUp()
 {
-  v000 =  0*e1+0*e2+0*e3;
-  v100 =  1*e1+0*e2+0*e3;
-  v200 = -1*e1+0*e2+0*e3;
-  v010 =  0*e1+1*e2+0*e3;
-  v110 =  1*e1+1*e2+0*e3;
-  v210 = -1*e1+1*e2+0*e3;
-  v020 =  0*e1-1*e2+0*e3;
-  v120 =  1*e1-1*e2+0*e3;
-  v220 = -1*e1-1*e2+0*e3;
-  v001 =  0*e1+0*e2+1*e3;
-  v101 =  1*e1+0*e2+1*e3;
-  v201 = -1*e1+0*e2+1*e3;
-  v011 =  0*e1+1*e2+1*e3;
-  v111 =  1*e1+1*e2+1*e3;
-  v211 = -1*e1+1*e2+1*e3;
-  v021 =  0*e1-1*e2+1*e3;
-  v121 =  1*e1-1*e2+1*e3;
-  v221 = -1*e1-1*e2+1*e3;
-  v002 =  0*e1+0*e2-1*e3;
-  v102 =  1*e1+0*e2-1*e3;
-  v202 = -1*e1+0*e2-1*e3;
-  v012 =  0*e1+1*e2-1*e3;
-  v112 =  1*e1+1*e2-1*e3;
-  v212 = -1*e1+1*e2-1*e3;
-  v022 =  0*e1-1*e2-1*e3;
-  v122 =  1*e1-1*e2-1*e3;
-  v222 = -1*e1-1*e2-1*e3;
+  v000 =  0*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v100 =  1*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v200 = -1*unitVec[0]+0*unitVec[1]+0*unitVec[2];
+  v010 =  0*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v110 =  1*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v210 = -1*unitVec[0]+1*unitVec[1]+0*unitVec[2];
+  v020 =  0*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v120 =  1*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v220 = -1*unitVec[0]-1*unitVec[1]+0*unitVec[2];
+  v001 =  0*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v101 =  1*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v201 = -1*unitVec[0]+0*unitVec[1]+1*unitVec[2];
+  v011 =  0*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v111 =  1*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v211 = -1*unitVec[0]+1*unitVec[1]+1*unitVec[2];
+  v021 =  0*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v121 =  1*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v221 = -1*unitVec[0]-1*unitVec[1]+1*unitVec[2];
+  v002 =  0*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v102 =  1*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v202 = -1*unitVec[0]+0*unitVec[1]-1*unitVec[2];
+  v012 =  0*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v112 =  1*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v212 = -1*unitVec[0]+1*unitVec[1]-1*unitVec[2];
+  v022 =  0*unitVec[0]-1*unitVec[1]-1*unitVec[2];
+  v122 =  1*unitVec[0]-1*unitVec[1]-1*unitVec[2];
+  v222 = -1*unitVec[0]-1*unitVec[1]-1*unitVec[2];
 }
 
@@ -113,4 +114,146 @@
   CPPUNIT_ASSERT(Everywhere().isInside(v122));
   CPPUNIT_ASSERT(Everywhere().isInside(v222));
+
+  CPPUNIT_ASSERT(Cuboid().isInside(v000));
+  CPPUNIT_ASSERT(Cuboid().isInside(v100));
+  CPPUNIT_ASSERT(Cuboid().isInside(v200));
+  CPPUNIT_ASSERT(Cuboid().isInside(v010));
+  CPPUNIT_ASSERT(Cuboid().isInside(v110));
+  CPPUNIT_ASSERT(Cuboid().isInside(v210));
+  CPPUNIT_ASSERT(Cuboid().isInside(v020));
+  CPPUNIT_ASSERT(Cuboid().isInside(v120));
+  CPPUNIT_ASSERT(Cuboid().isInside(v220));
+  CPPUNIT_ASSERT(Cuboid().isInside(v001));
+  CPPUNIT_ASSERT(Cuboid().isInside(v101));
+  CPPUNIT_ASSERT(Cuboid().isInside(v201));
+  CPPUNIT_ASSERT(Cuboid().isInside(v011));
+  CPPUNIT_ASSERT(Cuboid().isInside(v111));
+  CPPUNIT_ASSERT(Cuboid().isInside(v211));
+  CPPUNIT_ASSERT(Cuboid().isInside(v021));
+  CPPUNIT_ASSERT(Cuboid().isInside(v121));
+  CPPUNIT_ASSERT(Cuboid().isInside(v221));
+  CPPUNIT_ASSERT(Cuboid().isInside(v002));
+  CPPUNIT_ASSERT(Cuboid().isInside(v102));
+  CPPUNIT_ASSERT(Cuboid().isInside(v202));
+  CPPUNIT_ASSERT(Cuboid().isInside(v012));
+  CPPUNIT_ASSERT(Cuboid().isInside(v112));
+  CPPUNIT_ASSERT(Cuboid().isInside(v212));
+  CPPUNIT_ASSERT(Cuboid().isInside(v022));
+  CPPUNIT_ASSERT(Cuboid().isInside(v122));
+  CPPUNIT_ASSERT(Cuboid().isInside(v222));
+
+  CPPUNIT_ASSERT(Sphere().isInside(v000));
+  CPPUNIT_ASSERT(Sphere().isInside(v100));
+  CPPUNIT_ASSERT(Sphere().isInside(v200));
+  CPPUNIT_ASSERT(Sphere().isInside(v010));
+  CPPUNIT_ASSERT(!Sphere().isInside(v110));
+  CPPUNIT_ASSERT(!Sphere().isInside(v210));
+  CPPUNIT_ASSERT(Sphere().isInside(v020));
+  CPPUNIT_ASSERT(!Sphere().isInside(v120));
+  CPPUNIT_ASSERT(!Sphere().isInside(v220));
+  CPPUNIT_ASSERT(Sphere().isInside(v001));
+  CPPUNIT_ASSERT(!Sphere().isInside(v101));
+  CPPUNIT_ASSERT(!Sphere().isInside(v201));
+  CPPUNIT_ASSERT(!Sphere().isInside(v011));
+  CPPUNIT_ASSERT(!Sphere().isInside(v111));
+  CPPUNIT_ASSERT(!Sphere().isInside(v211));
+  CPPUNIT_ASSERT(!Sphere().isInside(v021));
+  CPPUNIT_ASSERT(!Sphere().isInside(v121));
+  CPPUNIT_ASSERT(!Sphere().isInside(v221));
+  CPPUNIT_ASSERT(Sphere().isInside(v002));
+  CPPUNIT_ASSERT(!Sphere().isInside(v102));
+  CPPUNIT_ASSERT(!Sphere().isInside(v202));
+  CPPUNIT_ASSERT(!Sphere().isInside(v012));
+  CPPUNIT_ASSERT(!Sphere().isInside(v112));
+  CPPUNIT_ASSERT(!Sphere().isInside(v212));
+  CPPUNIT_ASSERT(!Sphere().isInside(v022));
+  CPPUNIT_ASSERT(!Sphere().isInside(v122));
+  CPPUNIT_ASSERT(!Sphere().isInside(v222));
+}
+
+void ShapeUnittest::surfaceTest(){
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v000));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v100));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v200));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v010));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v110));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v210));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v020));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v120));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v220));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v001));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v101));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v201));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v011));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v111));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v211));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v021));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v121));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v221));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v002));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v102));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v202));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v012));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v112));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v212));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v022));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v122));
+  CPPUNIT_ASSERT(!Nowhere().isOnSurface(v222));
+
+  CPPUNIT_ASSERT(!Cuboid().isOnSurface(v000));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v100));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v200));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v010));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v110));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v210));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v020));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v120));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v220));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v001));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v101));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v201));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v011));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v111));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v211));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v021));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v121));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v221));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v002));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v102));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v202));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v012));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v112));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v212));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v022));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v122));
+  CPPUNIT_ASSERT(Cuboid().isOnSurface(v222));
+
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v000));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v100));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v200));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v010));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v110));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v210));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v020));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v120));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v220));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v001));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v101));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v201));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v011));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v111));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v211));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v021));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v121));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v221));
+  CPPUNIT_ASSERT(Sphere().isOnSurface(v002));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v102));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v202));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v012));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v112));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v212));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v022));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v122));
+  CPPUNIT_ASSERT(!Sphere().isOnSurface(v222));
 }
 
@@ -360,3 +503,122 @@
   }
 
+  {
+    Shape s1 = Sphere() && Cuboid(); // should be the same as Sphere
+
+    CPPUNIT_ASSERT(s1.isInside(v000));
+    CPPUNIT_ASSERT(s1.isInside(v100));
+    CPPUNIT_ASSERT(s1.isInside(v200));
+    CPPUNIT_ASSERT(s1.isInside(v010));
+    CPPUNIT_ASSERT(!s1.isInside(v110));
+    CPPUNIT_ASSERT(!s1.isInside(v210));
+    CPPUNIT_ASSERT(s1.isInside(v020));
+    CPPUNIT_ASSERT(!s1.isInside(v120));
+    CPPUNIT_ASSERT(!s1.isInside(v220));
+    CPPUNIT_ASSERT(s1.isInside(v001));
+    CPPUNIT_ASSERT(!s1.isInside(v101));
+    CPPUNIT_ASSERT(!s1.isInside(v201));
+    CPPUNIT_ASSERT(!s1.isInside(v011));
+    CPPUNIT_ASSERT(!s1.isInside(v111));
+    CPPUNIT_ASSERT(!s1.isInside(v211));
+    CPPUNIT_ASSERT(!s1.isInside(v021));
+    CPPUNIT_ASSERT(!s1.isInside(v121));
+    CPPUNIT_ASSERT(!s1.isInside(v221));
+    CPPUNIT_ASSERT(s1.isInside(v002));
+    CPPUNIT_ASSERT(!s1.isInside(v102));
+    CPPUNIT_ASSERT(!s1.isInside(v202));
+    CPPUNIT_ASSERT(!s1.isInside(v012));
+    CPPUNIT_ASSERT(!s1.isInside(v112));
+    CPPUNIT_ASSERT(!s1.isInside(v212));
+    CPPUNIT_ASSERT(!s1.isInside(v022));
+    CPPUNIT_ASSERT(!s1.isInside(v122));
+    CPPUNIT_ASSERT(!s1.isInside(v222));
+
+    CPPUNIT_ASSERT(!s1.isOnSurface(v000));
+    CPPUNIT_ASSERT(s1.isOnSurface(v100));
+    CPPUNIT_ASSERT(s1.isOnSurface(v200));
+    CPPUNIT_ASSERT(s1.isOnSurface(v010));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v110));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v210));
+    CPPUNIT_ASSERT(s1.isOnSurface(v020));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v120));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v220));
+    CPPUNIT_ASSERT(s1.isOnSurface(v001));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v101));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v201));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v011));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v111));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v211));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v021));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v121));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v221));
+    CPPUNIT_ASSERT(s1.isOnSurface(v002));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v102));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v202));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v012));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v112));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v212));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v022));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v122));
+    CPPUNIT_ASSERT(!s1.isOnSurface(v222));
+
+
+    Shape s2 = Sphere() || Cuboid(); // Should be same as Cuboid
+
+    CPPUNIT_ASSERT(s2.isInside(v000));
+    CPPUNIT_ASSERT(s2.isInside(v100));
+    CPPUNIT_ASSERT(s2.isInside(v200));
+    CPPUNIT_ASSERT(s2.isInside(v010));
+    CPPUNIT_ASSERT(s2.isInside(v110));
+    CPPUNIT_ASSERT(s2.isInside(v210));
+    CPPUNIT_ASSERT(s2.isInside(v020));
+    CPPUNIT_ASSERT(s2.isInside(v120));
+    CPPUNIT_ASSERT(s2.isInside(v220));
+    CPPUNIT_ASSERT(s2.isInside(v001));
+    CPPUNIT_ASSERT(s2.isInside(v101));
+    CPPUNIT_ASSERT(s2.isInside(v201));
+    CPPUNIT_ASSERT(s2.isInside(v011));
+    CPPUNIT_ASSERT(s2.isInside(v111));
+    CPPUNIT_ASSERT(s2.isInside(v211));
+    CPPUNIT_ASSERT(s2.isInside(v021));
+    CPPUNIT_ASSERT(s2.isInside(v121));
+    CPPUNIT_ASSERT(s2.isInside(v221));
+    CPPUNIT_ASSERT(s2.isInside(v002));
+    CPPUNIT_ASSERT(s2.isInside(v102));
+    CPPUNIT_ASSERT(s2.isInside(v202));
+    CPPUNIT_ASSERT(s2.isInside(v012));
+    CPPUNIT_ASSERT(s2.isInside(v112));
+    CPPUNIT_ASSERT(s2.isInside(v212));
+    CPPUNIT_ASSERT(s2.isInside(v022));
+    CPPUNIT_ASSERT(s2.isInside(v122));
+    CPPUNIT_ASSERT(s2.isInside(v222));
+
+    CPPUNIT_ASSERT(!s2.isOnSurface(v000));
+    CPPUNIT_ASSERT(s2.isOnSurface(v100));
+    CPPUNIT_ASSERT(s2.isOnSurface(v200));
+    CPPUNIT_ASSERT(s2.isOnSurface(v010));
+    CPPUNIT_ASSERT(s2.isOnSurface(v110));
+    CPPUNIT_ASSERT(s2.isOnSurface(v210));
+    CPPUNIT_ASSERT(s2.isOnSurface(v020));
+    CPPUNIT_ASSERT(s2.isOnSurface(v120));
+    CPPUNIT_ASSERT(s2.isOnSurface(v220));
+    CPPUNIT_ASSERT(s2.isOnSurface(v001));
+    CPPUNIT_ASSERT(s2.isOnSurface(v101));
+    CPPUNIT_ASSERT(s2.isOnSurface(v201));
+    CPPUNIT_ASSERT(s2.isOnSurface(v011));
+    CPPUNIT_ASSERT(s2.isOnSurface(v111));
+    CPPUNIT_ASSERT(s2.isOnSurface(v211));
+    CPPUNIT_ASSERT(s2.isOnSurface(v021));
+    CPPUNIT_ASSERT(s2.isOnSurface(v121));
+    CPPUNIT_ASSERT(s2.isOnSurface(v221));
+    CPPUNIT_ASSERT(s2.isOnSurface(v002));
+    CPPUNIT_ASSERT(s2.isOnSurface(v102));
+    CPPUNIT_ASSERT(s2.isOnSurface(v202));
+    CPPUNIT_ASSERT(s2.isOnSurface(v012));
+    CPPUNIT_ASSERT(s2.isOnSurface(v112));
+    CPPUNIT_ASSERT(s2.isOnSurface(v212));
+    CPPUNIT_ASSERT(s2.isOnSurface(v022));
+    CPPUNIT_ASSERT(s2.isOnSurface(v122));
+    CPPUNIT_ASSERT(s2.isOnSurface(v222));
+  }
+
 }
Index: src/unittests/ShapeUnittest.hpp
===================================================================
--- src/unittests/ShapeUnittest.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/unittests/ShapeUnittest.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -20,4 +20,5 @@
   CPPUNIT_TEST_SUITE( ShapeUnittest) ;
   CPPUNIT_TEST ( baseShapesTest );
+  CPPUNIT_TEST ( surfaceTest );
   CPPUNIT_TEST ( assignmentTest );
   CPPUNIT_TEST ( operatorTest );
@@ -29,4 +30,5 @@
 
   void baseShapesTest();
+  void surfaceTest();
   void assignmentTest();
   void operatorTest();
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/vector.cpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -193,4 +193,10 @@
   (*this) *= 1/factor;
 };
+
+Vector Vector::getNormalized() const{
+  Vector res= *this;
+  res.Normalize();
+  return res;
+}
 
 /** Zeros all components of this vector.
@@ -534,5 +540,3 @@
 // some comonly used vectors
 const Vector zeroVec(0,0,0);
-const Vector e1(1,0,0);
-const Vector e2(0,1,0);
-const Vector e3(0,0,1);
+const Vector unitVec[NDIM]={Vector(1,0,0),Vector(0,1,0),Vector(0,0,1)};
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision a7b777c63b3b50eee5da3df0edbd51b6e7eaad91)
+++ src/vector.hpp	(revision 2ad48247b6f9de9eea8821024e501debd95f776e)
@@ -80,4 +80,5 @@
   double NormSquared() const;
   void Normalize();
+  Vector getNormalized() const;
   void Zero();
   void One(const double one);
@@ -106,7 +107,5 @@
 // some commonly used and fixed vectors
 const extern Vector zeroVec;
-const extern Vector e1;
-const extern Vector e2;
-const extern Vector e3;
+const extern Vector unitVec[NDIM];
 
 ostream & operator << (ostream& ost, const Vector &m);
