Changeset e4ea46

Timestamp:

Dec 16, 2008, 6:39:28 PM (16 years ago)

Author:

Christian Neuen <neuen@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

10af0d

Parents:

caf5d6

Message:

Tesselation starts to look good, minor discrepancies are still there and a segmentation fault.

Files:

• configure.ac (modified) (1 diff)
• src/boundary.cpp (modified) (30 diffs)
• src/boundary.hpp (modified) (3 diffs)
• src/vector.cpp (modified) (2 diffs)
• src/vector.hpp (modified) (4 diffs)

configure.ac

@@ -14 +14 @@
 AC_PROG_CC
 AM_MISSING_PROG([DOXYGEN], [doxygen])
-
-AC_ARG_ENABLE([debug],AS_HELP_STRING([--enable-debug],[debugging level of compiler. Argument is yes or debugging level. (default is no)]),
-              [enable_debugging=$enableval], [enable_debugging=no])
-AC_ARG_ENABLE([optimization],AS_HELP_STRING([--enable-optimization],[Optimization level of compiler. Argument is yes or optimization. (default is 2)]),
-              [enable_optimization=$enableval], [enable_optimization=2])
-AC_ARG_ENABLE([warnings], AS_HELP_STRING([--enable-warnings],[Output compiler warnings, argument is none, some or full (default is some).]),
-              [enable_warnings=$enableval], [enable_warnings=some])
-AC_SET_COMPILER_FLAGS([$enable_optimization], [$enable_debugging], [$enable_warnings])
 
 # Checks for libraries.

src/boundary.cpp

@@ -2 +2 @@
 #include "boundary.hpp"
 
-
-
+#define DEBUG 0
 
 // ======================================== Points on Boundary =================================
@@ -11 +10 @@
   LinesCount = 0;
   Nr = -1;
-};
+}
+;
 
 BoundaryPointSet::BoundaryPointSet(atom *Walker)
@@ -18 +18 @@
   LinesCount = 0;
   Nr = Walker->nr;
-};
+}
+;
 
 BoundaryPointSet::~BoundaryPointSet()
@@ -24 +25 @@
   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   node = NULL;
-};
-
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
-{
-  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
-  if (line->endpoints[0] == this) {
-    lines.insert ( LinePair( line->endpoints[1]->Nr, line) );
-  } else {
-    lines.insert ( LinePair( line->endpoints[0]->Nr, line) );
-  }
+}
+;
+
+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++;
-};
-
-ostream & operator << (ostream &ost, BoundaryPointSet &a)
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryPointSet &a)
 {
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
   return ost;
-};
+}
+;
 
 // ======================================== Lines on Boundary =================================
@@ -47 +57 @@
 BoundaryLineSet::BoundaryLineSet()
 {
-  for (int i=0;i<2;i++)
+  for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
   TrianglesCount = 0;
   Nr = -1;
-};
+}
+;
 
 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
@@ -60 +71 @@
   SetEndpointsOrdered(endpoints, Point[0], Point[1]);
   // add this line to the hash maps of both endpoints
-  Point[0]->AddLine(this);
-  Point[1]->AddLine(this);
+  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()
 {
-  for (int i=0;i<2;i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-    endpoints[i]->lines.erase(Nr);
-    LineMap::iterator tester = endpoints[i]->lines.begin();
-    tester++;
-    if (tester == endpoints[i]->lines.end()) {
-      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-      delete(endpoints[i]);
-    } else
-      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
-  }
-};
-
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
-{
-  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
-  triangles.insert ( TrianglePair( TrianglesCount, triangle) );
+  for (int i = 0; i < 2; i++)
+    {
+      cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point "
+          << *endpoints[i] << "." << endl;
+      endpoints[i]->lines.erase(Nr);
+      LineMap::iterator tester = endpoints[i]->lines.begin();
+      tester++;
+      if (tester == endpoints[i]->lines.end())
+        {
+          cout << Verbose(5) << *endpoints[i]
+              << " has no more lines it's attached to, erasing." << endl;
+          //delete(endpoints[i]);
+        }
+      else
+        cout << Verbose(5) << *endpoints[i]
+            << " has still lines it's attached to." << endl;
+    }
+}
+;
+
+void
+BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+      << endl;
+  triangles.insert(TrianglePair(TrianglesCount, triangle));
   TrianglesCount++;
-};
-
-ostream & operator << (ostream &ost, BoundaryLineSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryLineSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "]";
   return ost;
-};
+}
+;
 
 // ======================================== Triangles on Boundary =================================
@@ -100 +125 @@
 BoundaryTriangleSet::BoundaryTriangleSet()
 {
-  for (int i=0;i<3;i++) {
-    endpoints[i] = NULL;
-    lines[i] = NULL;
-  }
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
   Nr = -1;
-};
-
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+}
+;
+
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3],
+    int number)
 {
   // set number
@@ -113 +141 @@
   // set lines
   cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
-  for (int i=0;i<3;i++) {
-    lines[i] = line[i];
-    lines[i]->AddTriangle(this);
-  }
+  for (int i = 0; i < 3; i++)
+    {
+      lines[i] = line[i];
+      lines[i]->AddTriangle(this);
+    }
   // get ascending order of endpoints
-  map <int, class BoundaryPointSet * > OrderMap;
-  for(int i=0;i<3;i++)  // for all three lines
-    for (int j=0;j<2;j++) { // for both endpoints
-      OrderMap.insert ( pair <int, class BoundaryPointSet * >( line[i]->endpoints[j]->Nr, line[i]->endpoints[j]) );
-      // and we don't care whether insertion fails
-    }
+  map<int, class BoundaryPointSet *> OrderMap;
+  for (int i = 0; i < 3; i++)
+    // for all three lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+        // and we don't care whether insertion fails
+      }
   // set endpoints
   int Counter = 0;
   cout << Verbose(6) << " with end points ";
-  for (map <int, class BoundaryPointSet * >::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) {
-    endpoints[Counter] = runner->second;
-    cout << " " << *endpoints[Counter];
-    Counter++;
-  }
-  if (Counter < 3) {
-    cerr << "ERROR! We have a triangle with only two distinct endpoints!" << endl;
-    //exit(1);
-  }
+  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+      != OrderMap.end(); runner++)
+    {
+      endpoints[Counter] = runner->second;
+      cout << " " << *endpoints[Counter];
+      Counter++;
+    }
+  if (Counter < 3)
+    {
+      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+          << endl;
+      //exit(1);
+    }
   cout << "." << endl;
-};
+}
+;
 
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-  for (int i=0;i<3;i++) {
-    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-    lines[i]->triangles.erase(Nr);
-    TriangleMap::iterator tester = lines[i]->triangles.begin();
-    tester++;
-    if (tester == lines[i]->triangles.end()) {
-      cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-      delete(lines[i]);
-    } else
-      cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
-  }
-};
-
-void BoundaryTriangleSet::GetNormalVector(Vector &NormalVector)
+  for (int i = 0; i < 3; i++)
+    {
+      cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+      lines[i]->triangles.erase(Nr);
+      TriangleMap::iterator tester = lines[i]->triangles.begin();
+      tester++;
+      if (tester == lines[i]->triangles.end())
+        {
+          cout << Verbose(5) << *lines[i]
+              << " is no more attached to any triangle, erasing." << endl;
+          delete (lines[i]);
+        }
+      else
+        cout << Verbose(5) << *lines[i] << " is still attached to a triangle."
+            << endl;
+    }
+}
+;
+
+void
+BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
   // get normal vector
-  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x, &endpoints[2]->node->x);
+  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
+      &endpoints[2]->node->x);
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (endpoints[0]->node->x.Projection(&NormalVector) > 0)
+  if (endpoints[0]->node->x.Projection(&OtherVector) > 0)
     NormalVector.Scale(-1.);
-};
-
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
   return ost;
-};
+}
+;
 
 // ========================================== F U N C T I O N S =================================
@@ -177 +226 @@
  * \return point which is shared or NULL if none
  */
-class BoundaryPointSet * GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
-{
-  class BoundaryLineSet * lines[2] = {line1, line2};
+class BoundaryPointSet *
+GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
+  class BoundaryLineSet * lines[2] =
+    { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map <int, class BoundaryPointSet * > OrderMap;
-  pair < map <int, class BoundaryPointSet * >::iterator, bool > OrderTest;
-  for(int i=0;i<2;i++)  // for both lines
-    for (int j=0;j<2;j++) { // for both endpoints
-      OrderTest = OrderMap.insert ( pair <int, class BoundaryPointSet * >( lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]) );
-      if (!OrderTest.second) { // if insertion fails, we have common endpoint
-        node = OrderTest.first->second;
-        cout << Verbose(5) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl;
-        j=2;
-        i=2;
-        break;
+  map<int, class BoundaryPointSet *> OrderMap;
+  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  for (int i = 0; i < 2; i++)
+    // for both lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+        if (!OrderTest.second)
+          { // if insertion fails, we have common endpoint
+            node = OrderTest.first->second;
+            cout << Verbose(5) << "Common endpoint of lines " << *line1
+                << " and " << *line2 << " is: " << *node << "." << endl;
+            j = 2;
+            i = 2;
+            break;
+          }
       }
-    }
   return node;
-};
+}
+;
 
 /** Determines the boundary points of a cluster.
@@ -204 +261 @@
  * \param *mol molecule structure representing the cluster
  */
