Changes in src/boundary.cpp [edb93c:8c54a3]

File:

• src/boundary.cpp (modified) (20 diffs)

src/boundary.cpp

@@ -1 @@
-/** \file boundary.hpp
- *
- * Implementations and super-function for envelopes
- */
-
-
 #include "boundary.hpp"
-
-#include<gsl/gsl_poly.h>
+#include "linkedcell.hpp"
+#include "molecules.hpp"
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_permutation.h>
+
+#define DEBUG 1
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 0
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
+#define HULLEPSILON 1e-7
+
+// ======================================== Points on Boundary =================================
+
+BoundaryPointSet::BoundaryPointSet()
+{
+  LinesCount = 0;
+  Nr = -1;
+}
+;
+
+BoundaryPointSet::BoundaryPointSet(atom *Walker)
+{
+  node = Walker;
+  LinesCount = 0;
+  Nr = Walker->nr;
+}
+;
+
+BoundaryPointSet::~BoundaryPointSet()
+{
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  if (!lines.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
+  node = NULL;
+}
+;
+
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
+  if (line->endpoints[0] == this)
+    {
+      lines.insert(LinePair(line->endpoints[1]->Nr, line));
+    }
+  else
+    {
+      lines.insert(LinePair(line->endpoints[0]->Nr, line));
+    }
+  LinesCount++;
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryPointSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  return ost;
+}
+;
+
+// ======================================== Lines on Boundary =================================
+
+BoundaryLineSet::BoundaryLineSet()
+{
+  for (int i = 0; i < 2; i++)
+    endpoints[i] = NULL;
+  TrianglesCount = 0;
+  Nr = -1;
+}
+;
+
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
+  // add this line to the hash maps of both endpoints
+  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+  Point[1]->AddLine(this); //
+  // clear triangles list
+  TrianglesCount = 0;
+  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+}
+;
+
+BoundaryLineSet::~BoundaryLineSet()
+{
+  int Numbers[2];
+  Numbers[0] = endpoints[1]->Nr;
+  Numbers[1] = endpoints[0]->Nr;
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
+    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
+    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
+      if ((*Runner).second == this) {
+        endpoints[i]->lines.erase(Runner);
+        break;
+      }
+    if (endpoints[i]->lines.empty()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+  if (!triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+}
+;
+
+void
+BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+      << endl;
+  triangles.insert(TrianglePair(triangle->Nr, triangle));
+  TrianglesCount++;
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryLineSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "]";
+  return ost;
+}
+;
+
+// ======================================== Triangles on Boundary =================================
+
+
+BoundaryTriangleSet::BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
+  Nr = -1;
+}
+;
+
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+{
+  // set number
+  Nr = number;
+  // set lines
+  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  for (int i = 0; i < 3; i++)
+    {
+      lines[i] = line[i];
+      lines[i]->AddTriangle(this);
+    }
+  // get ascending order of endpoints
+  map<int, class BoundaryPointSet *> OrderMap;
+  for (int i = 0; i < 3; i++)
+    // for all three lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+        // and we don't care whether insertion fails
+      }
+  // set endpoints
+  int Counter = 0;
+  cout << Verbose(6) << " with end points ";
+  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+      != OrderMap.end(); runner++)
+    {
+      endpoints[Counter] = runner->second;
+      cout << " " << *endpoints[Counter];
+      Counter++;
+    }
+  if (Counter < 3)
+    {
+      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+          << endl;
+      //exit(1);
+    }
+  cout << "." << endl;
+}
+;
+
+BoundaryTriangleSet::~BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    if (lines[i]->triangles.empty()) {
+      if (lines[i] != NULL) {
+        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+        delete (lines[i]);
+        lines[i] = NULL;
+      } else
+        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
+  }
+}
+;
+
+void
+BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+{
+  // get normal vector
+  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
+      &endpoints[2]->node->x);
+
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (NormalVector.Projection(&OtherVector) > 0)
+    NormalVector.Scale(-1.);
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  return ost;
+}
+;
+
+
+// ============================ CandidateForTesselation =============================
+
+CandidateForTesselation::CandidateForTesselation(
+        atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
+) {
+        point = candidate;
+        BaseLine = line;
+        OptCenter.CopyVector(&OptCandidateCenter);
+        OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
+}
+
+CandidateForTesselation::~CandidateForTesselation() {
+        point = NULL;
+        BaseLine = NULL;
+}
 
 // ========================================== F U N C T I O N S =================================
 
+/** Finds the endpoint two lines are sharing.
+ * \param *line1 first line
+ * \param *line2 second line
+ * \return point which is shared or NULL if none
+ */
+class BoundaryPointSet *
+GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
+  class BoundaryLineSet * lines[2] =
+    { line1, line2 };
+  class BoundaryPointSet *node = NULL;
+  map<int, class BoundaryPointSet *> OrderMap;
+  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  for (int i = 0; i < 2; i++)
+    // for both lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+        if (!OrderTest.second)
+          { // if insertion fails, we have common endpoint
+            node = OrderTest.first->second;
+            cout << Verbose(5) << "Common endpoint of lines " << *line1
+                << " and " << *line2 << " is: " << *node << "." << endl;
+            j = 2;
+            i = 2;
+            break;
+          }
+      }
+  return node;
+}
+;
+
+/** Determines the boundary points of a cluster.
+ * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
+ * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
+ * center and first and last point in the triple, it is thrown out.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure representing the cluster
+ */
+Boundaries *
+GetBoundaryPoints(ofstream *out, molecule *mol)
+{
+  atom *Walker = NULL;
+  PointMap PointsOnBoundary;
+  LineMap LinesOnBoundary;
+  TriangleMap TrianglesOnBoundary;
+
+  *out << Verbose(1) << "Finding all boundary points." << endl;
+  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+  BoundariesTestPair BoundaryTestPair;
+  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+  double radius, angle;
+  // 3a. Go through every axis
+  for (int axis = 0; axis < NDIM; axis++)
+    {
+      AxisVector.Zero();
+      AngleReferenceVector.Zero();
+      AngleReferenceNormalVector.Zero();
+      AxisVector.x[axis] = 1.;
+      AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+      AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+      //    *out << Verbose(1) << "Axisvector is ";
+      //    AxisVector.Output(out);
+      //    *out << " and AngleReferenceVector is ";
+      //    AngleReferenceVector.Output(out);
+      //    *out << "." << endl;
+      //    *out << " and AngleReferenceNormalVector is ";
+      //    AngleReferenceNormalVector.Output(out);
+      //    *out << "." << endl;
+      // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+      Walker = mol->start;
+      while (Walker->next != mol->end)
+        {
+          Walker = Walker->next;
+          Vector ProjectedVector;
+          ProjectedVector.CopyVector(&Walker->x);
+          ProjectedVector.ProjectOntoPlane(&AxisVector);
+          // correct for negative side
+          //if (Projection(y) < 0)
+          //angle = 2.*M_PI - angle;
+          radius = ProjectedVector.Norm();
+          if (fabs(radius) > MYEPSILON)
+            angle = ProjectedVector.Angle(&AngleReferenceVector);
+          else
+            angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+
+          //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+          if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+            {
+              angle = 2. * M_PI - angle;
+            }
+          //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+          //ProjectedVector.Output(out);
+          //*out << endl;
+          BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
+              DistancePair (radius, Walker)));
+          if (BoundaryTestPair.second)
+            { // successfully inserted
+            }
+          else
+            { // same point exists, check first r, then distance of original vectors to center of gravity
+              *out << Verbose(2)
+                  << "Encountered two vectors whose projection onto axis "
+                  << axis << " is equal: " << endl;
+              *out << Verbose(2) << "Present vector: ";
+              BoundaryTestPair.first->second.second->x.Output(out);
+              *out << endl;
+              *out << Verbose(2) << "New vector: ";
+              Walker->x.Output(out);
+              *out << endl;
+              double tmp = ProjectedVector.Norm();
+              if (tmp > BoundaryTestPair.first->second.first)
+                {
+                  BoundaryTestPair.first->second.first = tmp;
+                  BoundaryTestPair.first->second.second = Walker;
+                  *out << Verbose(2) << "Keeping new vector." << endl;
+                }
+              else if (tmp == BoundaryTestPair.first->second.first)
+                {
+                  if (BoundaryTestPair.first->second.second->x.ScalarProduct(
+                      &BoundaryTestPair.first->second.second->x)
+                      < Walker->x.ScalarProduct(&Walker->x))
+                    { // Norm() does a sqrt, which makes it a lot slower
+                      BoundaryTestPair.first->second.second = Walker;
+                      *out << Verbose(2) << "Keeping new vector." << endl;
+                    }
+                  else
+                    {
+                      *out << Verbose(2) << "Keeping present vector." << endl;
+                    }
+                }
+              else
+                {
+                  *out << Verbose(2) << "Keeping present vector." << endl;
+                }
+            }
+        }
+      // printing all inserted for debugging
+      //    {
+      //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      //      int i=0;
+      //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+      //        if (runner != BoundaryPoints[axis].begin())
+      //          *out << ", " << i << ": " << *runner->second.second;
+      //        else
+      //          *out << i << ": " << *runner->second.second;
+      //        i++;
+      //      }
+      //      *out << endl;
+      //    }
+      // 3c. throw out points whose distance is less than the mean of left and right neighbours
+      bool flag = false;
+      do
+        { // do as long as we still throw one out per round
+          *out << Verbose(1)
+              << "Looking for candidates to kick out by convex condition ... "
+              << endl;
+          flag = false;
+          Boundaries::iterator left = BoundaryPoints[axis].end();
+          Boundaries::iterator right = BoundaryPoints[axis].end();
+          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+              != BoundaryPoints[axis].end(); runner++)
+            {
+              // set neighbours correctly
+              if (runner == BoundaryPoints[axis].begin())
+                {
+                  left = BoundaryPoints[axis].end();
+                }
+              else
+                {
+                  left = runner;
+                }
+              left--;
+              right = runner;
+              right++;
+              if (right == BoundaryPoints[axis].end())
+                {
+                  right = BoundaryPoints[axis].begin();
+                }
+              // check distance
+
+              // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                {
+                  Vector SideA, SideB, SideC, SideH;
+                  SideA.CopyVector(&left->second.second->x);
+                  SideA.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideA: ";
+                  //          SideA.Output(out);
+                  //          *out << endl;
+
+                  SideB.CopyVector(&right->second.second->x);
+                  SideB.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideB: ";
+                  //          SideB.Output(out);
+                  //          *out << endl;
+
+                  SideC.CopyVector(&left->second.second->x);
+                  SideC.SubtractVector(&right->second.second->x);
+                  SideC.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideC: ";
+                  //          SideC.Output(out);
+                  //          *out << endl;
+
+                  SideH.CopyVector(&runner->second.second->x);
+                  SideH.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideH: ";
+                  //          SideH.Output(out);
+                  //          *out << endl;
+
+                  // calculate each length
+                  double a = SideA.Norm();
+                  //double b = SideB.Norm();
+                  //double c = SideC.Norm();
+                  double h = SideH.Norm();
+                  // calculate the angles
+                  double alpha = SideA.Angle(&SideH);
+                  double beta = SideA.Angle(&SideC);
+                  double gamma = SideB.Angle(&SideH);
+                  double delta = SideC.Angle(&SideH);
+                  double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
+                      < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+                  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+                  //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+                  if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
+                      < MYEPSILON) && (h < MinDistance))
+                    {
+                      // throw out point
+                      //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+                      BoundaryPoints[axis].erase(runner);
+                      flag = true;
+                    }
+                }
+            }
+        }
+      while (flag);
+    }
+  return BoundaryPoints;
+}
+;
 
 /** Determines greatest diameters of a cluster defined by its convex envelope.
@@ -27 +510 @@
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
-  if (BoundaryPoints == NULL) {
-    BoundaryFreeFlag = true;
-    BoundaryPoints = GetBoundaryPoints(out, mol);
-  } else {
-    *out << Verbose(1) << "Using given boundary points set." << endl;
-  }
+  if (BoundaryPoints == NULL)
+    {
+      BoundaryFreeFlag = true;
+      BoundaryPoints = GetBoundaryPoints(out, mol);
+    }
+  else
+    {
+      *out << Verbose(1) << "Using given boundary points set." << endl;
+    }
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -133 +619 @@
       *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
+        *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
       *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
     }
@@ -142 +628 @@
       *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
       for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
+        *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
       *vrmlfile << "\t0.03\t";
       for (i=0;i<NDIM;i++)
-        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
+        *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
       *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
     }
@@ -154 +640 @@
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
         *vrmlfile << "\t";
       }
@@ -187 +673 @@
       *rasterfile << "2" << endl << "  ";  // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
+        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
       *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
     }
@@ -196 +682 @@
       *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
       for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
+        *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
       *rasterfile << "\t0.03\t";
       for (i=0;i<NDIM;i++)
-        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
+        *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
       *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
     }
@@ -209 +695 @@
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
         *rasterfile << "\t";
       }
@@ -227 +713 @@
  * \param N arbitrary number to differentiate various zones in the tecplot format
  */
