Changes in src/boundary.cpp [cc2ee5:ed060e]

File:

• src/boundary.cpp (modified) (21 diffs, 1 prop)

src/boundary.cpp

@@ -2 +2 @@
 #include "boundary.hpp"
 
-#define DEBUG 1
-#define DoTecplotOutput 0
-#define DoRaster3DOutput 1
-#define TecplotSuffix ".dat"
-#define Raster3DSuffix ".r3d"
-
 // ======================================== Points on Boundary =================================
 
@@ -14 +8 @@
   LinesCount = 0;
   Nr = -1;
-}
-;
+};
 
 BoundaryPointSet::BoundaryPointSet(atom *Walker)
@@ -22 +15 @@
   LinesCount = 0;
   Nr = Walker->nr;
-}
-;
+};
 
 BoundaryPointSet::~BoundaryPointSet()
@@ -29 +21 @@
   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   node = NULL;
-  lines.clear();
-}
-;
-
-void
-BoundaryPointSet::AddLine(class BoundaryLineSet *line)
-{
-  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
-      << endl;
-  if (line->endpoints[0] == this)
-    {
-      lines.insert(LinePair(line->endpoints[1]->Nr, line));
-    }
-  else
-    {
-      lines.insert(LinePair(line->endpoints[0]->Nr, line));
-    }
+};
+
+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)
+};
+
+ostream & operator << (ostream &ost, BoundaryPointSet &a) 
 {
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
   return ost;
-}
-;
+};
 
 // ======================================== Lines on Boundary =================================
@@ -62 +44 @@
 BoundaryLineSet::BoundaryLineSet()
 {
-  for (int i = 0; i < 2; i++)
+  for (int i=0;i<2;i++)
     endpoints[i] = NULL;
   TrianglesCount = 0;
   Nr = -1;
-}
-;
+};
 
 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
@@ -76 +57 @@
   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); //
+  Point[0]->AddLine(this);
+  Point[1]->AddLine(this);
   // clear triangles list
   TrianglesCount = 0;
   cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
-}
-;
+};
 
 BoundaryLineSet::~BoundaryLineSet()
 {
-        for (int i = 0; i < 2; i++) {
-                cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-                endpoints[i]->lines.erase(Nr);
-                LineMap::iterator tester = endpoints[i]->lines.begin();
-                tester++;
-                if (tester == endpoints[i]->lines.end()) {
-                        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;
-        }
-}
-;
-
-void
-BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
-{
-  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
-      << endl;
-  triangles.insert(TrianglePair(TrianglesCount, triangle));
+  for (int i=0;i<2;i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 
+    endpoints[i]->lines.erase(Nr);
+    LineMap::iterator tester = endpoints[i]->lines.begin();
+    tester++;
+    if (tester == endpoints[i]->lines.end()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      delete(endpoints[i]);
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+};
+
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
+  triangles.insert ( TrianglePair( TrianglesCount, triangle) );
   TrianglesCount++;
-}
-;
-
-ostream &
-operator <<(ostream &ost, BoundaryLineSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
-      << a.endpoints[1]->node->Name << "]";
+};
+
+ostream & operator << (ostream &ost, BoundaryLineSet &a) 
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
   return ost;
-}
-;
+};
 
 // ======================================== Triangles on Boundary =================================
@@ -128 +97 @@
 BoundaryTriangleSet::BoundaryTriangleSet()
 {
-  for (int i = 0; i < 3; i++)
-    {
-      endpoints[i] = NULL;
-      lines[i] = NULL;
-    }
+  for (int i=0;i<3;i++) {
+    endpoints[i] = NULL;
+    lines[i] = NULL;
+  }
   Nr = -1;
-}
-;
-
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3],
-    int number)
+};
+
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
 {
   // set number
@@ -144 +110 @@
   // set lines
   cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
-  for (int i = 0; i < 3; i++)
-    {
-      lines[i] = line[i];
-      lines[i]->AddTriangle(this);
+  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 
     }
-  // 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);
-    }
+  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()) {
-                        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-                        if (lines[i] != NULL) {
-                                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 a triangle." << endl;
-        }
-}
-;
-
-void
-BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+  for (int i=0;i<3;i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    TriangleMap::iterator tester = lines[i]->triangles.begin();
+    tester++;
+    if (tester == lines[i]->triangles.end()) {
+      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+      delete(lines[i]);
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+  }
+};
+
+void BoundaryTriangleSet::GetNormalVector(Vector &NormalVector)
 {
   // get normal vector
-  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
-      &endpoints[2]->node->x);
-
+  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 (endpoints[0]->node->x.Projection(&OtherVector) > 0)
+  if (endpoints[0]->node->x.Projection(&NormalVector) > 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 << "]";
+};
+
+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;
-}
-;
+};
 
 // ========================================== F U N C T I O N S =================================
@@ -226 +174 @@
  * \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 * 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;
-          }
+  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.
@@ -261 +201 @@
  * \param *mol molecule structure representing the cluster
  */
-Boundaries *
-GetBoundaryPoints(ofstream *out, molecule *mol)
+Boundaries * GetBoundaryPoints(ofstream *out, molecule *mol)
 {
   atom *Walker = NULL;
@@ -268 +207 @@
   LineMap LinesOnBoundary;
   TriangleMap TrianglesOnBoundary;
-
+  
   *out << Verbose(1) << "Finding all boundary points." << endl;
-  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+  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)
+  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)
         {
-          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;
-                }
-            }
+          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;
+          }
         }
-      // 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);
-    }
+      }
+    } while (flag);
+  } 
   return BoundaryPoints;
