Context Navigation

-              r178f92
+              r437922
 BoundaryPointSet::BoundaryPointSet()
+{
         LinesCount = 0;
         Nr = -1;
+  LinesCount = 0;
+  Nr = -1;
+}
+;
 …
 BoundaryPointSet::BoundaryPointSet(atom *Walker)
+{
         node = Walker;
         LinesCount = 0;
         Nr = Walker->nr;
+  node = Walker;
+  LinesCount = 0;
+  Nr = Walker->nr;
+}
+;
 …
 BoundaryPointSet::~BoundaryPointSet()
+{
         cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
         if (!lines.empty())
                 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
         node = NULL;
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  if (!lines.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
+  node = NULL;
+}
+;
 …
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
         cout << Verbose(6) << "Adding " << *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));
+                }
         LinesCount++;
+  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));
+    }
+  LinesCount++;
+}
+;
 …
 operator <<(ostream &ost, BoundaryPointSet &a)
+{
         ost << "[" << a.Nr << "|" << a.node->Name << "]";
         return ost;
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  return ost;
+}
+;
 …
 BoundaryLineSet::BoundaryLineSet()
+{
         for (int i = 0; i < 2; i++)
                 endpoints[i] = NULL;
         TrianglesCount = 0;
         Nr = -1;
+  for (int i = 0; i < 2; i++)
+    endpoints[i] = NULL;
+  TrianglesCount = 0;
+  Nr = -1;
+}
+;
 …
 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
         // set number
         Nr = number;
         // set endpoints in ascending order
         SetEndpointsOrdered(endpoints, Point[0], Point[1]);
         // add this line to the hash maps of both endpoints
         Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
         Point[1]->AddLine(this); //
         // clear triangles list
         TrianglesCount = 0;
         cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
+  // add this line to the hash maps of both endpoints
+  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+  Point[1]->AddLine(this); //
+  // clear triangles list
+  TrianglesCount = 0;
+  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+}
+;
 …
 BoundaryLineSet::~BoundaryLineSet()
+{
         int Numbers[2];
         Numbers[0] = endpoints[1]->Nr;
         Numbers[1] = endpoints[0]->Nr;
         for (int i = 0; i < 2; i++) {
                 cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
                 endpoints[i]->lines.erase(Numbers[i]);
                 if (endpoints[i]->lines.empty()) {
                         cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
                         if (endpoints[i] != NULL) {
                                 delete(endpoints[i]);
                                 endpoints[i] = NULL;
                         } else
                                 cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
                 } else
                         cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+        }
         if (!triangles.empty())
                 cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+  int Numbers[2];
+  Numbers[0] = endpoints[1]->Nr;
+  Numbers[1] = endpoints[0]->Nr;
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    endpoints[i]->lines.erase(Numbers[i]);
+    if (endpoints[i]->lines.empty()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+  if (!triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+}
+;
 …
 BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
         cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
                         << endl;
         triangles.insert(TrianglePair(triangle->Nr, triangle));
         TrianglesCount++;
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+      << endl;
+  triangles.insert(TrianglePair(triangle->Nr, triangle));
+  TrianglesCount++;
+}
+;
 …
 operator <<(ostream &ost, BoundaryLineSet &a)
+{
         ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
                         << a.endpoints[1]->node->Name << "]";
         return ost;
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "]";
+  return ost;
+}
+;
 …
 BoundaryTriangleSet::BoundaryTriangleSet()
+{
         for (int i = 0; i < 3; i++)
+                {
                         endpoints[i] = NULL;
                         lines[i] = NULL;
+                }
         Nr = -1;
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
+  Nr = -1;
+}
+;
 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3],
                 int number)
+{
         // set number
         Nr = number;
         // set lines
         cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
         for (int i = 0; i < 3; i++)
+                {
                         lines[i] = line[i];
                         lines[i]->AddTriangle(this);
+                }
         // get ascending order of endpoints
         map<int, class BoundaryPointSet *> OrderMap;
         for (int i = 0; i < 3; i++)
                 // for all three lines
                 for (int j = 0; j < 2; j++)
                         { // for both endpoints
                                 OrderMap.insert(pair<int, class BoundaryPointSet *> (
                                                 line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
                                 // and we don't care whether insertion fails
+                        }
         // set endpoints
         int Counter = 0;
         cout << Verbose(6) << " with end points ";
         for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
                         != OrderMap.end(); runner++)
+                {
                         endpoints[Counter] = runner->second;
                         cout << " " << *endpoints[Counter];
                         Counter++;
+                }
         if (Counter < 3)
+                {
                         cerr << "ERROR! We have a triangle with only two distinct endpoints!"
                                         << endl;
                         //exit(1);
+                }
         cout << "." << endl;
+    int number)
+{
+  // set number
+  Nr = number;
+  // set lines
+  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  for (int i = 0; i < 3; i++)
+    {
+      lines[i] = line[i];
+      lines[i]->AddTriangle(this);
+    }
+  // get ascending order of endpoints
+  map<int, class BoundaryPointSet *> OrderMap;
+  for (int i = 0; i < 3; i++)
+    // for all three lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+        // and we don't care whether insertion fails
+      }
+  // set endpoints
+  int Counter = 0;
+  cout << Verbose(6) << " with end points ";
+  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+      != OrderMap.end(); runner++)
+    {
+      endpoints[Counter] = runner->second;
+      cout << " " << *endpoints[Counter];
+      Counter++;
+    }
+  if (Counter < 3)
+    {
+      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+          << endl;
+      //exit(1);
+    }
+  cout << "." << endl;
+}
+;
 …
 BoundaryTriangleSet::~BoundaryTriangleSet()
+{
         for (int i = 0; i < 3; i++) {
                 cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
                 lines[i]->triangles.erase(Nr);
                 if (lines[i]->triangles.empty()) {
                         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;
+        }
+  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;
+  }
+}
+;
 …
 BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+{
         // get normal vector
         NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
                         &endpoints[2]->node->x);
         // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
         if (NormalVector.Projection(&OtherVector) > 0)
                 NormalVector.Scale(-1.);
+  // get normal vector
+  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
+      &endpoints[2]->node->x);
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (NormalVector.Projection(&OtherVector) > 0)
+    NormalVector.Scale(-1.);
+}
+;
 …
 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;
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  return ost;
+}
+;
 …
 GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
         class BoundaryLineSet * lines[2] =
                 { line1, line2 };
         class BoundaryPointSet *node = NULL;
         map<int, class BoundaryPointSet *> OrderMap;
         pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
         for (int i = 0; i < 2; i++)
                 // for both lines
                 for (int j = 0; j < 2; j++)
                         { // for both endpoints
                                 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
                                                 lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
                                 if (!OrderTest.second)
                                         { // if insertion fails, we have common endpoint
                                                 node = OrderTest.first->second;
                                                 cout << Verbose(5) << "Common endpoint of lines " << *line1
                                                                 << " and " << *line2 << " is: " << *node << "." << endl;
                                                 j = 2;
                                                 i = 2;
                                                 break;
+                                        }
+                        }
         return node;
+  class BoundaryLineSet * lines[2] =
+    { line1, line2 };
+  class BoundaryPointSet *node = NULL;
+  map<int, class BoundaryPointSet *> OrderMap;
+  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  for (int i = 0; i < 2; i++)
+    // for both lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+        if (!OrderTest.second)
+          { // if insertion fails, we have common endpoint
+            node = OrderTest.first->second;
+            cout << Verbose(5) << "Common endpoint of lines " << *line1
+                << " and " << *line2 << " is: " << *node << "." << endl;
+            j = 2;
+            i = 2;
+            break;
+          }
+      }
+  return node;
+}
+;
 …
 GetBoundaryPoints(ofstream *out, molecule *mol)
+{
         atom *Walker = NULL;
         PointMap PointsOnBoundary;
         LineMap LinesOnBoundary;
         TriangleMap TrianglesOnBoundary;
         *out << Verbose(1) << "Finding all boundary points." << endl;
         Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
         BoundariesTestPair BoundaryTestPair;
         Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
         double radius, angle;
         // 3a. Go through every axis
         for (int axis = 0; axis < NDIM; axis++)
+                {
                         AxisVector.Zero();
                         AngleReferenceVector.Zero();
                         AngleReferenceNormalVector.Zero();
                         AxisVector.x[axis] = 1.;
                         AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
                         AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
                         //              *out << Verbose(1) << "Axisvector is ";
                         //              AxisVector.Output(out);
                         //              *out << " and AngleReferenceVector is ";
                         //              AngleReferenceVector.Output(out);
                         //              *out << "." << endl;
                         //              *out << " and AngleReferenceNormalVector is ";
                         //              AngleReferenceNormalVector.Output(out);
                         //              *out << "." << endl;
                         // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
                         Walker = mol->start;
                         while (Walker->next != mol->end)
+                                {
                                         Walker = Walker->next;
                                         Vector ProjectedVector;
                                         ProjectedVector.CopyVector(&Walker->x);
                                         ProjectedVector.ProjectOntoPlane(&AxisVector);
                                         // correct for negative side
                                         //if (Projection(y) < 0)
                                         //angle = 2.*M_PI - angle;
                                         radius = ProjectedVector.Norm();
                                         if (fabs(radius) > MYEPSILON)
                                                 angle = ProjectedVector.Angle(&AngleReferenceVector);
                                         else
                                                 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
                                         //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
                                         if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+                                                {
                                                         angle = 2. * M_PI - angle;
+                                                }
                                         //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
                                         //ProjectedVector.Output(out);
                                         //*out << endl;
                                         BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
                                                         DistancePair (radius, Walker)));
                                         if (BoundaryTestPair.second)
                                                 { // successfully inserted
+                                                }
                                         else
                                                 { // same point exists, check first r, then distance of original vectors to center of gravity
                                                         *out << Verbose(2)
                                                                         << "Encountered two vectors whose projection onto axis "
                                                                         << axis << " is equal: " << endl;
                                                         *out << Verbose(2) << "Present vector: ";
                                                         BoundaryTestPair.first->second.second->x.Output(out);
                                                         *out << endl;
                                                         *out << Verbose(2) << "New vector: ";
                                                         Walker->x.Output(out);
                                                         *out << endl;
                                                         double tmp = ProjectedVector.Norm();
                                                         if (tmp > BoundaryTestPair.first->second.first)
+                                                                {
                                                                         BoundaryTestPair.first->second.first = tmp;
                                                                         BoundaryTestPair.first->second.second = Walker;
                                                                         *out << Verbose(2) << "Keeping new vector." << endl;
+                                                                }
                                                         else if (tmp == BoundaryTestPair.first->second.first)
+                                                                {
                                                                         if (BoundaryTestPair.first->second.second->x.ScalarProduct(
                                                                                         &BoundaryTestPair.first->second.second->x)
                                                                                         < Walker->x.ScalarProduct(&Walker->x))
                                                                                 { // Norm() does a sqrt, which makes it a lot slower
                                                                                         BoundaryTestPair.first->second.second = Walker;
                                                                                         *out << Verbose(2) << "Keeping new vector." << endl;
+                                                                                }
                                                                         else
+                                                                                {
                                                                                         *out << Verbose(2) << "Keeping present vector." << endl;
+                                                                                }
+                                                                }
                                                         else
+                                                                {
                                                                         *out << Verbose(2) << "Keeping present vector." << endl;
+                                                                }
+                                                }
+                                }
                         // printing all inserted for debugging
                         //              {
                         //                      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
                         //                      int i=0;
                         //                      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
                         //                              if (runner != BoundaryPoints[axis].begin())
                         //                                      *out << ", " << i << ": " << *runner->second.second;
                         //                              else
                         //                                      *out << i << ": " << *runner->second.second;
                         //                              i++;
                         //                      }
                         //                      *out << endl;
                         //              }
                         // 3c. throw out points whose distance is less than the mean of left and right neighbours
                         bool flag = false;
                         do
                                 { // do as long as we still throw one out per round
                                         *out << Verbose(1)
                                                         << "Looking for candidates to kick out by convex condition ... "
                                                         << endl;
                                         flag = false;
                                         Boundaries::iterator left = BoundaryPoints[axis].end();
                                         Boundaries::iterator right = BoundaryPoints[axis].end();
                                         for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
                                                         != BoundaryPoints[axis].end(); runner++)