-Boundaries * GetBoundaryPoints(ofstream *out, molecule *mol)
+Boundaries *
+GetBoundaryPoints(ofstream *out, molecule *mol)
 {
   atom *Walker = NULL;
@@ -212 +270 @@
 
   *out << Verbose(1) << "Finding all boundary points." << endl;
-  Boundaries *BoundaryPoints = new Boundaries [NDIM];  // first is alpha, second is (r, nr)
+  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
   BoundariesTestPair BoundaryTestPair;
   Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
   double radius, angle;
   // 3a. Go through every axis
-  for (int axis=0; axis<NDIM; axis++)  {
-    AxisVector.Zero();
-    AngleReferenceVector.Zero();
-    AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis+1)%NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis+2)%NDIM] = 1.;
-  //    *out << Verbose(1) << "Axisvector is ";
-  //    AxisVector.Output(out);
-  //    *out << " and AngleReferenceVector is ";
-  //    AngleReferenceVector.Output(out);
-  //    *out << "." << endl;
-  //    *out << " and AngleReferenceNormalVector is ";
-  //    AngleReferenceNormalVector.Output(out);
-  //    *out << "." << endl;
-    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
-      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);
-  }
+  for (int axis = 0; axis < NDIM; axis++)
+    {
+      AxisVector.Zero();
+      AngleReferenceVector.Zero();
+      AngleReferenceNormalVector.Zero();
+      AxisVector.x[axis] = 1.;
+      AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+      AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+      //    *out << Verbose(1) << "Axisvector is ";
+      //    AxisVector.Output(out);
+      //    *out << " and AngleReferenceVector is ";
+      //    AngleReferenceVector.Output(out);
+      //    *out << "." << endl;
+      //    *out << " and AngleReferenceNormalVector is ";
+      //    AngleReferenceNormalVector.Output(out);
+      //    *out << "." << endl;
+      // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+      Walker = mol->start;
+      while (Walker->next != mol->end)
+        {
+          Walker = Walker->next;
+          Vector ProjectedVector;
+          ProjectedVector.CopyVector(&Walker->x);
+          ProjectedVector.ProjectOntoPlane(&AxisVector);
+          // correct for negative side
+          //if (Projection(y) < 0)
+          //angle = 2.*M_PI - angle;
+          radius = ProjectedVector.Norm();
+          if (fabs(radius) > MYEPSILON)
+            angle = ProjectedVector.Angle(&AngleReferenceVector);
+          else
+            angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+
+          //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+          if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+            {
+              angle = 2. * M_PI - angle;
+            }
+          //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+          //ProjectedVector.Output(out);
+          //*out << endl;
+          BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
+              DistancePair (radius, Walker)));
+          if (BoundaryTestPair.second)
+            { // successfully inserted
+            }
+          else
+            { // same point exists, check first r, then distance of original vectors to center of gravity
+              *out << Verbose(2)
+                  << "Encountered two vectors whose projection onto axis "
+                  << axis << " is equal: " << endl;
+              *out << Verbose(2) << "Present vector: ";
+              BoundaryTestPair.first->second.second->x.Output(out);
+              *out << endl;
+              *out << Verbose(2) << "New vector: ";
+              Walker->x.Output(out);
+              *out << endl;
+              double tmp = ProjectedVector.Norm();
+              if (tmp > BoundaryTestPair.first->second.first)
+                {
+                  BoundaryTestPair.first->second.first = tmp;
+                  BoundaryTestPair.first->second.second = Walker;
+                  *out << Verbose(2) << "Keeping new vector." << endl;
+                }
+              else if (tmp == BoundaryTestPair.first->second.first)
+                {
+                  if (BoundaryTestPair.first->second.second->x.ScalarProduct(
+                      &BoundaryTestPair.first->second.second->x)
+                      < Walker->x.ScalarProduct(&Walker->x))
+                    { // Norm() does a sqrt, which makes it a lot slower
+                      BoundaryTestPair.first->second.second = Walker;
+                      *out << Verbose(2) << "Keeping new vector." << endl;
+                    }
+                  else
+                    {
+                      *out << Verbose(2) << "Keeping present vector." << endl;
+                    }
+                }
+              else
+                {
+                  *out << Verbose(2) << "Keeping present vector." << endl;
+                }
+            }
+        }
+      // printing all inserted for debugging
+      //    {
+      //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      //      int i=0;
+      //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+      //        if (runner != BoundaryPoints[axis].begin())
+      //          *out << ", " << i << ": " << *runner->second.second;
+      //        else
+      //          *out << i << ": " << *runner->second.second;
+      //        i++;
+      //      }
+      //      *out << endl;
+      //    }
+      // 3c. throw out points whose distance is less than the mean of left and right neighbours
+      bool flag = false;
+      do
+        { // do as long as we still throw one out per round
+          *out << Verbose(1)
+              << "Looking for candidates to kick out by convex condition ... "
+              << endl;
+          flag = false;
+          Boundaries::iterator left = BoundaryPoints[axis].end();
+          Boundaries::iterator right = BoundaryPoints[axis].end();
+          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+              != BoundaryPoints[axis].end(); runner++)
+            {
+              // set neighbours correctly
+              if (runner == BoundaryPoints[axis].begin())
+                {
+                  left = BoundaryPoints[axis].end();
+                }
+              else
+                {
+                  left = runner;
+                }
+              left--;
+              right = runner;
+              right++;
+              if (right == BoundaryPoints[axis].end())
+                {
+                  right = BoundaryPoints[axis].begin();
+                }
+              // check distance
+
+              // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                {
+                  Vector SideA, SideB, SideC, SideH;
+                  SideA.CopyVector(&left->second.second->x);
+                  SideA.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideA: ";
+                  //          SideA.Output(out);
+                  //          *out << endl;
+
+                  SideB.CopyVector(&right->second.second->x);
+                  SideB.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideB: ";
+                  //          SideB.Output(out);
+                  //          *out << endl;
+
+                  SideC.CopyVector(&left->second.second->x);
+                  SideC.SubtractVector(&right->second.second->x);
+                  SideC.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideC: ";
+                  //          SideC.Output(out);
+                  //          *out << endl;
+
+                  SideH.CopyVector(&runner->second.second->x);
+                  SideH.ProjectOntoPlane(&AxisVector);
+                  //          *out << "SideH: ";
+                  //          SideH.Output(out);
+                  //          *out << endl;
+
+                  // calculate each length
+                  double a = SideA.Norm();
+                  //double b = SideB.Norm();
+                  //double c = SideC.Norm();
+                  double h = SideH.Norm();
+                  // calculate the angles
+                  double alpha = SideA.Angle(&SideH);
+                  double beta = SideA.Angle(&SideC);
+                  double gamma = SideB.Angle(&SideH);
+                  double delta = SideC.Angle(&SideH);
+                  double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
+                      < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+                  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+                  //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+                  if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
+                      < MYEPSILON) && (h < MinDistance))
+                    {
+                      // throw out point
+                      //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+                      BoundaryPoints[axis].erase(runner);
+                      flag = true;
+                    }
+                }
+            }
+        }
+      while (flag);
+    }
   return BoundaryPoints;
-};
+}
+;
 
 /** Determines greatest diameters of a cluster defined by its convex envelope.
@@ -379 +470 @@
  * \return NDIM array of the diameters
  */