-void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
-{
-  if (tecplot != NULL) {
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int *LookupList = new int[mol->AtomCount];
-    for (int i = 0; i < mol->AtomCount; i++)
-      LookupList[i] = -1;
-
-    // print atom coordinates
-    *out << Verbose(2) << "The following triangles were created:";
-    int Counter = 1;
-    TesselPoint *Walker = NULL;
-    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
-      Walker = target->second->node;
-      LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
-    }
-    *tecplot << endl;
-    // print connectivity
-    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
-      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
-    }
-    delete[] (LookupList);
-    *out << endl;
-  }
-}
-
-
-/** Determines the boundary points of a cluster.
- * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
- * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
- * center and first and last point in the triple, it is thrown out.
+void
+write_tecplot_file(ofstream *out, ofstream *tecplot,
+    class Tesselation *TesselStruct, class molecule *mol, int N)
+{
+  if (tecplot != NULL)
+    {
+      *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+      *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+      *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N="
+          << TesselStruct->PointsOnBoundaryCount << ", E="
+          << TesselStruct->TrianglesOnBoundaryCount
+          << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+      int *LookupList = new int[mol->AtomCount];
+      for (int i = 0; i < mol->AtomCount; i++)
+        LookupList[i] = -1;
+
+      // print atom coordinates
+      *out << Verbose(2) << "The following triangles were created:";
+      int Counter = 1;
+      atom *Walker = NULL;
+      for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target
+          != TesselStruct->PointsOnBoundary.end(); target++)
+        {
+          Walker = target->second->node;
+          LookupList[Walker->nr] = Counter++;
+          *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " "
+              << Walker->x.x[2] << " " << endl;
+        }
+      *tecplot << endl;
+      // print connectivity
+      for (TriangleMap::iterator runner =
+          TesselStruct->TrianglesOnBoundary.begin(); runner
+          != TesselStruct->TrianglesOnBoundary.end(); runner++)
+        {
+          *out << " " << runner->second->endpoints[0]->node->Name << "<->"
+              << runner->second->endpoints[1]->node->Name << "<->"
+              << runner->second->endpoints[2]->node->Name;
+          *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " "
+              << LookupList[runner->second->endpoints[1]->node->nr] << " "
+              << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+        }
+      delete[] (LookupList);
+      *out << endl;
+    }
+}
+
+/** Determines the volume of a cluster.
+ * Determines first the convex envelope, then tesselates it and calculates its volume.
  * \param *out output stream for debugging
+ * \param *filename filename prefix for output of vertex data
+ * \param *configuration needed for path to store convex envelope file
+ * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
  * \param *mol molecule structure representing the cluster
- */
-Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
-{
+ * \return determined volume of the cluster in cubed config:GetIsAngstroem()
+ */
+double
+VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
+    Boundaries *BoundaryPtr, molecule *mol)
+{
+  bool IsAngstroem = configuration->GetIsAngstroem();
   atom *Walker = NULL;
-  PointMap PointsOnBoundary;
-  LineMap LinesOnBoundary;
-  TriangleMap TrianglesOnBoundary;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  Vector helper;
-
-  *out << Verbose(1) << "Finding all boundary points." << endl;
-  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
-  BoundariesTestPair BoundaryTestPair;
-  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
-  double radius, angle;
-  // 3a. Go through every axis
-  for (int axis = 0; axis < NDIM; axis++) {
-    AxisVector.Zero();
-    AngleReferenceVector.Zero();
-    AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
-
-    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
-
-    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
+  struct Tesselation *TesselStruct = new Tesselation;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G, h;
+  Vector x, y;
+  double a, b, c;
+
+  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
+
+  // 1. calculate center of gravity
+  *out << endl;
+  Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
+
+  // 2. translate all points into CoG
+  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
       Walker = Walker->next;
-      Vector ProjectedVector;
-      ProjectedVector.CopyVector(&Walker->x);
-      ProjectedVector.SubtractVector(MolCenter);
-      ProjectedVector.ProjectOntoPlane(&AxisVector);
-
-      // correct for negative side
-      radius = ProjectedVector.NormSquared();
-      if (fabs(radius) > MYEPSILON)
-        angle = ProjectedVector.Angle(&AngleReferenceVector);
-      else
-        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
-
-      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
-        angle = 2. * M_PI - angle;
-      }
-      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
-      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
-        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
-        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
-        *out << Verbose(2) << "New vector: " << *Walker << endl;
-        double tmp = ProjectedVector.NormSquared();
-        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
-          BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker;
-          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
-        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
-          helper.SubtractVector(MolCenter);
-          tmp = helper.NormSquared();
-          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
-          helper.SubtractVector(MolCenter);
-          if (helper.NormSquared() < tmp) {
-            BoundaryTestPair.first->second.second = Walker;
-            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
-          } else {
-            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
-          }
-        } else {
-          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
-        }
-      }
-    }
-    // printing all inserted for debugging
-    //    {
-    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
-    //      int i=0;
-    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-    //        if (runner != BoundaryPoints[axis].begin())
-    //          *out << ", " << i << ": " << *runner->second.second;
-    //        else
-    //          *out << i << ": " << *runner->second.second;
-    //        i++;
-    //      }
-    //      *out << endl;
-    //    }
-    // 3c. throw out points whose distance is less than the mean of left and right neighbours
-    bool flag = false;
-    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
-    do { // do as long as we still throw one out per round
-      flag = false;
-      Boundaries::iterator left = BoundaryPoints[axis].end();
-      Boundaries::iterator right = BoundaryPoints[axis].end();
-      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        // set neighbours correctly
-        if (runner == BoundaryPoints[axis].begin()) {
-          left = BoundaryPoints[axis].end();
-        } else {
-          left = runner;
-        }
-        left--;
-        right = runner;
-        right++;
-        if (right == BoundaryPoints[axis].end()) {
-          right = BoundaryPoints[axis].begin();
-        }
-        // check distance
-
-        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
-        {
-          Vector SideA, SideB, SideC, SideH;
-          SideA.CopyVector(&left->second.second->x);
-          SideA.SubtractVector(MolCenter);
-          SideA.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideA: ";
-          //          SideA.Output(out);
-          //          *out << endl;
-
-          SideB.CopyVector(&right->second.second->x);
-          SideB.SubtractVector(MolCenter);
-          SideB.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideB: ";
-          //          SideB.Output(out);
-          //          *out << endl;
-
-          SideC.CopyVector(&left->second.second->x);
-          SideC.SubtractVector(&right->second.second->x);
-          SideC.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideC: ";
-          //          SideC.Output(out);
-          //          *out << endl;
-
-          SideH.CopyVector(&runner->second.second->x);
-          SideH.SubtractVector(MolCenter);
-          SideH.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideH: ";
-          //          SideH.Output(out);
-          //          *out << endl;
-
-          // calculate each length
-          double a = SideA.Norm();
-          //double b = SideB.Norm();
-          //double c = SideC.Norm();
-          double h = SideH.Norm();
-          // calculate the angles
-          double alpha = SideA.Angle(&SideH);
-          double beta = SideA.Angle(&SideC);
-          double gamma = SideB.Angle(&SideH);
-          double delta = SideC.Angle(&SideH);
-          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
-          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
-          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
-          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
-            // throw out point
-            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
-            BoundaryPoints[axis].erase(runner);
-            flag = true;
-          }
-        }
-      }
-    } while (flag);
-  }
-  delete(MolCenter);
-  return BoundaryPoints;
-};
-
-/** Tesselates the convex boundary by finding all boundary points.
- * \param *out output stream for debugging
- * \param *mol molecule structure with Atom's and Bond's
- * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
- * \param *LCList atoms in LinkedCell list
- * \param *filename filename prefix for output of vertex data
- * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
- */
-void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
-{
-  bool BoundaryFreeFlag = false;
-  Boundaries *BoundaryPoints = NULL;
-
-  cout << Verbose(1) << "Begin of find_convex_border" << endl;
-
-  if (TesselStruct != NULL) // free if allocated
-    delete(TesselStruct);
-  TesselStruct = new class Tesselation;
-
-  // 1. Find all points on the boundary
-  if (BoundaryPoints == NULL) {
+      Walker->x.Translate(CenterOfGravity);
+    }
+
+  // 3. Find all points on the boundary
+  if (BoundaryPoints == NULL)
+    {
       BoundaryFreeFlag = true;
       BoundaryPoints = GetBoundaryPoints(out, mol);
-  } else {
+    }
+  else
+    {
       *out << Verbose(1) << "Using given boundary points set." << endl;
-  }
-
-// printing all inserted for debugging
-  for (int axis=0; axis < NDIM; axis++)
-    {
-      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
-      int i=0;
-      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        if (runner != BoundaryPoints[axis].begin())
-          *out << ", " << i << ": " << *runner->second.second;
-        else
-          *out << i << ": " << *runner->second.second;
-        i++;
+    }
+
+  // 4. fill the boundary point list
+  for (int axis = 0; axis < NDIM; axis++)
+    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+        != BoundaryPoints[axis].end(); runner++)
+      {
+        TesselStruct->AddPoint(runner->second.second);
       }