+                                                {
                                                         // set neighbours correctly
                                                         if (runner == BoundaryPoints[axis].begin())
+                                                                {
                                                                         left = BoundaryPoints[axis].end();
+                                                                }
                                                         else
+                                                                {
                                                                         left = runner;
+                                                                }
                                                         left--;
                                                         right = runner;
                                                         right++;
                                                         if (right == BoundaryPoints[axis].end())
+                                                                {
                                                                         right = BoundaryPoints[axis].begin();
+                                                                }
                                                         // check distance
                                                         // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                                                                {
                                                                         Vector SideA, SideB, SideC, SideH;
                                                                         SideA.CopyVector(&left->second.second->x);
                                                                         SideA.ProjectOntoPlane(&AxisVector);
                                                                         //                                      *out << "SideA: ";
                                                                         //                                      SideA.Output(out);
                                                                         //                                      *out << endl;
                                                                         SideB.CopyVector(&right->second.second->x);
                                                                         SideB.ProjectOntoPlane(&AxisVector);
                                                                         //                                      *out << "SideB: ";
                                                                         //                                      SideB.Output(out);
                                                                         //                                      *out << endl;
                                                                         SideC.CopyVector(&left->second.second->x);
                                                                         SideC.SubtractVector(&right->second.second->x);
                                                                         SideC.ProjectOntoPlane(&AxisVector);
                                                                         //                                      *out << "SideC: ";
                                                                         //                                      SideC.Output(out);
                                                                         //                                      *out << endl;
                                                                         SideH.CopyVector(&runner->second.second->x);
                                                                         SideH.ProjectOntoPlane(&AxisVector);
                                                                         //                                      *out << "SideH: ";
                                                                         //                                      SideH.Output(out);
                                                                         //                                      *out << endl;
                                                                         // calculate each length
                                                                         double a = SideA.Norm();
                                                                         //double b = SideB.Norm();
                                                                         //double c = SideC.Norm();
                                                                         double h = SideH.Norm();
                                                                         // calculate the angles
                                                                         double alpha = SideA.Angle(&SideH);
                                                                         double beta = SideA.Angle(&SideC);
                                                                         double gamma = SideB.Angle(&SideH);
                                                                         double delta = SideC.Angle(&SideH);
                                                                         double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
                                                                                         < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
                                                                         //                                      *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
                                                                         //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
                                                                         if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
                                                                                         < MYEPSILON) && (h < MinDistance))
+                                                                                {
                                                                                         // throw out point
                                                                                         //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
                                                                                         BoundaryPoints[axis].erase(runner);
                                                                                         flag = true;
+                                                                                }
+                                                                }
+                                                }
+                                }
                         while (flag);
+                }
         return BoundaryPoints;
+  atom *Walker = NULL;
+  PointMap PointsOnBoundary;
+  LineMap LinesOnBoundary;
+  TriangleMap TrianglesOnBoundary;
+  *out << Verbose(1) << "Finding all boundary points." << endl;
+  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+  BoundariesTestPair BoundaryTestPair;
+  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+  double radius, angle;
+  // 3a. Go through every axis
+  for (int axis = 0; axis < NDIM; axis++)
+    {
+      AxisVector.Zero();
+      AngleReferenceVector.Zero();
+      AngleReferenceNormalVector.Zero();
+      AxisVector.x[axis] = 1.;
+      AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+      AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+      //    *out << Verbose(1) << "Axisvector is ";
+      //    AxisVector.Output(out);
+      //    *out << " and AngleReferenceVector is ";
+      //    AngleReferenceVector.Output(out);
+      //    *out << "." << endl;
+      //    *out << " and AngleReferenceNormalVector is ";
+      //    AngleReferenceNormalVector.Output(out);
+      //    *out << "." << endl;
+      // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+      Walker = mol->start;
+      while (Walker->next != mol->end)
+        {
+          Walker = Walker->next;
+          Vector ProjectedVector;
+          ProjectedVector.CopyVector(&Walker->x);
+          ProjectedVector.ProjectOntoPlane(&AxisVector);
+          // correct for negative side
+          //if (Projection(y) < 0)
+          //angle = 2.*M_PI - angle;
+          radius = ProjectedVector.Norm();
+          if (fabs(radius) > MYEPSILON)
+            angle = ProjectedVector.Angle(&AngleReferenceVector);
+          else
+            angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+          //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+          if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+            {
+              angle = 2. * M_PI - angle;
+            }
+          //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+          //ProjectedVector.Output(out);
+          //*out << endl;
+          BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
+              DistancePair (radius, Walker)));
+          if (BoundaryTestPair.second)
+            { // successfully inserted
+            }
+          else
+            { // same point exists, check first r, then distance of original vectors to center of gravity
+              *out << Verbose(2)
+                  << "Encountered two vectors whose projection onto axis "
+                  << axis << " is equal: " << endl;
+              *out << Verbose(2) << "Present vector: ";
+              BoundaryTestPair.first->second.second->x.Output(out);
+              *out << endl;
+              *out << Verbose(2) << "New vector: ";
+              Walker->x.Output(out);
+              *out << endl;
+              double tmp = ProjectedVector.Norm();
+              if (tmp > BoundaryTestPair.first->second.first)
+                {
+                  BoundaryTestPair.first->second.first = tmp;
+                  BoundaryTestPair.first->second.second = Walker;
+                  *out << Verbose(2) << "Keeping new vector." << endl;
+                }
+              else if (tmp == BoundaryTestPair.first->second.first)
+                {
+                  if (BoundaryTestPair.first->second.second->x.ScalarProduct(
+                      &BoundaryTestPair.first->second.second->x)
+                      < Walker->x.ScalarProduct(&Walker->x))
+                    { // Norm() does a sqrt, which makes it a lot slower
+                      BoundaryTestPair.first->second.second = Walker;
+                      *out << Verbose(2) << "Keeping new vector." << endl;
+                    }
+                  else
+                    {
+                      *out << Verbose(2) << "Keeping present vector." << endl;
+                    }
+                }
+              else
+                {
+                  *out << Verbose(2) << "Keeping present vector." << endl;
+                }
+            }
+        }
+      // printing all inserted for debugging
+      //    {
+      //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      //      int i=0;
+      //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+      //        if (runner != BoundaryPoints[axis].begin())
+      //          *out << ", " << i << ": " << *runner->second.second;
+      //        else
+      //          *out << i << ": " << *runner->second.second;
+      //        i++;
+      //      }
+      //      *out << endl;
+      //    }
+      // 3c. throw out points whose distance is less than the mean of left and right neighbours
+      bool flag = false;
+      do
+        { // do as long as we still throw one out per round
+          *out << Verbose(1)
+              << "Looking for candidates to kick out by convex condition ... "
+              << endl;
+          flag = false;
+          Boundaries::iterator left = BoundaryPoints[axis].end();
+          Boundaries::iterator right = BoundaryPoints[axis].end();
+          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+              != BoundaryPoints[axis].end(); runner++)
+            {
+              // set neighbours correctly
+              if (runner == BoundaryPoints[axis].begin())
+                {
+                  left = BoundaryPoints[axis].end();
+                }
+              else
+                {
+                  left = runner;
+                }
+              left--;
+              right = runner;
+              right++;
+              if (right == BoundaryPoints[axis].end())
+                {
+                  right = BoundaryPoints[axis].begin();
+                }
+              // check distance
+              // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                {
+                  Vector SideA, SideB, SideC, SideH;
+                  SideA.CopyVector(&left->second.second->x);
+                  SideA.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideA: ";
+                  //          SideA.Output(out);
+                  //          *out << endl;
+                  SideB.CopyVector(&right->second.second->x);
+                  SideB.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideB: ";
+                  //          SideB.Output(out);
+                  //          *out << endl;
+                  SideC.CopyVector(&left->second.second->x);
+                  SideC.SubtractVector(&right->second.second->x);
+                  SideC.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideC: ";
+                  //          SideC.Output(out);
+                  //          *out << endl;
+                  SideH.CopyVector(&runner->second.second->x);
+                  SideH.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideH: ";
+                  //          SideH.Output(out);
+                  //          *out << endl;
+                  // calculate each length
+                  double a = SideA.Norm();
+                  //double b = SideB.Norm();
+                  //double c = SideC.Norm();
+                  double h = SideH.Norm();
+                  // calculate the angles
+                  double alpha = SideA.Angle(&SideH);
+                  double beta = SideA.Angle(&SideC);
+                  double gamma = SideB.Angle(&SideH);
+                  double delta = SideC.Angle(&SideH);
+                  double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
+                      < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+                  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+                  //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+                  if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
+                      < MYEPSILON) && (h < MinDistance))
+                    {
+                      // throw out point
+                      //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+                      BoundaryPoints[axis].erase(runner);
+                      flag = true;
+                    }
+                }
+            }
+        }
+      while (flag);
+    }
+  return BoundaryPoints;
+}
+;
 …
 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;
+                }
         // determine biggest "diameter" of cluster for each axis
         Boundaries::iterator Neighbour, OtherNeighbour;
         double *GreatestDiameter = new double[NDIM];
         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;
+                                                }
+                                }
+                }
         *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);
         return GreatestDiameter;
+    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;
+    }
+  // determine biggest "diameter" of cluster for each axis
+  Boundaries::iterator Neighbour, OtherNeighbour;
+  double *GreatestDiameter = new double[NDIM];
+  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;
+            }
+        }
+    }
+  *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);
+  return GreatestDiameter;
+}
+;
 …
 void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
+{
         atom *Walker = mol->start;
         bond *Binder = mol->first;
         int i;
         Vector *center = mol->DetermineCenterOfAll(out);
         if (vrmlfile != NULL) {
                 //cout << Verbose(1) << "Writing Raster3D file ... ";
                 *vrmlfile << "#VRML V2.0 utf8" << endl;
                 *vrmlfile << "#Created by molecuilder" << endl;
                 *vrmlfile << "#All atoms as spheres" << endl;
                 while (Walker->next != mol->end) {
                         Walker = Walker->next;
                         *vrmlfile << "Sphere {" << endl << "    "; // 2 is sphere type
                         for (i=0;i<NDIM;i++)
                                 *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
                         *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+                }
                 *vrmlfile << "# All bonds as vertices" << endl;
                 while (Binder->next != mol->last) {
                         Binder = Binder->next;
                         *vrmlfile << "3" << endl << "   "; // 2 is round-ended cylinder type
                         for (i=0;i<NDIM;i++)
                                 *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
                         *vrmlfile << "\t0.03\t";
                         for (i=0;i<NDIM;i++)
                                 *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
                         *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+                }
                 *vrmlfile << "# All tesselation triangles" << endl;
                 for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
                         *vrmlfile << "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
                                         *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
                                 *vrmlfile << "\t";
+                        }
                         *vrmlfile << "1. 0. 0." << endl;        // red as colour
                         *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+                }
         } else {
                 cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+        }
         delete(center);
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (vrmlfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *vrmlfile << "#VRML V2.0 utf8" << endl;
+    *vrmlfile << "#Created by molecuilder" << endl;
+    *vrmlfile << "#All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+    *vrmlfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+    *vrmlfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *vrmlfile << "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
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *vrmlfile << "\t";
+      }
+      *vrmlfile << "1. 0. 0." << endl;  // red as colour
+      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+  }
+  delete(center);
 };
 …
 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;
                 *rasterfile << "8\n     25. -1.  1. 1. 1.        0.0            0 0 0 2\n       SOLID            1.0 0.0 0.0\n  BACKFACE        0.3 0.3 1.0      0 0\n";
                 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
+                }
                 *rasterfile << "9\n     terminating special property\n";
         } else {
                 cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+        }
         delete(center);
+  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;
+    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
+    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
+    }
+    *rasterfile << "9\n  terminating special property\n";
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  delete(center);
 };
 …
 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;
+                }
+    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;
+    }
+}
 …
 double
 VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
                 Boundaries *BoundaryPtr, molecule *mol)