-double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
+double *
+GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
+    bool IsAngstroem)
 {
   // get points on boundary of NULL was given as parameter
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
-  if (BoundaryPoints == NULL) {
-    BoundaryFreeFlag = true;
-    BoundaryPoints = GetBoundaryPoints(out, mol);
-  } else {
-    *out << Verbose(1) << "Using given boundary points set." << endl;
-  }
-
+  if (BoundaryPoints == NULL)
+    {
+      BoundaryFreeFlag = true;
+      BoundaryPoints = GetBoundaryPoints(out, mol);
+    }
+  else
+    {
+      *out << Verbose(1) << "Using given boundary points set." << endl;
+    }
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
   double *GreatestDiameter = new double[NDIM];
-  for(int i=0;i<NDIM;i++)
+  for (int i = 0; i < NDIM; i++)
     GreatestDiameter[i] = 0.;
   double OldComponent, tmp, w1, w2;
   Vector DistanceVector, OtherVector;
   int component, Othercomponent;
-  for(int axis=0;axis<NDIM;axis++) { // regard each projected plane
-    //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
-    for (int j=0;j<2;j++) { // and for both axis on the current plane
-      component = (axis+j+1)%NDIM;
-      Othercomponent = (axis+1+((j+1) & 1))%NDIM;
-      //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
-      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
-        // seek for the neighbours pair where the Othercomponent sign flips
-        Neighbour = runner;
-        Neighbour++;
-        if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
-          Neighbour = BoundaryPoints[axis].begin();
-        DistanceVector.CopyVector(&runner->second.second->x);
-        DistanceVector.SubtractVector(&Neighbour->second.second->x);
-        do {  // seek for neighbour pair where it flips
-          OldComponent = DistanceVector.x[Othercomponent];
-          Neighbour++;
-          if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
-            Neighbour = BoundaryPoints[axis].begin();
-          DistanceVector.CopyVector(&runner->second.second->x);
-          DistanceVector.SubtractVector(&Neighbour->second.second->x);
-          //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
-        } while ((runner != Neighbour) && ( fabs( OldComponent/fabs(OldComponent) - DistanceVector.x[Othercomponent]/fabs(DistanceVector.x[Othercomponent]) ) < MYEPSILON)); // as long as sign does not flip
-        if (runner != Neighbour) {
-          OtherNeighbour = Neighbour;
-          if (OtherNeighbour == BoundaryPoints[axis].begin())  // make it wrap around
-            OtherNeighbour = BoundaryPoints[axis].end();
-          OtherNeighbour--;
-          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
-          // now we have found the pair: Neighbour and OtherNeighbour
-          OtherVector.CopyVector(&runner->second.second->x);
-          OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
-          //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
-          //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
-          // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
-          w1 = fabs(OtherVector.x[Othercomponent]);
-          w2 = fabs(DistanceVector.x[Othercomponent]);
-          tmp = fabs((w1*DistanceVector.x[component] + w2*OtherVector.x[component])/(w1+w2));
-          // mark if it has greater diameter
-          //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
-          GreatestDiameter[component] = (GreatestDiameter[component] > tmp) ? GreatestDiameter[component] : tmp;
-        } //else
-          //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
-      }
-    }
-  }
-  *out << Verbose(0) << "RESULT: The biggest diameters are " << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and " << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "." << endl;
+  for (int axis = 0; axis < NDIM; axis++)
+    { // regard each projected plane
+      //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
+      for (int j = 0; j < 2; j++)
+        { // and for both axis on the current plane
+          component = (axis + j + 1) % NDIM;
+          Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
+          //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
+          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+              != BoundaryPoints[axis].end(); runner++)
+            {
+              //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
+              // seek for the neighbours pair where the Othercomponent sign flips
+              Neighbour = runner;
+              Neighbour++;
+              if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
+                Neighbour = BoundaryPoints[axis].begin();
+              DistanceVector.CopyVector(&runner->second.second->x);
+              DistanceVector.SubtractVector(&Neighbour->second.second->x);
+              do
+                { // seek for neighbour pair where it flips
+                  OldComponent = DistanceVector.x[Othercomponent];
+                  Neighbour++;
+                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
+                    Neighbour = BoundaryPoints[axis].begin();
+                  DistanceVector.CopyVector(&runner->second.second->x);
+                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
+                  //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+                }
+              while ((runner != Neighbour) && (fabs(OldComponent / fabs(
+                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
+                  DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
+              if (runner != Neighbour)
+                {
+                  OtherNeighbour = Neighbour;
+                  if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
+                    OtherNeighbour = BoundaryPoints[axis].end();
+                  OtherNeighbour--;
+                  //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
+                  // now we have found the pair: Neighbour and OtherNeighbour
+                  OtherVector.CopyVector(&runner->second.second->x);
+                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
+                  //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
+                  //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
+                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
+                  w1 = fabs(OtherVector.x[Othercomponent]);
+                  w2 = fabs(DistanceVector.x[Othercomponent]);
+                  tmp = fabs((w1 * DistanceVector.x[component] + w2
+                      * OtherVector.x[component]) / (w1 + w2));
+                  // mark if it has greater diameter
+                  //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+                  GreatestDiameter[component] = (GreatestDiameter[component]
+                      > tmp) ? GreatestDiameter[component] : tmp;
+                } //else
+              //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+            }
+        }
+    }
+  *out << Verbose(0) << "RESULT: The biggest diameters are "
+      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
+      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
+      : "atomiclength") << "." << endl;
 
   // free reference lists
   if (BoundaryFreeFlag)
-    delete[](BoundaryPoints);
+    delete[] (BoundaryPoints);
 
   return GreatestDiameter;
-};
-
+}
+;
 
 /*
@@ -460 +569 @@
  * \param *out output stream for debugging
  * \param *tecplot output stream for tecplot data
+ * \param N arbitrary number to differentiate various zones in the tecplot format
  */
-void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol)
-{
-  // 8. Store triangles in tecplot file
-  if (tecplot != NULL) {
-    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
-    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
-    *tecplot << "ZONE T=\"TRIANGLES\", 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;
+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;
-  }
+      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;
+    }
 }
 
@@ -501 +626 @@
  * \return determined volume of the cluster in cubed config:GetIsAngstroem()
  */
-double VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration, Boundaries *BoundaryPtr, molecule *mol)
+double
+VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration,
+    Boundaries *BoundaryPtr, molecule *mol)
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -510 +637 @@
   double volume = 0.;
   double PyramidVolume = 0.;
-  double G,h;
-  Vector x,y;
-  double a,b,c;
+  double G, h;
+  Vector x, y;
+  double a, b, c;
 
   //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
@@ -523 +650 @@
   *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
   Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    Walker->x.Translate(CenterOfGravity);
-  }
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
 
   // 3. Find all points on the boundary
-  if (BoundaryPoints == NULL) {
-    BoundaryFreeFlag = true;
-    BoundaryPoints = GetBoundaryPoints(out, mol);
-  } else {
-    *out << Verbose(1) << "Using given boundary points set." << endl;
-  }
+  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;
+  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;
+  //  *out << Verbose(1) << "Listing PointsOnBoundary:";
+  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+  //    *out << " " << *runner->second;
+  //  }
+  //  *out << endl;
 
   // 5a. guess starting triangle
@@ -558 +692 @@
   TesselStruct->TesselateOnBoundary(out, configuration, mol);
 
-  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
-
-
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
+      << " triangles with " << TesselStruct->LinesOnBoundaryCount
+      << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
+      << endl;
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
-  *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
-  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
-    x.CopyVector(&runner->second->endpoints[0]->node->x);
-    x.SubtractVector(&runner->second->endpoints[1]->node->x);
-    y.CopyVector(&runner->second->endpoints[0]->node->x);
-    y.SubtractVector(&runner->second->endpoints[2]->node->x);
-    a = sqrt(runner->second->endpoints[0]->node->x.Distance(&runner->second->endpoints[1]->node->x));
-    b = sqrt(runner->second->endpoints[0]->node->x.Distance(&runner->second->endpoints[2]->node->x));
-    c = sqrt(runner->second->endpoints[2]->node->x.Distance(&runner->second->endpoints[1]->node->x));
-    G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
-    x.MakeNormalVector(&runner->second->endpoints[0]->node->x, &runner->second->endpoints[1]->node->x, &runner->second->endpoints[2]->node->x);
-    x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
-    h = x.Norm(); // distance of CoG to triangle
-    PyramidVolume = (1./3.) * G * h;    // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
-    *out << Verbose(2) << "Area of triangle is " << G << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-    volume += PyramidVolume;
-  }
-  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-
+  *out << Verbose(1)
+      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
+      << endl;
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
+      != TesselStruct->TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+      x.CopyVector(&runner->second->endpoints[0]->node->x);
+      x.SubtractVector(&runner->second->endpoints[1]->node->x);
+      y.CopyVector(&runner->second->endpoints[0]->node->x);
+      y.SubtractVector(&runner->second->endpoints[2]->node->x);
+      a = sqrt(runner->second->endpoints[0]->node->x.Distance(
+          &runner->second->endpoints[1]->node->x));
+      b = sqrt(runner->second->endpoints[0]->node->x.Distance(
+          &runner->second->endpoints[2]->node->x));
+      c = sqrt(runner->second->endpoints[2]->node->x.Distance(
+          &runner->second->endpoints[1]->node->x));
+      G = sqrt(((a * a + b * b + c * c) * (a * a + b * b + c * c) - 2 * (a * a
+          * a * a + b * b * b * b + c * c * c * c)) / 16.); // area of tesselated triangle
+      x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
+          &runner->second->endpoints[1]->node->x,
+          &runner->second->endpoints[2]->node->x);
+      x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
+      h = x.Norm(); // distance of CoG to triangle
+      PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+      *out << Verbose(2) << "Area of triangle is " << G << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
+          << h << " and the volume is " << PyramidVolume << " "
+          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+      volume += PyramidVolume;
+    }
+  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
+      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
+      << endl;
 
   // 7. translate all points back from CoG
