Index: src/LinearAlgebra/Line.cpp
===================================================================
--- src/LinearAlgebra/Line.cpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Line.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -50,4 +50,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;
+}
 
 
@@ -224,4 +232,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);
@@ -247,4 +282,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){
@@ -252,4 +305,49 @@
   }
   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/LinearAlgebra/Line.hpp
===================================================================
--- src/LinearAlgebra/Line.hpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Line.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -17,11 +17,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;
@@ -36,4 +41,6 @@
 
   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;
@@ -41,8 +48,16 @@
   std::vector<Vector> getSphereIntersections() const;
 
+  LinePoint getLinePoint(const Vector&) const;
+  LinePoint posEndpoint() const;
+  LinePoint negEndpoint() const;
+
+
+
 private:
   std::auto_ptr<Vector> origin;
   std::auto_ptr<Vector> direction;
 };
+
+bool operator==(const Line&,const Line&);
 
 std::ostream & operator << (std::ostream& ost, const Line &m);
@@ -53,3 +68,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/LinearAlgebra/LineSegment.cpp
===================================================================
--- src/LinearAlgebra/LineSegment.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
+++ src/LinearAlgebra/LineSegment.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -0,0 +1,90 @@
+/*
+ * LineSegment.cpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#include "LinearAlgebra/LineSegment.hpp"
+#include "Helpers/Assert.hpp"
+
+#include "Line.hpp"
+#include "LinearAlgebra/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/LinearAlgebra/LineSegment.hpp
===================================================================
--- src/LinearAlgebra/LineSegment.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
+++ src/LinearAlgebra/LineSegment.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -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/LinearAlgebra/LineSegmentSet.cpp
===================================================================
--- src/LinearAlgebra/LineSegmentSet.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
+++ src/LinearAlgebra/LineSegmentSet.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -0,0 +1,223 @@
+/*
+ * LineSegmentSet.cpp
+ *
+ *  Created on: Jul 22, 2010
+ *      Author: crueger
+ */
+
+#include "LinearAlgebra/LineSegmentSet.hpp"
+
+#include "Helpers/Assert.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/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/LinearAlgebra/LineSegmentSet.hpp
===================================================================
--- src/LinearAlgebra/LineSegmentSet.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
+++ src/LinearAlgebra/LineSegmentSet.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -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/LinearAlgebra/Makefile.am
===================================================================
--- src/LinearAlgebra/Makefile.am	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Makefile.am	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -14,4 +14,6 @@
   gslvector.cpp \
   Line.cpp \
+  LineSegment.cpp \
+  LineSegmentSet.cpp \
   linearsystemofequations.cpp \
   Matrix.cpp \
@@ -27,4 +29,6 @@
   gslvector.hpp \
   Line.hpp \
+  LineSegment.hpp \
+  LineSegmentSet.hpp \
   linearsystemofequations.hpp \
   Matrix.hpp \
Index: src/LinearAlgebra/Matrix.cpp
===================================================================
--- src/LinearAlgebra/Matrix.cpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Matrix.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -230,24 +230,4 @@
        - at(2,1)*at(1,2)*at(0,0)
        - at(2,2)*at(1,0)*at(0,1);
-}
-
-Matrix Matrix::transpose() const
-{
-  Matrix copy(*this);
-  copy.transpose();
-  return copy;
-}
-
-
-Matrix &Matrix::transpose()
-{
-  double tmp;
-  for (int i=0;i<NDIM;i++)
-    for (int j=i+1;j<NDIM;j++) {
-      tmp = at(j,i);
-      at(j,i) = at(i,j);
-      at(i,j) = tmp;
-    }
-  return *this;
 }
 
@@ -273,4 +253,25 @@
 }
 
+Matrix Matrix::transpose() const{
+  MatrixContent *newContent = new MatrixContent();
+  newContent->content = gsl_matrix_alloc(NDIM, NDIM);
+  gsl_matrix_transpose_memcpy(newContent->content, content->content);
+  Matrix res = Matrix(newContent);
+  return res;
+}
+
+Matrix &Matrix::transpose()
+{
+  double tmp;
+  for (int i=0;i<NDIM;i++)
+    for (int j=i+1;j<NDIM;j++) {
+      tmp = at(j,i);
+      at(j,i) = at(i,j);
+      at(i,j) = tmp;
+    }
+  return *this;
+}
+
+
 Vector Matrix::transformToEigenbasis()
 {
Index: src/LinearAlgebra/Plane.cpp
===================================================================
--- src/LinearAlgebra/Plane.cpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Plane.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -27,4 +27,5 @@
 #include "Helpers/Verbose.hpp"
 #include "Helpers/Assert.hpp"
+#include "Helpers/helpers.hpp"
 #include <cmath>
 #include "LinearAlgebra/Line.hpp"
@@ -118,4 +119,12 @@
 Plane::~Plane()
 {}
+
+Plane &Plane::operator=(const Plane &rhs){
+  if(&rhs!=this){
+    normalVector.reset(new Vector(*rhs.normalVector));
+    offset = rhs.offset;
+  }
+  return *this;
+}
 
 
@@ -200,4 +209,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 ****************/
 
@@ -219,4 +233,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/LinearAlgebra/Plane.hpp
===================================================================
--- src/LinearAlgebra/Plane.hpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Plane.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -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/LinearAlgebra/Vector.cpp
===================================================================
--- src/LinearAlgebra/Vector.cpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Vector.cpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -237,4 +237,10 @@
   (*this) *= 1/factor;
 };
+
+Vector Vector::getNormalized() const{
+  Vector res= *this;
+  res.Normalize();
+  return res;
+}
 
 /** Zeros all components of this vector.
@@ -578,5 +584,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/LinearAlgebra/Vector.hpp
===================================================================
--- src/LinearAlgebra/Vector.hpp	(revision fd49053e27109f166b1c159007a0e33103499934)
+++ src/LinearAlgebra/Vector.hpp	(revision 6c438f95957e5e36d17eaf3a715bb46a41596300)
@@ -83,4 +83,5 @@
   double NormSquared() const;
   void Normalize();
+  Vector getNormalized() const;
   void Zero();
   void One(const double one);
@@ -109,7 +110,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);