+{
         bool IsAngstroem = configuration->GetIsAngstroem();
         atom *Walker = NULL;
         struct Tesselation *TesselStruct = new Tesselation;
         bool BoundaryFreeFlag = false;
         Boundaries *BoundaryPoints = BoundaryPtr;
         double volume = 0.;
         double PyramidVolume = 0.;
         double G, h;
         Vector x, y;
         double a, b, c;
         //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
         // 1. calculate center of gravity
         *out << endl;
         Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
         // 2. translate all points into CoG
         *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
         Walker = mol->start;
         while (Walker->next != mol->end)
+                {
                         Walker = Walker->next;
                         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);
+                        }
         *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;
         // 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;
         // 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.DistanceSquared(
                                         &runner->second->endpoints[1]->node->x));
                         b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
                                         &runner->second->endpoints[2]->node->x));
                         c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
                                         &runner->second->endpoints[1]->node->x));
                         G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
                         x.MakeNormalVector(&runner->second->endpoints[0]->node->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;
         CenterOfGravity->Scale(-1);
         Walker = mol->start;
         while (Walker->next != mol->end)
+                {
                         Walker = Walker->next;
                         Walker->x.Translate(CenterOfGravity);
+                }
         // 8. Store triangles in tecplot file
         string OutputName(filename);
         OutputName.append(TecplotSuffix);
         ofstream *tecplot = new ofstream(OutputName.c_str());
         write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
         tecplot->close();
         delete(tecplot);
         // free reference lists
         if (BoundaryFreeFlag)
                 delete[] (BoundaryPoints);
         return volume;
+    Boundaries *BoundaryPtr, molecule *mol)
+{
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  atom *Walker = NULL;
+  struct Tesselation *TesselStruct = new Tesselation;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G, h;
+  Vector x, y;
+  double a, b, c;
+  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
+  // 1. calculate center of gravity
+  *out << endl;
+  Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
+  // 2. translate all points into CoG
+  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      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);
+      }
+  *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;
+  // 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;
+  // 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.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
+      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+          &runner->second->endpoints[2]->node->x));
+      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
+          &runner->second->endpoints[1]->node->x));
+      G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
+      x.MakeNormalVector(&runner->second->endpoints[0]->node->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;
+  CenterOfGravity->Scale(-1);
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
+  // 8. Store triangles in tecplot file
+  string OutputName(filename);
+  OutputName.append(TecplotSuffix);
+  ofstream *tecplot = new ofstream(OutputName.c_str());
+  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+  tecplot->close();
+  delete(tecplot);
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[] (BoundaryPoints);
+  return volume;
+}
+;
 …
 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();
         Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
         double clustervolume;
         if (ClusterVolume == 0)
                 clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
                                 BoundaryPoints, mol);
         else
                 clustervolume = ClusterVolume;
         double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
                         IsAngstroem);
         Vector BoxLengths;
         int repetition[NDIM] =
                 { 1, 1, 1 };
         int TotalNoClusters = 1;
         for (int i = 0; i < NDIM; i++)
                 TotalNoClusters *= repetition[i];
         // sum up the atomic masses
         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;
         // solve cubic polynomial
         *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
                         << endl;
         double cellvolume;
         if (IsAngstroem)
                 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);
+                }
         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);
+                }
         // 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;
+    double ClusterVolume, double celldensity)
+{
+  // transform to PAS
+  mol->PrincipalAxisSystem(out, true);
+  // some preparations beforehand
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
+  double clustervolume;
+  if (ClusterVolume == 0)
+    clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
+        BoundaryPoints, mol);
+  else
+    clustervolume = ClusterVolume;
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
+      IsAngstroem);
+  Vector BoxLengths;
+  int repetition[NDIM] =
+    { 1, 1, 1 };
+  int TotalNoClusters = 1;
+  for (int i = 0; i < NDIM; i++)
+    TotalNoClusters *= repetition[i];
+  // sum up the atomic masses
+  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;
+  // solve cubic polynomial
+  *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
+      << endl;
+  double cellvolume;
+  if (IsAngstroem)
+    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);
+    }
+  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);
+    }
+  // 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;
+}
+;
 …
 Tesselation::Tesselation()
+{
         PointsOnBoundaryCount = 0;
         LinesOnBoundaryCount = 0;
         TrianglesOnBoundaryCount = 0;
         TriangleFilesWritten = 0;
+  PointsOnBoundaryCount = 0;
+  LinesOnBoundaryCount = 0;
+  TrianglesOnBoundaryCount = 0;
+  TriangleFilesWritten = 0;
+}
+;
 …
 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;
+        }
+  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;
+  }
+}
+;
 …
 Tesselation::GuessStartingTriangle(ofstream *out)
+{
         // 4b. create a starting triangle
         // 4b1. create all distances
         DistanceMultiMap DistanceMMap;
         double distance, tmp;
         Vector PlaneVector, TrialVector;
         PointMap::iterator A, B, C; // three nodes of the first triangle
         A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
         // with A chosen, take each pair B,C and sort
         if (A != PointsOnBoundary.end())
+                {
                         B = A;
                         B++;
                         for (; B != PointsOnBoundary.end(); B++)
+                                {
                                         C = B;
                                         C++;
                                         for (; C != PointsOnBoundary.end(); C++)
+                                                {
                                                         tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
                                                         distance = tmp * tmp;
                                                         tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
                                                         distance += tmp * tmp;
                                                         tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
                                                         distance += tmp * tmp;
                                                         DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
                                                                         PointMap::iterator, PointMap::iterator> (B, C)));
+                                                }
+                                }
+                }
         //              // listing distances
         //              *out << Verbose(1) << "Listing DistanceMMap:";
         //              for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
         //                      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
         //              }
         //              *out << endl;
         // 4b2. pick three baselines forming a triangle
         // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
         DistanceMultiMap::iterator baseline = DistanceMMap.begin();
         for (; baseline != DistanceMMap.end(); baseline++)
+                {
                         // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
                         // 2. next, we have to check whether all points reside on only one side of the triangle
                         // 3. construct plane vector
                         PlaneVector.MakeNormalVector(&A->second->node->x,
                                         &baseline->second.first->second->node->x,
                                         &baseline->second.second->second->node->x);
                         *out << Verbose(2) << "Plane vector of candidate triangle is ";
                         PlaneVector.Output(out);
                         *out << endl;
                         // 4. loop over all points
                         double sign = 0.;
                         PointMap::iterator checker = PointsOnBoundary.begin();
                         for (; checker != PointsOnBoundary.end(); checker++)
+                                {
                                         // (neglecting A,B,C)
                                         if ((checker == A) || (checker == baseline->second.first) || (checker
                                                         == baseline->second.second))
                                                 continue;
                                         // 4a. project onto plane vector
                                         TrialVector.CopyVector(&checker->second->node->x);
                                         TrialVector.SubtractVector(&A->second->node->x);
                                         distance = TrialVector.Projection(&PlaneVector);
                                         if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
                                                 continue;
                                         *out << Verbose(3) << "Projection of " << checker->second->node->Name
                                                         << " yields distance of " << distance << "." << endl;
                                         tmp = distance / fabs(distance);
                                         // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
                                         if ((sign != 0) && (tmp != sign))
+                                                {
                                                         // 4c. If so, break 4. loop and continue with next candidate in 1. loop
                                                         *out << Verbose(2) << "Current candidates: "
                                                                         << A->second->node->Name << ","
                                                                         << baseline->second.first->second->node->Name << ","
                                                                         << baseline->second.second->second->node->Name << " leave "
                                                                         << checker->second->node->Name << " outside the convex hull."
                                                                         << endl;
                                                         break;
+                                                }
                                         else
                                                 { // note the sign for later
                                                         *out << Verbose(2) << "Current candidates: "
                                                                         << A->second->node->Name << ","
                                                                         << baseline->second.first->second->node->Name << ","
                                                                         << baseline->second.second->second->node->Name << " leave "
                                                                         << checker->second->node->Name << " inside the convex hull."
                                                                         << endl;
                                                         sign = tmp;
+                                                }
                                         // 4d. Check whether the point is inside the triangle (check distance to each node
                                         tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
                                         int innerpoint = 0;
                                         if ((tmp < A->second->node->x.DistanceSquared(
                                                         &baseline->second.first->second->node->x)) && (tmp
                                                         < A->second->node->x.DistanceSquared(
                                                                         &baseline->second.second->second->node->x)))
                                                 innerpoint++;
                                         tmp = checker->second->node->x.DistanceSquared(
                                                         &baseline->second.first->second->node->x);
                                         if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
                                                         &A->second->node->x)) && (tmp
                                                         < baseline->second.first->second->node->x.DistanceSquared(
                                                                         &baseline->second.second->second->node->x)))
                                                 innerpoint++;
                                         tmp = checker->second->node->x.DistanceSquared(
                                                         &baseline->second.second->second->node->x);
                                         if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
                                                         &baseline->second.first->second->node->x)) && (tmp
                                                         < baseline->second.second->second->node->x.DistanceSquared(
                                                                         &A->second->node->x)))
                                                 innerpoint++;
                                         // 4e. If so, break 4. loop and continue with next candidate in 1. loop
                                         if (innerpoint == 3)
                                                 break;
+                                }
                         // 5. come this far, all on same side? Then break 1. loop and construct triangle
                         if (checker == PointsOnBoundary.end())
+                                {
                                         *out << "Looks like we have a candidate!" << endl;
                                         break;
+                                }
+                }
         if (baseline != DistanceMMap.end())
+                {
                         BPS[0] = baseline->second.first->second;
                         BPS[1] = baseline->second.second->second;
                         BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
                         BPS[0] = A->second;
                         BPS[1] = baseline->second.second->second;
                         BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
                         BPS[0] = baseline->second.first->second;
                         BPS[1] = A->second;
                         BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
                         // 4b3. insert created triangle
                         BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
                         TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
                         TrianglesOnBoundaryCount++;
                         for (int i = 0; i < NDIM; i++)
+                                {
                                         LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
                                         LinesOnBoundaryCount++;
+                                }
                         *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+                }
         else
+                {
                         *out << Verbose(1) << "No starting triangle found." << endl;
                         exit(255);
+                }
+  // 4b. create a starting triangle
+  // 4b1. create all distances
+  DistanceMultiMap DistanceMMap;
+  double distance, tmp;
+  Vector PlaneVector, TrialVector;
+  PointMap::iterator A, B, C; // three nodes of the first triangle
+  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
+  // with A chosen, take each pair B,C and sort
+  if (A != PointsOnBoundary.end())
+    {
+      B = A;
+      B++;
+      for (; B != PointsOnBoundary.end(); B++)
+        {
+          C = B;
+          C++;
+          for (; C != PointsOnBoundary.end(); C++)
+            {
+              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
+              distance = tmp * tmp;
+              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
+              distance += tmp * tmp;
+              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
+              distance += tmp * tmp;
+              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
+                  PointMap::iterator, PointMap::iterator> (B, C)));
+            }
+        }
+    }
+  //    // listing distances
+  //    *out << Verbose(1) << "Listing DistanceMMap:";
+  //    for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
+  //      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
+  //    }
+  //    *out << endl;
+  // 4b2. pick three baselines forming a triangle
+  // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+  DistanceMultiMap::iterator baseline = DistanceMMap.begin();
+  for (; baseline != DistanceMMap.end(); baseline++)
+    {
+      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+      // 2. next, we have to check whether all points reside on only one side of the triangle
+      // 3. construct plane vector
+      PlaneVector.MakeNormalVector(&A->second->node->x,
+          &baseline->second.first->second->node->x,
+          &baseline->second.second->second->node->x);
+      *out << Verbose(2) << "Plane vector of candidate triangle is ";
+      PlaneVector.Output(out);
+      *out << endl;
+      // 4. loop over all points
+      double sign = 0.;
+      PointMap::iterator checker = PointsOnBoundary.begin();
+      for (; checker != PointsOnBoundary.end(); checker++)
+        {
+          // (neglecting A,B,C)
+          if ((checker == A) || (checker == baseline->second.first) || (checker
+              == baseline->second.second))
+            continue;
+          // 4a. project onto plane vector
+          TrialVector.CopyVector(&checker->second->node->x);
+          TrialVector.SubtractVector(&A->second->node->x);
+          distance = TrialVector.Projection(&PlaneVector);
+          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+            continue;
+          *out << Verbose(3) << "Projection of " << checker->second->node->Name
+              << " yields distance of " << distance << "." << endl;
+          tmp = distance / fabs(distance);
+          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+          if ((sign != 0) && (tmp != sign))
+            {
+              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " outside the convex hull."
+                  << endl;
+              break;
+            }
+          else
+            { // note the sign for later
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " inside the convex hull."
+                  << endl;
+              sign = tmp;
+            }
+          // 4d. Check whether the point is inside the triangle (check distance to each node
+          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
+          int innerpoint = 0;
+          if ((tmp < A->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x)) && (tmp
+              < A->second->node->x.DistanceSquared(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x);
+          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
+              &A->second->node->x)) && (tmp
+              < baseline->second.first->second->node->x.DistanceSquared(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.DistanceSquared(
+              &baseline->second.second->second->node->x);
+          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
+              &baseline->second.first->second->node->x)) && (tmp
+              < baseline->second.second->second->node->x.DistanceSquared(
+                  &A->second->node->x)))
+            innerpoint++;
+          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+          if (innerpoint == 3)
+            break;
+        }
+      // 5. come this far, all on same side? Then break 1. loop and construct triangle
+      if (checker == PointsOnBoundary.end())
+        {
+          *out << "Looks like we have a candidate!" << endl;
+          break;
+        }
+    }
+  if (baseline != DistanceMMap.end())
+    {
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = A->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = A->second;
+      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      // 4b3. insert created triangle
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+      TrianglesOnBoundaryCount++;
+      for (int i = 0; i < NDIM; i++)
+        {
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+          LinesOnBoundaryCount++;
+        }
+      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
+  else
+    {
+      *out << Verbose(1) << "No starting triangle found." << endl;
+      exit(255);
+    }
+}
+;
 …
  * -# if the lines contains to only one triangle
  * -# We search all points in the boundary
  *              -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
  *              -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
  *                       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
  *              -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