-}
-;
+};
 
 /** Determines greatest diameters of a cluster defined by its convex envelope.
@@ -469 +375 @@
  * \param IsAngstroem whether we have angstroem or atomic units
  * \return NDIM array of the diameters
- */
-double *
-GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
-    bool IsAngstroem)
+ */ 
+double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
 {
   // get points on boundary of NULL was given as parameter
   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;
   double *GreatestDiameter = new double[NDIM];
-  for (int i = 0; i < NDIM; i++)
+  for(int i=0;i<NDIM;i++)
     GreatestDiameter[i] = 0.;
   double OldComponent, tmp, w1, w2;
   Vector DistanceVector, OtherVector;
   int component, Othercomponent;
-  for (int axis = 0; axis < NDIM; axis++)
-    { // regard each projected plane
-      //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
-      for (int j = 0; j < 2; j++)
-        { // and for both axis on the current plane
-          component = (axis + j + 1) % NDIM;
-          Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
-          //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
-          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
-              != BoundaryPoints[axis].end(); runner++)
-            {
-              //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
-              // seek for the neighbours pair where the Othercomponent sign flips
-              Neighbour = runner;
-              Neighbour++;
-              if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
-                Neighbour = BoundaryPoints[axis].begin();
-              DistanceVector.CopyVector(&runner->second.second->x);
-              DistanceVector.SubtractVector(&Neighbour->second.second->x);
-              do
-                { // seek for neighbour pair where it flips
-                  OldComponent = DistanceVector.x[Othercomponent];
-                  Neighbour++;
-                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
-                    Neighbour = BoundaryPoints[axis].begin();
-                  DistanceVector.CopyVector(&runner->second.second->x);
-                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
-                  //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
-                }
-              while ((runner != Neighbour) && (fabs(OldComponent / fabs(
-                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
-                  DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
-              if (runner != Neighbour)
-                {
-                  OtherNeighbour = Neighbour;
-                  if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
-                    OtherNeighbour = BoundaryPoints[axis].end();
-                  OtherNeighbour--;
-                  //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
-                  // now we have found the pair: Neighbour and OtherNeighbour
-                  OtherVector.CopyVector(&runner->second.second->x);
-                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
-                  //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
-                  //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
-                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
-                  w1 = fabs(OtherVector.x[Othercomponent]);
-                  w2 = fabs(DistanceVector.x[Othercomponent]);
-                  tmp = fabs((w1 * DistanceVector.x[component] + w2
-                      * OtherVector.x[component]) / (w1 + w2));
-                  // mark if it has greater diameter
-                  //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
-                  GreatestDiameter[component] = (GreatestDiameter[component]
-                      > tmp) ? GreatestDiameter[component] : tmp;
-                } //else
-              //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
-            }
-        }
+  for(int axis=0;axis<NDIM;axis++) { // regard each projected plane
+    //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
+    for (int j=0;j<2;j++) { // and for both axis on the current plane
+      component = (axis+j+1)%NDIM;
+      Othercomponent = (axis+1+((j+1) & 1))%NDIM;
+      //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
+        // seek for the neighbours pair where the Othercomponent sign flips
+        Neighbour = runner;
+        Neighbour++;
+        if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
+          Neighbour = BoundaryPoints[axis].begin();
+        DistanceVector.CopyVector(&runner->second.second->x);
+        DistanceVector.SubtractVector(&Neighbour->second.second->x);
+        do {  // seek for neighbour pair where it flips
+          OldComponent = DistanceVector.x[Othercomponent]; 
+          Neighbour++;
+          if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
+            Neighbour = BoundaryPoints[axis].begin();
+          DistanceVector.CopyVector(&runner->second.second->x);
+          DistanceVector.SubtractVector(&Neighbour->second.second->x);
+          //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+        } while ((runner != Neighbour) && ( fabs( OldComponent/fabs(OldComponent) - DistanceVector.x[Othercomponent]/fabs(DistanceVector.x[Othercomponent]) ) < MYEPSILON)); // as long as sign does not flip
+        if (runner != Neighbour) {
+          OtherNeighbour = Neighbour;
+          if (OtherNeighbour == BoundaryPoints[axis].begin())  // make it wrap around
+            OtherNeighbour = BoundaryPoints[axis].end();
+          OtherNeighbour--;
+          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 
+          // now we have found the pair: Neighbour and OtherNeighbour
+          OtherVector.CopyVector(&runner->second.second->x);
+          OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
+          //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
+          //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
+          // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
+          w1 = fabs(OtherVector.x[Othercomponent]);
+          w2 = fabs(DistanceVector.x[Othercomponent]);
+          tmp = fabs((w1*DistanceVector.x[component] + w2*OtherVector.x[component])/(w1+w2));
+          // mark if it has greater diameter
+          //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+          GreatestDiameter[component] = (GreatestDiameter[component] > tmp) ? GreatestDiameter[component] : tmp;
+        } //else
+          //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+      }
     }
-  *out << Verbose(0) << "RESULT: The biggest diameters are "
-      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
-      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
-      : "atomiclength") << "." << endl;
+  }
+  *out << Verbose(0) << "RESULT: The biggest diameters are " << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "." << endl;
 
   // free reference lists
   if (BoundaryFreeFlag)
-    delete[] (BoundaryPoints);
+    delete[](BoundaryPoints);
 
   return GreatestDiameter;
-}
-;
-
-/** Creates the objects in a raster3d file (renderable with a header.r3d)
- * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
- * \param *Tess Tesselation structure with constructed triangles
- * \param *mol molecule structure with atom positions
- */
-void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
-{
-        atom *Walker = mol->start;
-        bond *Binder = mol->first;
-        int i;
-        Vector *center = mol->DetermineCenterOfAll(out);
-        if (rasterfile != NULL) {
-                //cout << Verbose(1) << "Writing Raster3D file ... ";
-                *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
-                *rasterfile << "@header.r3d" << endl;
-                *rasterfile << "# All atoms as spheres" << endl;
-                while (Walker->next != mol->end) {
-                        Walker = Walker->next;
-                        *rasterfile << "2" << endl << "  ";     // 2 is sphere type
-                        for (i=0;i<NDIM;i++)
-                                *rasterfile << Walker->x.x[i]+center->x[i] << " ";
-                        *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
-                }
-
-                *rasterfile << "# All bonds as vertices" << endl;
-                while (Binder->next != mol->last) {
-                        Binder = Binder->next;
-                        *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 << "\t0.03\t";
-                        for (i=0;i<NDIM;i++)
-                                *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
-                        *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
-                }
-
-                *rasterfile << "# All tesselation triangles" << endl;
-                for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
-                        *rasterfile << "1" << endl << "  ";     // 1 is triangle type
-                        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->x.x[j]+center->x[j] << " ";
-                                *rasterfile << "\t";
-                        }
-                        *rasterfile << "1. 0. 0." << endl;      // red as colour
-                        *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
-                }
-        } else {
-                cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
-        }
-        delete(center);
-};
-
-/** This function creates the tecplot file, displaying the tesselation of the hull.
- * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
- * \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;
-      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 *tecplot output stream for tecplot 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
- * \return determined volume of the cluster in cubed config:GetIsAngstroem()
  */
-double
-VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration,
-    Boundaries *BoundaryPtr, molecule *mol)
+double VolumeOfConvexEnvelope(ofstream *out, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -690 +469 @@
   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.
+  double G,h;
+  Vector x,y;
+  double a,b,c;
 
   // 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;
-      Walker->x.Translate(CenterOfGravity);
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
+  
+  // 3. Find all points on the boundary
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+  
+  // 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);
     }
 
-  // 3. Find all points on the boundary
-  if (BoundaryPoints == NULL)
-    {
-      BoundaryFreeFlag = true;
-      BoundaryPoints = GetBoundaryPoints(out, mol);
-    }
-  else
-    {
-      *out << Verbose(1) << "Using given boundary points set." << endl;
-    }
-
-  // 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 << 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
 
+
   // listing for debugging
-  //  *out << Verbose(1) << "Listing PointsOnBoundary:";
-  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-  //    *out << " " << *runner->second;
-  //  }
-  //  *out << endl;
-
+//  *out << Verbose(1) << "Listing PointsOnBoundary:";
+//  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+//    *out << " " << *runner->second;
+//  }
+//  *out << endl;
+  
   // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
-
+  
   // 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;
+  *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
-  *out << Verbose(1)
-      << "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++)
-    { // go through every triangle, calculate volume of its pyramid with CoG as peak
-      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.Distance(
-          &runner->second->endpoints[1]->node->x));
-      b = sqrt(runner->second->endpoints[0]->node->x.Distance(
-          &runner->second->endpoints[2]->node->x));
-      c = sqrt(runner->second->endpoints[2]->node->x.Distance(
-          &runner->second->endpoints[1]->node->x));
-      G = sqrt(((a * a + b * b + c * c) * (a * a + b * b + c * c) - 2 * (a * a
-          * a * a + b * b * b * b + c * c * c * c)) / 16.); // area of tesselated triangle
-      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)
-      *out << Verbose(2) << "Area of triangle is " << G << " "
-          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
-          << h << " and the volume is " << PyramidVolume << " "
-          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-      volume += PyramidVolume;
-    }
-  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
-      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
-      << endl;
+  *out << Verbose(1) << "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++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
+    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.Distance(&runner->second->endpoints[1]->node->x));
+    b = sqrt(runner->second->endpoints[0]->node->x.Distance(&runner->second->endpoints[2]->node->x));
+    c = sqrt(runner->second->endpoints[2]->node->x.Distance(&runner->second->endpoints[1]->node->x));
+    G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
+    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)
+    *out << Verbose(2) << "Area of triangle is " << G << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+    volume += PyramidVolume; 
+  }
+  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+
 
   // 7. translate all points back from CoG
