Changes in / [a7c344:b32dbb]

Location:

src

Files:

• Exceptions/MultipleSolutionsException.hpp (deleted)
• Exceptions/SkewException.cpp (deleted)
• Exceptions/SkewException.hpp (deleted)
• Legacy/oldmenu.cpp (modified) (3 diffs)
• Line.cpp (modified) (2 diffs)
• Line.hpp (modified) (2 diffs)
• Makefile.am (modified) (2 diffs)
• Plane.cpp (modified) (4 diffs)
• Plane.hpp (modified) (2 diffs)
• Space.cpp (modified) (1 diff)
• Space.hpp (modified) (1 diff)
• atom.cpp (modified) (1 diff)
• boundary.cpp (modified) (1 diff)
• molecule_geometry.cpp (modified) (2 diffs)
• periodentafel.cpp (modified) (6 diffs)
• tesselation.cpp (modified) (9 diffs)
• tesselation.hpp (modified) (2 diffs)
• tesselationhelpers.cpp (modified) (2 diffs)
• unittests/LineUnittest.cpp (deleted)
• unittests/LineUnittest.hpp (deleted)
• unittests/Makefile.am (modified) (4 diffs)
• unittests/PlaneUnittest.cpp (modified) (2 diffs)
• unittests/PlaneUnittest.hpp (modified) (2 diffs)
• unittests/vectorunittest.cpp (modified) (1 diff)
• unittests/vectorunittest.hpp (modified) (1 diff)
• vector.cpp (modified) (5 diffs)
• vector.hpp (modified) (5 diffs)
• vector_ops.cpp (modified) (2 diffs)
• vector_ops.hpp (modified) (1 diff)

src/Legacy/oldmenu.cpp

@@ -27 +27 @@
 #include "vector_ops.hpp"
 #include "Plane.hpp"
-#include "Line.hpp"
 
 #include "UIElements/UIFactory.hpp"
@@ -186 +185 @@
           // rotate vector around first angle
           first->x = x;