-  *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
+  *out << Verbose(1) << "Translating system back from Center of Gravity."
+      << endl;
   CenterOfGravity->Scale(-1);
   Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    Walker->x.Translate(CenterOfGravity);
-  }
-
-
-
-
-
-  write_tecplot_file(out, tecplot, TesselStruct, mol);
-
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
+
+  // 8. Store triangles in tecplot file
+  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
 
   // free reference lists
   if (BoundaryFreeFlag)
-    delete[](BoundaryPoints);
+    delete[] (BoundaryPoints);
 
   return volume;
-};
-
+}
+;
 
 /** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
@@ -615 +762 @@
  * \param celldensity desired average density in final cell
  */
-void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity)
+void
+PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
+    double ClusterVolume, double celldensity)
 {
   // transform to PAS
@@ -625 +774 @@
   double clustervolume;
   if (ClusterVolume == 0)
-    clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration, BoundaryPoints, mol);
+    clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
+        BoundaryPoints, mol);
   else
     clustervolume = ClusterVolume;
-  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
+      IsAngstroem);
   Vector BoxLengths;
-  int repetition[NDIM] = {1, 1, 1};
+  int repetition[NDIM] =
+    { 1, 1, 1 };
   int TotalNoClusters = 1;
-  for (int i=0;i<NDIM;i++)
+  for (int i = 0; i < NDIM; i++)
     TotalNoClusters *= repetition[i];
 
@@ -638 +790 @@
   double totalmass = 0.;
   atom *Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
-    totalmass += Walker->type->mass;
-  }
-  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
-
-  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass/clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      totalmass += Walker->type->mass;
+    }
+  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
+      << totalmass << " atomicmassunit." << endl;
+
+  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
+      << totalmass / clustervolume << " atomicmassunit/"
+      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 
   // solve cubic polynomial
-  *out << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
+  *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
+      << endl;
   double cellvolume;
   if (IsAngstroem)
-    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_A - (totalmass/clustervolume))/(celldensity-1);
+    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
+        / clustervolume)) / (celldensity - 1);
   else
-    cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
-  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-
-  double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
-  *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  if (minimumvolume > cellvolume) {
-    cerr << Verbose(0) << "ERROR: the containing box already has a greater volume than the envisaged cell volume!" << endl;
-    cout << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
-    for(int i=0;i<NDIM;i++)
-      BoxLengths.x[i] = GreatestDiameter[i];
-    mol->CenterEdge(out, &BoxLengths);
-  } else {
-    BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
-    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1]
-              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2]
-              + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
-    BoxLengths.x[2] = minimumvolume - cellvolume;
-    double x0 = 0.,x1 = 0.,x2 = 0.;
-    if (gsl_poly_solve_cubic(BoxLengths.x[0],BoxLengths.x[1],BoxLengths.x[2],&x0,&x1,&x2) == 1) // either 1 or 3 on return
-      *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
-    else {
-      *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " ." << endl;
-      x0 = x2;  // sorted in ascending order
-    }
-
-    cellvolume = 1;
-    for(int i=0;i<NDIM;i++) {
-      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
-      cellvolume *= BoxLengths.x[i];
-    }
-
-    // set new box dimensions
-    *out << Verbose(0) << "Translating to box with these boundaries." << endl;
-    mol->CenterInBox((ofstream *)&cout, &BoxLengths);
-  }
+    cellvolume = (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;
-};
-
+  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
+      << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
+      << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
+      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+}
+;
 
 // =========================================================== class TESSELATION ===========================================
@@ -702 +884 @@
   LinesOnBoundaryCount = 0;
   TrianglesOnBoundaryCount = 0;
-};
+  TriangleFilesWritten = 0;
+}
+;
 
 /** Constructor of class Tesselation.
@@ -709 +893 @@
 Tesselation::~Tesselation()
 {
-        cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-    delete(runner->second);
-  }
-};
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner
+      != TrianglesOnBoundary.end(); runner++)
+    {
+      delete (runner->second);
+    }
+}
+;
 
 /** Gueses first starting triangle of the convex envelope.
@@ -720 +907 @@
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
  */