-      *out << endl;
-    }
-
-  // 2. fill the boundary point list
-  for (int axis = 0; axis < NDIM; axis++)
-    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
-        TesselStruct->AddPoint(runner->second.second);
-
-  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
+
+  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
+      << " points on the convex boundary." << endl;
   // now we have the whole set of edge points in the BoundaryList
 
@@ -488 +828 @@
   //  *out << endl;
 
-  // 3a. guess starting triangle
+  // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
 
-  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
-  TesselStruct->TesselateOnBoundary(out, mol);
-
-  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
-  if (!TesselStruct->InsertStraddlingPoints(out, mol))
-    *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
-
-  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
-  if (!TesselStruct->CorrectConcaveBaselines(out))
-    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
-
-  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
-
-  // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      write_raster3d_file(out, rasterplot, TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
-
-  // free reference lists
-  if (BoundaryFreeFlag)
-    delete[] (BoundaryPoints);
-
-  cout << Verbose(1) << "End of find_convex_border" << endl;
-};
-
-
-/** Determines the volume of a cluster.
- * Determines first the convex envelope, then tesselates it and calculates its volume.
- * \param *out output stream for debugging
- * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
- * \param *configuration needed for path to store convex envelope file
- * \return determined volume of the cluster in cubed config:GetIsAngstroem()
- */
-double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
-{
-  bool IsAngstroem = configuration->GetIsAngstroem();
-  double volume = 0.;
-  double PyramidVolume = 0.;
-  double G, h;
-  Vector x, y;
-  double a, b, c;
+  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
+  TesselStruct->TesselateOnBoundary(out, configuration, mol);
+
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
+      << " triangles with " << TesselStruct->LinesOnBoundaryCount
+      << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
+      << endl;
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
@@ -552 +843 @@
       << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
       << endl;
-  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
+  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.CopyVector(runner->second->endpoints[0]->node->node);
-      x.SubtractVector(runner->second->endpoints[1]->node->node);
-      y.CopyVector(runner->second->endpoints[0]->node->node);
-      y.SubtractVector(runner->second->endpoints[2]->node->node);
-      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
-      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
-      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
+      x.CopyVector(&runner->second->endpoints[0]->node->x);
+      x.SubtractVector(&runner->second->endpoints[1]->node->x);
+      y.CopyVector(&runner->second->endpoints[0]->node->x);
+      y.SubtractVector(&runner->second->endpoints[2]->node->x);
+      a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
+      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+          &runner->second->endpoints[2]->node->x));
+      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
       G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
-      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
-      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
+      x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
+          &runner->second->endpoints[1]->node->x,
+          &runner->second->endpoints[2]->node->x);
+      x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
       h = x.Norm(); // distance of CoG to triangle
       PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
@@ -576 +873 @@
       << endl;
 
+  // 7. translate all points back from CoG
+  *out << Verbose(1) << "Translating system back from Center of Gravity."
+      << endl;
+  CenterOfGravity->Scale(-1);
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
+
+  // 8. Store triangles in tecplot file
+  string OutputName(filename);
+  OutputName.append(TecplotSuffix);
+  ofstream *tecplot = new ofstream(OutputName.c_str());
+  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+  tecplot->close();
+  delete(tecplot);
+
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[] (BoundaryPoints);
+
   return volume;
 }
@@ -598 +918 @@
   bool IsAngstroem = configuration->GetIsAngstroem();
   Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
-  class Tesselation *TesselStruct = NULL;
-  LinkedCell LCList(mol, 10.);
-  Find_convex_border(out, mol, TesselStruct, &LCList, NULL);
   double clustervolume;
   if (ClusterVolume == 0)
-    clustervolume = VolumeOfConvexEnvelope(out, TesselStruct, configuration);
+    clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
+        BoundaryPoints, mol);
   else
     clustervolume = ClusterVolume;
-  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
+      IsAngstroem);
   Vector BoxLengths;
   int repetition[NDIM] =
@@ -691 +1010 @@
       // set new box dimensions
       *out << Verbose(0) << "Translating to box with these boundaries." << endl;
-      mol->SetBoxDimension(&BoxLengths);
-      mol->CenterInBox((ofstream *) &cout);
+      mol->CenterInBox((ofstream *) &cout, &BoxLengths);
     }
   // update Box of atoms by boundary
@@ -703 +1021 @@
 ;
 
