Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/analysis_correlation.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -267,5 +267,5 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
+          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
           if (triangle != NULL) {
             distance = DistanceToTrianglePlane(Walker->node, triangle);
@@ -327,5 +327,5 @@
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
-                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                triangle = Surface->FindClosestTriangleToVector(&checkX, LC );
                 if (triangle != NULL) {
                   distance = DistanceToTrianglePlane(&checkX, triangle);
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/tesselation.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -1423,5 +1423,5 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1438,5 +1438,5 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -2339,13 +2339,14 @@
 
   // fill the set of neighbours
-  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  Center.SubtractVector(TurningPoint->node);
-  set<TesselPoint*> SetOfNeighbours;
+  TesselPointSet SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
@@ -2357,5 +2358,4 @@
     AddTesselationPoint((*Sprinter), 2);
 
-
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2374,4 +2374,7 @@
     Runner = Sprinter;
     Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
   }
   delete(connectedClosestPoints);
@@ -2958,6 +2961,6 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -2979,12 +2982,12 @@
 };
 
-/** Finds the triangle that is closest to a given Vector \a *x.
+/** Finds the boundary points that are closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
- */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
-{
-	Info FunctionInfo(__func__);
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
   PointMap::const_iterator FindPoint;
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
@@ -3001,5 +3004,5 @@
   Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
 
-  set<BoundaryPointSet *> points;
+  DistanceMap * points = new DistanceMap;
   LC->GetNeighbourBounds(Nlower, Nupper);
   //Log() << Verbose(1) << endl;
@@ -3013,5 +3016,5 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end())
-              points.insert(FindPoint->second);
+              points->insert(DistancePair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
           }
         } else {
@@ -3021,6 +3024,25 @@
 
   // check whether we found some points
-  if (points.empty()) {
-    Log() << Verbose(1) << "There is no nearest triangle, too far away from the surface." << endl;
+  if (points->empty()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+
+  // get closest points
+  DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
     return NULL;
   }
@@ -3033,6 +3055,6 @@
   Vector Center;
   Vector BaseLine;
-  for (set<BoundaryPointSet *>::iterator Runner = points.begin(); Runner != points.end(); Runner++) {
-    for (LineMap::iterator LineRunner = (*Runner)->lines.begin(); LineRunner != (*Runner)->lines.end(); LineRunner++) {
+  for (DistanceMap::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
       helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
@@ -3067,4 +3089,5 @@
     }
   }
+  delete(points);
   // check whether closest line is "too close" :), then it's inside
   if (ClosestLine == NULL) {
@@ -3072,9 +3095,79 @@
     return NULL;
   }
-  list<BoundaryTriangleSet*> *triangles = new list<BoundaryTriangleSet*>;
+  return ClosestLine;
+};
+
+
+/** Finds the triangle that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+	Info FunctionInfo(__func__);
+
+	// get closest points
+	DistanceMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (DistanceMap::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++) {
+        // 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);
+        Center.ProjectOntoPlane(&BaseLine);
+        const double distance = Center.NormSquared();
+        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;
+          }
+        } else {
+          Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+        }
+      }
+    }
+  }
+  delete(points);
+
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
+    return NULL;
+  }
+  TriangleList * candidates = new TriangleList;
   for (TriangleMap::iterator Runner = ClosestLine->triangles.begin(); Runner != ClosestLine->triangles.end(); Runner++) {
-    triangles->push_back(Runner->second);
-  }
-  return triangles;
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 
@@ -3085,10 +3178,10 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
 {
 	Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
-  list<BoundaryTriangleSet*> candidates;
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList candidates;
   Vector Center;
   Vector helper;
@@ -3099,5 +3192,5 @@
   // reject all which are not closest
   double MinDistance = -1.;
-  for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
     (*Runner)->GetCenter(&Center);
     helper.CopyVector(x);
@@ -3149,5 +3242,5 @@
 {
 	Info FunctionInfo(__func__);
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+  class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
   Vector helper;
@@ -3197,128 +3290,4 @@
 }
 
-/**
- * Finds the point which is second closest to the provided one.
- *
- * @param Point to which to find the second closest other point
- * @param linked cell structure
- *
- * @return point which is second closest to the provided one
- * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
- */
-TesselPoint* Tesselation::FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TesselPoint* closestPoint = NULL;
-  TesselPoint* secondClosestPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //Log() << Verbose(1) << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
-            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
-              if (currentNorm < distance) {
-                // remember second point
-                secondDistance = distance;
-                secondClosestPoint = closestPoint;
-                // mark down new closest point
-                distance = currentNorm;
-                closestPoint = (*Runner);
-                //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
-              }
-            }
-          }
-        } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return secondClosestPoint;
-};
-
-/**
- * Finds the point which is closest to the provided one.
- *
- * @param Point to which to find the closest other point
- * @param SecondPoint the second closest other point on return, NULL if none found
- * @param linked cell structure
- *
- * @return point which is closest to the provided one, NULL if none found
- * @TODO Optimize me! These finds are unnecessary if we have a LC list of BoundaryPoints only!
- */
-TesselPoint* Tesselation::FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const
-{
-  Info FunctionInfo(__func__);
-  TesselPoint* closestPoint = NULL;
-  SecondPoint = NULL;
-  double distance = 1e16;
-  double secondDistance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //Log() << Verbose(1) << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        const LinkedNodes *List = LC->GetCurrentCell();
-        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper.NormSquared();
-            if (PointsOnBoundary.find((*Runner)->nr) != PointsOnBoundary.end()) {
-              if (currentNorm < distance) {
-                secondDistance = distance;
-                SecondPoint = closestPoint;
-                distance = currentNorm;
-                closestPoint = (*Runner);
-                //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
-              } else if (currentNorm < secondDistance) {
-                secondDistance = currentNorm;
-                SecondPoint = (*Runner);
-                //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
-              }
-            }
-          }
-        } else {
-          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-  // output
-  if ((closestPoint != NULL) && (SecondPoint != NULL)) {
-    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << *SecondPoint << "." << endl;
-  } else if (closestPoint != NULL) {
-    Log() << Verbose(1) << "Closest point is " << *closestPoint << " and second closest is " << SecondPoint << "." << endl;
-  } else {
-    Log() << Verbose(1) << "No closest point has been found, only " << closestPoint << " and " << SecondPoint << "." << endl;
-  }
-  return closestPoint;
-};
-
-
 /** Gets all points connected to the provided point by triangulation lines.
  *
@@ -3327,8 +3296,8 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-	Info FunctionInfo(__func__);
-  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+	Info FunctionInfo(__func__);
+	TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3392,9 +3361,9 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
 	Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector PlaneNormal;
   Vector AngleZero;
@@ -3402,5 +3371,5 @@
   Vector helper;
   const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3413,5 +3382,5 @@
   triangles = FindTriangles(TrianglePoints);
   if ((triangles != NULL) && (!triangles->empty())) {
-    for (list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
       PlaneNormal.AddVector(&(*Runner)->NormalVector);
   } else {
@@ -3447,5 +3416,5 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3469,5 +3438,4 @@
  * triangle we are looking for.
  *
- * @param *out output stream for debugging
  * @param *SetOfNeighbours all points for which the angle should be calculated
  * @param *Point of which get all connected points
@@ -3475,9 +3443,9 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
   Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector center;
   Vector PlaneNormal;
@@ -3493,5 +3461,5 @@
 
   // calculate central point
-  for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+  for (TesselPointSet::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
     center.AddVector((*TesselRunner)->node);
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -3530,11 +3498,16 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
     helper.ProjectOntoPlane(&PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
-    Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
-    anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
   }
 
@@ -3552,10 +3525,10 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list< TesselPointList *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
 	Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3594,5 +3567,5 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3668,18 +3641,18 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-	Info FunctionInfo(__func__);
-  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
-  list<TesselPoint*> *newPath = NULL;
+list<TesselPointList *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+	Info FunctionInfo(__func__);
+  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
+  TesselPointList *newPath = NULL;
   int count = 0;
 
 
-  list<TesselPoint*>::iterator CircleRunner;
-  list<TesselPoint*>::iterator CircleStart;
-
-  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3690,11 +3663,11 @@
 
     // go through list, look for reappearance of starting Point and create list
-    list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::iterator Marker = CircleStart;
     for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
-        newPath = new list<TesselPoint*>;
-        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3730,8 +3703,8 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-	Info FunctionInfo(__func__);
-  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+	Info FunctionInfo(__func__);
+	TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
@@ -3782,16 +3755,16 @@
   }
 
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
     }
   }
@@ -3800,16 +3773,16 @@
   count=0;
   NormalVector.Zero();
-  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
-    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
-  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
-  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
@@ -3825,5 +3798,5 @@
 
       // re-create all triangles by going through connected points list
-      list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -3931,10 +3904,10 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -3983,8 +3956,8 @@
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-	Info FunctionInfo(__func__);
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+	Info FunctionInfo(__func__);
+	TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
@@ -4114,9 +4087,9 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedLines()
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
 {
 	Info FunctionInfo(__func__);
 	UniqueLines AllLines;
-  map<int, int> * DegeneratedLines = new map<int, int>;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
@@ -4139,5 +4112,5 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -4158,9 +4131,9 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-	Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+	Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -4168,5 +4141,5 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -4189,5 +4162,5 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -4203,10 +4176,10 @@
 {
 	Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -4463,5 +4436,5 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
-  map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4615,7 +4588,7 @@
   }
 
-  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/tesselation.hpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -52,4 +52,6 @@
 // ======================================================= some template functions =========================================
 
+#define IndexToIndex map <int, int>
+
 #define PointMap map < int, class BoundaryPointSet * >
 #define PointSet set < class BoundaryPointSet * >
@@ -76,4 +78,7 @@
 #define PolygonSet set < class BoundaryPolygonSet * >
 #define PolygonList list < class BoundaryPolygonSet * >
+
+#define DistanceMap multimap <double, class BoundaryPointSet * >
+#define DistancePair pair <double, class BoundaryPointSet * >
 
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
@@ -310,6 +315,8 @@
     bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
-    TesselPoint* FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const;
-    TesselPoint* FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) const;
+    DistanceMap * 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;
+    BoundaryTriangleSet* FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
 
     // print for debugging
Index: src/unittests/Makefile.am
===================================================================
--- src/unittests/Makefile.am	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/unittests/Makefile.am	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -4,5 +4,5 @@
 AM_CXXFLAGS = $(CPPUNIT_CFLAGS)
 
-TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest VectorUnitTest
+TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest TesselationUnitTest TesselationInOutsideUnitTest VectorUnitTest
 check_PROGRAMS = $(TESTS)
 noinst_PROGRAMS = $(TESTS)
@@ -44,4 +44,10 @@
 StackClassUnitTest_LDADD = ../libmolecuilder.a
 
+TesselationUnitTest_SOURCES = tesselationunittest.cpp tesselationunittest.hpp
+TesselationUnitTest_LDADD = ../libmolecuilder.a
+
+TesselationInOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
+TesselationInOutsideUnitTest_LDADD = ../libmolecuilder.a
+
 VectorUnitTest_SOURCES = vectorunittest.cpp vectorunittest.hpp 
 VectorUnitTest_LDADD = ../libmolecuilder.a
Index: src/unittests/tesselation_insideoutsideunittest.cpp
===================================================================
--- src/unittests/tesselation_insideoutsideunittest.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/unittests/tesselation_insideoutsideunittest.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -0,0 +1,181 @@
+/*
+ * tesselation_insideoutsideunittest.cpp
+ *
+ *  Created on: Dec 28, 2009
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "defs.hpp"
+#include "tesselation.hpp"
+#include "tesselation_insideoutsideunittest.hpp"
+
+#define SPHERERADIUS 2.
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( TesselationInOutsideTest );
+
+
+void TesselationInOutsideTest::setUp()
+{
+  setVerbosity(3);
+
+  // create corners
+  class TesselPoint *Walker;
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "1");
+  Walker->nr = 1;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 1., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "2");
+  Walker->nr = 2;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "3");
+  Walker->nr = 3;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 1., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "4");
+  Walker->nr = 4;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 0., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "5");
+  Walker->nr = 5;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 1., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "6");
+  Walker->nr = 6;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 0., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "7");
+  Walker->nr = 7;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 1., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp");
+  strcpy(Walker->Name, "8");
+  Walker->nr = 8;
+  Corners.push_back(Walker);
+
+  // create linkedcell
+  LinkedList = new LinkedCell(&Corners, 2.*SPHERERADIUS);
+
+  // create tesselation
+  TesselStruct = new Tesselation;
+  TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
+
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
+      Log() << Verbose(2) << *(Runner->second) << endl;
+
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
+    }
+
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
+      Log() << Verbose(2) << *(Runner->second) << endl;
+
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+        //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
+    }
+  }
+};
+
+
+void TesselationInOutsideTest::tearDown()
+{
+  delete(LinkedList);
+  delete(TesselStruct);
+  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
+    delete((*Runner)->node);
+    delete(*Runner);
+  }
+  Corners.clear();
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
+};
+
+/** UnitTest for Tesselation::IsInnerPoint()
+ */
+void TesselationInOutsideTest::IsInnerPointTest()
+{
+  double n[3];
+  const double boundary = 5.;
+  const double step = 1.;
+
+  // go through the mesh and check each point
+  for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
+    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.))
+          CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
+        else
+          CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
+      }
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/tesselation_insideoutsideunittest.hpp
===================================================================
--- src/unittests/tesselation_insideoutsideunittest.hpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
+++ src/unittests/tesselation_insideoutsideunittest.hpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -0,0 +1,38 @@
+/*
+ * tesselation_insideoutsideunittest.hpp
+ *
+ *  Created on: Dec 28, 2009
+ *      Author: heber
+ */
+
+#ifndef TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_
+#define TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_
+
+/*********************************************** includes ***********************************/
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "linkedcell.hpp"
+#include "tesselation.hpp"
+
+/********************************************** Test classes **************************************/
+
+class TesselationInOutsideTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( TesselationInOutsideTest) ;
+    CPPUNIT_TEST ( IsInnerPointTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void IsInnerPointTest();
+
+private:
+      class Tesselation *TesselStruct;
+      LinkedNodes Corners;
+      class LinkedCell *LinkedList;
+};
+
+
+#endif /* TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_ */
Index: src/unittests/tesselationunittest.cpp
===================================================================
--- src/unittests/tesselationunittest.cpp	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/unittests/tesselationunittest.cpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -30,18 +30,18 @@
   class TesselPoint *Walker;
   Walker = new TesselPoint;