-  *out << Verbose(1) << "Translating system back from Center of Gravity."
-      << endl;
+  *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
-  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    Walker->x.Translate(CenterOfGravity);
+  }
 
   // free reference lists
   if (BoundaryFreeFlag)
-    delete[] (BoundaryPoints);
-
+    delete[](BoundaryPoints);
+  
   return volume;
-}
-;
+};
+
 
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
@@ -815 +564 @@
  * \param celldensity desired average density in final cell
  */
-void
-PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
-    double ClusterVolume, double celldensity)
+void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity)
 {
   // transform to PAS
   mol->PrincipalAxisSystem(out, true);
-
+  
   // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -827 +574 @@
   double clustervolume;
   if (ClusterVolume == 0)
-    clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
-        BoundaryPoints, mol);
-  else
+    clustervolume = VolumeOfConvexEnvelope(out, 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] =
-    { 1, 1, 1 };
+  int repetition[NDIM] = {1, 1, 1};
   int TotalNoClusters = 1;
-  for (int i = 0; i < NDIM; i++)
+  for (int i=0;i<NDIM;i++)
     TotalNoClusters *= repetition[i];
 
@@ -843 +587 @@
   double totalmass = 0.;
   atom *Walker = mol->start;
-  while (Walker->next != mol->end)
-    {
-      Walker = Walker->next;
-      totalmass += Walker->type->mass;
-    }
-  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
-      << totalmass << " atomicmassunit." << endl;
-
-  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
-      << totalmass / clustervolume << " atomicmassunit/"
-      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    totalmass += Walker->type->mass;
+  }
+  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
+  
+  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass/clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
   // solve cubic polynomial
-  *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
-      << endl;
+  *out << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
   double cellvolume;
   if (IsAngstroem)
-    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
-        / clustervolume)) / (celldensity - 1);
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
   else