-
-/** Fills the empty space of the simulation box with water/
+// =========================================================== class TESSELATION ===========================================
+
+/** Constructor of class Tesselation.
+ */
+Tesselation::Tesselation()
+{
+  PointsOnBoundaryCount = 0;
+  LinesOnBoundaryCount = 0;
+  TrianglesOnBoundaryCount = 0;
+  TriangleFilesWritten = 0;
+}
+;
+
+/** Constructor of class Tesselation.
+ * We have to free all points, lines and triangles.
+ */
+Tesselation::~Tesselation()
+{
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+}
+;
+
+/** Gueses first starting triangle of the convex envelope.
+ * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
  * \param *out output stream for debugging
- * \param *List list of molecules already present in box
- * \param *TesselStruct contains tesselated surface
- * \param *filler molecule which the box is to be filled with
- * \param configuration contains box dimensions
- * \param distance[NDIM] distance between filling molecules in each direction
- * \param RandAtomDisplacement maximum distance for random displacement per atom
- * \param RandMolDisplacement maximum distance for random displacement per filler molecule
- * \param DoRandomRotation true - do random rotiations, false - don't
- * \return *mol pointer to new molecule with filled atoms
- */
-molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
-{
-  molecule *Filling = new molecule(filler->elemente);
-  Vector CurrentPosition;
-  int N[NDIM];
-  int n[NDIM];
-  double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
-  double Rotations[NDIM*NDIM];
-  Vector AtomTranslations;
-  Vector FillerTranslations;
-  Vector FillerDistance;
-  double FillIt = false;
-  atom *Walker = NULL, *Runner = NULL;
-  bond *Binder = NULL;
-
-  // Center filler at origin
-  filler->CenterOrigin(out);
-  filler->Center.Zero();
-
-  // calculate filler grid in [0,1]^3
-  FillerDistance.Init(distance[0], distance[1], distance[2]);
-  FillerDistance.InverseMatrixMultiplication(M);
-  for(int i=0;i<NDIM;i++)
-    N[i] = ceil(1./FillerDistance.x[i]);
-
-  // go over [0,1]^3 filler grid
-  for (n[0] = 0; n[0] < N[0]; n[0]++)
-    for (n[1] = 0; n[1] < N[1]; n[1]++)
-      for (n[2] = 0; n[2] < N[2]; n[2]++) {
-        // calculate position of current grid vector in untransformed box
-        CurrentPosition.Init(n[0], n[1], n[2]);
-        CurrentPosition.MatrixMultiplication(M);
-        // Check whether point is in- or outside
-        FillIt = true;
-        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
-          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+ * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
+ */
+void
+Tesselation::GuessStartingTriangle(ofstream *out)
+{
+  // 4b. create a starting triangle
+  // 4b1. create all distances
+  DistanceMultiMap DistanceMMap;
+  double distance, tmp;
+  Vector PlaneVector, TrialVector;
+  PointMap::iterator A, B, C; // three nodes of the first triangle
+  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
+
+  // with A chosen, take each pair B,C and sort
+  if (A != PointsOnBoundary.end())
+    {
+      B = A;
+      B++;
+      for (; B != PointsOnBoundary.end(); B++)
+        {
+          C = B;
+          C++;
+          for (; C != PointsOnBoundary.end(); C++)
+            {
+              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
+              distance = tmp * tmp;
+              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
+              distance += tmp * tmp;
+              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
+              distance += tmp * tmp;
+              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
+                  PointMap::iterator, PointMap::iterator> (B, C)));
+            }
         }
-
-        if (FillIt) {
-          // fill in Filler
-          *out << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
-
-          // create molecule random translation vector ...
-          for (int i=0;i<NDIM;i++)
-            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-
-          // go through all atoms
-          Walker = filler->start;
-          while (Walker != filler->end) {
-            Walker = Walker->next;
-            // copy atom ...
-            *Runner = new atom(Walker);
-
-            // create atomic random translation vector ...
-            for (int i=0;i<NDIM;i++)
-              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
-
-            // ... and rotation matrix
-            if (DoRandomRotation) {
-              double phi[NDIM];
-              for (int i=0;i<NDIM;i++) {
-                phi[i] = rand()/(RAND_MAX/(2.*M_PI));
+    }
+  //    // listing distances
+  //    *out << Verbose(1) << "Listing DistanceMMap:";
+  //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
+  //      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
+  //    }
+  //    *out << endl;
+  // 4b2. pick three baselines forming a triangle
+  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
+  for (; baseline != DistanceMMap.end(); baseline++)
+    {
+      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+      // 2. next, we have to check whether all points reside on only one side of the triangle
+      // 3. construct plane vector
+      PlaneVector.MakeNormalVector(&A->second->node->x,
+          &baseline->second.first->second->node->x,
+          &baseline->second.second->second->node->x);
+      *out << Verbose(2) << "Plane vector of candidate triangle is ";
+      PlaneVector.Output(out);
+      *out << endl;
+      // 4. loop over all points
+      double sign = 0.;
+      PointMap::iterator checker = PointsOnBoundary.begin();
+      for (; checker != PointsOnBoundary.end(); checker++)
+        {
+          // (neglecting A,B,C)
+          if ((checker == A) || (checker == baseline->second.first) || (checker
+              == baseline->second.second))
+            continue;
+          // 4a. project onto plane vector
+          TrialVector.CopyVector(&checker->second->node->x);
+          TrialVector.SubtractVector(&A->second->node->x);
+          distance = TrialVector.Projection(&PlaneVector);
+          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+            continue;
+          *out << Verbose(3) << "Projection of " << checker->second->node->Name
+              << " yields distance of " << distance << "." << endl;
+          tmp = distance / fabs(distance);
+          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+          if ((sign != 0) && (tmp != sign))
+            {
+              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " outside the convex hull."
+                  << endl;
+              break;
+            }
+          else
+            { // note the sign for later
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " inside the convex hull."
+                  << endl;
+              sign = tmp;
+            }
+          // 4d. Check whether the point is inside the triangle (check distance to each node
+          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
+          int innerpoint = 0;
+          if ((tmp < A->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x)) && (tmp
+              < A->second->node->x.DistanceSquared(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x);
+          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
+              &A->second->node->x)) && (tmp
+              < baseline->second.first->second->node->x.DistanceSquared(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.DistanceSquared(
+              &baseline->second.second->second->node->x);
+          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x)) && (tmp
+              < baseline->second.second->second->node->x.DistanceSquared(
+                  &A->second->node->x)))
+            innerpoint++;
+          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+          if (innerpoint == 3)
+            break;
+        }
+      // 5. come this far, all on same side? Then break 1. loop and construct triangle
+      if (checker == PointsOnBoundary.end())
+        {
+          *out << "Looks like we have a candidate!" << endl;
+          break;
+        }
+    }
+  if (baseline != DistanceMMap.end())
+    {
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = A->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = A->second;
+      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+      // 4b3. insert created triangle
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+      TrianglesOnBoundaryCount++;
+      for (int i = 0; i < NDIM; i++)
+        {
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+          LinesOnBoundaryCount++;
+        }
+
+      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
+  else
+    {
+      *out << Verbose(1) << "No starting triangle found." << endl;
+      exit(255);
+    }
+}
+;
+
+/** Tesselates the convex envelope of a cluster from a single starting triangle.
+ * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
+ * 2 triangles. Hence, we go through all current lines:
+ * -# if the lines contains to only one triangle
+ * -# We search all points in the boundary
+ *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
+ *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
+ *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
+ *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
+ * \param *out output stream for debugging
+ * \param *configuration for IsAngstroem
+ * \param *mol the cluster as a molecule structure
+ */
+void
+Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
+    molecule *mol)
+{
+  bool flag;
+  PointMap::iterator winner;
+  class BoundaryPointSet *peak = NULL;
+  double SmallestAngle, TempAngle;
+  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
+      PropagationVector;
+  LineMap::iterator LineChecker[2];
+  do
+    {
+      flag = false;
+      for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
+          != LinesOnBoundary.end(); baseline++)
+        if (baseline->second->TrianglesCount == 1)
+          {
+            *out << Verbose(2) << "Current baseline is between "
+                << *(baseline->second) << "." << endl;
+            // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
+            SmallestAngle = M_PI;
+            BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+            // get peak point with respect to this base line's only triangle
+            for (int i = 0; i < 3; i++)
+              if ((BTS->endpoints[i] != baseline->second->endpoints[0])
+                  && (BTS->endpoints[i] != baseline->second->endpoints[1]))
+                peak = BTS->endpoints[i];
+            *out << Verbose(3) << " and has peak " << *peak << "." << endl;
+            // normal vector of triangle
+            BTS->GetNormalVector(NormalVector);
+            *out << Verbose(4) << "NormalVector of base triangle is ";
+            NormalVector.Output(out);
+            *out << endl;
+            // offset to center of triangle
+            CenterVector.Zero();
+            for (int i = 0; i < 3; i++)
+              CenterVector.AddVector(&BTS->endpoints[i]->node->x);
+            CenterVector.Scale(1. / 3.);
+            *out << Verbose(4) << "CenterVector of base triangle is ";
+            CenterVector.Output(out);
+            *out << endl;
+            // vector in propagation direction (out of triangle)
+            // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+            TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
+            PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
+            TempVector.CopyVector(&CenterVector);
+            TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+            //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
+            if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+              PropagationVector.Scale(-1.);
+            *out << Verbose(4) << "PropagationVector of base triangle is ";
+            PropagationVector.Output(out);
+            *out << endl;
+            winner = PointsOnBoundary.end();
+            for (PointMap::iterator target = PointsOnBoundary.begin(); target
+                != PointsOnBoundary.end(); target++)
+              if ((target->second != baseline->second->endpoints[0])
+                  && (target->second != baseline->second->endpoints[1]))
+                { // don't take the same endpoints
+                  *out << Verbose(3) << "Target point is " << *(target->second)
+                      << ":";
+                  bool continueflag = true;
+
+                  VirtualNormalVector.CopyVector(
+                      &baseline->second->endpoints[0]->node->x);
+                  VirtualNormalVector.AddVector(
+                      &baseline->second->endpoints[0]->node->x);
+                  VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
+                  VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
+                  TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+                  continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4)
+                          << "Angle between propagation direction and base line to "
+                          << *(target->second) << " is " << TempAngle
+                          << ", bad direction!" << endl;
+                      continue;
+                    }
+                  else
+                    *out << Verbose(4)
+                        << "Angle between propagation direction and base line to "
+                        << *(target->second) << " is " << TempAngle
+                        << ", good direction!" << endl;
+                  LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                      target->first);
+                  LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                      target->first);
+                  //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+                  //            else
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                  //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+                  //            else
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                  // check first endpoint (if any connecting line goes to target or at least not more than 1)
+                  continueflag = continueflag && (((LineChecker[0]
+                      == baseline->second->endpoints[0]->lines.end())
+                      || (LineChecker[0]->second->TrianglesCount == 1)));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << *(baseline->second->endpoints[0])
+                          << " has line " << *(LineChecker[0]->second)
+                          << " to " << *(target->second)
+                          << " as endpoint with "
+                          << LineChecker[0]->second->TrianglesCount
+                          << " triangles." << endl;
+                      continue;
+                    }
+                  // check second endpoint (if any connecting line goes to target or at least not more than 1)
+                  continueflag = continueflag && (((LineChecker[1]
+                      == baseline->second->endpoints[1]->lines.end())
+                      || (LineChecker[1]->second->TrianglesCount == 1)));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << *(baseline->second->endpoints[1])
+                          << " has line " << *(LineChecker[1]->second)
+                          << " to " << *(target->second)
+                          << " as endpoint with "
+                          << LineChecker[1]->second->TrianglesCount
+                          << " triangles." << endl;
+                      continue;
+                    }
+                  // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
+                  continueflag = continueflag && (!(((LineChecker[0]
+                      != baseline->second->endpoints[0]->lines.end())
+                      && (LineChecker[1]
+                          != baseline->second->endpoints[1]->lines.end())
+                      && (GetCommonEndpoint(LineChecker[0]->second,
+                          LineChecker[1]->second) == peak))));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << "Current target is peak!" << endl;
+                      continue;
+                    }
+                  // in case NOT both were found
+                  if (continueflag)
+                    { // create virtually this triangle, get its normal vector, calculate angle
+                      flag = true;
+                      VirtualNormalVector.MakeNormalVector(
+                          &baseline->second->endpoints[0]->node->x,
+                          &baseline->second->endpoints[1]->node->x,
+                          &target->second->node->x);
+                      // make it always point inward
+                      if (baseline->second->endpoints[0]->node->x.Projection(
+                          &VirtualNormalVector) > 0)
+                        VirtualNormalVector.Scale(-1.);
+                      // calculate angle
+                      TempAngle = NormalVector.Angle(&VirtualNormalVector);
+                      *out << Verbose(4) << "NormalVector is ";
+                      VirtualNormalVector.Output(out);
+                      *out << " and the angle is " << TempAngle << "." << endl;
+                      if (SmallestAngle > TempAngle)
+                        { // set to new possible winner
+                          SmallestAngle = TempAngle;
+                          winner = target;
+                        }
+                    }
+                }
+            // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
+            if (winner != PointsOnBoundary.end())
+              {
+                *out << Verbose(2) << "Winning target point is "
+                    << *(winner->second) << " with angle " << SmallestAngle
+                    << "." << endl;
+                // create the lins of not yet present
+                BLS[0] = baseline->second;
+                // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
+                LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                    winner->first);
+                LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                    winner->first);
+                if (LineChecker[0]
+                    == baseline->second->endpoints[0]->lines.end())
+                  { // create
+                    BPS[0] = baseline->second->endpoints[0];
+                    BPS[1] = winner->second;
+                    BLS[1] = new class BoundaryLineSet(BPS,
+                        LinesOnBoundaryCount);
+                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                        BLS[1]));
+                    LinesOnBoundaryCount++;
+                  }
+                else
+                  BLS[1] = LineChecker[0]->second;
+                if (LineChecker[1]
+                    == baseline->second->endpoints[1]->lines.end())
+                  { // create
+                    BPS[0] = baseline->second->endpoints[1];
+                    BPS[1] = winner->second;
+                    BLS[2] = new class BoundaryLineSet(BPS,
+                        LinesOnBoundaryCount);
+                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                        BLS[2]));
+                    LinesOnBoundaryCount++;
+                  }
+                else
+                  BLS[2] = LineChecker[1]->second;
+                BTS = new class BoundaryTriangleSet(BLS,
+                    TrianglesOnBoundaryCount);
+                TrianglesOnBoundary.insert(TrianglePair(
+                    TrianglesOnBoundaryCount, BTS));
+                TrianglesOnBoundaryCount++;
               }
