Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/tesselation.cpp	(revision 97498a1c712a846028cfe0e2f484ccda52680f60)
@@ -388,5 +388,5 @@
 };
 
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
  * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
@@ -411,19 +411,23 @@
   }
 
+  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
+  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
+
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
   do {
     if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
+      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
+      break;
     } else 
       i++;
-    if (i>2)
-      break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
+  } while (i<3);
   if (i==3) {
     eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
   }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+  Log() << Verbose(1) << "INFO: Crosspoint is " << CrossPoint << "." << endl;
 
   // check whether intersection is inside or not by comparing length of intersection and length of cross point
@@ -2987,5 +2991,5 @@
  * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
  */
-DistanceMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
 {
   Info FunctionInfo(__func__);
@@ -3004,5 +3008,5 @@
   Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
 
-  DistanceMap * points = new DistanceMap;
+  DistanceToPointMap * points = new DistanceToPointMap;
   LC->GetNeighbourBounds(Nlower, Nupper);
   //Log() << Verbose(1) << endl;
@@ -3015,6 +3019,8 @@
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
-            if (FindPoint != PointsOnBoundary.end())
-              points->insert(DistancePair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+            if (FindPoint != PointsOnBoundary.end()) {
+              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
+            }
           }
         } else {
@@ -3042,5 +3048,5 @@
 
   // get closest points
-  DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
   if (points == NULL) {
     eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
@@ -3055,5 +3061,5 @@
   Vector Center;
   Vector BaseLine;
-  for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
       // calculate closest point on line to desired point
@@ -3109,5 +3115,5 @@
 
 	// get closest points
-	DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+	DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
   if (points == NULL) {
     eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
@@ -3118,40 +3124,53 @@
   Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
   BoundaryLineSet *ClosestLine = NULL;
-  double MinDistance = -1.;
-  Vector helper;
+  double MinDistance = 1e+16;
+  Vector BaseLineIntersection;
   Vector Center;
   Vector BaseLine;
-  for (DistanceMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+  Vector BaseLineCenter;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
     for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
-      for (TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++) {
+
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthBase = BaseLine.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      const double lengthEndA = BaseLineIntersection.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthEndB = BaseLineIntersection.NormSquared();
+
+      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+        if ((lengthEndA < MinDistance) && (lengthEndA < lengthEndB)) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndA;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEndA << "." << endl;
+        } else if (lengthEndB < MinDistance) {
+          ClosestLine = LineRunner->second;
+          MinDistance = lengthEndB;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is closer with " << lengthEndB << "." << endl;
+        } else {
+          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
+        }
+      } else { // intersection is closer, calculate
         // calculate closest point on line to desired point
-        helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-        helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
-        helper.Scale(0.5);
-        Center.CopyVector(x);
-        Center.SubtractVector(&helper);
-        BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
-        BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        BaseLineIntersection.CopyVector(x);
+        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.CopyVector(&BaseLineIntersection);
         Center.ProjectOntoPlane(&BaseLine);
-        const double distance = Center.NormSquared();
+        BaseLineIntersection.SubtractVector(&Center);
+        const double distance = BaseLineIntersection.NormSquared();
+        if (Center.NormSquared() > BaseLine.NormSquared()) {
+          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+        }
         if ((ClosestLine == NULL) || (distance < MinDistance)) {
-          // additionally calculate intersection on line (whether it's on the line section or not)
-          helper.CopyVector(x);
-          helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
-          helper.SubtractVector(&Center);
-          const double lengthA = helper.ScalarProduct(&BaseLine);
-          helper.CopyVector(x);
-          helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
-          helper.SubtractVector(&Center);
-          const double lengthB = helper.ScalarProduct(&BaseLine);
-          if (lengthB*lengthA < 0) {  // if have different sign
-            ClosestLine = LineRunner->second;
-            MinDistance = distance;
-            Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
-          } else {
-            Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
-          }
+          ClosestLine = LineRunner->second;
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
         } else {
-          Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
         }
       }
@@ -3190,46 +3209,27 @@
     return NULL;
 
-  // reject all which are not closest
-  double MinDistance = -1.;
+  // go through all and pick the one with the best alignment to x
+  double MinAlignment = 2.*M_PI;
   for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
     helper.CopyVector(x);
     helper.SubtractVector(&Center);
-    helper.ProjectOntoPlane(&(*Runner)->NormalVector);
-    const double distance = helper.NormSquared();
-    if (fabs(distance - MinDistance) < MYEPSILON) {  // degenerate case, keep both
-      candidates.push_back(*Runner);
-    } else if (distance < MinDistance) {
-      candidates.clear();
-      candidates.push_back(*Runner);
-      MinDistance = distance;
+    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    if (Alignment < MinAlignment) {
+      result = *Runner;
+      MinAlignment = Alignment;
+      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
+    } else {
+      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
     }
   }
   delete(triangles);
-  if (!candidates.empty()) {
-    if (candidates.size() == 1) { // there is no degenerate case
-      result = candidates.front();
-      Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
-    } else {
-      result = candidates.front();
-      result->GetCenter(&Center);
-      Center.SubtractVector(x);
-      Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
-      if (Center.ScalarProduct(&result->NormalVector) < 0) {
-        result = candidates.back();
-        Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
-        if (Center.ScalarProduct(&result->NormalVector) < 0) {
-          eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
-        }
-      }
-    }
-  } else {
-    result = NULL;
-    Log() << Verbose(0) << "No closest triangle found." << endl;
-  }
+
   return result;
 };
 
 /** Checks whether the provided Vector is within the tesselation structure.
+ * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
+ * towards or away from the given \a &Point.
  *
  * @param point of which to check the position
@@ -3245,4 +3245,6 @@
   Vector Center;
   Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
 
   if (result == NULL) {// is boundary point or only point in point cloud?
@@ -3255,37 +3257,34 @@
   result->GetCenter(&Center);
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  Center.SubtractVector(&Point);
+  DistanceToCenter.CopyVector(&Center);
+  DistanceToCenter.SubtractVector(&Point);
   Log() << Verbose(2) << "INFO: Vector from point to test to center is " << Center << "." << endl;
-  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
+
+  // check whether we are on boundary
+  if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
+    // calculate whether inside of triangle
+    DistanceToCenter.ProjectOntoPlane(&result->NormalVector);
+    DistanceToCenter.AddVector(&Center);
+    Center.CopyVector(&DistanceToCenter);
+    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
+    DistanceToCenter.AddVector(&result->NormalVector); // points outside
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
+      return true;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
+      return false;
+    }
+  }
+
+  // then check direction to boundary
+  if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
     Log() << Verbose(1) << Point << " is an inner point." << endl;
     return true;
   } else {
-    for (int i=0;i<3;i++) {
-      helper.CopyVector(result->endpoints[i]->node->node);
-      helper.SubtractVector(&Point);
-      if (Center.NormSquared() > helper.NormSquared())
-        Center.CopyVector(&helper);
-    }
-    if (Center.NormSquared() < epsilon*epsilon) {
-      Log() << Verbose(1) << Point << " is inner point, being sufficiently close (wrt " << epsilon << ") to boundary." << endl;
-      return true;
-    }
     Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
     return false;
   }
-}
-
-/** Checks whether the provided TesselPoint is within the tesselation structure.
- *
- * @param *Point of which to check the position
- * @param *LC Linked Cell structure
- * @param epsilon Distance of \a &Point to Center of triangle is tested greater against this value (standard: -MYEPSILON)
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon) const
-{
-	Info FunctionInfo(__func__);
-  return IsInnerPoint(*(Point->node), LC, epsilon);
 }
 
@@ -3460,16 +3459,39 @@
   }
 
+  // check whether there's something to do
+  if (SetOfNeighbours->size() < 3) {
+    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+      connectedCircle->push_back(*TesselRunner);
+    return connectedCircle;
+  }
+
+  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
   // calculate central point
-  for (TesselPointSet::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
-    center.AddVector((*TesselRunner)->node);
+
+  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
+  TesselB++;
+  TesselC++;
+  TesselC++;
+  int counter = 0;
+  while (TesselC != SetOfNeighbours->end()) {
+    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
+    counter++;
+    TesselA++;
+    TesselB++;
+    TesselC++;
+    PlaneNormal.AddVector(&helper);
+  }
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
-  //  << "; scale factor " << 1.0/connectedPoints.size();
-  center.Scale(1.0/SetOfNeighbours->size());
-  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
-
-  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
-  PlaneNormal.CopyVector(Point->node);
-  PlaneNormal.SubtractVector(&center);
-  PlaneNormal.Normalize();
+  //  << "; scale factor " << counter;
+  PlaneNormal.Scale(1.0/(double)counter);
+//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+//
+//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+//  PlaneNormal.CopyVector(Point->node);
+//  PlaneNormal.SubtractVector(&center);
+//  PlaneNormal.Normalize();
   Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
@@ -3504,4 +3526,6 @@
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
+    if (angle > M_PI) // the correction is of no use here (and not desired)
+      angle = 2.*M_PI - angle;
     Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
     InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/tesselation.hpp	(revision 97498a1c712a846028cfe0e2f484ccda52680f60)
@@ -79,6 +79,6 @@
 #define PolygonList list < class BoundaryPolygonSet * >
 
-#define DistanceMap multimap <double, class BoundaryPointSet * >
-#define DistancePair pair <double, class BoundaryPointSet * >
+#define DistanceToPointMap multimap <double, class BoundaryPointSet * >
+#define DistanceToPointPair pair <double, class BoundaryPointSet * >
 
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
@@ -313,7 +313,6 @@
     class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
-    bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
-    DistanceMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
+    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
     BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;
     TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
Index: src/unittests/tesselation_insideoutsideunittest.cpp
===================================================================
--- src/unittests/tesselation_insideoutsideunittest.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/unittests/tesselation_insideoutsideunittest.cpp	(revision 97498a1c712a846028cfe0e2f484ccda52680f60)
@@ -26,5 +26,5 @@
 void TesselationInOutsideTest::setUp()
 {
-  setVerbosity(3);
+  setVerbosity(5);
 
   // create corners
@@ -32,5 +32,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 0., 0.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "1");
   Walker->nr = 1;
@@ -38,5 +38,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 1., 0.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "2");
   Walker->nr = 2;
@@ -44,5 +44,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., 0.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "3");
   Walker->nr = 3;
@@ -50,5 +50,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 1., 0.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "4");
   Walker->nr = 4;
@@ -56,5 +56,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 0., 1.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "5");
   Walker->nr = 5;
@@ -62,5 +62,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(0., 1., 1.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "6");
   Walker->nr = 6;
@@ -68,5 +68,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 0., 1.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "7");
   Walker->nr = 7;
@@ -74,5 +74,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(1., 1., 1.);
-  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
   strcpy(Walker->Name, "8");
   Walker->nr = 8;
@@ -146,5 +146,5 @@
 {
   double n[3];
-  const double boundary = 5.;
+  const double boundary = 2.;
   const double step = 1.;
 
@@ -153,5 +153,5 @@
     for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
       for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
-        if ( ((n[0] > 0) && (n[0] > 0) && (n[0] > 0)) && (n[0]+n[1]+n[2] <= 1.))
+        if ( ((n[0] >= 0.) && (n[1] >= 0.) && (n[2] >= 0.)) && (fabs(n[0])+fabs(n[1])+fabs(n[2]) <= 1.))
           CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
         else
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/vector.cpp	(revision 97498a1c712a846028cfe0e2f484ccda52680f60)
@@ -8,9 +8,15 @@
 #include "defs.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
+#include "info.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
+#include "memoryallocator.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
 
 /************************************ Functions for class vector ************************************/
@@ -219,4 +225,5 @@
 bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
 {
+  Info FunctionInfo(__func__);
   double factor;
   Vector Direction, helper;
@@ -226,5 +233,5 @@
   Direction.SubtractVector(Origin);
   Direction.Normalize();
-  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
   if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
@@ -236,5 +243,5 @@
   factor = helper.ScalarProduct(PlaneNormal)/factor;
   if (factor < MYEPSILON) { // Origin is in-plane
-    //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
+    Log() << Verbose(1) << "Origin of line is in-plane, simple." << endl;
     CopyVector(Origin);
     return true;
@@ -243,5 +250,5 @@
   Direction.Scale(factor);
   CopyVector(Origin);
-  //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 
@@ -250,5 +257,5 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
   } else {
@@ -299,14 +306,13 @@
 bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
 {
-  bool result = true;
+  Info FunctionInfo(__func__);
   Vector Direction, OtherDirection;
-  Vector AuxiliaryNormal;
-  Vector Distance;
-  const Vector *Normal = NULL;
-  Vector *ConstructedNormal = NULL;
-  bool FreeNormal = false;
+  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
+  gsl_vector *b = gsl_vector_alloc(NDIM);
+  gsl_vector *x = gsl_vector_alloc(NDIM);
+  gsl_permutation *perm = NULL;
+  int signum;
 
   // construct both direction vectors
-  Zero();
   Direction.CopyVector(Line1b);
   Direction.SubtractVector(Line1a);
@@ -318,61 +324,48 @@
     return false;
 
-  Direction.Normalize();
-  OtherDirection.Normalize();
-
-  //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
-
-  if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
-    if ((Line1a == Line2a) || (Line1a == Line2b))
-      CopyVector(Line1a);
-    else if ((Line1b == Line2b) || (Line1b == Line2b))
-        CopyVector(Line1b);
-    else
-      return false;
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
+  // set vector
+  for (int i=0;i<NDIM;i++)
+    gsl_vector_set(b, i, Line1a->x[i]-Line2a->x[i]);
+  Log() << Verbose(1) << "b = " << endl;
+  gsl_vector_fprintf(stdout, b, "%g");
+
+  // set matrix
+  for (int i=0;i<NDIM;i++)
+    gsl_matrix_set(A, 0, i, -Direction.x[i]);
+  for (int i=0;i<NDIM;i++)
+    gsl_matrix_set(A, 1, i, OtherDirection.x[i]);
+  for (int i=0;i<NDIM;i++)
+    gsl_matrix_set(A, 2, i, 1.);
+  Log() << Verbose(1) << "A = " << endl;
+  gsl_matrix_fprintf(stdout, A, "%g");
+
+  // Solve A x=b for x
+  perm = gsl_permutation_alloc(NDIM);
+  gsl_linalg_LU_decomp(A, perm, &signum);
+  gsl_linalg_LU_solve(A, perm, b, x);
+  gsl_permutation_free(perm);
+  gsl_vector_free(b);
+  gsl_matrix_free(A);
+
+  Log() << Verbose(1) << "Solution is " << gsl_vector_get(x,0) << ", " << gsl_vector_get(x,1) << "." << endl;
+
+  CopyVector(&Direction);
+  Scale(gsl_vector_get(x,0));
+  AddVector(Line1a);
+  Log() << Verbose(1) << "INFO: First intersection is " << *this << "." << endl;
+
+  Vector test;
+  test.CopyVector(&OtherDirection);
+  test.Scale(gsl_vector_get(x,1));
+  test.AddVector(Line2a);
+  Log() << Verbose(1) << "INFO: Second intersection is " << test << "." << endl;
+  test.SubtractVector(this);
+
+  gsl_vector_free(x);
+
+  if (test.IsZero())
     return true;
-  } else {
-    // check whether we have a plane normal vector
-    if (PlaneNormal == NULL) {
-      ConstructedNormal = new Vector;
-      ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
-      Normal = ConstructedNormal;
-      FreeNormal = true;
-    } else
-      Normal = PlaneNormal;
-
-    AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
-    //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
-
-    Distance.CopyVector(Line2a);
-    Distance.SubtractVector(Line1a);
-    //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
-    if (Distance.IsZero()) {
-      // offsets are equal, match found
-      CopyVector(Line1a);
-      result = true;
-    } else {
-      CopyVector(Distance.Projection(&AuxiliaryNormal));
-      //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
-      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
-      //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
-      Scale(1./(factor*factor));
-      //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
-      CopyVector(Projection(&Direction));
-      //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
-      if (this->IsZero())
-        result = false;
-      else
-        result = true;
-      AddVector(Line1a);
-    }
-
-    if (FreeNormal)
-      delete(ConstructedNormal);
-  }
-  if (result)
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
-
-  return result;
+  else
+    return false;
 };
 