+ *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
+ *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
+ *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
+ *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
  * \param *out output stream for debugging
  * \param *configuration for IsAngstroem
 …
 void
 Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
                 molecule *mol)
+{
         bool flag;
         PointMap::iterator winner;
         class BoundaryPointSet *peak = NULL;
         double SmallestAngle, TempAngle;
         Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
                         PropagationVector;
         LineMap::iterator LineChecker[2];
         do
+                {
                         flag = false;
                         for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
                                         != LinesOnBoundary.end(); baseline++)
                                 if (baseline->second->TrianglesCount == 1)
+                                        {
                                                 *out << Verbose(2) << "Current baseline is between "
                                                                 << *(baseline->second) << "." << endl;
                                                 // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
                                                 SmallestAngle = M_PI;
                                                 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
                                                 // get peak point with respect to this base line's only triangle
                                                 for (int i = 0; i < 3; i++)
                                                         if ((BTS->endpoints[i] != baseline->second->endpoints[0])
                                                                         && (BTS->endpoints[i] != baseline->second->endpoints[1]))
                                                                 peak = BTS->endpoints[i];
                                                 *out << Verbose(3) << " and has peak " << *peak << "." << endl;
                                                 // normal vector of triangle
                                                 BTS->GetNormalVector(NormalVector);
                                                 *out << Verbose(4) << "NormalVector of base triangle is ";
                                                 NormalVector.Output(out);
                                                 *out << endl;
                                                 // offset to center of triangle
                                                 CenterVector.Zero();
                                                 for (int i = 0; i < 3; i++)
                                                         CenterVector.AddVector(&BTS->endpoints[i]->node->x);
                                                 CenterVector.Scale(1. / 3.);
                                                 *out << Verbose(4) << "CenterVector of base triangle is ";
                                                 CenterVector.Output(out);
                                                 *out << endl;
                                                 // vector in propagation direction (out of triangle)
                                                 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
                                                 TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
                                                 TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
                                                 PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
                                                 TempVector.CopyVector(&CenterVector);
                                                 TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
                                                 //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
                                                 if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
                                                         PropagationVector.Scale(-1.);
                                                 *out << Verbose(4) << "PropagationVector of base triangle is ";
                                                 PropagationVector.Output(out);
                                                 *out << endl;
                                                 winner = PointsOnBoundary.end();
                                                 for (PointMap::iterator target = PointsOnBoundary.begin(); target
                                                                 != PointsOnBoundary.end(); target++)
                                                         if ((target->second != baseline->second->endpoints[0])
                                                                         && (target->second != baseline->second->endpoints[1]))
                                                                 { // don't take the same endpoints
                                                                         *out << Verbose(3) << "Target point is " << *(target->second)
                                                                                         << ":";
                                                                         bool continueflag = true;
                                                                         VirtualNormalVector.CopyVector(
                                                                                         &baseline->second->endpoints[0]->node->x);
                                                                         VirtualNormalVector.AddVector(
                                                                                         &baseline->second->endpoints[0]->node->x);
                                                                         VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
                                                                         VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
                                                                         TempAngle = VirtualNormalVector.Angle(&PropagationVector);
                                                                         continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
                                                                         if (!continueflag)
+                                                                                {
                                                                                         *out << Verbose(4)
                                                                                                         << "Angle between propagation direction and base line to "
                                                                                                         << *(target->second) << " is " << TempAngle
                                                                                                         << ", bad direction!" << endl;
                                                                                         continue;
+                                                                                }
                                                                         else
                                                                                 *out << Verbose(4)
                                                                                                 << "Angle between propagation direction and base line to "
                                                                                                 << *(target->second) << " is " << TempAngle
                                                                                                 << ", good direction!" << endl;
                                                                         LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                                                                                         target->first);
                                                                         LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                                                                                         target->first);
                                                                         //                                              if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
                                                                         //                                                      *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
                                                                         //                                              else
                                                                         //                                                      *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                                                                         //                                              if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
                                                                         //                                                      *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
                                                                         //                                              else
                                                                         //                                                      *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                                                                         // check first endpoint (if any connecting line goes to target or at least not more than 1)
                                                                         continueflag = continueflag && (((LineChecker[0]
                                                                                         == baseline->second->endpoints[0]->lines.end())
                                                                                         || (LineChecker[0]->second->TrianglesCount == 1)));
                                                                         if (!continueflag)
+                                                                                {
                                                                                         *out << Verbose(4) << *(baseline->second->endpoints[0])
                                                                                                         << " has line " << *(LineChecker[0]->second)
                                                                                                         << " to " << *(target->second)
                                                                                                         << " as endpoint with "
                                                                                                         << LineChecker[0]->second->TrianglesCount
                                                                                                         << " triangles." << endl;
                                                                                         continue;
+                                                                                }
                                                                         // check second endpoint (if any connecting line goes to target or at least not more than 1)
                                                                         continueflag = continueflag && (((LineChecker[1]
                                                                                         == baseline->second->endpoints[1]->lines.end())
                                                                                         || (LineChecker[1]->second->TrianglesCount == 1)));
                                                                         if (!continueflag)
+                                                                                {
                                                                                         *out << Verbose(4) << *(baseline->second->endpoints[1])
                                                                                                         << " has line " << *(LineChecker[1]->second)
                                                                                                         << " to " << *(target->second)
                                                                                                         << " as endpoint with "
                                                                                                         << LineChecker[1]->second->TrianglesCount
                                                                                                         << " triangles." << endl;
                                                                                         continue;
+                                                                                }
                                                                         // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
                                                                         continueflag = continueflag && (!(((LineChecker[0]
                                                                                         != baseline->second->endpoints[0]->lines.end())
                                                                                         && (LineChecker[1]
                                                                                                         != baseline->second->endpoints[1]->lines.end())
                                                                                         && (GetCommonEndpoint(LineChecker[0]->second,
                                                                                                         LineChecker[1]->second) == peak))));
                                                                         if (!continueflag)
+                                                                                {
                                                                                         *out << Verbose(4) << "Current target is peak!" << endl;
                                                                                         continue;
+                                                                                }
                                                                         // in case NOT both were found
                                                                         if (continueflag)
                                                                                 { // create virtually this triangle, get its normal vector, calculate angle
                                                                                         flag = true;
                                                                                         VirtualNormalVector.MakeNormalVector(
                                                                                                         &baseline->second->endpoints[0]->node->x,
                                                                                                         &baseline->second->endpoints[1]->node->x,
                                                                                                         &target->second->node->x);
                                                                                         // make it always point inward
                                                                                         if (baseline->second->endpoints[0]->node->x.Projection(
                                                                                                         &VirtualNormalVector) > 0)
                                                                                                 VirtualNormalVector.Scale(-1.);
                                                                                         // calculate angle
                                                                                         TempAngle = NormalVector.Angle(&VirtualNormalVector);
                                                                                         *out << Verbose(4) << "NormalVector is ";
                                                                                         VirtualNormalVector.Output(out);
                                                                                         *out << " and the angle is " << TempAngle << "." << endl;
                                                                                         if (SmallestAngle > TempAngle)
                                                                                                 { // set to new possible winner
                                                                                                         SmallestAngle = TempAngle;
                                                                                                         winner = target;
+                                                                                                }
+                                                                                }
+                                                                }
                                                 // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
                                                 if (winner != PointsOnBoundary.end())
+                                                        {
                                                                 *out << Verbose(2) << "Winning target point is "
                                                                                 << *(winner->second) << " with angle " << SmallestAngle
                                                                                 << "." << endl;
                                                                 // create the lins of not yet present
                                                                 BLS[0] = baseline->second;
                                                                 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
                                                                 LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                                                                                 winner->first);
                                                                 LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                                                                                 winner->first);
                                                                 if (LineChecker[0]
                                                                                 == baseline->second->endpoints[0]->lines.end())
                                                                         { // create
                                                                                 BPS[0] = baseline->second->endpoints[0];
                                                                                 BPS[1] = winner->second;
                                                                                 BLS[1] = new class BoundaryLineSet(BPS,
                                                                                                 LinesOnBoundaryCount);
                                                                                 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                                                                                                 BLS[1]));
                                                                                 LinesOnBoundaryCount++;
+                                                                        }
                                                                 else
                                                                         BLS[1] = LineChecker[0]->second;
                                                                 if (LineChecker[1]
                                                                                 == baseline->second->endpoints[1]->lines.end())
                                                                         { // create
                                                                                 BPS[0] = baseline->second->endpoints[1];
                                                                                 BPS[1] = winner->second;
                                                                                 BLS[2] = new class BoundaryLineSet(BPS,
                                                                                                 LinesOnBoundaryCount);
                                                                                 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                                                                                                 BLS[2]));
                                                                                 LinesOnBoundaryCount++;
+                                                                        }
                                                                 else
                                                                         BLS[2] = LineChecker[1]->second;
                                                                 BTS = new class BoundaryTriangleSet(BLS,
                                                                                 TrianglesOnBoundaryCount);
                                                                 TrianglesOnBoundary.insert(TrianglePair(
                                                                                 TrianglesOnBoundaryCount, BTS));
                                                                 TrianglesOnBoundaryCount++;
+                                                        }
                                                 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);