-void Tesselation::GuessStartingTriangle(ofstream *out)
+void
+Tesselation::GuessStartingTriangle(ofstream *out)
 {
   // 4b. create a starting triangle
@@ -731 +919 @@
 
   // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end()) {
-    B = A;
-    B++;
-    for (; B != PointsOnBoundary.end(); B++) {
-      C = B;
-      C++;
-      for (; C != PointsOnBoundary.end(); C++) {
-        tmp = A->second->node->x.Distance(&B->second->node->x);
-        distance = tmp*tmp;
-        tmp = A->second->node->x.Distance(&C->second->node->x);
-        distance += tmp*tmp;
-        tmp = B->second->node->x.Distance(&C->second->node->x);
-        distance += tmp*tmp;
-        DistanceMMap.insert( DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(B,C) ) );
-      }
-    }
-  }
-//    // 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;
-
+  if (A != PointsOnBoundary.end())
+    {
+      B = A;
+      B++;
+      for (; B != PointsOnBoundary.end(); B++)
+        {
+          C = B;
+          C++;
+          for (; C != PointsOnBoundary.end(); C++)
+            {
+              tmp = A->second->node->x.Distance(&B->second->node->x);
+              distance = tmp * tmp;
+              tmp = A->second->node->x.Distance(&C->second->node->x);
+              distance += tmp * tmp;
+              tmp = B->second->node->x.Distance(&C->second->node->x);
+              distance += tmp * tmp;
+              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
+                  PointMap::iterator, PointMap::iterator> (B, C)));
+            }
+        }
+    }
+  //    // 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.Distance(&A->second->node->x);
-      int innerpoint = 0;
-      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x))
-          && (tmp < A->second->node->x.Distance(&baseline->second.second->second->node->x)))
-        innerpoint++;
-      tmp = checker->second->node->x.Distance(&baseline->second.first->second->node->x);
-      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x))
-          && (tmp < baseline->second.first->second->node->x.Distance(&baseline->second.second->second->node->x)))
-        innerpoint++;
-      tmp = checker->second->node->x.Distance(&baseline->second.second->second->node->x);
-      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x))
-          && (tmp < baseline->second.second->second->node->x.Distance(&A->second->node->x)))
-        innerpoint++;
-      // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-      if (innerpoint == 3)
-        break;
-    }
-    // 5. come this far, all on same side? Then break 1. loop and construct triangle
-    if (checker == PointsOnBoundary.end()) {
-      *out << "Looks like we have a candidate!" << endl;
-      break;
-    }
-  }
-  if (baseline != DistanceMMap.end()) {
-    BPS[0] = baseline->second.first->second;
-    BPS[1] = baseline->second.second->second;
-    BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-    BPS[0] = A->second;
-    BPS[1] = baseline->second.second->second;
-    BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-    BPS[0] = baseline->second.first->second;
-    BPS[1] = A->second;
-    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-
-    // 4b3. insert created triangle
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
-    TrianglesOnBoundaryCount++;
-    for(int i=0;i<NDIM;i++) {
-      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-      LinesOnBoundaryCount++;
-    }
-
-    *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
-  } else {
-    *out << Verbose(1) << "No starting triangle found." << endl;
-    exit(255);
-  }
-};
-
+  for (; baseline != DistanceMMap.end(); baseline++)
+    {
+      // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+      // 2. next, we have to check whether all points reside on only one side of the triangle
+      // 3. construct plane vector
+      PlaneVector.MakeNormalVector(&A->second->node->x,
+          &baseline->second.first->second->node->x,
+          &baseline->second.second->second->node->x);
+      *out << Verbose(2) << "Plane vector of candidate triangle is ";
+      PlaneVector.Output(out);
+      *out << endl;
+      // 4. loop over all points
+      double sign = 0.;
+      PointMap::iterator checker = PointsOnBoundary.begin();
+      for (; checker != PointsOnBoundary.end(); checker++)
+        {
+          // (neglecting A,B,C)
+          if ((checker == A) || (checker == baseline->second.first) || (checker
+              == baseline->second.second))
+            continue;
+          // 4a. project onto plane vector
+          TrialVector.CopyVector(&checker->second->node->x);
+          TrialVector.SubtractVector(&A->second->node->x);
+          distance = TrialVector.Projection(&PlaneVector);
+          if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+            continue;
+          *out << Verbose(3) << "Projection of " << checker->second->node->Name
+              << " yields distance of " << distance << "." << endl;
+          tmp = distance / fabs(distance);
+          // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+          if ((sign != 0) && (tmp != sign))
+            {
+              // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " outside the convex hull."
+                  << endl;
+              break;
+            }
+          else
+            { // note the sign for later
+              *out << Verbose(2) << "Current candidates: "
+                  << A->second->node->Name << ","
+                  << baseline->second.first->second->node->Name << ","
+                  << baseline->second.second->second->node->Name << " leave "
+                  << checker->second->node->Name << " inside the convex hull."
+                  << endl;
+              sign = tmp;
+            }
+          // 4d. Check whether the point is inside the triangle (check distance to each node
+          tmp = checker->second->node->x.Distance(&A->second->node->x);
+          int innerpoint = 0;
+          if ((tmp < A->second->node->x.Distance(
+              &baseline->second.first->second->node->x)) && (tmp
+              < A->second->node->x.Distance(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.Distance(
+              &baseline->second.first->second->node->x);
+          if ((tmp < baseline->second.first->second->node->x.Distance(
+              &A->second->node->x)) && (tmp
+              < baseline->second.first->second->node->x.Distance(
+                  &baseline->second.second->second->node->x)))
+            innerpoint++;
+          tmp = checker->second->node->x.Distance(
+              &baseline->second.second->second->node->x);
+          if ((tmp < baseline->second.second->second->node->x.Distance(
+              &baseline->second.first->second->node->x)) && (tmp
+              < baseline->second.second->second->node->x.Distance(
+                  &A->second->node->x)))
+            innerpoint++;
+          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+          if (innerpoint == 3)
+            break;
+        }
+      // 5. come this far, all on same side? Then break 1. loop and construct triangle
+      if (checker == PointsOnBoundary.end())
+        {
+          *out << "Looks like we have a candidate!" << endl;
+          break;
+        }
+    }
+  if (baseline != DistanceMMap.end())
+    {
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = A->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = A->second;
+      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+      // 4b3. insert created triangle
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+      TrianglesOnBoundaryCount++;
+      for (int i = 0; i < NDIM; i++)
+        {
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+          LinesOnBoundaryCount++;
+        }
+
+      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
+  else
+    {
+      *out << Verbose(1) << "No starting triangle found." << endl;
+      exit(255);
+    }
+}
+;
 
 /** Tesselates the convex envelope of a cluster from a single starting triangle.
@@ -855 +1078 @@
  * \param *mol the cluster as a molecule structure
  */
-void Tesselation::TesselateOnBoundary(ofstream *out, config *configuration, molecule *mol)
+void
+Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
+    molecule *mol)
 {
   bool flag;
@@ -861 +1086 @@
   class BoundaryPointSet *peak = NULL;
   double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, PropagationVector;
+  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
+      PropagationVector;
   LineMap::iterator LineChecker[2];
-  do {
-    flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
-      if (baseline->second->TrianglesCount == 1) {
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
-        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
-        SmallestAngle = M_PI;
-        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
-        // get peak point with respect to this base line's only triangle
-        for(int i=0;i<3;i++)
-          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
-            peak = BTS->endpoints[i];
-        *out << Verbose(3) << " and has peak " << *peak << "." << endl;
-        // normal vector of triangle
-        BTS->GetNormalVector(NormalVector);
-        *out << Verbose(4) << "NormalVector of base triangle is ";
-        NormalVector.Output(out);
-        *out << endl;
-        // offset to center of triangle
-        CenterVector.Zero();
-        for(int i=0;i<3;i++)
-          CenterVector.AddVector(&BTS->endpoints[i]->node->x);
-        CenterVector.Scale(1./3.);
-        *out << Verbose(4) << "CenterVector of base triangle is ";
-        CenterVector.Output(out);
-        *out << endl;
-        // vector in propagation direction (out of triangle)
-        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-        TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
-        PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
-        TempVector.CopyVector(&CenterVector);
-        TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x);  // TempVector is vector on triangle plane pointing from one baseline egde towards center!
-        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
-        if (PropagationVector.Projection(&TempVector) > 0)  // make sure normal propagation vector points outward from baseline
-          PropagationVector.Scale(-1.);
-        *out << Verbose(4) << "PropagationVector of base triangle is ";
-        PropagationVector.Output(out);
-        *out << endl;
-        winner = PointsOnBoundary.end();
-        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++)
-          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
-            *out << Verbose(3) << "Target point is " << *(target->second) << ":";
-            bool continueflag = true;
-
-            VirtualNormalVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            VirtualNormalVector.AddVector(&baseline->second->endpoints[0]->node->x);
-            VirtualNormalVector.Scale(-1./2.);   // points now to center of base line
-            VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
-            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
-            if (!continueflag) {
-              *out << Verbose(4) << "Angle between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
-              continue;
-            } else
-              *out << Verbose(4) << "Angle between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
-            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
-            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
-  //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
-  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
-  //            else
-  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
-  //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
-  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
-  //            else
-  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
-            // check first endpoint (if any connecting line goes to target or at least not more than 1)
-            continueflag = continueflag && (( (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) || (LineChecker[0]->second->TrianglesCount == 1)));
-            if (!continueflag) {
-              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
-              continue;
-            }
-            // check second endpoint (if any connecting line goes to target or at least not more than 1)
-            continueflag = continueflag && (( (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) || (LineChecker[1]->second->TrianglesCount == 1)));
-            if (!continueflag) {
-              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
-              continue;
-            }
-            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
-            continueflag = continueflag && (!(
-                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
-                && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak))
-               ));
-            if (!continueflag) {
-              *out << Verbose(4) << "Current target is peak!" << endl;
-              continue;
-            }
-            // in case NOT both were found
-            if (continueflag) {  // create virtually this triangle, get its normal vector, calculate angle
-              flag = true;
-              VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x);
-              // make it always point inward
-              if (baseline->second->endpoints[0]->node->x.Projection(&VirtualNormalVector) > 0)
-                VirtualNormalVector.Scale(-1.);
-              // calculate angle
-              TempAngle = NormalVector.Angle(&VirtualNormalVector);
-              *out << Verbose(4) << "NormalVector is ";
-              VirtualNormalVector.Output(out);
-              *out << " and the angle is " << TempAngle << "." << endl;
-              if (SmallestAngle > TempAngle) {  // set to new possible winner
-                SmallestAngle = TempAngle;
-                winner = target;
+  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
           }
-        // 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);
-
-};
+        else
+          *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
+              << " has a triangle count of "
+              << baseline->second->TrianglesCount << "." << endl;
+    }
+  while (flag);
+
+}
+;
 
 /** Adds an atom to the tesselation::PointsOnBoundary list.
  * \param *Walker atom to add
  */
-void Tesselation::AddPoint(atom *Walker)
+void
+Tesselation::AddPoint(atom *Walker)
 {
   PointTestPair InsertUnique;
   BPS[0] = new class BoundaryPointSet(Walker);
-  InsertUnique = PointsOnBoundary.insert( PointPair(Walker->nr, BPS[0]) );
-  if (InsertUnique.second)  // if new point was not present before, increase counter
+  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
+  if (InsertUnique.second) // if new point was not present before, increase counter
     PointsOnBoundaryCount++;
-};
-
-void Tesselation::AddTrianglePoint(atom* Candidate, int n)
-{
-          PointTestPair InsertUnique;
-          TPS[n] = new class BoundaryPointSet(Candidate);
-          InsertUnique = PointsOnBoundary.insert( PointPair(Candidate->nr, TPS[n]) );
-          if (InsertUnique.second)  // if new point was not present before, increase counter
-          {
-                  PointsOnBoundaryCount++;
-          }
-          else
-          {
-                  delete TPS[n];
-                  cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node) << " gibt's schon in der PointMap." << endl;
-                  TPS[n] = (InsertUnique.first)->second;
-          }
-};
-
-/*
- * Function tries to add line from current Points in BPS to BoundaryLineSet.
- * If succesfull it raises the line count and iserts the new line into the BLS, 
- * if unsuccesfull, it writes the line which had been present into the BLS, deleting the new constructed one.  
+}
+;
+
+/** Adds point to Tesselation::PointsOnBoundary if not yet present.
+ * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
+ * @param Candidate point to add
+ * @param n index for this point in Tesselation::TPS array
  */