-
-              Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
-              Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
-              Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
-              Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
-              Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
-              Rotations[7] =               sin(phi[1])                                                ;
-              Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
-              Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
-              Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            else
+              {
+                *out << Verbose(1)
+                    << "I could not determine a winner for this baseline "
+                    << *(baseline->second) << "." << endl;
+              }
+
+            // 5d. If the set of lines is not yet empty, go to 5. and continue
+          }
+        else
+          *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
+              << " has a triangle count of "
+              << baseline->second->TrianglesCount << "." << endl;
+    }
+  while (flag);
+
+}
+;
+
+/** Adds an atom to the tesselation::PointsOnBoundary list.
+ * \param *Walker atom to add
+ */
+void
+Tesselation::AddPoint(atom *Walker)
+{
+  PointTestPair InsertUnique;
+  BPS[0] = new class BoundaryPointSet(Walker);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
+  if (InsertUnique.second) // if new point was not present before, increase counter
+    PointsOnBoundaryCount++;
+}
+;
+
+/** Adds point to Tesselation::PointsOnBoundary if not yet present.
+ * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
+ * @param Candidate point to add
+ * @param n index for this point in Tesselation::TPS array
+ */
+void
+Tesselation::AddTrianglePoint(atom* Candidate, int n)
+{
+  PointTestPair InsertUnique;
+  TPS[n] = new class BoundaryPointSet(Candidate);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
+  if (InsertUnique.second) { // if new point was not present before, increase counter
+    PointsOnBoundaryCount++;
+  } else {
+    delete TPS[n];
+    cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    TPS[n] = (InsertUnique.first)->second;
+  }
+}
+;
+
+/** Function tries to add line from current Points in BPS to BoundaryLineSet.
+ * If successful it raises the line count and inserts the new line into the BLS,
+ * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
+ * @param *a first endpoint
+ * @param *b second endpoint
+ * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
+  bool insertNewLine = true;
+
+  if (a->lines.find(b->node->nr) != a->lines.end()) {
+    LineMap::iterator FindLine;
+    pair<LineMap::iterator,LineMap::iterator> FindPair;
+    FindPair = a->lines.equal_range(b->node->nr);
+
+    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+      // If there is a line with less than two attached triangles, we don't need a new line.
+      if (FindLine->second->TrianglesCount < 2) {
+        insertNewLine = false;
+        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
+
+        BPS[0] = FindLine->second->endpoints[0];
+        BPS[1] = FindLine->second->endpoints[1];
+        BLS[n] = FindLine->second;
+
+        break;
+      }
+    }
+  }
+
+  if (insertNewLine) {
+    AlwaysAddTriangleLine(a, b, n);
+  }
+}
+;
+
+/**
+ * Adds lines from each of the current points in the BPS to BoundaryLineSet.
+ * Raises the line count and inserts the new line into the BLS.
+ *
+ * @param *a first endpoint
+ * @param *b second endpoint
+ * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+{
+  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  BPS[0] = a;
+  BPS[1] = b;
+  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
+  // add line to global map
+  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
+  // increase counter
+  LinesOnBoundaryCount++;
+};
+
+/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
+ * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
+ */
+void
+Tesselation::AddTriangle()
+{
+  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
+
+  // add triangle to global map
+  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+  TrianglesOnBoundaryCount++;
+
+  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
+}
+;
+
+
+double det_get(gsl_matrix *A, int inPlace) {
+  /*
+  inPlace = 1 => A is replaced with the LU decomposed copy.
+  inPlace = 0 => A is retained, and a copy is used for LU.
+  */
+
+  double det;
+  int signum;
+  gsl_permutation *p = gsl_permutation_alloc(A->size1);
+  gsl_matrix *tmpA;
+
+  if (inPlace)
+  tmpA = A;
+  else {
+  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+  gsl_matrix_memcpy(tmpA , A);
+  }
+
+
+  gsl_linalg_LU_decomp(tmpA , p , &signum);
+  det = gsl_linalg_LU_det(tmpA , signum);
+  gsl_permutation_free(p);
+  if (! inPlace)
+  gsl_matrix_free(tmpA);
+
+  return det;
+};
+
+void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
+  gsl_matrix *A = gsl_matrix_calloc(3,3);
+  double m11, m12, m13, m14;
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m11 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m12 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m13 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, b.x[i]);
+  }
+  m14 = det_get(A, 1);
+
+  if (fabs(m11) < MYEPSILON)
+    cerr << "ERROR: three points are colinear." << endl;
+
+  center->x[0] =  0.5 * m12/ m11;
+  center->x[1] = -0.5 * m13/ m11;
+  center->x[2] =  0.5 * m14/ m11;
+
+  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
+    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+
+  gsl_matrix_free(A);
+};
+
+
+
+/**
+ * Function returns center of sphere with RADIUS, which rests on points a, b, c
+ * @param Center this vector will be used for return
+ * @param a vector first point of triangle
+ * @param b vector second point of triangle
+ * @param c vector third point of triangle
+ * @param *Umkreismittelpunkt new cneter point of circumference
+ * @param Direction vector indicates up/down
+ * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
+ * @param Halfplaneindicator double indicates whether Direction is up or down
+ * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
+ * @param alpha double angle at a
+ * @param beta double, angle at b
+ * @param gamma, double, angle at c
+ * @param Radius, double
+ * @param Umkreisradius double radius of circumscribing circle
+ */
+void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
+    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
+  Vector TempNormal, helper;
+  double Restradius;
+  Vector OtherCenter;
+  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+  Center->Zero();
+  helper.CopyVector(&a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(&b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(&c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+  NewUmkreismittelpunkt->CopyVector(Center);
+  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
+  // Here we calculated center of circumscribing circle, using barycentric coordinates
+  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
+
+  TempNormal.CopyVector(&a);
+  TempNormal.SubtractVector(&b);
+  helper.CopyVector(&a);
+  helper.SubtractVector(&c);
+  TempNormal.VectorProduct(&helper);
+  if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  else
+    {
+      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+
+  TempNormal.Normalize();
+  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
+  TempNormal.Scale(Restradius);
+  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
+
+  Center->AddVector(&TempNormal);
+  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
+  get_sphere(&OtherCenter, a, b, c, RADIUS);
+  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
+  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
+};
+
+
+/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
+ * \param *Center new center on return
+ * \param *a first point
+ * \param *b second point
+ * \param *c third point
+ */
+void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
+  Vector helper;
+  double alpha, beta, gamma;
+  Vector SideA, SideB, SideC;
+  SideA.CopyVector(b);
+  SideA.SubtractVector(c);
+  SideB.CopyVector(c);
+  SideB.SubtractVector(a);
+  SideC.CopyVector(a);
+  SideC.SubtractVector(b);
+  alpha = M_PI - SideB.Angle(&SideC);
+  beta = M_PI - SideC.Angle(&SideA);
+  gamma = M_PI - SideA.Angle(&SideB);
+  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
+  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
+    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+
+  Center->Zero();
+  helper.CopyVector(a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+};
+
+/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
+ * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
+ * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
+ * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
+ * \param CircleCenter Center of the parameter circle
+ * \param CirclePlaneNormal normal vector to plane of the parameter circle
+ * \param CircleRadius radius of the parameter circle
+ * \param NewSphereCenter new center of a circumcircle
+ * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
+ * \param NormalVector normal vector
+ * \param SearchDirection search direction to make angle unique on return.
+ * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
+ */
+double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
+  Vector helper;
+  double radius, alpha;
+
+  helper.CopyVector(&NewSphereCenter);
+  // test whether new center is on the parameter circle's plane
+  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
+  radius = helper.ScalarProduct(&helper);
+  // test whether the new center vector has length of CircleRadius
+  if (fabs(radius - CircleRadius) > HULLEPSILON)
+    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+  alpha = helper.Angle(&OldSphereCenter);
+  // make the angle unique by checking the halfplanes/search direction
+  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+    alpha = 2.*M_PI - alpha;
+  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&OldSphereCenter);
+  helper.ProjectOntoPlane(&NormalVector);
+  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
+  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
+    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    return alpha;
+  } else {
+    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    return 2.*M_PI;
+  }
+};
+
+
+/** Checks whether the triangle consisting of the three atoms is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
+ *                 triangles exist which is the maximum for three points
+ */
+int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
+// TODO: use findTriangles here and return result.size();
+  LineMap::iterator FindLine;
+  PointMap::iterator FindPoint;
+  TriangleMap::iterator FindTriangle;
+  int adjacentTriangleCount = 0;
+  class BoundaryPointSet *Points[3];
+
+  //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+  // builds a triangle point set (Points) of the end points
+  for (int i = 0; i < 3; i++) {
+    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      Points[i] = FindPoint->second;
+    } else {
+      Points[i] = NULL;
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (Points[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (Points[j] != NULL) {
+          FindLine = Points[i]->lines.find(Points[j]->node->nr);
+          if (FindLine != Points[i]->lines.end()) {
+            for (; FindLine->first == Points[j]->node->nr; FindLine++) {
+              FindTriangle = FindLine->second->triangles.begin();
+              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
+                if ((
+                  (FindTriangle->second->endpoints[0] == Points[0])
+                    || (FindTriangle->second->endpoints[0] == Points[1])
+                    || (FindTriangle->second->endpoints[0] == Points[2])
+                  ) && (
+                    (FindTriangle->second->endpoints[1] == Points[0])
+                    || (FindTriangle->second->endpoints[1] == Points[1])
+                    || (FindTriangle->second->endpoints[1] == Points[2])
+                  ) && (
+                    (FindTriangle->second->endpoints[2] == Points[0])
+                    || (FindTriangle->second->endpoints[2] == Points[1])
+                    || (FindTriangle->second->endpoints[2] == Points[2])
+                  )
+                ) {
+                  adjacentTriangleCount++;
+                }
+              }
             }
-
-            // ... and put at new position
-            if (DoRandomRotation)
-              Runner->x.MatrixMultiplication(Rotations);
-            Runner->x.AddVector(&AtomTranslations);
-            Runner->x.AddVector(&FillerTranslations);
-            Runner->x.AddVector(&CurrentPosition);
-            // insert into Filling
-            Filling->AddAtom(Runner);
+            // Only one of the triangle lines must be considered for the triangle count.
+            *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+            return adjacentTriangleCount;
+
           }
-
-          // go through all bonds and add as well
-          Binder = filler->first;
-          while(Binder != filler->last) {
-            Binder = Binder->next;
+        }
+      }
+    }
+  }
+
+  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+  return adjacentTriangleCount;
+};
+
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * Note that this function is for the starting triangle.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
+ * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
+ * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
+ * @param BaseLine BoundaryLineSet with the current base line
+ * @param ThirdNode third atom to avoid in search
+ * @param candidates list of equally good candidates to return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_third_point_for_Tesselation(
+  Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
+  class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
+  double *ShortestAngle, const double RADIUS, LinkedCell *LC
+) {
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NewSphereCenter;        // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  Vector NewNormalVector;   // normal vector of the Candidate's triangle
+  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  LinkedAtoms *List = NULL;
+  double CircleRadius; // radius of this circle
+  double radius;
+  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  atom *Candidate = NULL;
+  CandidateForTesselation *optCandidate = NULL;
+
+  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
+
+  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+
+  // construct center of circle
+  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+
+  // calculate squared radius atom *ThirdNode,f circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    // test whether old center is on the band's plane
+    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(radius - CircleRadius) < HULLEPSILON) {
+
+      // check SearchDirection
+      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+      }
+
+      // get cell for the starting atom
+      if (LC->SetIndexToVector(&CircleCenter)) {
+        for(int i=0;i<NDIM;i++) // store indices of this cell
+        N[i] = LC->n[i];
+        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+      } else {
+        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+        return;
+      }
+      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+      //cout << Verbose(2) << "LC Intervals:";
+      for (int i=0;i<NDIM;i++) {
+        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+      }
+      //cout << 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]++) {
+            List = LC->GetCurrentCell();
+            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+            if (List != NULL) {
+              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                Candidate = (*Runner);
+
+                // check for three unique points
+                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
+
+                  // construct both new centers
+                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+
+                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
+                        && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  ) {
+                    helper.CopyVector(&NewNormalVector);
+                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+                    if (radius < RADIUS*RADIUS) {
+                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      NewSphereCenter.AddVector(&helper);
+                      NewSphereCenter.SubtractVector(&CircleCenter);
+                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+
+                      // OtherNewSphereCenter is created by the same vector just in the other direction
+                      helper.Scale(-1.);
+                      OtherNewSphereCenter.AddVector(&helper);
+                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
+                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+
+                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      alpha = min(alpha, Otheralpha);
+                      // if there is a better candidate, drop the current list and add the new candidate
+                      // otherwise ignore the new candidate and keep the list
+                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
+                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
+                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
+                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
+                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                        } else {
+                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
+                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
+                        }
+                        // if there is an equal candidate, add it to the list without clearing the list
+                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
+                          candidates->push_back(optCandidate);
+                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
+                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        } else {
+                          // remove all candidates from the list and then the list itself
+                          class CandidateForTesselation *remover = NULL;
+                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
+                            remover = *it;
+                            delete(remover);
+                          }
+                          candidates->clear();
+                          candidates->push_back(optCandidate);
+                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
+                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        }
+                        *ShortestAngle = alpha;
+                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
+                      } else {
+                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
+                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+                        } else {
+                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                        }
+                      }
+
+                    } else {
+                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                    }
+                  } else {
+                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                  }
+                } else {
+                  if (ThirdNode != NULL) {
+                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                  } else {
+                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                  }
+                }
+              }
+            }
+          }
+    } else {
+      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+    }
+  } else {
+    if (ThirdNode != NULL)
+      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+    else
+      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+  }
+
+  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
+  if (candidates->size() > 1) {
+          candidates->unique();
+          candidates->sort(sortCandidates);
+  }
+
+  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
+};
+
+struct Intersection {
+        Vector x1;
+        Vector x2;
+        Vector x3;
+        Vector x4;
+};
+
+/**
+ * Intersection calculation function.
+ *
+ * @param x to find the result for
+ * @param function parameter
+ */
+double MinIntersectDistance(const gsl_vector * x, void *params) {
+        double retval = 0;
+        struct Intersection *I = (struct Intersection *)params;
+        Vector intersection;
+        Vector SideA,SideB,HeightA, HeightB;
+        for (int i=0;i<NDIM;i++)
+                intersection.x[i] = gsl_vector_get(x, i);
+
+        SideA.CopyVector(&(I->x1));
+        SideA.SubtractVector(&I->x2);
+        HeightA.CopyVector(&intersection);
+        HeightA.SubtractVector(&I->x1);
+        HeightA.ProjectOntoPlane(&SideA);
+
+        SideB.CopyVector(&I->x3);
+        SideB.SubtractVector(&I->x4);
+        HeightB.CopyVector(&intersection);
+        HeightB.SubtractVector(&I->x3);
+        HeightB.ProjectOntoPlane(&SideB);
+
+        retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
+        //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
+
+        return retval;
+};
+
+
+/**
+ * Calculates whether there is an intersection between two lines. The first line
+ * always goes through point 1 and point 2 and the second line is given by the
+ * connection between point 4 and point 5.
+ *
+ * @param point 1 of line 1
+ * @param point 2 of line 1
+ * @param point 1 of line 2
+ * @param point 2 of line 2
+ *
+ * @return true if there is an intersection between the given lines, false otherwise
+ */
+bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
+        bool result;
+
+        struct Intersection par;
+                par.x1.CopyVector(&point1);
+                par.x2.CopyVector(&point2);
+                par.x3.CopyVector(&point3);
+                par.x4.CopyVector(&point4);
+
+    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
+    gsl_multimin_fminimizer *s = NULL;
+    gsl_vector *ss, *x;
+    gsl_multimin_function minex_func;
+
+    size_t iter = 0;
+    int status;
+    double size;
+
+    /* Starting point */
+    x = gsl_vector_alloc(NDIM);
+    gsl_vector_set(x, 0, point1.x[0]);
+    gsl_vector_set(x, 1, point1.x[1]);
+    gsl_vector_set(x, 2, point1.x[2]);
+
+    /* Set initial step sizes to 1 */
+    ss = gsl_vector_alloc(NDIM);
+    gsl_vector_set_all(ss, 1.0);
+
+    /* Initialize method and iterate */
+    minex_func.n = NDIM;
+    minex_func.f = &MinIntersectDistance;
+    minex_func.params = (void *)&par;
+
+    s = gsl_multimin_fminimizer_alloc(T, NDIM);
+    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
+
+    do {
+        iter++;
+        status = gsl_multimin_fminimizer_iterate(s);
+
+        if (status) {
+          break;
+        }
+
+        size = gsl_multimin_fminimizer_size(s);
+        status = gsl_multimin_test_size(size, 1e-2);
+
+        if (status == GSL_SUCCESS) {
+                cout << Verbose(2) << "converged to minimum" <<  endl;
+        }
+    } while (status == GSL_CONTINUE && iter < 100);
+
+    // check whether intersection is in between or not
+        Vector intersection, SideA, SideB, HeightA, HeightB;
+        double t1, t2;
+        for (int i = 0; i < NDIM; i++) {
+                intersection.x[i] = gsl_vector_get(s->x, i);
+        }
+
+        SideA.CopyVector(&par.x2);
+        SideA.SubtractVector(&par.x1);
+        HeightA.CopyVector(&intersection);
+        HeightA.SubtractVector(&par.x1);
+
+        t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
+
+        SideB.CopyVector(&par.x4);
+        SideB.SubtractVector(&par.x3);
+        HeightB.CopyVector(&intersection);
+        HeightB.SubtractVector(&par.x3);
+
+        t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
+
+        cout << Verbose(2) << "Intersection " << intersection << " is at "
+                << t1 << " for (" << point1 << "," << point2 << ") and at "
+                << t2 << " for (" << point3 << "," << point4 << "): ";
+
+        if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
+                cout << "true intersection." << endl;
+                result = true;
+        } else {
+                cout << "intersection out of region of interest." << endl;
+                result = false;
+        }
+
+        // free minimizer stuff
+    gsl_vector_free(x);
+    gsl_vector_free(ss);
+    gsl_multimin_fminimizer_free(s);
+
+        return result;
+}
+
+/** Finds the second point of starting triangle.
+ * \param *a first atom
+ * \param *Candidate pointer to candidate atom on return
+ * \param Oben vector indicating the outside
+ * \param Opt_Candidate reference to recommended candidate on return
+ * \param Storage[3] array storing angles and other candidate information
+ * \param RADIUS radius of virtual sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
+  Vector AngleCheck;
+  double norm = -1., angle;
+  LinkedAtoms *List = NULL;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LC->SetIndexToAtom(*a)) {  // get cell for the starting atom
+    for(int i=0;i<NDIM;i++) // store indices of this cell
+      N[i] = LC->n[i];
+  } else {
+    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+    return;
+  }
+  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+  cout << Verbose(3) << "LC Intervals from [";
+  for (int i=0;i<NDIM;i++) {
+  cout << " " << N[i] << "<->" << LC->N[i];
+  }
+  cout << "] :";
+  for (int i=0;i<NDIM;i++) {
+    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+  }
+  cout << 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]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            Candidate = (*Runner);
+            // check if we only have one unique point yet ...
+            if (a != Candidate) {
+              // Calculate center of the circle with radius RADIUS through points a and Candidate
+              Vector OrthogonalizedOben, a_Candidate, Center;
+              double distance, scaleFactor;
+
+              OrthogonalizedOben.CopyVector(&Oben);
+              a_Candidate.CopyVector(&(a->x));
+              a_Candidate.SubtractVector(&(Candidate->x));
+              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
+              OrthogonalizedOben.Normalize();
+              distance = 0.5 * a_Candidate.Norm();
+              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
+              OrthogonalizedOben.Scale(scaleFactor);
+
+              Center.CopyVector(&(Candidate->x));
+              Center.AddVector(&(a->x));
+              Center.Scale(0.5);
+              Center.AddVector(&OrthogonalizedOben);
+
+              AngleCheck.CopyVector(&Center);
+              AngleCheck.SubtractVector(&(a->x));
+              norm = a_Candidate.Norm();
+              // second point shall have smallest angle with respect to Oben vector
+              if (norm < RADIUS*2.) {
+                angle = AngleCheck.Angle(&Oben);
+                if (angle < Storage[0]) {
+                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
+                  Opt_Candidate = Candidate;
+                  Storage[0] = angle;
+                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+                } else {
+                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+                }
+              } else {
+                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
+              }
+            } else {
+              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
+            }
           }
         } else {
-          // leave empty
-          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
+          cout << Verbose(3) << "Linked cell list is empty." << endl;
         }
       }
-  return Filling;
+  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
 };
+
+/** Finds the starting triangle for find_non_convex_border().
+ * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
+ * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
+ * point are called.
+ * \param RADIUS radius of virtual rolling sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
+{
+  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  int i = 0;
+  LinkedAtoms *List = NULL;
+  atom* FirstPoint = NULL;
+  atom* SecondPoint = NULL;
+  atom* MaxAtom[NDIM];
+  double max_coordinate[NDIM];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector SearchDirection;
+
+  Oben.Zero();
+
+  for (i = 0; i < 3; i++) {
+    MaxAtom[i] = NULL;
+    max_coordinate[i] = -1;
+  }
+
+  // 1. searching topmost atom with respect to each axis
+  for (int i=0;i<NDIM;i++) { // each axis
+    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+            if ((*Runner)->x.x[i] > max_coordinate[i]) {
+              cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
+              max_coordinate[i] = (*Runner)->x.x[i];
+              MaxAtom[i] = (*Runner);
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  }
+
+  cout << Verbose(2) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    cout << i << ": " << *MaxAtom[i] << "\t";
+  cout << endl;
+
+  BTS = NULL;
+  CandidateList *Opt_Candidates = new CandidateList();
+  for (int k=0;k<NDIM;k++) {
+    Oben.x[k] = 1.;
+    FirstPoint = MaxAtom[k];
+    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
+
+    double ShortestAngle;
+    atom* Opt_Candidate = NULL;
+    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+
+    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    SecondPoint = Opt_Candidate;
+    if (SecondPoint == NULL)  // have we found a second point?
+      continue;
+    else
+      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+
+    helper.CopyVector(&(FirstPoint->x));
+    helper.SubtractVector(&(SecondPoint->x));
+    helper.Normalize();
+    Oben.ProjectOntoPlane(&helper);
+    Oben.Normalize();
+    helper.VectorProduct(&Oben);
+    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+
+    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+    Chord.SubtractVector(&(SecondPoint->x));
+    double radius = Chord.ScalarProduct(&Chord);
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+    helper.CopyVector(&Oben);
+    helper.Scale(CircleRadius);
+    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+
+    // look in one direction of baseline for initial candidate
+    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+
+    // adding point 1 and point 2 and the line between them
+    AddTrianglePoint(FirstPoint, 0);
+    AddTrianglePoint(SecondPoint, 1);
+    AddTriangleLine(TPS[0], TPS[1], 0);
+
+    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
+    Find_third_point_for_Tesselation(
+      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
+    );
+    cout << Verbose(1) << "List of third Points is ";
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        cout << " " << *(*it)->point;
+    }
+    cout << endl;
+
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+      // add third triangle point
+      AddTrianglePoint((*it)->point, 2);
+      // add the second and third line
+      AddTriangleLine(TPS[1], TPS[2], 1);
+      AddTriangleLine(TPS[0], TPS[2], 2);
+      // ... and triangles to the Maps of the Tesselation class
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      AddTriangle();
+      // ... and calculate its normal vector (with correct orientation)
+      (*it)->OptCenter.Scale(-1.);
+      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
+      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
+      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
+      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
+
+      // if we do not reach the end with the next step of iteration, we need to setup a new first line
+      if (it != Opt_Candidates->end()--) {
+        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
+        SecondPoint = (*it)->point;
+        // adding point 1 and point 2 and the line between them
+        AddTrianglePoint(FirstPoint, 0);
+        AddTrianglePoint(SecondPoint, 1);
+        AddTriangleLine(TPS[0], TPS[1], 0);
+      }
+      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+    }
+    if (BTS != NULL) // we have created one starting triangle
+      break;
+    else {
+      // remove all candidates from the list and then the list itself
+      class CandidateForTesselation *remover = NULL;
+      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        remover = *it;
+        delete(remover);
+      }
+      Opt_Candidates->clear();
+    }
+  }
+
+  // remove all candidates from the list and then the list itself
+  class CandidateForTesselation *remover = NULL;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    remover = *it;
+    delete(remover);
+  }
+  delete(Opt_Candidates);
+  cout << Verbose(1) << "End of Find_starting_triangle\n";
+};
+
+/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
+ * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
+ * make it bigger (i.e. closing one (the baseline) and opening two new ones).
+ * \param TPS[3] nodes of the triangle
+ * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
+ */
+bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
+{
+  bool result = false;
+  int counter = 0;
+
+  // check all three points
+  for (int i=0;i<3;i++)
+    for (int j=i+1; j<3; j++) {
+      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
+        LineMap::iterator FindLine;
+        pair<LineMap::iterator,LineMap::iterator> FindPair;
+        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
+        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+          // If there is a line with less than two attached triangles, we don't need a new line.
+          if (FindLine->second->TrianglesCount < 2) {
+            counter++;
+            break;  // increase counter only once per edge
+          }
+        }
+      } else { // no line
+        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
+        result = true;
+      }
+    }
+  if (counter > 1) {
+    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
+    result = true;
+  }
+  return result;
+};
+
+
+/** This function finds a triangle to a line, adjacent to an existing one.
+ * @param out output stream for debugging
+ * @param *mol molecule with Atom's and Bond's
+ * @param Line current baseline to search from
+ * @param T current triangle which \a Line is edge of
+ * @param RADIUS radius of the rolling ball
+ * @param N number of found triangles
+ * @param *filename filename base for intermediate envelopes
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+bool Tesselation::Find_next_suitable_triangle(ofstream *out,
+    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
+  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  bool result = true;
+  CandidateList *Opt_Candidates = new CandidateList();
+
+  Vector CircleCenter;
+  Vector CirclePlaneNormal;
+  Vector OldSphereCenter;
+  Vector SearchDirection;
+  Vector helper;
+  atom *ThirdNode = NULL;
+  LineMap::iterator testline;
+  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+  double radius, CircleRadius;
+
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+  for (int i=0;i<3;i++)
+    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+      ThirdNode = T.endpoints[i]->node;
+
+  // construct center of circle
+  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
+  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
+
+  // calculate squared radius of circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    // construct old center
+    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
+    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+    helper.Scale(sqrt(RADIUS*RADIUS - radius));
+    OldSphereCenter.AddVector(&helper);
+    OldSphereCenter.SubtractVector(&CircleCenter);
+    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+
+    // construct SearchDirection
+    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
+    helper.CopyVector(&Line.endpoints[0]->node->x);
+    helper.SubtractVector(&ThirdNode->x);
+    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+      SearchDirection.Scale(-1.);
+    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+    SearchDirection.Normalize();
+    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+      // rotated the wrong way!
+      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+    }
+
+    // add third point
+    Find_third_point_for_Tesselation(
+      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
+      &ShortestAngle, RADIUS, LC
+    );
+
+  } else {
+    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+  }
+
+  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << "Third Points are ";
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+          cout << " " << *(*it)->point;
+  }
+  cout << endl;
+
+  BoundaryLineSet *BaseRay = &Line;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+          cout << Verbose(1) << " Third point candidate is " << *(*it)->point
+                << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
+          cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
+
+          // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+          atom *AtomCandidates[3];
+          AtomCandidates[0] = (*it)->point;
+          AtomCandidates[1] = BaseRay->endpoints[0]->node;
+          AtomCandidates[2] = BaseRay->endpoints[1]->node;
+          int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
+
+          BTS = NULL;
+          // If there is no triangle, add it regularly.
+          if (existentTrianglesCount == 0) {
+      AddTrianglePoint((*it)->point, 0);
+      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
+      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+
+      AddTriangleLine(TPS[0], TPS[1], 0);
+      AddTriangleLine(TPS[0], TPS[2], 1);
+      AddTriangleLine(TPS[1], TPS[2], 2);
+
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      AddTriangle();
+      (*it)->OptCenter.Scale(-1.);
+      BTS->GetNormalVector((*it)->OptCenter);
+      (*it)->OptCenter.Scale(-1.);
+
+      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
+        << " for this triangle ... " << endl;
+      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
+          } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
+        AddTrianglePoint((*it)->point, 0);
+        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
+        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+
+        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
+        // i.e. at least one of the three lines must be present with TriangleCount <= 1
+        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
+          AddTriangleLine(TPS[0], TPS[1], 0);
+          AddTriangleLine(TPS[0], TPS[2], 1);
+          AddTriangleLine(TPS[1], TPS[2], 2);
+
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          AddTriangle();
+
+          (*it)->OtherOptCenter.Scale(-1.);
+          BTS->GetNormalVector((*it)->OtherOptCenter);
+          (*it)->OtherOptCenter.Scale(-1.);
+
+          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
+          << " for this triangle ... " << endl;
+          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
+        } else {
+          cout << Verbose(1) << "WARNING: This triangle consisting of ";
+          cout << *(*it)->point << ", ";
+          cout << *BaseRay->endpoints[0]->node << " and ";
+          cout << *BaseRay->endpoints[1]->node << " ";
+          cout << "exists and is not added, as it does not seem helpful!" << endl;
+          result = false;
+        }
+    } else {
+      cout << Verbose(1) << "This triangle consisting of ";
+      cout << *(*it)->point << ", ";
+      cout << *BaseRay->endpoints[0]->node << " and ";
+      cout << *BaseRay->endpoints[1]->node << " ";
+      cout << "is invalid!" << endl;
+      result = false;
+    }
+
+          if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
+      if (DoTecplotOutput) {
+        string NameofTempFile(tempbasename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(TecplotSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+
+      if (DoRaster3DOutput) {
+        string NameofTempFile(tempbasename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(Raster3DSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_raster3d_file(out, tempstream, this, mol);
+        // include the current position of the virtual sphere in the temporary raster3d file
+        // make the circumsphere's center absolute again
+        helper.CopyVector(&BaseRay->endpoints[0]->node->x);
+        helper.AddVector(&BaseRay->endpoints[1]->node->x);
+        helper.Scale(0.5);
+        (*it)->OptCenter.AddVector(&helper);
+        Vector *center = mol->DetermineCenterOfAll(out);
+        (*it)->OptCenter.AddVector(center);
+        delete(center);
+        // and add to file plus translucency object
+        *tempstream << "# current virtual sphere\n";
+        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+          << "\t" << RADIUS << "\t1 0 0\n";
+        *tempstream << "9\n  terminating special property\n";
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+      if (DoTecplotOutput || DoRaster3DOutput)
+        TriangleFilesWritten++;
+          }
+
+    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
+          BaseRay = BLS[0];
+  }
+
+  // remove all candidates from the list and then the list itself
+  class CandidateForTesselation *remover = NULL;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    remover = *it;
+    delete(remover);
+  }
+  delete(Opt_Candidates);
+  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
+  return result;
+};
+
+/**
+ * Sort function for the candidate list.
+ */
+bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
+  // TODO: use get angle and remove duplicate code
+  Vector BaseLineVector, OrthogonalVector, helper;
+        if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
+          cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
+          //return false;
+          exit(1);
+        }
+        // create baseline vector
+        BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
+        BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
+        BaseLineVector.Normalize();
+
+  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
+        helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
+        helper.SubtractVector(&(candidate1->point->x));
+        OrthogonalVector.CopyVector(&helper);
+        helper.VectorProduct(&BaseLineVector);
+        OrthogonalVector.SubtractVector(&helper);
+        OrthogonalVector.Normalize();
+
+  // calculate both angles and correct with in-plane vector
+        helper.CopyVector(&(candidate1->point->x));
+        helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
+        double phi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  helper.CopyVector(&(candidate2->point->x));
+        helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
+        double psi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+                psi = 2.*M_PI - psi;
+        }
+
+        cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
+        cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
+
+        // return comparison
+        return phi < psi;
+}
+
+/**
+ * Checks whether the provided atom is within the tesselation structure.
+ *
+ * @param atom of which to check the position
+ * @param tesselation structure
+ *
+ * @return true if the atom is inside the tesselation structure, false otherwise
+ */
+bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC) {
+  if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
+    cout << Verbose(0) << "Error: There is no tesselation structure to compare the atom with, "
+        << "please create one first.";
+    exit(1);
+  }
+
+  class atom *trianglePoints[3];
+  trianglePoints[0] = findClosestAtom(Atom, LC);
+  // check whether closest atom is "too close" :), then it's inside
+  if (trianglePoints[0]->x.DistanceSquared(&Atom->x) < MYEPSILON)
+    return true;
+  list<atom*> *connectedClosestAtoms = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
+  trianglePoints[1] = connectedClosestAtoms->front();
+  trianglePoints[2] = connectedClosestAtoms->back();
+  for (int i=0;i<3;i++) {
+    if (trianglePoints[i] == NULL) {
+      cout << Verbose(1) << "IsInnerAtom encounters serious error, point " << i << " not found." << endl;
+    }
+
+    cout << Verbose(1) << "List of possible atoms:" << endl;
+    cout << *trianglePoints[i] << endl;
+
+//    for(list<atom*>::iterator runner = connectedClosestAtoms->begin(); runner != connectedClosestAtoms->end(); runner++)
+//      cout << Verbose(2) << *(*runner) << endl;
+  }
+  delete(connectedClosestAtoms);
+
+  list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);
+
+  if (triangles->empty()) {
+    cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
+    exit(1);
+  }
+
+  Vector helper;
+  helper.CopyVector(&Atom->x);
+
+  // Only in case of degeneration, there will be two different scalar products.
+  // If at least one scalar product is positive, the point is considered to be outside.
+  bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
+      && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
+  delete(triangles);
+  return status;
+}
+
+/**
+ * Finds the atom which is closest to the provided one.
+ *
+ * @param atom to which to find the closest other atom
+ * @param linked cell structure
+ *
+ * @return atom which is closest to the provided one
+ */
+atom* findClosestAtom(const atom* Atom, LinkedCell* LC) {
+  LinkedAtoms *List = NULL;
+  atom* closestAtom = NULL;
+  double distance = 1e16;
+  Vector helper;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  LC->SetIndexToVector(&Atom->x); // 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];
+  //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //cout << 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]++) {
+        List = LC->GetCurrentCell();
+        //cout << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            helper.CopyVector(&Atom->x);
+            helper.SubtractVector(&(*Runner)->x);
+            double currentNorm = helper. Norm();
+            if (currentNorm < distance) {
+              distance = currentNorm;
+              closestAtom = (*Runner);
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
+            << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  return closestAtom;
+}
+
+/**
+ * Gets all atoms connected to the provided atom by triangulation lines.
+ *
+ * @param atom of which get all connected atoms
+ * @param atom to be checked whether it is an inner atom
+ *
+ * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
+ */
+list<atom*> * Tesselation::getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck) {
+  list<atom*> connectedAtoms;
+  map<double, atom*> anglesOfAtoms;
+  map<double, atom*>::iterator runner;
+  LineMap::iterator findLines = LinesOnBoundary.begin();
+  list<atom*>::iterator listRunner;
+  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
+  atom* current;
+  bool takeAtom = false;
+
+  planeNorm.CopyVector(&Atom->x);
+  planeNorm.SubtractVector(&AtomToCheck->x);
+  planeNorm.Normalize();
+
+  while (findLines != LinesOnBoundary.end()) {
+    takeAtom = false;
+
+    if (findLines->second->endpoints[0]->Nr == Atom->nr) {
+      takeAtom = true;
+      current = findLines->second->endpoints[1]->node;
+    } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
+      takeAtom = true;
+      current = findLines->second->endpoints[0]->node;
+    }
+
+    if (takeAtom) {
+      connectedAtoms.push_back(current);
+      currentPoint.CopyVector(&current->x);
+      currentPoint.ProjectOntoPlane(&planeNorm);
+      center.AddVector(&currentPoint);
+    }
+
+    findLines++;
+  }
+
+  cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
+    << "; scale factor " << 1.0/connectedAtoms.size();
+
+  center.Scale(1.0/connectedAtoms.size());
+  listRunner = connectedAtoms.begin();
+
+  cout << " calculated center " << center <<  endl;
+
+  // construct one orthogonal vector
+  helper.CopyVector(&AtomToCheck->x);
+  helper.ProjectOntoPlane(&planeNorm);
+  OrthogonalVector.MakeNormalVector(&center, &helper, &(*listRunner)->x);
+  while (listRunner != connectedAtoms.end()) {
+    double angle = getAngle((*listRunner)->x, AtomToCheck->x, center, OrthogonalVector);
+    cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
+    anglesOfAtoms.insert(pair<double, atom*>(angle, (*listRunner)));
+    listRunner++;
+  }
+
+  list<atom*> *result = new list<atom*>;
+  runner = anglesOfAtoms.begin();
+  cout << "First value is " << *runner->second << endl;
+  result->push_back(runner->second);
+  runner = anglesOfAtoms.end();
+  runner--;
+  cout << "Second value is " << *runner->second << endl;
+  result->push_back(runner->second);
+
+  cout << "List of closest atoms has " << result->size() << " elements, which are "
+    << *(result->front()) << " and " << *(result->back()) << endl;
+
+  return result;
+}
+
+/**
+ * Finds triangles belonging to the three provided atoms.
+ *
+ * @param atom list, is expected to contain three atoms
+ *
+ * @return triangles which belong to the provided atoms, will be empty if there are none,
+ *         will usually be one, in case of degeneration, there will be two
+ */
+list<BoundaryTriangleSet*> *Tesselation::FindTriangles(atom* Points[3]) {
+  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+  LineMap::iterator FindLine;
+  PointMap::iterator FindPoint;
+  TriangleMap::iterator FindTriangle;
+  class BoundaryPointSet *TrianglePoints[3];
+
+  for (int i = 0; i < 3; i++) {
+    FindPoint = PointsOnBoundary.find(Points[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      TrianglePoints[i] = FindPoint->second;
+    } else {
+      TrianglePoints[i] = NULL;
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (TrianglePoints[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (TrianglePoints[j] != NULL) {
+          FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
+          if (FindLine != TrianglePoints[i]->lines.end()) {
+            for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
+              FindTriangle = FindLine->second->triangles.begin();
+              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
+                if ((
+                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
+                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
+                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
+                  ) && (
+                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
+                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
+                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
+                  ) && (
+                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
+                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
+                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
+                  )
+                ) {
+                  result->push_back(FindTriangle->second);
+                }
+              }
+            }
+            // Is it sufficient to consider one of the triangle lines for this.
+            return result;
+
+          }
+        }
+      }
+    }
+  }
+
+  return result;
+}
+
+/**
+ * Gets the angle between a point and a reference relative to the provided center.
+ *
+ * @param point to calculate the angle for
+ * @param reference to which to calculate the angle
+ * @param center for which to calculate the angle between the vectors
+ * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
+ *
+ * @return angle between point and reference
+ */
+double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) {
+  Vector ReferenceVector, helper;
+  cout << Verbose(2) << reference << " is our reference " << endl;
+  cout << Verbose(2) << center << " is our center " << endl;
+  // create baseline vector
+  ReferenceVector.CopyVector(&reference);
+  ReferenceVector.SubtractVector(&center);
+  if (ReferenceVector.IsNull())
+    return M_PI;
+
+  // calculate both angles and correct with in-plane vector
+  helper.CopyVector(&point);
+  helper.SubtractVector(&center);
+  if (helper.IsNull())
+    return M_PI;
+  double phi = ReferenceVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+
+  cout << Verbose(2) << point << " has angle " << phi << endl;
+
+  return phi;
+}
+
+/**
+ * Checks whether the provided point is within the tesselation structure.
+ *
+ * This is a wrapper function for IsInnerAtom, so it can be used with a simple
+ * vector, too.
+ *
+ * @param point of which to check the position
+ * @param tesselation structure
+ *
+ * @return true if the point is inside the tesselation structure, false otherwise
+ */
+bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) {
+  atom *temp_atom = new atom;
+  bool status = false;
+  temp_atom->x.CopyVector(&Point);
+
+  status = IsInnerAtom(temp_atom, Tess, LC);
+  delete(temp_atom);
+
+  return status;
+}
 
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
@@ -818 +3090 @@
  * \param *mol molecule structure with Atom's and Bond's
  * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