+    molecule *mol)
+{
+  bool flag;
+  PointMap::iterator winner;
+  class BoundaryPointSet *peak = NULL;
+  double SmallestAngle, TempAngle;
+  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
+      PropagationVector;
+  LineMap::iterator LineChecker[2];
+  do
+    {
+      flag = false;
+      for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
+          != LinesOnBoundary.end(); baseline++)
+        if (baseline->second->TrianglesCount == 1)
+          {
+            *out << Verbose(2) << "Current baseline is between "
+                << *(baseline->second) << "." << endl;
+            // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
+            SmallestAngle = M_PI;
+            BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+            // get peak point with respect to this base line's only triangle
+            for (int i = 0; i < 3; i++)
+              if ((BTS->endpoints[i] != baseline->second->endpoints[0])
+                  && (BTS->endpoints[i] != baseline->second->endpoints[1]))
+                peak = BTS->endpoints[i];
+            *out << Verbose(3) << " and has peak " << *peak << "." << endl;
+            // normal vector of triangle
+            BTS->GetNormalVector(NormalVector);
+            *out << Verbose(4) << "NormalVector of base triangle is ";
+            NormalVector.Output(out);
+            *out << endl;
+            // offset to center of triangle
+            CenterVector.Zero();
+            for (int i = 0; i < 3; i++)
+              CenterVector.AddVector(&BTS->endpoints[i]->node->x);
+            CenterVector.Scale(1. / 3.);
+            *out << Verbose(4) << "CenterVector of base triangle is ";
+            CenterVector.Output(out);
+            *out << endl;
+            // vector in propagation direction (out of triangle)
+            // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+            TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
+            PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
+            TempVector.CopyVector(&CenterVector);
+            TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+            //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
+            if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+              PropagationVector.Scale(-1.);
+            *out << Verbose(4) << "PropagationVector of base triangle is ";
+            PropagationVector.Output(out);
+            *out << endl;
+            winner = PointsOnBoundary.end();
+            for (PointMap::iterator target = PointsOnBoundary.begin(); target
+                != PointsOnBoundary.end(); target++)
+              if ((target->second != baseline->second->endpoints[0])
+                  && (target->second != baseline->second->endpoints[1]))
+                { // don't take the same endpoints
+                  *out << Verbose(3) << "Target point is " << *(target->second)
+                      << ":";
+                  bool continueflag = true;
+                  VirtualNormalVector.CopyVector(
+                      &baseline->second->endpoints[0]->node->x);
+                  VirtualNormalVector.AddVector(
+                      &baseline->second->endpoints[0]->node->x);
+                  VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
+                  VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
+                  TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+                  continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4)
+                          << "Angle between propagation direction and base line to "
+                          << *(target->second) << " is " << TempAngle
+                          << ", bad direction!" << endl;
+                      continue;
+                    }
+                  else
+                    *out << Verbose(4)
+                        << "Angle between propagation direction and base line to "
+                        << *(target->second) << " is " << TempAngle
+                        << ", good direction!" << endl;
+                  LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                      target->first);
+                  LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                      target->first);
+                  //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+                  //            else
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                  //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+                  //            else
+                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                  // check first endpoint (if any connecting line goes to target or at least not more than 1)
+                  continueflag = continueflag && (((LineChecker[0]
+                      == baseline->second->endpoints[0]->lines.end())
+                      || (LineChecker[0]->second->TrianglesCount == 1)));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << *(baseline->second->endpoints[0])
+                          << " has line " << *(LineChecker[0]->second)
+                          << " to " << *(target->second)
+                          << " as endpoint with "
+                          << LineChecker[0]->second->TrianglesCount
+                          << " triangles." << endl;
+                      continue;
+                    }
+                  // check second endpoint (if any connecting line goes to target or at least not more than 1)
+                  continueflag = continueflag && (((LineChecker[1]
+                      == baseline->second->endpoints[1]->lines.end())
+                      || (LineChecker[1]->second->TrianglesCount == 1)));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << *(baseline->second->endpoints[1])
+                          << " has line " << *(LineChecker[1]->second)
+                          << " to " << *(target->second)
+                          << " as endpoint with "
+                          << LineChecker[1]->second->TrianglesCount
+                          << " triangles." << endl;
+                      continue;
+                    }
+                  // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
+                  continueflag = continueflag && (!(((LineChecker[0]
+                      != baseline->second->endpoints[0]->lines.end())
+                      && (LineChecker[1]
+                          != baseline->second->endpoints[1]->lines.end())
+                      && (GetCommonEndpoint(LineChecker[0]->second,
+                          LineChecker[1]->second) == peak))));
+                  if (!continueflag)
+                    {
+                      *out << Verbose(4) << "Current target is peak!" << endl;
+                      continue;
+                    }
+                  // in case NOT both were found
+                  if (continueflag)
+                    { // create virtually this triangle, get its normal vector, calculate angle
+                      flag = true;
+                      VirtualNormalVector.MakeNormalVector(
+                          &baseline->second->endpoints[0]->node->x,
+                          &baseline->second->endpoints[1]->node->x,
+                          &target->second->node->x);
+                      // make it always point inward
+                      if (baseline->second->endpoints[0]->node->x.Projection(
+                          &VirtualNormalVector) > 0)
+                        VirtualNormalVector.Scale(-1.);
+                      // calculate angle
+                      TempAngle = NormalVector.Angle(&VirtualNormalVector);
+                      *out << Verbose(4) << "NormalVector is ";
+                      VirtualNormalVector.Output(out);
+                      *out << " and the angle is " << TempAngle << "." << endl;
+                      if (SmallestAngle > TempAngle)
+                        { // set to new possible winner
+                          SmallestAngle = TempAngle;
+                          winner = target;
+                        }
+                    }
+                }
+            // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
+            if (winner != PointsOnBoundary.end())
+              {
+                *out << Verbose(2) << "Winning target point is "
+                    << *(winner->second) << " with angle " << SmallestAngle
+                    << "." << endl;
+                // create the lins of not yet present
+                BLS[0] = baseline->second;
+                // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
+                LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                    winner->first);
+                LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                    winner->first);
+                if (LineChecker[0]
+                    == baseline->second->endpoints[0]->lines.end())
+                  { // create
+                    BPS[0] = baseline->second->endpoints[0];
+                    BPS[1] = winner->second;
+                    BLS[1] = new class BoundaryLineSet(BPS,
+                        LinesOnBoundaryCount);
+                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                        BLS[1]));
+                    LinesOnBoundaryCount++;
+                  }
+                else
+                  BLS[1] = LineChecker[0]->second;
+                if (LineChecker[1]
+                    == baseline->second->endpoints[1]->lines.end())
+                  { // create
+                    BPS[0] = baseline->second->endpoints[1];
+                    BPS[1] = winner->second;
+                    BLS[2] = new class BoundaryLineSet(BPS,
+                        LinesOnBoundaryCount);
+                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                        BLS[2]));
+                    LinesOnBoundaryCount++;
+                  }
+                else
+                  BLS[2] = LineChecker[1]->second;
+                BTS = new class BoundaryTriangleSet(BLS,
+                    TrianglesOnBoundaryCount);
+                TrianglesOnBoundary.insert(TrianglePair(
+                    TrianglesOnBoundaryCount, BTS));
+                TrianglesOnBoundaryCount++;
+              }
+            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);
+}
 …
 Tesselation::AddPoint(atom *Walker)
+{
         PointTestPair InsertUnique;
         BPS[0] = new class BoundaryPointSet(Walker);
         InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
         if (InsertUnique.second) // if new point was not present before, increase counter
                 PointsOnBoundaryCount++;
+  PointTestPair InsertUnique;
+  BPS[0] = new class BoundaryPointSet(Walker);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
+  if (InsertUnique.second) // if new point was not present before, increase counter
+    PointsOnBoundaryCount++;
+}
+;
 …
 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;
+                }
+  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;
+    }
+}
+;
 …
 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++;
+                }
+    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++;
+    }
+}
+;
 …
+{
         cout << Verbose(1) << "Adding triangle to its lines" << endl;
         TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
         TrianglesOnBoundaryCount++;
         /*
          * this is apparently done when constructing triangle
          for (i=0; i<3; i++)
+         {
          BLS[i]->AddTriangle(BTS);
+         }
          */
+  cout << Verbose(1) << "Adding triangle to its lines" << endl;
+  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+  TrianglesOnBoundaryCount++;
+  /*
+   * this is apparently done when constructing triangle
+   for (i=0; i<3; i++)
+   {
+   BLS[i]->AddTriangle(BTS);
+   }
+   */
+}
+;
 …
 double det_get(gsl_matrix *A, int inPlace) {
         /*
         inPlace = 1 => A is replaced with the LU decomposed copy.
         inPlace = 0 => A is retained, and a copy is used for LU.
         */
         double det;
         int signum;
         gsl_permutation *p = gsl_permutation_alloc(A->size1);
         gsl_matrix *tmpA;
         if (inPlace)
         tmpA = A;
         else {
         gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
         gsl_matrix_memcpy(tmpA , A);
+        }
         gsl_linalg_LU_decomp(tmpA , p , &signum);
         det = gsl_linalg_LU_det(tmpA , signum);
         gsl_permutation_free(p);
         if (! inPlace)
         gsl_matrix_free(tmpA);
         return det;
+  /*
+  inPlace = 1 => A is replaced with the LU decomposed copy.
+  inPlace = 0 => A is retained, and a copy is used for LU.
+  */
+  double det;
+  int signum;
+  gsl_permutation *p = gsl_permutation_alloc(A->size1);
+  gsl_matrix *tmpA;
+  if (inPlace)
+  tmpA = A;
+  else {
+  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+  gsl_matrix_memcpy(tmpA , A);
+  }
+  gsl_linalg_LU_decomp(tmpA , p , &signum);
+  det = gsl_linalg_LU_det(tmpA , signum);
+  gsl_permutation_free(p);
+  if (! inPlace)
+  gsl_matrix_free(tmpA);
+  return det;
 };
 void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
         gsl_matrix *A = gsl_matrix_calloc(3,3);
         double m11, m12, m13, m14;
         for(int i=0;i<3;i++) {
                 gsl_matrix_set(A, i, 0, a.x[i]);
                 gsl_matrix_set(A, i, 1, b.x[i]);
                 gsl_matrix_set(A, i, 2, c.x[i]);
+        }
         m11 = det_get(A, 1);
         for(int i=0;i<3;i++) {
                 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
                 gsl_matrix_set(A, i, 1, b.x[i]);
                 gsl_matrix_set(A, i, 2, c.x[i]);
+        }
         m12 = det_get(A, 1);
         for(int i=0;i<3;i++) {
                 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
                 gsl_matrix_set(A, i, 1, a.x[i]);
                 gsl_matrix_set(A, i, 2, c.x[i]);
+        }
         m13 = det_get(A, 1);
         for(int i=0;i<3;i++) {
                 gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
                 gsl_matrix_set(A, i, 1, a.x[i]);
                 gsl_matrix_set(A, i, 2, b.x[i]);
+        }
         m14 = det_get(A, 1);
         if (fabs(m11) < MYEPSILON)
                 cerr << "ERROR: three points are colinear." << endl;
         center->x[0] =  0.5 * m12/ m11;
         center->x[1] = -0.5 * m13/ m11;
         center->x[2] =  0.5 * m14/ m11;
         if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
                 cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
         gsl_matrix_free(A);
+  gsl_matrix *A = gsl_matrix_calloc(3,3);
+  double m11, m12, m13, m14;
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m11 = det_get(A, 1);
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m12 = det_get(A, 1);
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m13 = det_get(A, 1);
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, b.x[i]);
+  }
+  m14 = det_get(A, 1);
+  if (fabs(m11) < MYEPSILON)
+    cerr << "ERROR: three points are colinear." << endl;
+  center->x[0] =  0.5 * m12/ m11;
+  center->x[1] = -0.5 * m13/ m11;
+  center->x[2] =  0.5 * m14/ m11;
+  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
+    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+  gsl_matrix_free(A);
 };
 …
  */
 void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
                 double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
         Vector TempNormal, helper;
         double Restradius;
         Vector OtherCenter;
         cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
         Center->Zero();
         helper.CopyVector(&a);
         helper.Scale(sin(2.*alpha));
         Center->AddVector(&helper);
         helper.CopyVector(&b);
         helper.Scale(sin(2.*beta));
         Center->AddVector(&helper);
         helper.CopyVector(&c);
         helper.Scale(sin(2.*gamma));
         Center->AddVector(&helper);
         //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
         Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
         NewUmkreismittelpunkt->CopyVector(Center);
         cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
         // Here we calculated center of circumscribing circle, using barycentric coordinates
         cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
         TempNormal.CopyVector(&a);
         TempNormal.SubtractVector(&b);
         helper.CopyVector(&a);
         helper.SubtractVector(&c);
         TempNormal.VectorProduct(&helper);
         if (fabs(HalfplaneIndicator) < MYEPSILON)
+                {
                         if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+                                {
                                         TempNormal.Scale(-1);
+                                }
+                }
         else
+                {
                         if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+                                {
                                         TempNormal.Scale(-1);
+                                }
+                }
         TempNormal.Normalize();
         Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
         cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
         TempNormal.Scale(Restradius);
         cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
         Center->AddVector(&TempNormal);
         cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
         get_sphere(&OtherCenter, a, b, c, RADIUS);
         cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
         cout << Verbose(3) << "End of Get_center_of_sphere.\n";