-
-void Tesselation::AddTriangleLine(int n)
-{
-          LineMap::iterator LineWalker;
-          BLS[n] = new class BoundaryLineSet(BPS, BPS[1]->node->nr);
-          if ((BPS[0]->lines.find(BPS[1]->node->nr))->second->endpoints[0] == BLS[n]->endpoints[0]
-              and (BPS[0]->lines.find(BPS[1]->node->nr))->second->endpoints[1] == BLS[n]->endpoints[1]) 
-                  //If a line is there, how do I recognize that beyond a shadow of a doubt? 
-          {
-                  delete BLS[n];
-                  cout << Verbose(2) << "Tried to add an existing line, handled it." << endl;
-                  LineWalker = LinesOnBoundary.end();
-                  
-                  while(LineWalker->second->endpoints[0] != BLS[n]->endpoints[0] or LineWalker->second->endpoints[1] != BLS[n]->endpoints[1])
-                  {
-                          cout << Verbose(1) << "Looking for line which already exists"<< endl; 
-                          LineWalker--;
-                  }
-                  BLS[n]=LineWalker->second;
-          }
-          else
-          {
-                  cout << Verbose(2) << "Adding line which has not been used before." << endl;
-                  BPS[1]->lines.insert( LinePair(BPS[0]->node->nr, BLS[n]) );
-                  BPS[2]->lines.insert( LinePair(BPS[1]->node->nr, BLS[n]) );
-                  
-                  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BLS[n]));
-                  LinesOnBoundaryCount++;
-                  
-          }
-};
-
-/*
- * Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm). 
+void
+Tesselation::AddTrianglePoint(atom* Candidate, int n)
+{
+  PointTestPair InsertUnique;
+  TPS[n] = new class BoundaryPointSet(Candidate);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
+  if (InsertUnique.second) // if new point was not present before, increase counter
+    {
+      PointsOnBoundaryCount++;
+    }
+  else
+    {
+      delete TPS[n];
+      cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)
+          << " gibt's schon in der PointMap." << endl;
+      TPS[n] = (InsertUnique.first)->second;
+    }
+}
+;
+
+/** Function tries to add line from current Points in BPS to BoundaryLineSet.
+ * If succesful it raises the line count and inserts the new line into the BLS,
+ * if unsuccesful, it writes the line which had been present into the BLS, deleting the new constructed one.
+ * @param *a first endpoint
+ * @param *b second endpoint
+ * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void
+Tesselation::AddTriangleLine(class BoundaryPointSet *a,
+    class BoundaryPointSet *b, int n)
+{
+  LineMap::iterator LineWalker;
+  //cout << "Manually checking endpoints for line." << endl;
+  if ((a->lines.find(b->node->nr))->first == b->node->nr)
+  //If a line is there, how do I recognize that beyond a shadow of a doubt?
+    {
+      //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl;
+
+      LineWalker = LinesOnBoundary.end();
+      LineWalker--;
+
+      while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr,
+          b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max(
+          a->node->nr, b->node->nr))
+        {
+          //cout << Verbose(1) << "Looking for line which already exists"<< endl;
+          LineWalker--;
+        }
+      BPS[0] = LineWalker->second->endpoints[0];
+      BPS[1] = LineWalker->second->endpoints[1];
+      BLS[n] = LineWalker->second;
+
+    }
+  else
+    {
+      cout << Verbose(2)
+          << "Adding line which has not been used before between "
+          << *(a->node) << " and " << *(b->node) << "." << endl;
+      BPS[0] = a;
+      BPS[1] = b;
+      BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+      LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
+      LinesOnBoundaryCount++;
+
+    }
+}
+;
+
+/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
  * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
  */
-
-void Tesselation::AddTriangleToLines()
-{
-        cout <<Verbose(1) << "Adding triangle to its lines" <<endl;
-        int i=0;
-        TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-        TrianglesOnBoundaryCount++;
-        for (i=0; i<3; i++)
-        {
-                BLS[i]->AddTriangle(BTS);
-        }
-};
-
-
-
-/*!
- * This recursive function finds a third point, to form a triangle with two given ones.
+void
+Tesselation::AddTriangleToLines()
+{
+
+  cout << Verbose(1) << "Adding triangle to its lines" << endl;
+  int i = 0;
+  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+  TrianglesOnBoundaryCount++;
+
+  /*
+   * this is apparently done when constructing triangle
+
+   for (i=0; i<3; i++)
+   {
+   BLS[i]->AddTriangle(BTS);
+   }
+   */
+}
+;
+
+/** 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 id, its general \
+ *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
  *  direction and angle into Storage.
  *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
  *  with all neighbours of the candidate.
+ * @param a first point
+ * @param b second point
+ * @param Candidate base point along whose bonds to start looking from
+ * @param Parent point to avoid during search as its wrong direction
+ * @param RecursionLevel contains current recursion depth
+ * @param Chord baseline vector of first and second point
+ * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
+ * @param OldNormal normal of the triangle which the baseline belongs to
+ * @param Opt_Candidate candidate reference to return
+ * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
+ * @param Storage array containing two angles of current Opt_Candidate
+ * @param RADIUS radius of ball
+ * @param mol molecule structure with atoms and bonds
  */
-void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent, int n, Vector *Chord, Vector *d1, Vector *OldNormal, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
+void
+Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
+    atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
 {
   /* OldNormal is normal vector on the old triangle
    * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from a to b.
    */
   Vector dif_a; //Vector from a to candidate
   Vector dif_b; //Vector from b to candidate
   Vector AngleCheck; // Projection of a difference vector on plane orthogonal on triangle side.
+  Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
+
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
   atom *Walker; // variable atom point
+
 
   dif_a.CopyVector(&(a->x));
@@ -1116 +1461 @@
   dif_b.SubtractVector(&(Candidate->x));
   AngleCheck.CopyVector(&dif_a);
+  AngleCheck.Scale(-1);
   AngleCheck.ProjectOntoPlane(Chord);
 
-  if (problem)
-  {
-          cout << "Atom number" << Candidate->nr << endl;
-          Candidate->x.Output((ofstream *)&cout);
-          cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
-  }
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+
+
+  if (DEBUG)
+    {
+      cout << "Atom number" << Candidate->nr << endl;
+      Candidate->x.Output((ofstream *) &cout);
+      cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr]
+          << endl;
+    }
 
   if (a != Candidate and b != Candidate)
-  {
-
-          if (Chord->Norm()/(2*sin(0.5*dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find if Ball will touch atom
-          {
-                  if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
-                  {
-                          Opt_Candidate = Candidate;
-                          Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-                          Storage[1]=AngleCheck.Angle(OldNormal);
-                  }
-                  else
-                  {
-                          if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]> AngleCheck.Angle(OldNormal)) or \
-                                          (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]< AngleCheck.Angle(OldNormal)))
-                                  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-                          {
-                                  Opt_Candidate = Candidate;
-                                  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
-                                  Storage[1]=AngleCheck.Angle(OldNormal);
-                          }
-                          else
-                          {
-                                  if (problem)
-                                          cout << "Looses to better candidate" << endl;
-                          }
-                  }
-          }
-          else
-          {
-                  if (problem)
-                          cout << "erfuellt Dreiecksbedingung fuer sehne nicht" <<endl;
-          }
-  }
+    {
+//      alpha = dif_a.Angle(&dif_b) / 2.;
+//      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+//      SideB = dif_a.Norm();
+//      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+//          alpha); // note this is squared of center line length
+// Those are remains from Freddie. Needed?
+
+//      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+//      Old Chordtest. Replaced by test for Radius of Circumscribing circle
+
+      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)
+        {
+
+          sign = AngleCheck.ScalarProduct(d1);
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+
+          // intermediate calculations:
+
+          Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+          Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+
+          cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+          cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+          cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+          cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+
+          TempNormal.CopyVector(&dif_a);
+          TempNormal.VectorProduct(&dif_b);
+          if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0)
+            {
+              TempNormal.Scale(-1);
+            }
+          TempNormal.Normalize();
+          Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+          cout << "Umkreisradius ist " << Umkreisradius << endl;
+          cout << "Restradius ist " << Restradius << endl;
+          TempNormal.Scale(Restradius);
+          cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+          Mittelpunkt.CopyVector(&Umkreismittelpunkt);
+          Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
+          cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+          cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+          cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+          cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+          cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+          cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+          cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+
+          if (Storage[0]< -1.5) // first Candidate at all
+            {
+
+              cout << "Next better candidate is " << *Candidate << " with ";
+              Opt_Candidate = Candidate;
+              Storage[0] = sign;
+              Storage[1] = AngleCheck.Angle(OldNormal);
+              Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+              cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                  << Storage[0] << endl;
+            }
+          else
+            {
+              //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
+              //within the ball.
+
+              Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
+              cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+
+
+              if (Distance >RADIUS) // We have no interference and may now check whether the new point is better.
+                {
+                  cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+
+                  if ((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))// || //This will give absolute preference to those in "right-hand" quadrants
+                      //sqrt(Candidate->x.Distance(Opt_Mittelpunkt)) < RADIUS)    //and those where Candidate would be within old Sphere.
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                          << Storage[0] << endl;
+                    }
+                  else
+                    {
+                      if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0
+                          && Storage[1] > AngleCheck.Angle(OldNormal)) ||
+                          (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0
+                          && Storage[1] < AngleCheck.Angle(OldNormal)))
+                      //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+                        {
+                          cout << "Next better candidate is " << *Candidate << " with ";
+                          Opt_Candidate = Candidate;
+                          Storage[0] = sign;
+                          Storage[1] = AngleCheck.Angle(OldNormal);
+                          Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                          cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                              << Storage[0] << endl;
+                        }
+                    }
+                }
+              else
+                {
+                  if (DEBUG)
+                    cout << "Looses to better candidate" << endl;
+                }
+            }
+        }
+      else
+        {
+          if (DEBUG)
+            {
+              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+            }
+        }
+    }
+
   else