-          first->x = Line(zeroVec,z).rotateVector(first->x,b - M_PI);
+          first->x = RotateVector(first->x,z,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
           // remove the projection onto the rotation plane of the second angle
@@ -202 +201 @@
           // rotate another vector around second angle
           n = y;
-          n = Line(zeroVec,x).rotateVector(n,c - M_PI);
+          n = RotateVector(n,x,c - M_PI);
           Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 

src/Line.cpp

@@ -11 +11 @@
 
 #include "vector.hpp"
-#include "log.hpp"
-#include "verbose.hpp"
-#include "gslmatrix.hpp"
-#include "info.hpp"
-#include "Exceptions/LinearDependenceException.hpp"
-#include "Exceptions/SkewException.hpp"
-#include "Plane.hpp"
 
-using namespace std;
-
-Line::Line(const Vector &_origin, const Vector &_direction) :
+Line::Line(Vector &_origin, Vector &_direction) :
+  origin(new Vector(_origin)),
   direction(new Vector(_direction))
 {
   direction->Normalize();
-  origin.reset(new Vector(_origin.partition(*direction).second));
 }
-
-Line::Line(const Line &src) :
-  origin(new Vector(*src.origin)),
-  direction(new Vector(*src.direction))
-{}
 
 Line::~Line()
@@ -38 +24 @@
 
 double Line::distance(const Vector &point) const{
-  // get any vector from line to point
-  Vector helper = point - *origin;
-  // partition this vector along direction
-  // the residue points from the line to the point
-  return helper.partition(*direction).second.Norm();
+  return fabs(point.ScalarProduct(*direction) - origin->ScalarProduct(*direction));
 }
 
 Vector Line::getClosestPoint(const Vector &point) const{
-  // get any vector from line to point
-  Vector helper = point - *origin;
-  // partition this vector along direction
-  // add only the part along the direction
-  return *origin + helper.partition(*direction).first;
+  double factor = point.ScalarProduct(*direction) - origin->ScalarProduct(*direction);
+  return (point - (factor * (*direction)));
 }
-
-Vector Line::getDirection() const{
-  return *direction;
-}
-
-Vector Line::getOrigin() const{
-  return *origin;
-}
-
-vector<Vector> Line::getPointsOnLine() const{
-  vector<Vector> res;
-  res.reserve(2);
-  res.push_back(*origin);
-  res.push_back(*origin+*direction);
-  return res;
-}
-
-/** Calculates the intersection of the two lines that are both on the same plane.
- * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
- * \param *out output stream for debugging
- * \param *Line1a first vector of first line
- * \param *Line1b second vector of first line
- * \param *Line2a first vector of second line
- * \param *Line2b second vector of second line
- * \return true - \a this will contain the intersection on return, false - lines are parallel
- */
-Vector Line::getIntersection(const Line& otherLine) const{
-  Info FunctionInfo(__func__);
-
-  pointset line1Points = getPointsOnLine();
-
-  Vector Line1a = line1Points[0];
-  Vector Line1b = line1Points[1];
-
-  pointset line2Points = otherLine.getPointsOnLine();
-
-  Vector Line2a = line2Points[0];
-  Vector Line2b = line2Points[1];
-
-  Vector res;
-
-  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
-
-  M->SetAll(1.);
-  for (int i=0;i<3;i++) {
-    M->Set(0, i, Line1a[i]);
-    M->Set(1, i, Line1b[i]);
-    M->Set(2, i, Line2a[i]);
-    M->Set(3, i, Line2b[i]);
-  }
-
-  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    cout << "\t" << M->Get(i,j);
-  //  cout << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
-    throw SkewException(__FILE__,__LINE__);
-  }
-
-  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
-
-
-  // constuct a,b,c
-  Vector a = Line1b - Line1a;
-  Vector b = Line2b - Line2a;
-  Vector c = Line2a - Line1a;
-  Vector d = Line2b - Line1b;
-  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   res.Zero();
-   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
-   throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // check for parallelity
-  Vector parallel;
-  double factor = 0.;
-  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
-    parallel = Line1a - Line2a;
-    factor = parallel.ScalarProduct(a)/a.Norm();
-    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-      res = Line2a;
-      Log() << Verbose(1) << "Lines conincide." << endl;
-      return res;
-    } else {
-      parallel = Line1a - Line2b;
-      factor = parallel.ScalarProduct(a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        res = Line2b;
-        Log() << Verbose(1) << "Lines conincide." << endl;
-        return res;
-      }
-    }
-    Log() << Verbose(1) << "Lines are parallel." << endl;
-    res.Zero();
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-
-  // obtain s
-  double s;
-  Vector temp1, temp2;
-  temp1 = c;
-  temp1.VectorProduct(b);
-  temp2 = a;
-  temp2.VectorProduct(b);
-  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
-
-  // construct intersection
-  res = a;
-  res.Scale(s);
-  res += Line1a;
-  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
-
-  return res;
-}
-
-/** Rotates the vector by an angle of \a alpha around this line.
- * \param rhs Vector to rotate
- * \param alpha rotation angle in radian
- */
-Vector Line::rotateVector(const Vector &rhs, double alpha) const{
-  Vector helper = rhs;
-
-  // translate the coordinate system so that the line goes through (0,0,0)
-  helper -= *origin;
-
-  // partition the vector into a part that gets rotated and a part that lies along the line
-  pair<Vector,Vector> parts = helper.partition(*direction);
-
-  // we just keep anything that is along the axis
-  Vector res = parts.first;
-
-  // the rest has to be rotated
-  Vector a = parts.second;
-  // we only have to do the rest, if we actually could partition the vector
-  if(!a.IsZero()){
-    // construct a vector that is orthogonal to a and direction and has length |a|
-    Vector y = a;
-    // direction is normalized, so the result has length |a|
-    y.VectorProduct(*direction);
-
-    res += cos(alpha) * a + sin(alpha) * y;
-  }
-
-  // translate the coordinate system back
-  res += *origin;
-  return res;
-}
-
-Plane Line::getOrthogonalPlane(const Vector &origin) const{
-  return Plane(getDirection(),origin);
-}
-
-Line makeLineThrough(const Vector &x1, const Vector &x2){
-  if(x1==x2){
-    throw LinearDependenceException(__FILE__,__LINE__);
-  }
-  return Line(x1,x1-x2);
-}

src/Line.hpp

@@ -12 +12 @@
 
 #include <memory>
-#include <vector>
 
 class Vector;
-class Plane;
 
 class Line : public Space
 {
 public:
-  Line(const Vector &_origin, const Vector &_direction);
-  Line(const Line& _src);
+  Line(Vector &_origin, Vector &_direction);
   virtual ~Line();
 
-  virtual double distance(const Vector &point) const;
-  virtual Vector getClosestPoint(const Vector &point) const;
-
-  Vector getDirection() const;
-  Vector getOrigin() const;
-
-  std::vector<Vector> getPointsOnLine() const;
-
-  Vector getIntersection(const Line& otherLine) const;
-
-  Vector rotateVector(const Vector &rhs, double alpha) const;
-
-  Plane getOrthogonalPlane(const Vector &origin) const;
+  virtual double distance(const Vector &point) const=0;
+  virtual Vector getClosestPoint(const Vector &point) const=0;
 
 private:
@@ -43 +29 @@
 };
 