+    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
+  Vector TempNormal, helper;
+  double Restradius;
+  Vector OtherCenter;
+  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+  Center->Zero();
+  helper.CopyVector(&a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(&b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(&c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+  NewUmkreismittelpunkt->CopyVector(Center);
+  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
+  // Here we calculated center of circumscribing circle, using barycentric coordinates
+  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
+  TempNormal.CopyVector(&a);
+  TempNormal.SubtractVector(&b);
+  helper.CopyVector(&a);
+  helper.SubtractVector(&c);
+  TempNormal.VectorProduct(&helper);
+  if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  else
+    {
+      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  TempNormal.Normalize();
+  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
+  TempNormal.Scale(Restradius);
+  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
+  Center->AddVector(&TempNormal);
+  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
+  get_sphere(&OtherCenter, a, b, c, RADIUS);
+  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
+  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
 };
 /** 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.
+ *  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
 …
  */
 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 alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
         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*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
                         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);
+                                }
                                 if (sign >= 0) {
                                         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);
                                         Mittelpunkt.SubtractVector(&ReferencePoint);
                                         cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
                                         BallAngle = Mittelpunkt.Angle(OldNormal);
                                         cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
                                         //cout << "direction1 is " << *direction1 << "." << endl;
                                         //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
                                         //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
                                         NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
                                         if ((Mittelpunkt.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 << " is not in search direction." << 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";
+    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 alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
+  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*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
+      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);
+        }
+        if (sign >= 0) {
+          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);
+          Mittelpunkt.SubtractVector(&ReferencePoint);
+          cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
+          BallAngle = Mittelpunkt.Angle(OldNormal);
+          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+          //cout << "direction1 is " << *direction1 << "." << endl;
+          //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
+          //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+          NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+          if ((Mittelpunkt.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 << " is not in search direction." << 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";
+}
+;
 …
 void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
         Vector helper;
         double alpha, beta, gamma;
         Vector SideA, SideB, SideC;
         SideA.CopyVector(b);
         SideA.SubtractVector(c);
         SideB.CopyVector(c);
         SideB.SubtractVector(a);
         SideC.CopyVector(a);
         SideC.SubtractVector(b);
         alpha = M_PI - SideB.Angle(&SideC);
         beta = M_PI - SideC.Angle(&SideA);
         gamma = M_PI - SideA.Angle(&SideB);
         cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
         if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
                 cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
         Center->Zero();
         helper.CopyVector(a);
         helper.Scale(sin(2.*alpha));
         Center->AddVector(&helper);
         helper.CopyVector(b);
         helper.Scale(sin(2.*beta));
         Center->AddVector(&helper);
         helper.CopyVector(c);
         helper.Scale(sin(2.*gamma));
         Center->AddVector(&helper);
         Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+  Vector helper;
+  double alpha, beta, gamma;
+  Vector SideA, SideB, SideC;
+  SideA.CopyVector(b);
+  SideA.SubtractVector(c);
+  SideB.CopyVector(c);
+  SideB.SubtractVector(a);
+  SideC.CopyVector(a);
+  SideC.SubtractVector(b);
+  alpha = M_PI - SideB.Angle(&SideC);
+  beta = M_PI - SideC.Angle(&SideA);
+  gamma = M_PI - SideA.Angle(&SideB);
+  cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
+  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
+    cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+  Center->Zero();
+  helper.CopyVector(a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
 };
 …
 double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
         Vector helper;
         double radius, alpha;
         helper.CopyVector(&NewSphereCenter);
         // test whether new center is on the parameter circle's plane
         if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
                 cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))        << "!" << endl;
                 helper.ProjectOntoPlane(&CirclePlaneNormal);
+        }
         radius = helper.ScalarProduct(&helper);
         // test whether the new center vector has length of CircleRadius
         if (fabs(radius - CircleRadius) > HULLEPSILON)
                 cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
         alpha = helper.Angle(&OldSphereCenter);
         // make the angle unique by checking the halfplanes/search direction
         if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)      // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
                 alpha = 2.*M_PI - alpha;
         cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
         radius = helper.Distance(&OldSphereCenter);
         helper.ProjectOntoPlane(&NormalVector);
         // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
         if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
                 cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
                 return alpha;
         } else {
                 cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
                 return 2.*M_PI;
+        }
+  Vector helper;
+  double radius, alpha;
+  helper.CopyVector(&NewSphereCenter);
+  // test whether new center is on the parameter circle's plane
+  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
+  radius = helper.ScalarProduct(&helper);
+  // test whether the new center vector has length of CircleRadius
+  if (fabs(radius - CircleRadius) > HULLEPSILON)
+    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+  alpha = helper.Angle(&OldSphereCenter);
+  // make the angle unique by checking the halfplanes/search direction
+  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+    alpha = 2.*M_PI - alpha;
+  cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&OldSphereCenter);
+  helper.ProjectOntoPlane(&NormalVector);
+  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
+  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
+    cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    return alpha;
+  } else {
+    cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    return 2.*M_PI;
+  }
 };
 …
 // void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
 // {
 //       atom *Walker = NULL;
 //       Vector CircleCenter;   // center of the circle, i.e. of the band of sphere's centers
 //       Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
 //       Vector OldSphereCenter;         // center of the sphere defined by the three points of BaseTriangle
 //       Vector NewSphereCenter;         // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
 //       Vector OtherNewSphereCenter;    // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
 //       Vector NewNormalVector;         // normal vector of the Candidate's triangle
 //       Vector SearchDirection;         // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
 //       Vector helper;
 //       LinkedAtoms *List = NULL;
 //       double CircleRadius; // radius of this circle
 //       double radius;
 //       double alpha, Otheralpha; // angles (i.e. parameter for the circle).
 //       double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
 //       int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 //       atom *Candidate = NULL;
+//   atom *Walker = NULL;
+//   Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+//   Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+//   Vector OldSphereCenter;   // center of the sphere defined by the three points of BaseTriangle
+//   Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+//   Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+//   Vector NewNormalVector;   // normal vector of the Candidate's triangle
+//   Vector SearchDirection;   // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
+//   Vector helper;
+//   LinkedAtoms *List = NULL;
+//   double CircleRadius; // radius of this circle
+//   double radius;
+//   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+//   double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
+//   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+//   atom *Candidate = NULL;
 //
 //       cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
+//   cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
 //
 //       cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
+//   cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
 //
 //       // construct center of circle
 //       CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
 //       CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
 //       CircleCenter.Scale(0.5);
+//   // construct center of circle
+//   CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+//   CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+//   CircleCenter.Scale(0.5);
 //
 //       // construct normal vector of circle
 //       CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
 //       CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+//   // construct normal vector of circle
+//   CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+//   CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
 //
 //       // calculate squared radius of circle
 //       radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
 //       if (radius/4. < RADIUS*RADIUS) {
 //               CircleRadius = RADIUS*RADIUS - radius/4.;
 //               CirclePlaneNormal.Normalize();
 //               cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+//   // calculate squared radius of circle
+//   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+//   if (radius/4. < RADIUS*RADIUS) {
+//     CircleRadius = RADIUS*RADIUS - radius/4.;
+//     CirclePlaneNormal.Normalize();
+//     cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 //
 //               // construct old center
 //               GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
 //               helper.CopyVector(&BaseTriangle->NormalVector);        // normal vector ensures that this is correct center of the two possible ones
 //               radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
 //               helper.Scale(sqrt(RADIUS*RADIUS - radius));
 //               OldSphereCenter.AddVector(&helper);
 //               OldSphereCenter.SubtractVector(&CircleCenter);
 //               cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+//     // construct old center
+//     GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
+//     helper.CopyVector(&BaseTriangle->NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+//     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+//     helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//     OldSphereCenter.AddVector(&helper);
+//     OldSphereCenter.SubtractVector(&CircleCenter);
+//     cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
 //
 //               // test whether old center is on the band's plane
 //               if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
 //                       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
 //                       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
 //               }
 //               radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
 //               if (fabs(radius - CircleRadius) < HULLEPSILON) {
+//     // test whether old center is on the band's plane
+//     if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+//       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+//       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+//     }
+//     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+//     if (fabs(radius - CircleRadius) < HULLEPSILON) {
 //
 //                       // construct SearchDirection
 //                       SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
 //                       helper.CopyVector(&BaseLine->endpoints[0]->node->x);
 //                       for(int i=0;i<3;i++)   // just take next different endpoint
 //                               if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
 //                                       helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
 //                               }
 //                       if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)     // ohoh, SearchDirection points inwards!
 //                               SearchDirection.Scale(-1.);
 //                       SearchDirection.ProjectOntoPlane(&OldSphereCenter);
 //                       SearchDirection.Normalize();
 //                       cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
 //                       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {     // rotated the wrong way!
 //                               cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
 //                       }
+//       // construct SearchDirection
+//       SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
+//       helper.CopyVector(&BaseLine->endpoints[0]->node->x);
+//       for(int i=0;i<3;i++)  // just take next different endpoint
+//         if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
+//           helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
+//         }
+//       if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // ohoh, SearchDirection points inwards!
+//         SearchDirection.Scale(-1.);
+//       SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+//       SearchDirection.Normalize();
+//       cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+//       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+//         cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+//       }
 //
 //                       if (LC->SetIndexToVector(&CircleCenter)) {     // get cell for the starting atom
 //                               for(int i=0;i<NDIM;i++) // store indices of this cell
 //                                       N[i] = LC->n[i];
 //                               cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
 //                       } else {
 //                               cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
 //                               return;
 //                       }
 //                       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
 //                       cout << Verbose(2) << "LC Intervals:";
 //                       for (int i=0;i<NDIM;i++) {
 //                               Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
 //                               Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
 //                               cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
 //                       }
 //                       cout << endl;
 //                       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
 //                               for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
 //                                       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
 //                                               List = LC->GetCurrentCell();
 //                                               cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
 //                                               if (List != NULL) {
 //                                                       for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
 //                                                               Candidate = (*Runner);
+//       if (LC->SetIndexToVector(&CircleCenter)) {  // get cell for the starting atom
+//         for(int i=0;i<NDIM;i++) // store indices of this cell
+//           N[i] = LC->n[i];
+//         cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+//       } else {
+//         cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+//         return;
+//       }
+//       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+//       cout << Verbose(2) << "LC Intervals:";
+//       for (int i=0;i<NDIM;i++) {
+//         Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+//         Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+//         cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+//       }
+//       cout << endl;
+//       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+//         for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+//           for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+//             List = LC->GetCurrentCell();
+//             cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+//             if (List != NULL) {
+//               for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+//                 Candidate = (*Runner);
 //
 //                                                               // check for three unique points
 //                                                               if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
 //                                                                       cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+//                 // check for three unique points
+//                 if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
+//                   cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
 //
 //                                                                       // construct both new centers
 //                                                                       GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
 //                                                                       OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+//                   // construct both new centers
+//                   GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+//                   OtherNewSphereCenter.CopyVector(&NewSphereCenter);
 //
 //                                                                       if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
 //                                                                               helper.CopyVector(&NewNormalVector);
 //                                                                               cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
 //                                                                               radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
 //                                                                               if (radius < RADIUS*RADIUS) {
 //                                                                                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
 //                                                                                       cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
 //                                                                                       NewSphereCenter.AddVector(&helper);
 //                                                                                       NewSphereCenter.SubtractVector(&CircleCenter);
 //                                                                                       cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+//                   if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+//                     helper.CopyVector(&NewNormalVector);
+//                     cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+//                     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+//                     if (radius < RADIUS*RADIUS) {
+//                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//                       cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+//                       NewSphereCenter.AddVector(&helper);
+//                       NewSphereCenter.SubtractVector(&CircleCenter);
+//                       cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
 //
 //                                                                                       helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
 //                                                                                       OtherNewSphereCenter.AddVector(&helper);
 //                                                                                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
 //                                                                                       cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+//                       helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+//                       OtherNewSphereCenter.AddVector(&helper);
+//                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
+//                       cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 //
 //                                                                                       // check both possible centers
 //                                                                                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
 //                                                                                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
 //                                                                                       alpha = min(alpha, Otheralpha);
 //                                                                                       if (*ShortestAngle > alpha) {
 //                                                                                                       OptCandidate = Candidate;
 //                                                                                                       *ShortestAngle = alpha;
 //                                                                                                       if (alpha != Otheralpha)
 //                                                                                                               OptCandidateCenter->CopyVector(&NewSphereCenter);
 //                                                                                                       else
 //                                                                                                               OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
 //                                                                                                       cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
 //                                                                                       } else {
 //                                                                                               if (OptCandidate != NULL)
 //                                                                                                       cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
 //                                                                                               else
 //                                                                                                       cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
 //                                                                                       }
+//                       // check both possible centers
+//                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                       alpha = min(alpha, Otheralpha);
+//                       if (*ShortestAngle > alpha) {
+//                           OptCandidate = Candidate;
+//                           *ShortestAngle = alpha;
+//                           if (alpha != Otheralpha)
+//                             OptCandidateCenter->CopyVector(&NewSphereCenter);
+//                           else
+//                             OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+//                           cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+//                       } else {
+//                         if (OptCandidate != NULL)
+//                           cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+//                         else
+//                           cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+//                       }
 //
 //                                                                               } else {
 //                                                                                       cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
 //                                                                               }
 //                                                                       } else {
 //                                                                               cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
 //                                                                       }
 //                                                               } else {
 //                                                                       cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
 //                                                               }
 //                                                       }
 //                                               }
 //                                       }
 //               } else {
 //                       cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
 //               }
 //       } else {
 //               cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
 //       }