-    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
-        / clustervolume)) / (celldensity - 1);
-  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
-      << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
-      : "atomiclength") << "^3." << endl;
-
-  double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
-      * GreatestDiameter[1] * GreatestDiameter[2]);
-  *out << Verbose(1)
-      << "Minimum volume of the convex envelope contained in a rectangular box is "
-      << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
-      : "atomiclength") << "^3." << endl;
-  if (minimumvolume > cellvolume)
-    {
-      cerr << Verbose(0)
-          << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
-          << endl;
-      cout << Verbose(0)
-          << "Setting Box dimensions to minimum possible, the greatest diameters."
-          << endl;
-      for (int i = 0; i < NDIM; i++)
-        BoxLengths.x[i] = GreatestDiameter[i];
-      mol->CenterEdge(out, &BoxLengths);
+    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
+  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  
+  double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
+  *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  if (minimumvolume > cellvolume) {
+    cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
+    cout << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
+    for(int i=0;i<NDIM;i++)
+      BoxLengths.x[i] = GreatestDiameter[i];
+    mol->CenterEdge(out, &BoxLengths);
+  } else {
+    BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
+    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
+              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
+              + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
+    BoxLengths.x[2] = minimumvolume - cellvolume;
+    double x0 = 0.,x1 = 0.,x2 = 0.;
+    if (gsl_poly_solve_cubic(BoxLengths.x[0],BoxLengths.x[1],BoxLengths.x[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
+      *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
+    else {
+      *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
+      x0 = x2;  // sorted in ascending order
     }
-  else
-    {
-      BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
-          * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
-      BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
-          * GreatestDiameter[1] + repetition[0] * repetition[2]
-          * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
-          * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
-      BoxLengths.x[2] = minimumvolume - cellvolume;
-      double x0 = 0., x1 = 0., x2 = 0.;
-      if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
-          BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
-        *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
-            << " ." << endl;
-      else
-        {
-          *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
-              << " and " << x1 << " and " << x2 << " ." << endl;
-          x0 = x2; // sorted in ascending order
-        }
-
-      cellvolume = 1;
-      for (int i = 0; i < NDIM; i++)
-        {
-          BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
-          cellvolume *= BoxLengths.x[i];
-        }
-
-      // set new box dimensions
-      *out << Verbose(0) << "Translating to box with these boundaries." << endl;
-      mol->CenterInBox((ofstream *) &cout, &BoxLengths);
+  
+    cellvolume = 1;
+    for(int i=0;i<NDIM;i++) {
+      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths.x[i];
     }
+  
+    // set new box dimensions
+    *out << Verbose(0) << "Translating to box with these boundaries." << endl;
+    mol->CenterInBox((ofstream *)&cout, &BoxLengths);
+  }
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
-  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
-      << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
-      << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
-      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-}
-;
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+};
+
 
 // =========================================================== class TESSELATION ===========================================
@@ -934 +648 @@
 Tesselation::Tesselation()
 {
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
+  PointsOnBoundaryCount = 0; 
+  LinesOnBoundaryCount = 0; 
   TrianglesOnBoundaryCount = 0;
-  TriangleFilesWritten = 0;
-}
-;
+};
 
 /** Constructor of class Tesselation.
@@ -946 +658 @@
 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;
-        }
-        for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
-                if (runner->second != NULL) {
-                        delete (runner->second);
-                        runner->second = NULL;
-                } else
-                        cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
-        }
-        for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
-                if (runner->second != NULL) {
-                        delete (runner->second);
-                        runner->second = NULL;
-                } else
-                        cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
-        }
-}
-;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    delete(runner->second);
+  }
+};
 
 /** Gueses first starting triangle of the convex envelope.
@@ -975 +667 @@
  * \param *out output stream for debugging
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
- */
-void
-Tesselation::GuessStartingTriangle(ofstream *out)
+ */ 
+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.Distance(&B->second->node->x);
-              distance = tmp * tmp;
-              tmp = A->second->node->x.Distance(&C->second->node->x);
-              distance += tmp * tmp;
-              tmp = B->second->node->x.Distance(&C->second->node->x);
-              distance += tmp * tmp;
-              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
-                  PointMap::iterator, PointMap::iterator> (B, C)));
-            }
-        }
+  double distance;
+  for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+    for(PointMap::iterator sprinter = PointsOnBoundary.begin(); sprinter != PointsOnBoundary.end(); sprinter++) {
+      if (runner->first < sprinter->first) {
+        distance = runner->second->node->x.Distance(&sprinter->second->node->x);
+        DistanceMMap.insert( DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(runner,sprinter) ) );
+      }
     }
-  //    // 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
+  }
+
+//    // 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. take three smallest distance that form a triangle
+  // we take the smallest distance as the base line
   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.Distance(&A->second->node->x);
-          int innerpoint = 0;
-          if ((tmp < A->second->node->x.Distance(
-              &baseline->second.first->second->node->x)) && (tmp
-              < A->second->node->x.Distance(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.Distance(
-              &baseline->second.first->second->node->x);
-          if ((tmp < baseline->second.first->second->node->x.Distance(
-              &A->second->node->x)) && (tmp
-              < baseline->second.first->second->node->x.Distance(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.Distance(
-              &baseline->second.second->second->node->x);
-          if ((tmp < baseline->second.second->second->node->x.Distance(
-              &baseline->second.first->second->node->x)) && (tmp
-              < baseline->second.second->second->node->x.Distance(
-                  &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);
-    }
-}
-;
+  BPS[0] = baseline->second.first->second;
+  BPS[1] = baseline->second.second->second;
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
+  // take the second smallest as the second base line
+  DistanceMultiMap::iterator secondline = DistanceMMap.begin();
+  do {
+    secondline++;
+  } while (!(
+      ((BPS[0] == secondline->second.first->second) && (BPS[1] != secondline->second.second->second)) ||
+      ((BPS[0] == secondline->second.second->second) && (BPS[1] != secondline->second.first->second)) ||
+      ((BPS[1] == secondline->second.first->second) && (BPS[0] != secondline->second.second->second)) ||
+      ((BPS[1] == secondline->second.second->second) && (BPS[0] != secondline->second.first->second))
+    ));
+  BPS[0] = secondline->second.first->second;
+  BPS[1] = secondline->second.second->second;
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
+  // connection yields the third line (note: first and second endpoint are sorted!)
+  if (baseline->second.first->second == secondline->second.first->second) {
+    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.second->second);
+  } else if (baseline->second.first->second == secondline->second.second->second) {
+    SetEndpointsOrdered(BPS, baseline->second.second->second, secondline->second.first->second);
+  } else if (baseline->second.second->second == secondline->second.first->second) {
+    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.second->second);
+  } else if (baseline->second.second->second == secondline->second.second->second) {
+    SetEndpointsOrdered(BPS, baseline->second.first->second, baseline->second.first->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;
+};
+
 
 /** Tesselates the convex envelope of a cluster from a single starting triangle.
@@ -1147 +749 @@
  * \param *mol the cluster as a molecule structure
  */
-void
-Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
-    molecule *mol)
+void Tesselation::TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol)
 {
   bool flag;
@@ -1155 +755 @@
   class BoundaryPointSet *peak = NULL;
   double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
-      PropagationVector;
+  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++;
+  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;
               }
-            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);
-
-}
-;
+        // 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++;
+        } 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);
+  
+  stringstream line;
+  line << configuration->configpath << "/" << CONVEXENVELOPE;
+  *out << Verbose(1) << "Storing convex envelope in tecplot data file " << line.str() << "." << endl;
+  ofstream output(line.str().c_str());
+  output << "TITLE = \"3D CONVEX SHELL\"" << endl;
+  output << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+  output << "ZONE T=\"TRIANGLES\", N=" << PointsOnBoundaryCount << ", E=" << 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 = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
+    Walker = target->second->node;
+    LookupList[Walker->nr] = Counter++;
+    output << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+  }
+  output << endl;
+    // print connectivity
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
+    output << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 
+  }
+  output.close();
+  delete[](LookupList);
+  *out << endl;
+};
 
 /** Adds an atom to the tesselation::PointsOnBoundary list.
  * \param *Walker atom to add
  */