-  {
-          if (problem)
-                          cout << "identisch mit Ursprungslinie" << endl;
-  }
-
-  if (n<5) // Five is the recursion level threshold.
-  {
-          for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
-      {
-                  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-                  if (Walker->nr == Parent->nr){
-                          continue;
-                  }
-                  else{
-                          Find_next_suitable_point(a, b, Walker, Candidate, n+1, Chord, d1, OldNormal, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
-                  }
-      }
-  }
-};
-
-/*!
- * this function fins a triangle to a line, adjacent to an existing one.
+    {
+      if (DEBUG)
+        cout << "identisch mit Ursprungslinie" << endl;
+    }
+
+  if (RecursionLevel < 7) // Five is the recursion level threshold.
+    {
+      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
+        {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
+              Candidate);
+          if (Walker == Parent)
+            { // don't go back the same bond
+              continue;
+            }
+          else
+            {
+              Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel
+                  + 1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS,
+                  mol); //call function again
+            }
+        }
+    }
+}
+;
+
+/** This function finds a triangle to a line, adjacent to an existing one.
+ * @param out   output stream for debugging
+ * @param tecplot output stream for writing found triangles in TecPlot format
+ * @param mol molecule structure with all atoms and bonds
+ * @param Line current baseline to search from
+ * @param T current triangle which \a Line is edge of
+ * @param RADIUS radius of the rolling ball
+ * @param N number of found triangles
  */
-void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N)
-{
-        cout << Verbose(1) << "Looking for next suitable triangle \n";
+void
+Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
+    molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+    const double& RADIUS, int N)
+{
+  cout << Verbose(1) << "Looking for next suitable triangle \n";
   Vector direction1;
   Vector helper;
   Vector Chord;
+  ofstream *tempstream = NULL;
+  char filename[255];
   //atom* Walker;
 
   double Storage[2];
-  Storage[0]=-2.;    // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[1]=9999999.;  // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1] = 9999999.; // This is also lower then any value produced by an eligible atom, which are all positive
   atom* Opt_Candidate = NULL;
-
+  Vector Opt_Mittelpunkt;
 
   cout << Verbose(1) << "Constructing helpful vectors ... " << endl;
   helper.CopyVector(&(Line.endpoints[0]->node->x));
-  for (int i =0; i<3; i++)
-    {
-      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
-        {
-          helper.SubtractVector(&T.endpoints[i]->node->x);
-          break;
-        }
-    }
-
+  for (int i = 0; i < 3; i++)
+    {
+      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr
+          && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr)
+        {
+          helper.SubtractVector(&T.endpoints[i]->node->x);
+          break;
+        }
+    }
 
   direction1.CopyVector(&Line.endpoints[0]->node->x);
@@ -1213 +1677 @@
   direction1.VectorProduct(&(T.NormalVector));
 
-  if (direction1.ScalarProduct(&helper)<0)
+  if (direction1.ScalarProduct(&helper) < 0)
     {
       direction1.Scale(-1);
@@ -1221 +1685 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
-  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
-
-  if (N>7)
-  {
-          cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
-          write_tecplot_file(out, tecplot, this, mol);
-          tecplot->flush();
-          Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol, 1);
-          Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol, 1);
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... "
+      << endl;
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
+      Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
+      &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
+      Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
+      &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+
+  if ((TrianglesOnBoundaryCount % 10) == 0)
+    {
+    sprintf(filename, "testEnvelope-%d.dat", TriangleFilesWritten);
+    tempstream = new ofstream(filename, ios::trunc);
+    write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten++);
+    tempstream->close();
+    tempstream->flush();
+    delete(tempstream);
   }
+
+  if ((Opt_Candidate == NULL) || (N<12))
+    {
+      cout << Verbose(1)
+          << "No new Atom found, triangle construction will crash" << endl;
+      write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
+      Find_next_suitable_point(Line.endpoints[0]->node,
+          Line.endpoints[1]->node, Line.endpoints[0]->node,
+          Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector),
+          Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+      Find_next_suitable_point(Line.endpoints[0]->node,
+          Line.endpoints[1]->node, Line.endpoints[1]->node,
+          Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector),
+          Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+    }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
 
-  cout << Verbose(1) << "Adding exactly one Walker for reasons completely unknown to me ... " << endl;
   cout << " Optimal candidate is " << *Opt_Candidate << endl;
 
@@ -1241 +1725 @@
   AddTrianglePoint(Line.endpoints[0]->node, 1);
   AddTrianglePoint(Line.endpoints[1]->node, 2);
-  
-  BPS[0] = TPS[0];
-  BPS[1] = TPS[1];
-  AddTriangleLine(0);
-  BPS[0] = TPS[0];
-  BPS[1] = TPS[2];
-  AddTriangleLine(1);
-  BPS[0] = TPS[1];
-  BPS[1] = TPS[2];
-  AddTriangleLine(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();
-  
-  
-  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
+  cout << "New triangle with " << *BTS << endl;
+
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... "
+      << endl;
 
   BTS->GetNormalVector(BTS->NormalVector);
 
-  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))<0 && Storage[0]>0) ||
-                  ((BTS->NormalVector.ScalarProduct(&(T.NormalVector))>0 && Storage[0]<0)) )
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)
+      || ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0]
+          < 0)))
     {
       BTS->NormalVector.Scale(-1);
     };
 
-};
-
-
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent, int n, Vector Oben, atom*& Opt_Candidate, double Storage[2], molecule* mol, double RADIUS)
-{
-        cout << Verbose(1) << "Looking for second point of starting triangle, recursive level "<< n <<endl;;
-        int i;
-        Vector AngleCheck;
-        atom* Walker;
-
-        if (a->nr !=Candidate->nr)
-        {
-                AngleCheck.CopyVector(&(Candidate->x));
-                AngleCheck.SubtractVector(&(a->x));
-                if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
-                {
-                        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                        Opt_Candidate=Candidate;
-                        Storage[0]=AngleCheck.Angle(&Oben);
-                        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
-                }
-                else{
-                        if (AngleCheck.Norm() > RADIUS)
-                        {
-                                cout << Verbose(1) << "refused due to Radius" << AngleCheck.Norm() << endl;
-                        }
-                        else{
-                                cout << Verbose(1) << "Supposedly looses to a better candidate" << Opt_Candidate->nr << endl;
-                        }
-                }
-        }
-
-        if (n<5)
-    {
-      for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
-      {
-          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
-          if (Walker->nr == Parent->nr)
-                  continue;
-          else
-                  Find_second_point_for_Tesselation(a, Walker, Candidate, n+1, Oben, Opt_Candidate, Storage, mol, RADIUS);
-      };
+}
+;
+
+void
+Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[2],
+    molecule* mol, double RADIUS)
+{
+  cout << Verbose(1)
+      << "Looking for second point of starting triangle, recursive level "
+      << RecursionLevel << endl;;
+  int i;
+  Vector AngleCheck;
+  atom* Walker;
+  double norm = -1.;
+
+  // check if we only have one unique point yet ...
+  if (a != 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)
+        {
+          if (AngleCheck.Angle(&Oben) < Storage[0])
+            {
+              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+              cout << "Next better candidate is " << *Candidate
+                  << " with distance " << norm << ".\n";
+              Opt_Candidate = Candidate;
+              Storage[0] = AngleCheck.Angle(&Oben);
+              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
+            }
+          else
+            {
+              cout << Verbose(1) << "Supposedly looses to a better candidate "
+                  << *Opt_Candidate << endl;
+            }
+        }
+      else
+        {
+          cout << Verbose(1) << *Candidate << " refused due to Radius " << norm
+              << endl;
+        }
+    }
+
+  // if not recursed to deeply, look at all its bonds
+  if (RecursionLevel < 7)
+    {
+      for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
+        {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
+              Candidate);
+          if (Walker == Parent) // don't go back along the bond we came from
+            continue;
+          else
+            Find_second_point_for_Tesselation(a, Walker, Candidate,
+                RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+        };
     };