-  Walker->node = new Vector(1., 0., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(1., 0., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., 1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., 1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., -1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., -1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "3");
   Walker->nr = 3;
@@ -49,5 +49,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 0., 1.);
-  Walker->Name = new char[3];
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "4");
   Walker->nr = 4;
@@ -59,19 +59,35 @@
   // create tesselation
   TesselStruct = new Tesselation;
-  TesselStruct->PointsOnBoundary.clear();
-  TesselStruct->LinesOnBoundary.clear();
-  TesselStruct->TrianglesOnBoundary.clear();
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->PointsOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->LinesOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->TrianglesOnBoundary.empty() );
   TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
-  bool flag = false;
 
-  LineMap::iterator baseline = TesselStruct->LinesOnBoundary.begin();
-  while (baseline != TesselStruct->LinesOnBoundary.end()) {
-    if (baseline->second->triangles.size() == 1) {
-      flag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
     }
-    baseline++;
-    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (flag)) {
-      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-      flag = false;
+
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
     }
   }
@@ -84,33 +100,11 @@
   delete(TesselStruct);
   for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
-    delete[]((*Runner)->Name);
     delete((*Runner)->node);
     delete(*Runner);
   }
   Corners.clear();
-};
-
-/** UnitTest for Tesselation::IsInnerPoint()
- */
-void TesselationTest::IsInnerPointTest()
-{
-  // true inside points
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.5,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.5,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.5), LinkedList) );
-
-  // corners
-  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
-    CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint((*Runner), LinkedList) );
-
-  // true outside points
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,5.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,0.,5.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(1.,1.,1.), LinkedList) );
-
-  // tricky point, there are three equally close triangles
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(5.,0.,0.), LinkedList) );
-
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
 };
 
Index: src/unittests/tesselationunittest.hpp
===================================================================
--- src/unittests/tesselationunittest.hpp	(revision 71b20e8e52bf69205f641774cf6e270e8ac8e84c)
+++ src/unittests/tesselationunittest.hpp	(revision c15ca26c8aaeb433b6c80f98d5a7967ed4398a4e)
@@ -21,5 +21,4 @@
 {
     CPPUNIT_TEST_SUITE( TesselationTest) ;
-    CPPUNIT_TEST ( IsInnerPointTest );
     CPPUNIT_TEST ( GetAllTrianglesTest );
     CPPUNIT_TEST ( ContainmentTest );
@@ -29,5 +28,4 @@
       void setUp();
       void tearDown();
-      void IsInnerPointTest();
       void GetAllTrianglesTest();
       void ContainmentTest();