-void
-Tesselation::AddPoint(atom *Walker)
-{
-  PointTestPair InsertUnique;
+void Tesselation::AddPoint(atom *Walker)
+{
   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(2) << "Atom " << *((InsertUnique.first)->second->node)
-          << " gibt's schon in der PointMap." << endl;
-      TPS[n] = (InsertUnique.first)->second;
-    }
-}
-;
-
-/** Function tries to add line from current Points in BPS to BoundaryLineSet.
- * If succesful it raises the line count and inserts the new line into the BLS,
- * if unsuccesful, 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)
-{
-  LineMap::iterator LineWalker;
-  //cout << "Manually checking endpoints for line." << endl;
-  if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
-  //If a line is there, how do I recognize that beyond a shadow of a doubt?
-    {
-      //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl;
-
-      LineWalker = LinesOnBoundary.end();
-      LineWalker--;
-
-      while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr,
-          b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max(
-          a->node->nr, b->node->nr))
-        {
-          //cout << Verbose(1) << "Looking for line which already exists"<< endl;
-          LineWalker--;
-        }
-      BPS[0] = LineWalker->second->endpoints[0];
-      BPS[1] = LineWalker->second->endpoints[1];
-      BLS[n] = LineWalker->second;
-
-    }
-  else
-    {
-      cout << Verbose(2)
-          << "Adding line which has not been used before between "
-          << *(a->node) << " and " << *(b->node) << "." << endl;
-      BPS[0] = a;
-      BPS[1] = b;
-      BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-
-      LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
-      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::AddTriangleToLines()
-{
-
-  cout << Verbose(1) << "Adding triangle to its lines" << endl;
-  int i = 0;
-  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-  TrianglesOnBoundaryCount++;
-
-  /*
-   * this is apparently done when constructing triangle
-
-   for (i=0; i<3; i++)
-   {
-   BLS[i]->AddTriangle(BTS);
-   }
-   */
-}
-;
-
-/**
- * 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;
-    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(4) << "Center of sphere of circumference is " << *Center << ".\n";
-    cout << Verbose(3) << "End of Get_center_of_sphere.\n";
-  }
-  ;
-
-
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
- *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
- *  upon which we operate.
- *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
- *  direction and angle into Storage.
- *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
- *  with all neighbours of the candidate.
- * @param a first point
- * @param b second point
- * *param c atom old third point of old triangle
- * @param Candidate base point along whose bonds to start looking from
- * @param Parent point to avoid during search as its wrong direction
- * @param RecursionLevel contains current recursion depth
- * @param Chord baseline vector of first and second point
- * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
- * @param OldNormal normal of the triangle which the baseline belongs to
- * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere
- * @param Opt_Candidate candidate reference to return
- * @param Storage array containing two angles of current Opt_Candidate
- * @param RADIUS radius of ball
- * @param mol molecule structure with atoms and bonds
- */
-
-void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
-    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
-{
-        cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-        cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-        cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-        cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-        cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-        cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
-        /* OldNormal is normal vector on the old triangle
-         * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
-         * Chord points from b to a!!!
-         */
-        Vector dif_a; //Vector from a to candidate
-        Vector dif_b; //Vector from b to candidate
-        Vector AngleCheck;
-        Vector TempNormal, Umkreismittelpunkt;
-        Vector Mittelpunkt;
-
-        double CurrentEpsilon = 0.1;
-        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-        double BallAngle, AlternativeSign;
-        atom *Walker; // variable atom point
-
-        Vector NewUmkreismittelpunkt;
-
-        if (a != Candidate and b != Candidate and c != Candidate) {
-                cout << Verbose(3) << "We have a unique candidate!" << endl;
-                dif_a.CopyVector(&(a->x));
-                dif_a.SubtractVector(&(Candidate->x));
-                dif_b.CopyVector(&(b->x));
-                dif_b.SubtractVector(&(Candidate->x));
-                AngleCheck.CopyVector(&(Candidate->x));
-                AngleCheck.SubtractVector(&(a->x));
-                AngleCheck.ProjectOntoPlane(Chord);
-
-                SideA = dif_b.Norm();
-                SideB = dif_a.Norm();
-                SideC = Chord->Norm();
-                //Chord->Scale(-1);
-
-                alpha = Chord->Angle(&dif_a);
-                beta = M_PI - Chord->Angle(&dif_b);
-                gamma = dif_a.Angle(&dif_b);
-
-                cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-                if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
-                        cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-                        cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-                }
-
-                if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
-                        Umkreisradius = SideA / 2.0 / sin(alpha);
-                        //cout << Umkreisradius << endl;
-                        //cout << SideB / 2.0 / sin(beta) << endl;
-                        //cout << SideC / 2.0 / sin(gamma) << endl;
-
-                        if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
-                                cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
-                                cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
-                                sign = AngleCheck.ScalarProduct(direction1);
-                                if (fabs(sign)<MYEPSILON) {
-                                        if (AngleCheck.ScalarProduct(OldNormal)<0) {
-                                                sign =0;
-                                                AlternativeSign=1;
-                                        } else {
-                                                sign =0;
-                                                AlternativeSign=-1;
-                                        }
-                                } else {
-                                        sign /= fabs(sign);
-                                        cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-                                }
-
-                                Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-                                AngleCheck.CopyVector(&ReferencePoint);
-                                AngleCheck.Scale(-1);
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-                                AngleCheck.AddVector(&Mittelpunkt);
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-                                cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-                                BallAngle = AngleCheck.Angle(OldNormal);
-                                cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-
-                                //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-                                //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-
-                                cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
-
-                                NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
-
-                                if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
-                                        if (Storage[0]< -1.5) { // first Candidate at all
-                                                if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-                                                        cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
-                                                        Opt_Candidate = Candidate;
-                                                        Storage[0] = sign;
-                                                        Storage[1] = AlternativeSign;
-                                                        Storage[2] = BallAngle;
-                                                        cout << " angle " << Storage[2] << " and Up/Down "
-                                                        << Storage[0] << endl;
-                                                } else
-                                                        cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-                                        } else {
-                                                if ( Storage[2] > BallAngle) {
-                                                        if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
-                                                                cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
-                                                                Opt_Candidate = Candidate;
-                                                                Storage[0] = sign;
-                                                                Storage[1] = AlternativeSign;
-                                                                Storage[2] = BallAngle;
-                                                                cout << " angle " << Storage[2] << " and Up/Down "
-                                                                << Storage[0] << endl;
-                                                        } else
-                                                                cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
-                                                } else {
-                                                        if (DEBUG) {
-                                                                cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
-                                                        }
-                                                }
-                                        }
-                                } else {
-                                        if (DEBUG) {
-                                                cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
-                                        }
-                                }
-                        } else {
-                                if (DEBUG) {
-                                        cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
-                                }
-                        }
-                } else {
-                        if (DEBUG) {
-                                cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
-                        }
-                }
-        } else {
-                if (DEBUG) {
-                        cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
-                }
-        }
-
-        if (RecursionLevel < 5) { // Seven is the recursion level threshold.
-                for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) { // don't go back the same bond
-                                continue;
-                        } else {
-                                Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
-                        }
-                }
-        }
-        cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-}
-;
-
-
-  /** This recursive function finds a third point, to form a triangle with two given ones.
-   * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
-   *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
-   *  upon which we operate.
-   *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
-   *  direction and angle into Storage.
-   *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
-   *  with all neighbours of the candidate.
-   * @param a first point
-   * @param b second point
-   * @param Candidate base point along whose bonds to start looking from
-   * @param Parent point to avoid during search as its wrong direction
-   * @param RecursionLevel contains current recursion depth
-   * @param Chord baseline vector of first and second point
-   * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
-   * @param OldNormal normal of the triangle which the baseline belongs to
-   * @param Opt_Candidate candidate reference to return
-   * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
-   * @param Storage array containing two angles of current Opt_Candidate
-   * @param RADIUS radius of ball
-   * @param mol molecule structure with atoms and bonds
-   */
-void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
-    atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
-{
-        /* OldNormal is normal vector on the old triangle
-         * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
-         * Chord points from b to a!!!
-         */
-        Vector dif_a; //Vector from a to candidate
-        Vector dif_b; //Vector from b to candidate
-        Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
-        Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
-
-        double CurrentEpsilon = 0.1;
-        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-        double BallAngle;
-        atom *Walker; // variable atom point
-
-
-        dif_a.CopyVector(&(a->x));
-        dif_a.SubtractVector(&(Candidate->x));
-        dif_b.CopyVector(&(b->x));
-        dif_b.SubtractVector(&(Candidate->x));
-        DirectionCheckPoint.CopyVector(&dif_a);
-        DirectionCheckPoint.Scale(-1);
-        DirectionCheckPoint.ProjectOntoPlane(Chord);
-
-        SideA = dif_b.Norm();
-        SideB = dif_a.Norm();
-        SideC = Chord->Norm();
-        //Chord->Scale(-1);
-
-        alpha = Chord->Angle(&dif_a);
-        beta = M_PI - Chord->Angle(&dif_b);
-        gamma = dif_a.Angle(&dif_b);
-
-
-        if (DEBUG) {
-                cout << "Atom number" << Candidate->nr << endl;
-                Candidate->x.Output((ofstream *) &cout);
-                cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
-        }
-
-        if (a != Candidate and b != Candidate) {
-                //      alpha = dif_a.Angle(&dif_b) / 2.;
-                //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-                //      SideB = dif_a.Norm();
-                //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-                //          alpha); // note this is squared of center line length
-                //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-                // Those are remains from Freddie. Needed?
-
-                Umkreisradius = SideA / 2.0 / sin(alpha);
-                //cout << Umkreisradius << endl;
-                //cout << SideB / 2.0 / sin(beta) << endl;
-                //cout << SideC / 2.0 / sin(gamma) << endl;
-
-                if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
-                        // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
-                        Umkreismittelpunkt.Zero();
-                        helper.CopyVector(&a->x);
-                        helper.Scale(sin(2.*alpha));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        helper.CopyVector(&b->x);
-                        helper.Scale(sin(2.*beta));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        helper.CopyVector(&Candidate->x);
-                        helper.Scale(sin(2.*gamma));
-                        Umkreismittelpunkt.AddVector(&helper);
-                        //Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-                        Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
-                        TempNormal.CopyVector(&dif_a);
-                        TempNormal.VectorProduct(&dif_b);
-                        if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
-                                TempNormal.Scale(-1);
-                        }
-                        TempNormal.Normalize();
-                        Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-                        TempNormal.Scale(Restradius);
-
-                        Mittelpunkt.CopyVector(&Umkreismittelpunkt);
-                        Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-
-                        AngleCheck.CopyVector(Chord);
-                        AngleCheck.Scale(-0.5);
-                        AngleCheck.SubtractVector(&(b->x));
-                        AngleCheckReference.CopyVector(&AngleCheck);
-                        AngleCheckReference.AddVector(Opt_Mittelpunkt);
-                        AngleCheck.AddVector(&Mittelpunkt);
-
-                        BallAngle = AngleCheck.Angle(&AngleCheckReference);
-
-                        d1->ProjectOntoPlane(&AngleCheckReference);
-                        sign = AngleCheck.ScalarProduct(d1);
-                        if (fabs(sign) < MYEPSILON)
-                                sign = 0;
-                        else
-                                sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-
-                        if (Storage[0]< -1.5) { // first Candidate at all
-                                cout << "Next better candidate is " << *Candidate << " with ";
-                                Opt_Candidate = Candidate;
-                                Storage[0] = sign;
-                                Storage[1] = BallAngle;
-                                Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                        } else {
-                                /*
-                                 * removed due to change in criterium, now checking angle of ball to old normal.
-                                //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
-                                //within the ball.
-
-                                Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-                                //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
-
-
-                                if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
-                                 */
-                                        //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
-
-                                if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
-                                        //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
-                                        cout << "Next better candidate is " << *Candidate << " with ";
-                                        Opt_Candidate = Candidate;
-                                        Storage[0] = sign;
-                                        Storage[1] = BallAngle;
-                                        Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                        cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                                } else {
-                                        if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
-                                                //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-                                                cout << "Next better candidate is " << *Candidate << " with ";
-                                                Opt_Candidate = Candidate;
-                                                Storage[0] = sign;
-                                                Storage[1] = BallAngle;
-                                                Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                                                cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
-                                        }
-                                }
-                        }
-/*
-               * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
-               *
-                else
-                {
-                  if (sign>0 && BallAngle>0 && Storage[0]<0)
-                    {
-                      cout << "Next better candidate is " << *Candidate << " with ";
-                      Opt_Candidate = Candidate;
-                      Storage[0] = sign;
-                      Storage[1] = BallAngle;
-                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
-                      << Storage[0] << endl;
-
-//Debugging purposes only
-                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-                      cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-                      cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-                      cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
-                      cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-                      cout << "Umkreisradius ist " << Umkreisradius << endl;
-                      cout << "Restradius ist " << Restradius << endl;
-                      cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
-                      cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
-                      cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-                      cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-                      cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
-
-
-
-                    }
-                  else
-                    {
-                      if (DEBUG)
-                        cout << "Looses to better candidate" << endl;
-                    }
-                }
-                */
-                } else {
-                        if (DEBUG) {
-                                cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
-                        }
-                }
-        } else {
-                if (DEBUG) {
-                        cout << "identisch mit Ursprungslinie" << endl;
-                }
-        }
-
-        if (RecursionLevel < 9) { // Five is the recursion level threshold.
-                for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) { // don't go back the same bond
-                                continue;
-                        } else {
-                                Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
-                        }
-                }
-        }
-};
-
-/** This function finds a triangle to a line, adjacent to an existing one.
- * @param out   output stream for debugging
- * @param tecplot output stream for writing found triangles in TecPlot format
- * @param mol molecule structure with all atoms and bonds
- * @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
- */
-bool Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
-    molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N, const char *tempbasename)
-{
-        cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
-        Vector direction1;
-        Vector helper;
-        Vector Chord;
-        ofstream *tempstream = NULL;
-        char NumberName[255];
-        double tmp;
-        //atom* Walker;
-        atom* OldThirdPoint;
-
-        double Storage[3];
-        Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-        Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
-        Storage[2] = 9999999.;
-        atom* Opt_Candidate = NULL;
-        Vector Opt_Mittelpunkt;
-
-        cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-
-        helper.CopyVector(&(Line.endpoints[0]->node->x));
-        for (int i = 0; i < 3; i++) {
-                if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr) {
-                        OldThirdPoint = T.endpoints[i]->node;
-                        helper.SubtractVector(&T.endpoints[i]->node->x);
-                        break;
-                }
-        }
-
-        direction1.CopyVector(&Line.endpoints[0]->node->x);
-        direction1.SubtractVector(&Line.endpoints[1]->node->x);
-        direction1.VectorProduct(&(T.NormalVector));
-
-        if (direction1.ScalarProduct(&helper) < 0) {
-                direction1.Scale(-1);
-        }
-        cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
-
-        Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
-        Chord.SubtractVector(&(Line.endpoints[1]->node->x));
-        cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
-
-
-        Vector Umkreismittelpunkt, a, b, c;
-        double alpha, beta, gamma;
-        a.CopyVector(&(T.endpoints[0]->node->x));
-        b.CopyVector(&(T.endpoints[1]->node->x));
-        c.CopyVector(&(T.endpoints[2]->node->x));
-        a.SubtractVector(&(T.endpoints[1]->node->x));
-        b.SubtractVector(&(T.endpoints[2]->node->x));
-        c.SubtractVector(&(T.endpoints[0]->node->x));
-
-        alpha = M_PI - a.Angle(&c);
-        beta = M_PI - b.Angle(&a);
-        gamma = M_PI - c.Angle(&b);
-        if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-                cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-        Umkreismittelpunkt.Zero();
-        helper.CopyVector(&T.endpoints[0]->node->x);
-        helper.Scale(sin(2.*alpha));
-        Umkreismittelpunkt.AddVector(&helper);
-        helper.CopyVector(&T.endpoints[1]->node->x);
-        helper.Scale(sin(2.*beta));
-        Umkreismittelpunkt.AddVector(&helper);
-        helper.CopyVector(&T.endpoints[2]->node->x);
-        helper.Scale(sin(2.*gamma));
-        Umkreismittelpunkt.AddVector(&helper);
-        //Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-        //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-        Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-        //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-        cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
-        cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
-
-        cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-        cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
-        if (DEBUG)
-                cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
-        tmp = 0;
-        for (int i=0;i<NDIM;i++) {
-                helper.CopyVector(&T.endpoints[i]->node->x);
-                helper.SubtractVector(&Umkreismittelpunkt);
-                if (DEBUG)
-                        cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
-                if (tmp == 0) // set on first time for comparison against next ones
-                        tmp = helper.Norm();
-                if (fabs(helper.Norm() - tmp) > MYEPSILON)
-                        cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
-        }
-
-        cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
-
-        Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-        Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-        if (Opt_Candidate == NULL) {
-                cerr << "WARNING: Could not find a suitable candidate." << endl;
-                return false;
-        }
-        cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
-
-        // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-        bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
-
-        if (flag) { // if so, add
-                AddTrianglePoint(Opt_Candidate, 0);
-                AddTrianglePoint(Line.endpoints[0]->node, 1);
-                AddTrianglePoint(Line.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);
-                AddTriangleToLines();
-
-                BTS->GetNormalVector(BTS->NormalVector);
-
-                if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
-                        BTS->NormalVector.Scale(-1);
-                };
-                cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-                cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
-        } else { // else, yell and do nothing
-                cout << Verbose(1) << "This triangle consisting of ";
-                cout << *Opt_Candidate << ", ";
-                cout << *Line.endpoints[0]->node << " and ";
-                cout << *Line.endpoints[1]->node << " ";
-                cout << "is invalid!" << endl;
-                return false;
-        }
-
-        if ((TrianglesOnBoundaryCount % 10) == 0) {
-                sprintf(NumberName, "-%d", TriangleFilesWritten);
-                if (DoTecplotOutput) {
-                        string NameofTempFile(tempbasename);
-                        NameofTempFile.append(NumberName);
-                        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);
-                        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);
-                        tempstream->close();
-                        tempstream->flush();
-                        delete(tempstream);
-                }
-                if (DoTecplotOutput || DoRaster3DOutput)
-                        TriangleFilesWritten++;
-        }
-
-        cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
-        return true;
-};
-
-/** 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 *a first endpoint
- * \param *b second endpoint
- * \param *c third endpoint
- * \return false - triangle invalid due to edge criteria, true - triangle may be added.
- */
-bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c) {
-        LineMap::iterator TryAndError;
-        PointTestPair InsertPoint;
-        bool Present[3];
-
-        *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-        TPS[0] = new class BoundaryPointSet(a);
-        TPS[1] = new class BoundaryPointSet(b);
-        TPS[2] = new class BoundaryPointSet(c);
-        for (int i=0;i<3;i++) { // check through all endpoints
-                InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
-                Present[i] = !InsertPoint.second;
-                if (Present[i]) { // if new point was not present before, increase counter
-                        delete TPS[i];
-                        *out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
-                        TPS[i] = (InsertPoint.first)->second;
-                }
-        }
-
-        // check lines
-        for (int i=0;i<3;i++)
-                if (Present[i])
-                        for (int j=i;j<3;j++)
-                                if (Present[j]) {
-                                        TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
-                                        if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
-                                                *out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
-                                                *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-                                                return false;
-                                        }
-                                }
-        *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
-        return true;
-};
-
-
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[3],
-    molecule* mol, double RADIUS)
-{
-        cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
-        int i;
-        Vector AngleCheck;
-        atom* Walker;
-        double norm = -1., angle;
-
-        // check if we only have one unique point yet ...
-        if (a != Candidate) {
-                cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
-                AngleCheck.CopyVector(&(Candidate->x));
-                AngleCheck.SubtractVector(&(a->x));
-                norm = AngleCheck.Norm();
-                // second point shall have smallest angle with respect to Oben vector
-                if (norm < RADIUS) {
-                        angle = AngleCheck.Angle(&Oben);
-                        if (angle < Storage[0]) {
-                                //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                                cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
-                                Opt_Candidate = Candidate;
-                                Storage[0] = AngleCheck.Angle(&Oben);
-                                //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-                        } else {
-                                cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
-                        }
-                } else {
-                        cout << "Refused due to Radius " << norm << endl;
-                }
-        }
-
-        // if not recursed to deeply, look at all its bonds
-        if (RecursionLevel < 7) {
-                for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
-                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                        if (Walker == Parent) // don't go back along the bond we came from
-                                continue;
-                        else
-                                Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
-                }
-        }
-        cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
-};
-
-void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
-{
-        cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-        int i = 0;
-        atom* Walker;
-        atom* FirstPoint;
-        atom* SecondPoint;
-        atom* max_index[NDIM];
-        double max_coordinate[NDIM];
-        Vector Oben;
-        Vector helper;
-        Vector Chord;
-        Vector CenterOfFirstLine;
-
-        Oben.Zero();
-
-        for (i = 0; i < 3; i++) {
-                max_index[i] = NULL;
-                max_coordinate[i] = -1;
-        }
-        cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
-
-        // 1. searching topmost atom with respect to each axis
-        Walker = mol->start;
-        while (Walker->next != mol->end) {
-                Walker = Walker->next;
-                for (i = 0; i < 3; i++) {
-                        if (Walker->x.x[i] > max_coordinate[i]) {
-                                max_coordinate[i] = Walker->x.x[i];
-                                max_index[i] = Walker;
-                        }
-                }
-        }
-
-        cout << Verbose(2) << "Found maximum coordinates: ";
-        for (int i=0;i<NDIM;i++)
-                cout << i << ": " << *max_index[i] << "\t";
-        cout << endl;
-  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-        const int k = 1;
-        Oben.x[k] = 1.;
-        FirstPoint = max_index[k];
-
-        cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-        double Storage[3];
-        atom* Opt_Candidate = NULL;
-        Storage[0] = 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.
-        Storage[1] = 999999.; // This will be an angle looking for the third point.
-        Storage[2] = 999999.;
-
-        Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
-        SecondPoint = Opt_Candidate;
-        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);
-        Storage[0] = -2.; // This will indicate the quadrant.
-        Storage[1] = 9999999.; // This will be an angle looking for the third point.
-        Storage[2] = 9999999.;
-
-        Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-        Chord.SubtractVector(&(SecondPoint->x));
-        // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-        cout << Verbose(2) << "Looking for third point candidates \n";
-        // look in one direction of baseline for initial candidate
-        Opt_Candidate = NULL;
-        CenterOfFirstLine.CopyVector(&Chord);
-        CenterOfFirstLine.Scale(0.5);
-        CenterOfFirstLine.AddVector(&(SecondPoint->x));
-
-        cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-        Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
-        // look in other direction of baseline for possible better candidate
-        cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-        Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
-        cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
-
-        // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
-
-        // Finally, we only have to add the found points
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTrianglePoint(Opt_Candidate, 2);
-        // ... and respective lines
-        AddTriangleLine(TPS[0], TPS[1], 0);
-        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);
-        AddTriangleToLines();
-        // ... and calculate its normal vector (with correct orientation)
-        Oben.Scale(-1.);
-        BTS->GetNormalVector(Oben);
-        cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
-        cout << Verbose(1) << "End of Find_starting_triangle\n";
-};
-
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol, const char *filename, const double RADIUS)
-{
-        int N = 0;
-        struct Tesselation *Tess = new Tesselation;
-        cout << Verbose(1) << "Entering search for non convex hull. " << endl;
-        cout << flush;
-        LineMap::iterator baseline;
-        cout << Verbose(0) << "Begin of Find_non_convex_border\n";
-        bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-        bool failflag = false;
-        if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
-                mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
-
-        Tess->Find_starting_triangle(mol, RADIUS);
-
-        baseline = Tess->LinesOnBoundary.begin();
-        while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
-                if (baseline->second->TrianglesCount == 1) {
-                        failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
-                        flag = flag || failflag;
-                        if (!failflag)
-                                cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
-                } else {
-                        cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
-                }
-                N++;
-                baseline++;
-                if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-                        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-                        flag = false;
-                }
-        }
-        if (failflag) {
-                cout << Verbose(1) << "Writing final tecplot file\n";
-                if (DoTecplotOutput)
-                        write_tecplot_file(out, tecplot, Tess, mol, -1);
-                if (DoRaster3DOutput)
-                        write_raster3d_file(out, tecplot, Tess, mol);
-        } else {
-                cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
-        }
-
-        cout << Verbose(0) << "End of Find_non_convex_border\n";
-        delete(Tess);
-};
-
+  PointsOnBoundary.insert( PointPair(Walker->nr, BPS[0]) );
+  PointsOnBoundaryCount++;
+};