-/**
- * Named constructor to make a line through two points
- */
-Line makeLineThrough(const Vector &x1, const Vector &x2);
-
 #endif /* LINE_HPP_ */

src/Makefile.am

@@ -88 +88 @@
                                   Exceptions/LinearDependenceException.cpp \
                                   Exceptions/MathException.cpp \
-                                  Exceptions/SkewException.cpp \
                                   Exceptions/ZeroVectorException.cpp
                                   
@@ -94 +93 @@
                                   Exceptions/LinearDependenceException.hpp \
                                   Exceptions/MathException.hpp \
-                                  Exceptions/SkewException.hpp \
                                   Exceptions/ZeroVectorException.hpp
 

src/Plane.cpp

@@ -14 +14 @@
 #include "Helpers/Assert.hpp"
 #include <cmath>
-#include "Line.hpp"
-#include "Exceptions/MultipleSolutionsException.hpp"
 
 /**
@@ -44 +42 @@
 /**
  * Constructs a plane from two direction vectors and a offset.
+ * If no offset is given a plane through origin is assumed
  */
 Plane::Plane(const Vector &y1, const Vector &y2, double _offset) throw(ZeroVectorException,LinearDependenceException) :
@@ -114 +113 @@
 }
 
-Vector Plane::getOffsetVector() const {
+Vector Plane::getOffsetVector() {
   return getOffset()*getNormal();
 }
 
-vector<Vector> Plane::getPointsOnPlane() const{
+vector<Vector> Plane::getPointsOnPlane(){
   std::vector<Vector> res;
   res.reserve(3);
@@ -148 +147 @@
  * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
  */
-Vector Plane::GetIntersection(const Line& line) const
+Vector Plane::GetIntersection(const Vector &Origin, const Vector &LineVector)
 {
   Info FunctionInfo(__func__);
   Vector res;
 
-  double factor1 = getNormal().ScalarProduct(line.getDirection());
-  if(fabs(factor1)<MYEPSILON){
-    // the plane is parallel... under all circumstances this is bad luck
-    // we no have either no or infinite solutions
-    if(isContained(line.getOrigin())){
-      throw MultipleSolutionsException<Vector>(__FILE__,__LINE__,line.getOrigin());
-    }
-    else{
-      throw LinearDependenceException(__FILE__,__LINE__);
-    }
-  }
-
-  double factor2 = getNormal().ScalarProduct(line.getOrigin());
+  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
+  Vector Direction = LineVector - Origin;
+  Direction.Normalize();
+  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
+  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
+  double factor1 = Direction.ScalarProduct(*normalVector.get());
+  if (fabs(factor1) < MYEPSILON) { // Uniqueness: line parallel to plane?
+    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  double factor2 = Origin.ScalarProduct(*normalVector.get());
+  if (fabs(factor2-offset) < MYEPSILON) { // Origin is in-plane
+    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
+    res = Origin;
+    return res;
+  }
+
   double scaleFactor = (offset-factor2)/factor1;
 
-  res = line.getOrigin() + scaleFactor * line.getDirection();
-
-  // tests to make sure the resulting vector really is on plane and line
-  ASSERT(isContained(res),"Calculated line-Plane intersection does not lie on plane.");
-  ASSERT(line.isContained(res),"Calculated line-Plane intersection does not lie on line.");
+  //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
+  Direction.Scale(scaleFactor);
+  res = Origin + Direction;
+  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
+
+  // test whether resulting vector really is on plane
+  ASSERT(fabs(res.ScalarProduct(*normalVector) - offset) < MYEPSILON,
+         "Calculated line-Plane intersection does not lie on plane.");
   return res;
 };
-
-Vector Plane::mirrorVector(const Vector &rhs) const {
-  Vector helper = getVectorToPoint(rhs);
-  // substract twice the Vector to the plane
-  return rhs+2*helper;
-}
-
-Line Plane::getOrthogonalLine(const Vector &origin) const{
-  return Line(origin,getNormal());
-}
 
 /************ Methods inherited from Space ****************/

src/Plane.hpp

@@ -17 +17 @@
 
 class Vector;
-class Line;
 
 class Plane : public Space
@@ -43 +42 @@
    * Same as getOffset()*getNormal();
    */
-  Vector getOffsetVector() const;
+  Vector getOffsetVector();
 
   /**
    * returns three seperate points on this plane
    */
-  std::vector<Vector> getPointsOnPlane() const;
+  std::vector<Vector> getPointsOnPlane();
 
   // some calculations
-  Vector GetIntersection(const Line &Line) const;
-
-  Vector mirrorVector(const Vector &rhs) const;
-
-  /**
-   * get a Line that is orthogonal to this plane, going through a chosen
-   * point.
-   *
-   * The point does not have to lie on the plane itself.
-   */
-  Line getOrthogonalLine(const Vector &origin) const;
+  Vector GetIntersection(const Vector &Origin, const Vector &LineVector);
 
   /****** Methods inherited from Space ***********/

src/Space.cpp

@@ -19 +19 @@
 }
 
-Vector Space::getVectorToPoint(const Vector &origin) const{
-  Vector support = getClosestPoint(origin);
-  return support-origin;
-}
-
 bool Space::isContained(const Vector &point) const{
   return (distance(point)) < MYEPSILON;

src/Space.hpp

@@ -17 +17 @@
   virtual ~Space();
 
-  /**
-   * Calculates the distance between a Space and a Vector.
-   */
   virtual double distance(const Vector &point) const=0;
-
-  /**
-   * get the closest point inside the space to another point
-   */
   virtual Vector getClosestPoint(const Vector &point) const=0;
-
-  /**
-   * get the shortest Vector from a point to a Space.
-   *
-   * The Vector always points from the given Vector to the point in space
-   * returned by Plane::getClosestPoint().
-   */
-  virtual Vector getVectorToPoint(const Vector &point) const;
-
-  /**
-   * Test wether a point is contained in the space.
-   *
-   * returns true, when the point lies inside and false
-   * otherwise.
-   */
   virtual bool isContained(const Vector &point) const;
-
-  /**
-   * Tests if this space contains the center of the coordinate system.
-   */
   virtual bool hasZero() const;
 

src/atom.cpp

@@ -68 +68 @@
 atom *atom::GetTrueFather()
 {
-  if(father == this){ // top most father is the one that points on itself
-    return this;
-  }
-  else if(!father) {
-    return 0;
-  }
-  else {
-    return father->GetTrueFather();
-  }
+  atom *walker = this;
+  do {
+    if (walker == walker->father) // top most father is the one that points on itself
+      break;
+    walker = walker->father;
+  } while (walker != NULL);
+  return walker;
 };
 

src/boundary.cpp

@@ -620 +620 @@
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
-      x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1);
-      y = runner->second->getEndpoint(0) - runner->second->getEndpoint(2);
-      const double a = x.Norm();
-      const double b = y.Norm();
-      const double c = runner->second->getEndpoint(2).distance(runner->second->getEndpoint(1));
+      x = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[1]->node->node);
+      y = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[2]->node->node);
+      const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
+      const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[2]->node->node));
+      const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
       const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