- * \param *LCList atoms in LinkedCell list
  * \param *filename filename prefix for output of vertex data
  * \para RADIUS radius of the virtual sphere
@@ -827 +3098 @@
   bool freeTess = false;
   bool freeLC = false;
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  int TriangleFilesWritten = 0;
-
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
   if (Tess == NULL) {
@@ -848 +3115 @@
   }
 
-  Tess->Find_starting_triangle(out, RADIUS, LCList);
+  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
 
   baseline = Tess->LinesOnBoundary.begin();
   while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->TrianglesCount == 1) {
-      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
+      failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
@@ -869 +3136 @@
       flag = false;
     }
-
-    // write temporary envelope
-    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-
-      if (DoRaster3DOutput) {
-        string NameofTempFile(filename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, Tess, mol);
-//        // include the current position of the virtual sphere in the temporary raster3d file
-//        // make the circumsphere's center absolute again
-//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
-//        helper.AddVector(BaseRay->endpoints[1]->node->node);
-//        helper.Scale(0.5);
-//        (*it)->OptCenter.AddVector(&helper);
-//        Vector *center = mol->DetermineCenterOfAll(out);
-//        (*it)->OptCenter.SubtractVector(center);
-//        delete(center);
-//        // and add to file plus translucency object
-//        *tempstream << "# current virtual sphere\n";
-//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-//          << "\t" << RADIUS << "\t1 0 0\n";
-//        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
-    }
-  }
-
-  // write final envelope
-  if (filename != 0) {
+  }
+  if (1) { //failflag) {
     *out << Verbose(1) << "Writing final tecplot file\n";
     if (DoTecplotOutput) {
@@ -956 +3170 @@
     cout << Verbose(2) << "None." << endl;
 
+  // Tests the IsInnerAtom() function.
+  Vector x (0, 0, 0);
+  cout << Verbose(0) << "Point to check: " << x << endl;
+  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
+    << "for vector " << x << "." << endl;
+  atom* a = Tess->PointsOnBoundary.begin()->second->node;
+  cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
+  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
+    << "for atom " << a << " (on boundary)." << endl;
+  LinkedAtoms *List = NULL;
+  for (int i=0;i<NDIM;i++) { // each axis
+    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
+    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
+      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
+        List = LCList->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+            if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
+              a = *Runner;
+              i=3;
+              for (int j=0;j<NDIM;j++)
+                LCList->n[j] = LCList->N[j];
+              break;
+            }
+          }
+        }
+      }
+  }
+  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
+    << "for atom " << a << " (inside)." << endl;
+
+
   if (freeTess)
     delete(Tess);