-
-
-};
-
-
-void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
-{
-        cout << Verbose(1) << "Looking for starting triangle \n";
-  int i=0;
+}
+;
+
+void
+Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
+{
+  cout << Verbose(1) << "Looking for starting triangle \n";
+  int i = 0;
   atom* Walker;
   atom* FirstPoint;
   atom* SecondPoint;
-  int max_index[3];
+  atom* max_index[3];
   double max_coordinate[3];
   Vector Oben;
   Vector helper;
   Vector Chord;
+  Vector Opt_Mittelpunkt;
 
   Oben.Zero();
 
-
-  for(i =0; i<3; i++)
-    {
-      max_index[i] =-1;
-      max_coordinate[i] =-1;
-    }
-cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
+  for (i = 0; i < 3; i++)
+    {
+      max_index[i] = NULL;
+      max_coordinate[i] = -1;
+    }
+  cout << Verbose(1) << "Molecule mol is there and has " << mol->AtomCount
+      << " Atoms \n";
+
+  // 1. searching topmost atom with respect to each axis
   Walker = mol->start;
   while (Walker->next != mol->end)
     {
-        Walker = Walker->next;
-        for (i=0; i<3; i++)
-        {
-                if (Walker->x.x[i] > max_coordinate[i])
-                {
-                        max_coordinate[i]=Walker->x.x[i];
-                        max_index[i]=Walker->nr;
-                }
-        }
-    }
-
-  cout << Verbose(1) << "Found maximum coordinates. "<< endl;
+      Walker = Walker->next;
+      for (i = 0; i < 3; i++)
+        {
+          if (Walker->x.x[i] > max_coordinate[i])
+            {
+              max_coordinate[i] = Walker->x.x[i];
+              max_index[i] = Walker;
+            }
+        }
+    }
+
+  cout << Verbose(1) << "Found maximum coordinates. " << endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k=0;
-
-  Oben.x[k]=1.;
-  FirstPoint = mol->start;
-  FirstPoint = FirstPoint->next;
-  while (FirstPoint->nr != max_index[k])
-    {
-          FirstPoint = FirstPoint->next;
-    }
-  cout << Verbose(1) << "Coordinates of start atom: " << FirstPoint->x.x[0] << endl;
+  const int k = 1;
+  Oben.x[k] = 1.;
+  FirstPoint = max_index[k];
+
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
+      << FirstPoint->x.x[0] << endl;
   double Storage[2];
   atom* Opt_Candidate = NULL;
-  Storage[0]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[1]=999999.;  // This will be an angle looking for the third point.
-  cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
-
-  Find_second_point_for_Tesselation(FirstPoint, mol->ListOfBondsPerAtom[FirstPoint->nr][0]->GetOtherAtom(FirstPoint), FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
-
-
+  Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[1] = 999999.; // This will be an angle looking for the third point.
+  cout << Verbose(1) << "Number of Bonds: "
+      << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
+
+  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
+      Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
   SecondPoint = Opt_Candidate;
-  Opt_Candidate=NULL;
+  cout << Verbose(1) << "Found second point is " << *SecondPoint << ".\n";
+
   helper.CopyVector(&(FirstPoint->x));
   helper.SubtractVector(&(SecondPoint->x));
+  helper.Normalize();
   Oben.ProjectOntoPlane(&helper);
+  Oben.Normalize();
   helper.VectorProduct(&Oben);
-  Storage[0]=-2.;       // This will indicate the quadrant.
-  Storage[1]= 9999999.; // This will be an angle looking for the third point.
+  Storage[0] = -2.; // This will indicate the quadrant.
+  Storage[1] = 9999999.; // This will be an angle looking for the third point.
 
   Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
   Chord.SubtractVector(&(SecondPoint->x));
+  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
 
   cout << Verbose(1) << "Looking for third point candidates \n";
-  Find_next_suitable_point(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
-
-
-  //Starting Triangle is Walker, SecondPoint, Opt_Candidate
-
-    cout << Verbose(1) << "The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << "." << endl;
-
+  // look in one direction of baseline for initial candidate
+  Opt_Candidate = NULL;
+  Opt_Mittelpunkt;
+  Find_next_suitable_point(FirstPoint, SecondPoint, SecondPoint, FirstPoint, 0,
+      &Chord, &helper, &Oben, Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+  // look in other direction of baseline for possible better candidate
+  Find_next_suitable_point(FirstPoint, SecondPoint, FirstPoint, SecondPoint, 0,
+      &Chord, &helper, &Oben, Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+
+  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+
+  cout << Verbose(1) << "The found starting triangle consists of "
+      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
+      << "." << endl;
+
+  // Finally, we only have to add the found points
   AddTrianglePoint(FirstPoint, 0);
   AddTrianglePoint(SecondPoint, 1);
   AddTrianglePoint(Opt_Candidate, 2);
-  
-  BPS[0] = TPS[0];
-  BPS[1] = TPS[1];
-  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = TPS[1];
-  BPS[1] = TPS[2];
-  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = TPS[0];
-  BPS[1] = TPS[2];
-  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-
-  Tesselation::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++;
-    };
-
-       BTS->GetNormalVector(BTS->NormalVector);
-
-       if( BTS->NormalVector.ScalarProduct(&Oben)<0)
-         {
-           BTS->NormalVector.Scale(-1);
-         }
-};
-
-
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
-{
-        int N =0;
-        struct Tesselation *Tess = new Tesselation;
-        cout << Verbose(1) << "Entering search for non convex hull. " << endl;
-        cout << flush;
-        const double RADIUS =6.;
-        LineMap::iterator baseline;
-        Tess->Find_starting_triangle(mol, RADIUS);
-
-        baseline = Tess->LinesOnBoundary.begin();
-        while (baseline != Tess->LinesOnBoundary.end()) {
-                if (baseline->second->TrianglesCount == 1)
-                {
-                        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
-                        Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
-                        cout << Verbose(1) << "End of Tesselation ... " << endl;
-                }
-                else
-                {
-                        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
-                }
-                N++;
-                baseline++;
-        }
-
-};
+  // ... and respective lines
+  AddTriangleLine(TPS[0], TPS[1], 0);
+  AddTriangleLine(TPS[1], TPS[2], 1);
+  AddTriangleLine(TPS[0], TPS[2], 2);
+  // ... and triangles to the Maps of the Tesselation class
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTriangleToLines();
+  // ... and calculate its normal vector (with correct orientation)
+  Oben.Scale(-1.);
+  BTS->GetNormalVector(Oben);
+}
+;
+
+void
+Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
+{
+  int N = 0;
+  struct Tesselation *Tess = new Tesselation;
+  cout << Verbose(1) << "Entering search for non convex hull. " << endl;
+  cout << flush;
+  const double RADIUS = 6.;
+  LineMap::iterator baseline;
+  Tess->Find_starting_triangle(mol, RADIUS);
+
+  baseline = Tess->LinesOnBoundary.begin();
+  while (baseline != Tess->LinesOnBoundary.end())
+    {
+      if (baseline->second->TrianglesCount == 1)
+        {
+          cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+          Tess->Find_next_suitable_triangle(out, tecplot, mol,
+              *(baseline->second),
+              *(((baseline->second->triangles.begin()))->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
+          cout << Verbose(1) << "End of Tesselation ... " << endl;
+        }
+      else
+        {
+          cout << Verbose(1) << "There is a line with "
+              << baseline->second->TrianglesCount << " triangles adjacent"
+              << endl;
+        }
+      N++;
+      baseline++;
+    }
+  write_tecplot_file(out, tecplot, Tess, mol, -1);
+
+}
+;

src/boundary.hpp

@@ -83 +83 @@
     void AddPoint(atom * Walker);
     void AddTrianglePoint(atom* Candidate, int n);
-    void AddTriangleLine(int n);
+    void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
     void AddTriangleToLines();
     void Find_starting_triangle(molecule* mol, const double RADIUS);
@@ -91 +91 @@
     LineMap LinesOnBoundary;
     TriangleMap TrianglesOnBoundary;
-    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions 
+    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
     class BoundaryPointSet *BPS[2];
     class BoundaryLineSet *BLS[3];
@@ -98 +98 @@
     int LinesOnBoundaryCount;
     int TrianglesOnBoundaryCount;
+    int TriangleFilesWritten;
 };
 

src/vector.cpp

@@ -22 +22 @@
 Vector::~Vector() {};
 
+/** Calculates square of distance between this and another vector.
+ * \param *y array to second vector
+ * \return \f$| x - y |^2\f$
+ */
+double Vector::DistanceSquared(const Vector *y) const
+{
+  double res = 0.;
+  for (int i=NDIM;i--;)
+    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
+  return (res);
+};
+
 /** Calculates distance between this and another vector.
  * \param *y array to second vector
- * \return \f$| x - y |^2\f$
+ * \return \f$| x - y |\f$
  */
 double Vector::Distance(const Vector *y) const
@@ -31 +43 @@
   for (int i=NDIM;i--;)
     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (res);
+  return (sqrt(res));
 };
 

src/vector.hpp

@@ -7 +7 @@
  * basically, just a x[3] but with helpful functions
  */
-class Vector {   
+class Vector {
   public:
     double x[NDIM];
@@ -16 +16 @@
 
   double Distance(const Vector *y) const;
+  double DistanceSquared(const Vector *y) const;
   double PeriodicDistance(const Vector *y, const double *cell_size) const;
   double ScalarProduct(const Vector *y) const;
@@ -28 +29 @@
   void VectorProduct(const Vector *y);
   void ProjectOntoPlane(const Vector *y);
-  void Zero(); 
+  void Zero();
   void One(double one);
   void Init(double x1, double x2, double x3);
@@ -41 +42 @@
   void KeepPeriodic(ofstream *out, double *matrix);
   void LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors);
-  
+
   double CutsPlaneAt(Vector *A, Vector *B, Vector *C);
   bool GetOneNormalVector(const Vector *x1);