-      x = runner->second->getPlane().getNormal();
-      x.Scale(runner->second->getEndpoint(1).ScalarProduct(x));
+      x = Plane(*(runner->second->endpoints[0]->node->node),
+                *(runner->second->endpoints[1]->node->node),
+                *(runner->second->endpoints[2]->node->node)).getNormal();
+      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(x));
       const double h = x.Norm(); // distance of CoG to triangle
       const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)

src/molecule_geometry.cpp

@@ -16 +16 @@
 #include "molecule.hpp"
 #include "World.hpp"
-#include "Plane.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -252 +251 @@
 void molecule::Mirror(const Vector *n)
 {
-  Plane p(*n,0);
-  // TODO: replace with simpler construct (e.g. Boost::foreach)
-  // once the structure of the atom list is fully reworked
-  atom *Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    (*Walker->node) = p.mirrorVector(*Walker->node);
-  }
+  ActOnAllVectors( &Vector::Mirror, *n );
 };
 

src/periodentafel.cpp

@@ -212 +212 @@
 
   // fill elements DB
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDELEMENTSDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -256 +258 @@
 
   // fill valence DB per element
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDVALENCEDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDVALENCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -273 +277 @@
 
   // fill valence DB per element
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDORBITALDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDORBITALDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -290 +296 @@
 
   // fill H-BondDistance DB per element
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDHBONDDISTANCEDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDDISTANCEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -310 +318 @@
 
   // fill H-BondAngle DB per element
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDHBONDANGLEDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDHBONDANGLEDB, MAXSTRINGSIZE-strlen(filename));
   infile.open(filename);
   if (infile != NULL) {
@@ -353 +363 @@
   char filename[MAXSTRINGSIZE];
 
-  snprintf(filename,MAXSTRINGSIZE,"%s/%s",path,STANDARDELEMENTSDB);
+  strncpy(filename, path, MAXSTRINGSIZE);
+  strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));
+  strncat(filename, STANDARDELEMENTSDB, MAXSTRINGSIZE-strlen(filename));
   f.open(filename);
   if (f != NULL) {

src/tesselation.cpp

@@ -17 +17 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
-#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -470 +469 @@
 
   try {
-    Line centerLine = makeLineThrough(*MolCenter, *x);
-    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(centerLine);
-
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
-    DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
-
-    if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
-      return true;
-    }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
-      return true;
-    }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
-      DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
-      return true;
-    }
-    // Calculate cross point between one baseline and the line from the third endpoint to intersection
-    int i = 0;
-    do {
-      Line line1 = makeLineThrough(*(endpoints[i%3]->node->node),*(endpoints[(i+1)%3]->node->node));
-      Line line2 = makeLineThrough(*(endpoints[(i+2)%3]->node->node),*Intersection);
-      CrossPoint = line1.getIntersection(line2);
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
+  }
+  catch (LinearDependenceException &excp) {
+    Log() << Verbose(1) << excp;
+    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
+    return false;
+  }
+
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);
+
+  if (Intersection->DistanceSquared(*endpoints[0]->node->node) < MYEPSILON) {
+    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);
+    return true;
+  }   else if (Intersection->DistanceSquared(*endpoints[1]->node->node) < MYEPSILON) {
+    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);
+    return true;
+  }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
+    DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);
+    return true;
+  }
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  int i = 0;
+  do {
+    try {
+      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
+                                                    *(endpoints[(i+1)%3]->node->node),
+                                                    *(endpoints[(i+2)%3]->node->node),
+                                                    *Intersection);
       helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
       CrossPoint -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
@@ -499 +505 @@
       if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
         DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);