+//                     } else {
+//                       cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+//                     }
+//                   } else {
+//                     cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+//                   }
+//                 } else {
+//                   cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
+//                 }
+//               }
+//             }
+//           }
+//     } else {
+//       cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+//     }
+//   } else {
+//     cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
+//   }
 //
 //       cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
+//   cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
 // };
 …
  */
 bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
         LineMap::iterator FindLine;
         PointMap::iterator FindPoint;
         *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
         for (int i=0;i<3;i++) { // check through all endpoints
                 FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
                 if (FindPoint != PointsOnBoundary.end())
                         TPS[i] = FindPoint->second;
                 else
                         TPS[i] = NULL;
+        }
         // check lines
         for (int i=0;i<3;i++)
                 if (TPS[i] != NULL)
                         for (int j=i;j<3;j++)
                                 if (TPS[j] != NULL) {
                                         FindLine = TPS[i]->lines.find(TPS[j]->node->nr);
                                         if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) {
                                                 *out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->second->TrianglesCount << " triangles attached." << endl;
                                                 *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
                                                 return false;
+                                        }
+                                }
         *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
         return true;
+  LineMap::iterator FindLine;
+  PointMap::iterator FindPoint;
+  *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+  for (int i=0;i<3;i++) { // check through all endpoints
+    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end())
+      TPS[i] = FindPoint->second;
+    else
+      TPS[i] = NULL;
+  }
+  // check lines
+  for (int i=0;i<3;i++)
+    if (TPS[i] != NULL)
+      for (int j=i;j<3;j++)
+        if (TPS[j] != NULL) {
+          FindLine = TPS[i]->lines.find(TPS[j]->node->nr);
+          if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) {
+            *out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->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_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+{
         Vector CircleCenter;    // center of the circle, i.e. of the band of sphere's centers
         Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
         Vector NewSphereCenter;  // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
         Vector OtherNewSphereCenter;     // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
         Vector NewNormalVector;  // normal vector of the Candidate's triangle
         Vector helper;
         LinkedAtoms *List = NULL;
         double CircleRadius; // radius of this circle
         double radius;
         double alpha, Otheralpha; // angles (i.e. parameter for the circle).
         int N[NDIM], Nlower[NDIM], Nupper[NDIM];
         atom *Candidate = NULL;
         cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
         cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
         // construct center of circle
         CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
         CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
         CircleCenter.Scale(0.5);
         // construct normal vector of circle
         CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
         CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
         // calculate squared radius of circle
         radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
         if (radius/4. < RADIUS*RADIUS) {
                 CircleRadius = RADIUS*RADIUS - radius/4.;
                 CirclePlaneNormal.Normalize();
                 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
                 // test whether old center is on the band's plane
                 if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
                         cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
                         OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+                }
                 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
                 if (fabs(radius - CircleRadius) < HULLEPSILON) {
                         // check SearchDirection
                         cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
                         if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {      // rotated the wrong way!
                                 cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+                        }
                         // get cell for the starting atom
                         if (LC->SetIndexToVector(&CircleCenter)) {
                                         for(int i=0;i<NDIM;i++) // store indices of this cell
                                         N[i] = LC->n[i];
                                 cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
                         } else {
                                 cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
                                 return;
+                        }
                         // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
                         cout << Verbose(2) << "LC Intervals:";
                         for (int i=0;i<NDIM;i++) {
                                 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
                                 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
                                 cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+                        }
                         cout << endl;
                         for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
                                 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
                                         for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
                                                 List = LC->GetCurrentCell();
                                                 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
                                                 if (List != NULL) {
                                                         for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                                                                 Candidate = (*Runner);
                                                                 // check for three unique points
                                                                 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
                                                                         cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
                                                                         // construct both new centers
                                                                         GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
                                                                         OtherNewSphereCenter.CopyVector(&NewSphereCenter);
                                                                         if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
                                                                                 helper.CopyVector(&NewNormalVector);
                                                                                 cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                                                                                 radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
                                                                                 if (radius < RADIUS*RADIUS) {
                                                                                         helper.Scale(sqrt(RADIUS*RADIUS - radius));
                                                                                         cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
                                                                                         NewSphereCenter.AddVector(&helper);
                                                                                         NewSphereCenter.SubtractVector(&CircleCenter);
                                                                                         cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
                                                                                         helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
                                                                                         OtherNewSphereCenter.AddVector(&helper);
                                                                                         OtherNewSphereCenter.SubtractVector(&CircleCenter);
                                                                                         cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
                                                                                         // check both possible centers
                                                                                         alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                                                                                         Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                                                                                         alpha = min(alpha, Otheralpha);
                                                                                         if (*ShortestAngle > alpha) {
                                                                                                         OptCandidate = Candidate;
                                                                                                         *ShortestAngle = alpha;
                                                                                                         if (alpha != Otheralpha)
                                                                                                                 OptCandidateCenter->CopyVector(&NewSphereCenter);
                                                                                                         else
                                                                                                                 OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
                                                                                                         cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
                                                                                         } else {
                                                                                                 if (OptCandidate != NULL)
                                                                                                         cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
                                                                                                 else
                                                                                                         cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                                                                                        }
                                                                                 } else {
                                                                                         cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                                                                                }
                                                                         } else {
                                                                                 cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                                                                        }
                                                                 } else {
                                                                         if (ThirdNode != NULL)
                                                                                 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
                                                                         else
                                                                                 cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                                                                }
+                                                        }
+                                                }
+                                        }
                 } else {
                         cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+                }
         } else {
                 if (ThirdNode != NULL)
                         cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
                 else
                         cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+        }
         cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  Vector NewNormalVector;   // normal vector of the Candidate's triangle
+  Vector helper;
+  LinkedAtoms *List = NULL;
+  double CircleRadius; // radius of this circle
+  double radius;
+  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  atom *Candidate = NULL;
+  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
+  cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+  // construct center of circle
+  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+  // calculate squared radius of circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    // test whether old center is on the band's plane
+    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(radius - CircleRadius) < HULLEPSILON) {
+      // check SearchDirection
+      cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+      }
+      // get cell for the starting atom
+      if (LC->SetIndexToVector(&CircleCenter)) {
+          for(int i=0;i<NDIM;i++) // store indices of this cell
+          N[i] = LC->n[i];
+        cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+      } else {
+        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+        return;
+      }
+      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+      cout << Verbose(2) << "LC Intervals:";
+      for (int i=0;i<NDIM;i++) {
+        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+        cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+      }
+      cout << endl;
+      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+            List = LC->GetCurrentCell();
+            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+            if (List != NULL) {
+              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                Candidate = (*Runner);
+                // check for three unique points
+                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
+                  cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+                  // construct both new centers
+                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+                    helper.CopyVector(&NewNormalVector);
+                    cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+                    if (radius < RADIUS*RADIUS) {
+                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+                      NewSphereCenter.AddVector(&helper);
+                      NewSphereCenter.SubtractVector(&CircleCenter);
+                      cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+                      helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+                      OtherNewSphereCenter.AddVector(&helper);
+                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
+                      cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+                      // check both possible centers
+                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      alpha = min(alpha, Otheralpha);
+                      if (*ShortestAngle > alpha) {
+                          OptCandidate = Candidate;
+                          *ShortestAngle = alpha;
+                          if (alpha != Otheralpha)
+                            OptCandidateCenter->CopyVector(&NewSphereCenter);
+                          else
+                            OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+                          cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+                      } else {
+                        if (OptCandidate != NULL)
+                          cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+                        else
+                          cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                      }
+                    } else {
+                      cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                    }
+                  } else {
+                    cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                  }
+                } else {
+                  if (ThirdNode != NULL)
+                    cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                  else
+                    cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                }
+              }
+            }
+          }
+    } else {
+      cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+    }
+  } else {
+    if (ThirdNode != NULL)
+      cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+    else
+      cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+  }
+  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
 };
 …
 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
         cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
         Vector AngleCheck;
         double norm = -1., angle;
         LinkedAtoms *List = NULL;
         int N[NDIM], Nlower[NDIM], Nupper[NDIM];
         if (LC->SetIndexToAtom(a)) {    // get cell for the starting atom
                 for(int i=0;i<NDIM;i++) // store indices of this cell
                         N[i] = LC->n[i];
         } else {
                 cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
                 return;
+        }
         // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
         cout << Verbose(2) << "LC Intervals:";
         for (int i=0;i<NDIM;i++) {
                 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
                 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
                 cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+        }
         cout << endl;
         for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
                 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
                         for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
                                 List = LC->GetCurrentCell();
                                 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
                                 if (List != NULL) {
                                         for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                                                 Candidate = (*Runner);
                                                                 // check 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;
+                                                        }
+                                                }
+                                        }
+                                }
+                        }
         cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
+  Vector AngleCheck;
+  double norm = -1., angle;
+  LinkedAtoms *List = NULL;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
+    for(int i=0;i<NDIM;i++) // store indices of this cell
+      N[i] = LC->n[i];
+  } else {
+    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+    return;
+  }
+  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+  cout << Verbose(2) << "LC Intervals:";
+  for (int i=0;i<NDIM;i++) {
+    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+  }
+  cout << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            Candidate = (*Runner);
+                // check 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;
+              }
+            }
+          }
+        }
+      }
+  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
 };
 …
 void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
+{
         cout << Verbose(1) << "Begin of Find_starting_triangle\n";
         int i = 0;
         LinkedAtoms *List = NULL;
         atom* FirstPoint;
         atom* SecondPoint;
         atom* MaxAtom[NDIM];
         double max_coordinate[NDIM];
         Vector Oben;
         Vector helper;
         Vector Chord;
         Vector SearchDirection;
         Vector OptCandidateCenter;
         Oben.Zero();
         for (i = 0; i < 3; i++) {
                 MaxAtom[i] = NULL;
                 max_coordinate[i] = -1;
+        }
         // 1. searching topmost atom with respect to each axis
         for (int i=0;i<NDIM;i++) { // each axis
                 LC->n[i] = LC->N[i]-1; // current axis is topmost cell
                 for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
                         for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
                                 List = LC->GetCurrentCell();
                                 //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
                                 if (List != NULL) {
                                         for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
                                                 cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
                                                 if ((*Runner)->x.x[i] > max_coordinate[i]) {
                                                         max_coordinate[i] = (*Runner)->x.x[i];
                                                         MaxAtom[i] = (*Runner);
+                                                }
+                                        }
                                 } else {
                                         cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+                                }
+                        }
+        }
         cout << Verbose(2) << "Found maximum coordinates: ";
         for (int i=0;i<NDIM;i++)
                 cout << i << ": " << *MaxAtom[i] << "\t";
         cout << endl;
         const int k = 1;                // arbitrary choice
         Oben.x[k] = 1.;
         FirstPoint = MaxAtom[k];
         cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
         // add first point
         AddTrianglePoint(FirstPoint, 0);
         double ShortestAngle;
         atom* Opt_Candidate = NULL;
         ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
         Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
         SecondPoint = Opt_Candidate;
         cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
         // add second point and first baseline
         AddTrianglePoint(SecondPoint, 1);
         AddTriangleLine(TPS[0], TPS[1], 0);
         helper.CopyVector(&(FirstPoint->x));
         helper.SubtractVector(&(SecondPoint->x));
         helper.Normalize();
         Oben.ProjectOntoPlane(&helper);
         Oben.Normalize();
         helper.VectorProduct(&Oben);
         ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
         Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
         Chord.SubtractVector(&(SecondPoint->x));
         double radius = Chord.ScalarProduct(&Chord);
         double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
         helper.CopyVector(&Oben);
         helper.Scale(CircleRadius);
         // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
         cout << Verbose(2) << "Looking for third point candidates \n";
         // look in one direction of baseline for initial candidate
         Opt_Candidate = NULL;
         SearchDirection.MakeNormalVector(&Chord, &Oben);        // whether we look "left" first or "right" first is not important ...
         cout << Verbose(1) << "Looking for third point candidates ...\n";
         Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
         cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
         // add third point
         AddTrianglePoint(Opt_Candidate, 2);
         // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
         // Finally, we only have to add the found further lines
         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)
         OptCandidateCenter.Scale(-1.);
         cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
         BTS->GetNormalVector(OptCandidateCenter);
         cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
         cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
         cout << Verbose(1) << "End of Find_starting_triangle\n";