-        return false;
+        i=4;
+        break;
       }
       i++;
-    } while (i < 3);
+    } catch (LinearDependenceException &excp){
+      break;
+    }
+  } while (i < 3);
+  if (i == 3) {
     DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);
     return true;
-  }
-  catch (MathException &excp) {
-    Log() << Verbose(1) << excp;
-    DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);
+  } else {
+    DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);
     return false;
   }
@@ -535 +544 @@
   GetCenter(&Direction);
   try {
-    Line l = makeLineThrough(*x, Direction);
-    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(l);
-  }
-  catch (MathException &excp) {
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
+  }
+  catch (LinearDependenceException &excp) {
     (*ClosestPoint) = (*x);
   }
@@ -561 +569 @@
     Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    Line l = makeLineThrough(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
-    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(l);
+    CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(endpoints[i%3]->node->node), *(endpoints[(i+1)%3]->node->node));
     CrossDirection[i] = CrossPoint[i] - InPlane;
     CrossPoint[i] -= (*endpoints[i%3]->node->node);  // cross point was returned as absolute vector
@@ -733 +740 @@
 }
 
-Vector BoundaryTriangleSet::getEndpoint(int i) const{
-  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
-
-  return *endpoints[i]->node->node;
-}
-
-string BoundaryTriangleSet::getEndpointName(int i) const{
-  ASSERT(i>=0 && i<3,"Index of Endpoint out of Range");
-
-  return endpoints[i]->node->getName();
-}
-
 /** output operator for BoundaryTriangleSet.
  * \param &ost output stream
@@ -751 +746 @@
 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a)
 {
-  ost << "[" << a.Nr << "|" << a.getEndpointName(0) << "," << a.getEndpointName(1) << "," << a.getEndpointName(2) << "]";
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << "," << a.endpoints[1]->node->getName() << "," << a.endpoints[2]->node->getName() << "]";
   //  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","
   //      << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";
@@ -1519 +1514 @@
         CenterVector.Zero();
         for (int i = 0; i < 3; i++)
-          CenterVector += BTS->getEndpoint(i);
+          CenterVector += (*BTS->endpoints[i]->node->node);
         CenterVector.Scale(1. / 3.);
         DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);
@@ -4873 +4868 @@
   if (LastTriangle != NULL) {
     stringstream sstr;
-    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
+    sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->endpoints[0]->node->getName() << "_" << LastTriangle->endpoints[1]->node->getName() << "_" << LastTriangle->endpoints[2]->node->getName();
     NumberName = sstr.str();
     if (DoTecplotOutput) {

src/tesselation.hpp

@@ -171 +171 @@
 
     Plane getPlane() const;
-    Vector getEndpoint(int) const;
-    std::string getEndpointName(int) const;
 
     class BoundaryPointSet *endpoints[3];
@@ -179 +177 @@
     Vector SphereCenter;
     int Nr;
-
-  private:
-
 };
 

src/tesselationhelpers.cpp

@@ -15 +15 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
-#include "Line.hpp"
 #include "vector_ops.hpp"
 #include "verbose.hpp"
@@ -725 +724 @@
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  Line line1 = makeLineThrough(*(Base->endpoints[0]->node->node),*(Base->endpoints[1]->node->node));
-  Line line2 = makeLineThrough(NewOffset, NewDirection);
-  *Intersection = line1.getIntersection(line2);
+  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
+                                                   *(Base->endpoints[1]->node->node),
+                                                     NewOffset, NewDirection);
   Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
   DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);

src/unittests/Makefile.am

@@ -23 +23 @@
   InfoUnitTest \
   LinearSystemOfEquationsUnitTest \
-  LineUnittest \
   LinkedCellUnitTest \
   ListOfBondsUnitTest \
@@ -65 +64 @@
   infounittest.cpp \
   linearsystemofequationsunittest.cpp \
-  LineUnittest.cpp \
   LinkedCellUnitTest.cpp \
   listofbondsunittest.cpp \
@@ -99 +97 @@
   infounittest.hpp \
   linearsystemofequationsunittest.hpp \
-  LineUnittest.hpp \
   LinkedCellUnitTest.hpp \
   listofbondsunittest.hpp \
@@ -156 +153 @@
 LinearSystemOfEquationsUnitTest_LDADD = ${ALLLIBS}
 
-LineUnittest_SOURCES = UnitTestMain.cpp LineUnittest.cpp LineUnittest.hpp
-LineUnittest_LDADD = ${ALLLIBS}
-
 LinkedCellUnitTest_SOURCES = UnitTestMain.cpp LinkedCellUnitTest.cpp LinkedCellUnitTest.hpp 
 LinkedCellUnitTest_LDADD = ${ALLLIBS}

src/unittests/PlaneUnittest.cpp

@@ -17 +17 @@
 
 #include "vector.hpp"
-#include "Line.hpp"
 
 CPPUNIT_TEST_SUITE_REGISTRATION( PlaneUnittest );
@@ -154 +153 @@
   CPPUNIT_ASSERT(fabs(p4->distance(e1)-1) < MYEPSILON);
   CPPUNIT_ASSERT_EQUAL(zeroVec,p4->getClosestPoint(e1));
+
+
 }
-
-void PlaneUnittest::mirrorTest(){
-  Vector fixture;
-
-  // some Vectors that lie on the planes
-  fixture = p1->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p1->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = p1->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  fixture = p2->mirrorVector(zeroVec);
-  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p2->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p2->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-
-  fixture = p3->mirrorVector(zeroVec);
-  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p3->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,e1);
-  fixture = p3->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  fixture = p4->mirrorVector(zeroVec);
-  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
-  fixture = p4->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,e2);
-  fixture = p4->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,e3);
-
-  // some Vectors outside of the planes
-  {
-    Vector t = (2./3.)*(e1+e2+e3);
-    fixture = p1->mirrorVector(zeroVec);
-    CPPUNIT_ASSERT_EQUAL(fixture,t);
-  }
-
-  fixture = p2->mirrorVector(e3);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e3);
-  fixture = p3->mirrorVector(e2);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e2);
-  fixture = p4->mirrorVector(e1);
-  CPPUNIT_ASSERT_EQUAL(fixture,-1*e1);
-}
-
-void PlaneUnittest::LineIntersectionTest(){
-  Vector fixture;
-  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  Line l1 = makeLineThrough(zeroVec,Vector(2,1,0));
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e1, zeroVec).GetIntersection(l1) );
-  CPPUNIT_ASSERT_EQUAL( zeroVec, fixture );
-
-  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
-  Line l2 = makeLineThrough(e1,Vector(0,1,1));
-  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(e2, Vector(2,1,0)).GetIntersection(l2) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
-}

src/unittests/PlaneUnittest.hpp

@@ -20 +20 @@
   CPPUNIT_TEST ( pointsTest );
   CPPUNIT_TEST ( operationsTest );
-  CPPUNIT_TEST ( mirrorTest );
-  CPPUNIT_TEST ( LineIntersectionTest );
   CPPUNIT_TEST_SUITE_END();
 
@@ -32 +30 @@
   void pointsTest();
   void operationsTest();
-  void mirrorTest();
-  void LineIntersectionTest();
 
 private:

src/unittests/vectorunittest.cpp

@@ -216 +216 @@
 }
 
+/** UnitTest for line intersections.
+ */
+void VectorTest::LineIntersectionTest()
+{
+  // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unit, zero).GetIntersection(zero, two) );
+  CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(otherunit, two).GetIntersection( unit, notunit) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
+
+  // four vectors equal to zero
+  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(zero, zero, zero, zero), LinearDependenceException);
+  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // four vectors equal to unit
+  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, unit, unit, unit), LinearDependenceException);
+  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // two equal lines
+  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, two));
+  CPPUNIT_ASSERT_EQUAL( unit, fixture );
+
+  // line from (1,0,0) to (2,1,0) cuts line from (1,0,0) to (0,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, otherunit) );
+  CPPUNIT_ASSERT_EQUAL( unit, fixture );
+
+  // line from (1,0,0) to (0,0,0) cuts line from (0,0,0) to (2,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, zero, zero, two) );
+  CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // line from (1,0,0) to (2,1,0) cuts line from (0,0,0) to (0,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, zero, otherunit) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0., -1., 0.), fixture );
+};
+
+/** UnitTest for vector rotations.
+ */
+void VectorTest::VectorRotationTest()
+{
+  fixture = Vector(-1.,0.,0.);
+
+  // zero vector does not change
+  fixture = RotateVector(zero,unit, 1.);
+  CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  fixture = RotateVector(zero, two, 1.);
+  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
+
+  // vector on axis does not change
+  fixture = RotateVector(unit,unit, 1.);
+  CPPUNIT_ASSERT_EQUAL( unit, fixture );
+
+  // rotations
+  fixture = RotateVector(otherunit, unit, M_PI);
+  CPPUNIT_ASSERT_EQUAL( Vector(0.,-1.,0.), fixture );
+
+  fixture = RotateVector(otherunit, unit, 2. * M_PI);
+  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
+
+  fixture = RotateVector(otherunit,unit, 0);
+  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
+
+  fixture = RotateVector(Vector(0.,0.,1.), notunit, M_PI);
+  CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
+}
 
 /**

src/unittests/vectorunittest.hpp

@@ -27 +27 @@
     CPPUNIT_TEST ( ProjectionTest );
     CPPUNIT_TEST ( NormalsTest );
+    CPPUNIT_TEST ( LineIntersectionTest );
+    CPPUNIT_TEST ( VectorRotationTest );
     CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();

src/vector.cpp

@@ -213 +213 @@
 {
   Vector tmp;
-  tmp[0] = x[1]* y[2] - x[2]* y[1];
-  tmp[1] = x[2]* y[0] - x[0]* y[2];
-  tmp[2] = x[0]* y[1] - x[1]* y[0];
+  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
+  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
+  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
   (*this) = tmp;
 };
@@ -232 +232 @@
   *this -= tmp;
 };
+
+/** Calculates the minimum distance vector of this vector to the plane.
+ * \param *out output stream for debugging
+ * \param *PlaneNormal normal of plane
+ * \param *PlaneOffset offset of plane
+ * \return distance to plane
+ * \return distance vector onto to plane
+ */
+Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
+{
+  Vector temp = (*this) - PlaneOffset;
+  temp.MakeNormalTo(PlaneNormal);
+  temp.Scale(-1.);
+  // then add connecting vector from plane to point
+  temp += (*this)-PlaneOffset;
+  double sign = temp.ScalarProduct(PlaneNormal);
+  if (fabs(sign) > MYEPSILON)
+    sign /= fabs(sign);
+  else
+    sign = 0.;
+
+  temp.Normalize();
+  temp.Scale(sign);
+  return temp;
+};
+
 
 /** Calculates the minimum distance of this vector to the plane.
@@ -525 +551 @@
   MatrixMultiplication(M);
 };
-
-std::pair<Vector,Vector> Vector::partition(const Vector &rhs) const{
-  double factor = ScalarProduct(rhs)/rhs.NormSquared();
-  Vector res= factor * rhs;
-  return make_pair(res,(*this)-res);
-}
-
-std::pair<pointset,Vector> Vector::partition(const pointset &points) const{
-  Vector helper = *this;
-  pointset res;
-  for(pointset::const_iterator iter=points.begin();iter!=points.end();++iter){
-    pair<Vector,Vector> currPart = helper.partition(*iter);
-    res.push_back(currPart.first);
-    helper = currPart.second;
-  }
-  return make_pair(res,helper);
-}
 
 /** Do a matrix multiplication.
@@ -602 +611 @@
 };
 
+/** Mirrors atom against a given plane.
+ * \param n[] normal vector of mirror plane.
+ */
+void Vector::Mirror(const Vector &n)
+{
+  double projection;
+  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
+  // withdraw projected vector twice from original one
+  for (int i=NDIM;i--;)
+    at(i) -= 2.*projection*n[i];
+};
+
 /** Calculates orthonormal vector to one given vectors.
  * Just subtracts the projection onto the given vector from this vector.
@@ -612 +633 @@
   bool result = false;
   double factor = y1.ScalarProduct(*this)/y1.NormSquared();
-  Vector x1 = factor * y1;
+  Vector x1;
+  x1 = factor * y1;
   SubtractVector(x1);
   for (int i=NDIM;i--;)

src/vector.hpp

@@ -16 +16 @@
 
 #include <memory>
-#include <vector>
 
 #include "defs.hpp"
@@ -22 +21 @@
 
 /********************************************** declarations *******************************/
-
-class Vector;
-
-typedef std::vector<Vector> pointset;
 
 /** Single vector.
@@ -44 +39 @@
 
   double DistanceSquared(const Vector &y) const;
+  Vector GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
   double DistanceToSpace(const Space& space) const;
   double PeriodicDistance(const Vector &y, const double * const cell_size) const;
@@ -60 +56 @@
   void ProjectIt(const Vector &y);
   Vector Projection(const Vector &y) const;
+  void Mirror(const Vector &x);
   void ScaleAll(const double *factor);
   void Scale(const double factor);
@@ -69 +66 @@
   bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   void WrapPeriodically(const double * const M, const double * const Minv);
-  std::pair<Vector,Vector> partition(const Vector&) const;
-  std::pair<pointset,Vector> partition(const pointset&) const;
 
   // Accessors ussually come in pairs... and sometimes even more than that

src/vector_ops.cpp

@@ -15 +15 @@
 #include "Helpers/fast_functions.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
-#include "Exceptions/SkewException.hpp"
 
 #include <gsl/gsl_linalg.h>
@@ -111 +110 @@
   return true;
 };
+
+/** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
+ * \param *axis rotation axis
+ * \param alpha rotation angle in radian
+ */
+Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha)
+{
+  Vector a,y;
+  Vector res;
+  // normalise this vector with respect to axis
+  a = vec;
+  a.ProjectOntoPlane(axis);
+  // construct normal vector
+  try {
+    y = Plane(axis,a,0).getNormal();
+  }
+  catch (MathException &excp) {
+    // The normal vector cannot be created if there is linar dependency.
+    // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
+    return vec;
+  }
+  y.Scale(vec.Norm());
+  // scale normal vector by sine and this vector by cosine
+  y.Scale(sin(alpha));
+  a.Scale(cos(alpha));
+  res = vec.Projection(axis);
+  // add scaled normal vector onto this vector
+  res += y;
+  // add part in axis direction
+  res += a;
+  return res;
+};
+
+/** Calculates the intersection of the two lines that are both on the same plane.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
+ * \param *out output stream for debugging
+ * \param *Line1a first vector of first line
+ * \param *Line1b second vector of first line
+ * \param *Line2a first vector of second line
+ * \param *Line2b second vector of second line
+ * \return true - \a this will contain the intersection on return, false - lines are parallel
+ */
+Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b)
+{
+  Info FunctionInfo(__func__);
+
+  Vector res;
+
+  auto_ptr<GSLMatrix> M = auto_ptr<GSLMatrix>(new GSLMatrix(4,4));
+
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a[i]);
+    M->Set(1, i, Line1b[i]);
+    M->Set(2, i, Line2a[i]);
+    M->Set(3, i, Line2b[i]);
+  }
+
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  Log() << Verbose(1) << "INFO: Line1a = " << Line1a << ", Line1b = " << Line1b << ", Line2a = " << Line2a << ", Line2b = " << Line2b << "." << endl;
+
+
+  // constuct a,b,c
+  Vector a = Line1b - Line1a;
+  Vector b = Line2b - Line2a;
+  Vector c = Line2a - Line1a;
+  Vector d = Line2b - Line1b;
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   res.Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel = Line1a - Line2a;
+    factor = parallel.ScalarProduct(a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      res = Line2a;
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return res;
+    } else {
+      parallel = Line1a - Line2b;
+      factor = parallel.ScalarProduct(a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        res = Line2b;
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return res;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    res.Zero();
+    throw LinearDependenceException(__FILE__,__LINE__);
+  }
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1 = c;
+  temp1.VectorProduct(b);
+  temp2 = a;
+  temp2.VectorProduct(b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+
+  // construct intersection
+  res = a;
+  res.Scale(s);
+  res += Line1a;
+  Log() << Verbose(1) << "Intersection is at " << res << "." << endl;
+
+  return res;
+};

src/vector_ops.hpp

@@ -10 +10 @@
 
 bool LSQdistance(Vector &res,const Vector **vectors, int num);
+Vector RotateVector(const Vector &vec,const Vector &axis, const double alpha);
+Vector GetIntersectionOfTwoLinesOnPlane(const Vector &Line1a, const Vector &Line1b, const Vector &Line2a, const Vector &Line2b);
 
 #endif /* VECTOR_OPS_HPP_ */