+  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  int i = 0;
+  LinkedAtoms *List = NULL;
+  atom* FirstPoint;
+  atom* SecondPoint;
+  atom* MaxAtom[NDIM];
+  double max_coordinate[NDIM];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector SearchDirection;
+  Vector OptCandidateCenter;
+  Oben.Zero();
+  for (i = 0; i < 3; i++) {
+    MaxAtom[i] = NULL;
+    max_coordinate[i] = -1;
+  }
+  // 1. searching topmost atom with respect to each axis
+  for (int i=0;i<NDIM;i++) { // each axis
+    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+            cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
+            if ((*Runner)->x.x[i] > max_coordinate[i]) {
+              max_coordinate[i] = (*Runner)->x.x[i];
+              MaxAtom[i] = (*Runner);
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  }
+  cout << Verbose(2) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    cout << i << ": " << *MaxAtom[i] << "\t";
+  cout << endl;
+  const int k = 1;    // arbitrary choice
+  Oben.x[k] = 1.;
+  FirstPoint = MaxAtom[k];
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
+  // add first point
+  AddTrianglePoint(FirstPoint, 0);
+  double ShortestAngle;
+  atom* Opt_Candidate = NULL;
+  ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+  SecondPoint = Opt_Candidate;
+  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+  // add second point and first baseline
+  AddTrianglePoint(SecondPoint, 1);
+  AddTriangleLine(TPS[0], TPS[1], 0);
+  helper.CopyVector(&(FirstPoint->x));
+  helper.SubtractVector(&(SecondPoint->x));
+  helper.Normalize();
+  Oben.ProjectOntoPlane(&helper);
+  Oben.Normalize();
+  helper.VectorProduct(&Oben);
+  ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+  Chord.SubtractVector(&(SecondPoint->x));
+  double radius = Chord.ScalarProduct(&Chord);
+  double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+  helper.CopyVector(&Oben);
+  helper.Scale(CircleRadius);
+  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+  cout << Verbose(2) << "Looking for third point candidates \n";
+  // look in one direction of baseline for initial candidate
+  Opt_Candidate = NULL;
+  SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+  cout << Verbose(1) << "Looking for third point candidates ...\n";
+  Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
+  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+  // add third point
+  AddTrianglePoint(Opt_Candidate, 2);
+  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+  // Finally, we only have to add the found further lines
+  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)
+  OptCandidateCenter.Scale(-1.);
+  cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
+  BTS->GetNormalVector(OptCandidateCenter);
+  cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+  cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+  cout << Verbose(1) << "End of Find_starting_triangle\n";
 };
 …
  */
 bool Tesselation::Find_next_suitable_triangle(ofstream *out,
                 molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
                 const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
         cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
         ofstream *tempstream = NULL;
         char NumberName[255];
         atom* Opt_Candidate = NULL;
         Vector OptCandidateCenter;
         Vector CircleCenter;
         Vector CirclePlaneNormal;
         Vector OldSphereCenter;
         Vector SearchDirection;
         Vector helper;
         atom *ThirdNode = NULL;
         double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
         double radius, CircleRadius;
         cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
         for (int i=0;i<3;i++)
                 if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
                         ThirdNode = T.endpoints[i]->node;
         // construct center of circle
         CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
         CircleCenter.AddVector(&Line.endpoints[1]->node->x);
         CircleCenter.Scale(0.5);
         // construct normal vector of circle
         CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
         CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
         // calculate squared radius of circle
         radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
         if (radius/4. < RADIUS*RADIUS) {
                 CircleRadius = RADIUS*RADIUS - radius/4.;
                 CirclePlaneNormal.Normalize();
                 cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
                 // construct old center
                 GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
                 helper.CopyVector(&T.NormalVector);     // normal vector ensures that this is correct center of the two possible ones
                 radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
                 helper.Scale(sqrt(RADIUS*RADIUS - radius));
                 OldSphereCenter.AddVector(&helper);
                 OldSphereCenter.SubtractVector(&CircleCenter);
                 cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
                 // construct SearchDirection
                 SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
                 helper.CopyVector(&Line.endpoints[0]->node->x);
                 helper.SubtractVector(&ThirdNode->x);
                 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)      // ohoh, SearchDirection points inwards!
                         SearchDirection.Scale(-1.);
                 SearchDirection.ProjectOntoPlane(&OldSphereCenter);
                 SearchDirection.Normalize();
                 cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
                 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {      // rotated the wrong way!
                         cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+                }
                 // add third point
                 cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
                 Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
         } else {
                 cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+        }
         if (Opt_Candidate == NULL) {
                 cerr << "WARNING: Could not find a suitable candidate." << endl;
                 return false;
+        }
         cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl;
         // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
         atom *AtomCandidates[3];
         AtomCandidates[0] = Opt_Candidate;
         AtomCandidates[1] = Line.endpoints[0]->node;
         AtomCandidates[2] = Line.endpoints[1]->node;
         bool flag = CheckPresenceOfTriangle(out, AtomCandidates);
         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();
                 OptCandidateCenter.Scale(-1.);
                 BTS->GetNormalVector(OptCandidateCenter);
                 OptCandidateCenter.Scale(-1.);
                 cout << "--> 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 (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
                 sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
                 if (DoTecplotOutput) {
                         string NameofTempFile(tempbasename);
                         NameofTempFile.append(NumberName);
                         for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
                                 NameofTempFile.erase(npos, 1);
                         NameofTempFile.append(TecplotSuffix);
                         cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
                         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
                         write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
                         tempstream->close();
                         tempstream->flush();
                         delete(tempstream);
+                }
                 if (DoRaster3DOutput) {
                         string NameofTempFile(tempbasename);
                         NameofTempFile.append(NumberName);
                         for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
                                 NameofTempFile.erase(npos, 1);
                         NameofTempFile.append(Raster3DSuffix);
                         cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
                         tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
                         write_raster3d_file(out, tempstream, this, mol);
                         // include the current position of the virtual sphere in the temporary raster3d file
                         // make the circumsphere's center absolute again
                         helper.CopyVector(&Line.endpoints[0]->node->x);
                         helper.AddVector(&Line.endpoints[1]->node->x);
                         helper.Scale(0.5);
                         OptCandidateCenter.AddVector(&helper);
                         Vector *center = mol->DetermineCenterOfAll(out);
                         OptCandidateCenter.AddVector(center);
                         delete(center);
                         // and add to file plus translucency object
                         *tempstream << "# current virtual sphere\n";
                         *tempstream << "8\n     25.0            0.6              -1.0 -1.0 -1.0          0.2                            0 0 0 0\n";
                         *tempstream << "2\n     " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
                         *tempstream << "9\n     terminating special property\n";
                         tempstream->close();
                         tempstream->flush();
                         delete(tempstream);
+                }
                 if (DoTecplotOutput || DoRaster3DOutput)
                         TriangleFilesWritten++;
+        }
         cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
         return true;
+    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
+  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  atom* Opt_Candidate = NULL;
+  Vector OptCandidateCenter;
+  Vector CircleCenter;
+  Vector CirclePlaneNormal;
+  Vector OldSphereCenter;
+  Vector SearchDirection;
+  Vector helper;
+  atom *ThirdNode = NULL;
+  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+  double radius, CircleRadius;
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+  for (int i=0;i<3;i++)
+    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+      ThirdNode = T.endpoints[i]->node;
+  // construct center of circle
+  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
+  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
+  // calculate squared radius of circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    // construct old center
+    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
+    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+    helper.Scale(sqrt(RADIUS*RADIUS - radius));
+    OldSphereCenter.AddVector(&helper);
+    OldSphereCenter.SubtractVector(&CircleCenter);
+    cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+    // construct SearchDirection
+    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
+    helper.CopyVector(&Line.endpoints[0]->node->x);
+    helper.SubtractVector(&ThirdNode->x);
+    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // ohoh, SearchDirection points inwards!
+      SearchDirection.Scale(-1.);
+    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+    SearchDirection.Normalize();
+    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+    }
+    // add third point
+    cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+    Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
+  } else {
+    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+  }
+  if (Opt_Candidate == NULL) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl;
+  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+  atom *AtomCandidates[3];
+  AtomCandidates[0] = Opt_Candidate;
+  AtomCandidates[1] = Line.endpoints[0]->node;
+  AtomCandidates[2] = Line.endpoints[1]->node;
+  bool flag = CheckPresenceOfTriangle(out, AtomCandidates);
+  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();
+    OptCandidateCenter.Scale(-1.);
+    BTS->GetNormalVector(OptCandidateCenter);
+    OptCandidateCenter.Scale(-1.);
+    cout << "--> 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 (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+    sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(TecplotSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoRaster3DOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+      NameofTempFile.append(Raster3DSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_raster3d_file(out, tempstream, this, mol);
+      // include the current position of the virtual sphere in the temporary raster3d file
+      // make the circumsphere's center absolute again
+      helper.CopyVector(&Line.endpoints[0]->node->x);
+      helper.AddVector(&Line.endpoints[1]->node->x);
+      helper.Scale(0.5);
+      OptCandidateCenter.AddVector(&helper);
+      Vector *center = mol->DetermineCenterOfAll(out);
+      OptCandidateCenter.AddVector(center);
+      delete(center);
+      // and add to file plus translucency object
+      *tempstream << "# current virtual sphere\n";
+      *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+      *tempstream << "2\n  " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
+      *tempstream << "9\n  terminating special property\n";
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
+  }
+  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+  return true;
 };
 …
 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
+{
         int N = 0;
         bool freeTess = false;
         *out << Verbose(1) << "Entering search for non convex hull. " << endl;
         if (Tess == NULL) {
                 *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
                 Tess = new Tesselation;
                 freeTess = true;
+        }
         bool freeLC = false;
         LineMap::iterator baseline;
         *out << 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 (LCList == NULL) {
                 LCList = new LinkedCell(mol, 2.*RADIUS);
                 freeLC = true;
+        }
         Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
         baseline = Tess->LinesOnBoundary.begin();
         while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
                 if (baseline->second->TrianglesCount == 1) {
                         failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
                         flag = flag || failflag;
                         if (!failflag)
                                 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 (1) { //failflag) {
                 *out << Verbose(1) << "Writing final tecplot file\n";
                 if (DoTecplotOutput) {
                         string OutputName(filename);
                         OutputName.append(TecplotSuffix);
                         ofstream *tecplot = new ofstream(OutputName.c_str());
                         write_tecplot_file(out, tecplot, Tess, mol, -1);
                         tecplot->close();
                         delete(tecplot);
+                }
                 if (DoRaster3DOutput) {
                         string OutputName(filename);
                         OutputName.append(Raster3DSuffix);
                         ofstream *raster = new ofstream(OutputName.c_str());
                         write_raster3d_file(out, raster, Tess, mol);
                         raster->close();
                         delete(raster);
+                }
         } else {
                 cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+        }
         if (freeTess)
                 delete(Tess);
         if (freeLC)
                 delete(LCList);
         *out << Verbose(0) << "End of Find_non_convex_border\n";
+  int N = 0;
+  bool freeTess = false;
+  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
+  if (Tess == NULL) {
+    *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
+    Tess = new Tesselation;
+    freeTess = true;
+  }
+  bool freeLC = false;
+  LineMap::iterator baseline;
+  *out << 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 (LCList == NULL) {
+    LCList = new LinkedCell(mol, 2.*RADIUS);
+    freeLC = true;
+  }
+  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
+  baseline = Tess->LinesOnBoundary.begin();
+  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+    if (baseline->second->TrianglesCount == 1) {
+      failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
+      flag = flag || failflag;
+      if (!failflag)
+        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 (1) { //failflag) {
+    *out << Verbose(1) << "Writing final tecplot file\n";
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      write_tecplot_file(out, tecplot, Tess, mol, -1);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(Raster3DSuffix);
+      ofstream *raster = new ofstream(OutputName.c_str());
+      write_raster3d_file(out, raster, Tess, mol);
+      raster->close();
+      delete(raster);
+    }
+  } else {
+    cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+  }
+  if (freeTess)
+    delete(Tess);
+  if (freeLC)
+    delete(LCList);
+  *out << Verbose(0) << "End of Find_non_convex_border\n";
 };

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 437922 for src/boundary.cpp

Legend:

src/boundary.cpp

Download in other formats: