Context Navigation

← Previous Change
Next Change →

Changeset 6ac7ee for src/boundary.cpp

Timestamp:

Feb 9, 2009, 5:24:10 PM (16 years ago)

Author:

Frederik Heber <heber@…>

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:

d8b94a

Parents:

124df1

git-author:

Frederik Heber <heber@…> (02/09/09 15:55:37)

git-committer:

Frederik Heber <heber@…> (02/09/09 17:24:10)

Message:

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

Basically, this resulted from a lot of conversions two from spaces to one tab, which is my standard indentation. The mess was caused by eclipse auto-indenting. And in espack2:ConcaveHull was the new stuff, so all from ConcaveHull was replaced in case of doubt.
Additionally, vector had ofstream << operator instead ostream << ...

File:

: 1 edited

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

Legend:

: Unmodified
: Added
: Removed

src/boundary.cpp

Property mode changed from 100644 to 100755

-              r124df1
+              r6ac7ee
-#include "molecules.hpp"
 #include "boundary.hpp"
 #define DEBUG 1
+#define DoTecplotOutput 0
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
 #define DoRaster3DOutput 1
+#define DoVRMLOutput 1
 #define TecplotSuffix ".dat"
 #define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
+#define HULLEPSILON MYEPSILON
 // ======================================== Points on Boundary =================================
 …
 BoundaryPointSet::BoundaryPointSet()
+{
   LinesCount = 0;
   Nr = -1;
+        LinesCount = 0;
+        Nr = -1;
+}
+;
 …
 BoundaryPointSet::BoundaryPointSet(atom *Walker)
+{
   node = Walker;
   LinesCount = 0;
   Nr = Walker->nr;
+        node = Walker;
+        LinesCount = 0;
+        Nr = Walker->nr;
+}
+;
 …
 BoundaryPointSet::~BoundaryPointSet()
+{
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  node = NULL;
+}
+;
+void
+BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+      << endl;
+  if (line->endpoints[0] == this)
+    {
+      lines.insert(LinePair(line->endpoints[1]->Nr, line));
+    }
+  else
+    {
+      lines.insert(LinePair(line->endpoints[0]->Nr, line));
+    }
+  LinesCount++;
+        cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+        if (!lines.empty())
+                cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
+        node = NULL;
+}
+;
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+        cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+                        << endl;
+        if (line->endpoints[0] == this)
+                {
+                        lines.insert(LinePair(line->endpoints[1]->Nr, line));
+                }
+        else
+                {
+                        lines.insert(LinePair(line->endpoints[0]->Nr, line));
+                }
+        LinesCount++;
+}
+;
 …
 operator <<(ostream &ost, BoundaryPointSet &a)
+{
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
   return ost;
+        ost << "[" << a.Nr << "|" << a.node->Name << "]";
+        return ost;
+}
+;
 …
 BoundaryLineSet::BoundaryLineSet()
+{
   for (int i = 0; i < 2; i++)
     endpoints[i] = NULL;
   TrianglesCount = 0;
   Nr = -1;
+        for (int i = 0; i < 2; i++)
+                endpoints[i] = NULL;
+        TrianglesCount = 0;
+        Nr = -1;
+}
+;
 …
 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
   // set number
   Nr = number;
   // set endpoints in ascending order
   SetEndpointsOrdered(endpoints, Point[0], Point[1]);
   // add this line to the hash maps of both endpoints
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
   // clear triangles list
   TrianglesCount = 0;
   cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+        // set number
+        Nr = number;
+        // set endpoints in ascending order
+        SetEndpointsOrdered(endpoints, Point[0], Point[1]);
+        // add this line to the hash maps of both endpoints
+        Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+        Point[1]->AddLine(this); //
+        // clear triangles list
+        TrianglesCount = 0;
+        cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+}
+;
 …
 BoundaryLineSet::~BoundaryLineSet()
+{
+  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;
+    }
+        int Numbers[2];
+        Numbers[0] = endpoints[1]->Nr;
+        Numbers[1] = endpoints[0]->Nr;
+        for (int i = 0; i < 2; i++) {
+                cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+                endpoints[i]->lines.erase(Numbers[i]);
+                if (endpoints[i]->lines.empty()) {
+                        cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+                        if (endpoints[i] != NULL) {
+                                delete(endpoints[i]);
+                                endpoints[i] = NULL;
+                        } else
+                                cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+                } else
+                        cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+        }
+        if (!triangles.empty())
+                cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+}
+;
 …
 BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
       << endl;
   triangles.insert(TrianglePair(TrianglesCount, triangle));
   TrianglesCount++;
+        cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+                        << endl;
+        triangles.insert(TrianglePair(triangle->Nr, triangle));
+        TrianglesCount++;
+}
+;
 …
 operator <<(ostream &ost, BoundaryLineSet &a)
+{
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
       << a.endpoints[1]->node->Name << "]";
   return ost;
+        ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+                        << a.endpoints[1]->node->Name << "]";
+        return ost;
+}
+;
 …
 BoundaryTriangleSet::BoundaryTriangleSet()
+{
   for (int i = 0; i < 3; i++)
+    {
       endpoints[i] = NULL;
       lines[i] = NULL;
+    }
   Nr = -1;
+        for (int i = 0; i < 3; i++)
+                {
+                        endpoints[i] = NULL;
+                        lines[i] = NULL;
+                }
+        Nr = -1;
+}
+;
 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3],
     int number)
+{
   // set number
   Nr = number;
   // set lines
   cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
   for (int i = 0; i < 3; i++)
+    {
       lines[i] = line[i];
       lines[i]->AddTriangle(this);
+    }
   // get ascending order of endpoints
   map<int, class BoundaryPointSet *> OrderMap;
   for (int i = 0; i < 3; i++)
     // for all three lines
     for (int j = 0; j < 2; j++)
       { // for both endpoints
         OrderMap.insert(pair<int, class BoundaryPointSet *> (
             line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
         // and we don't care whether insertion fails
+      }
   // set endpoints
   int Counter = 0;
   cout << Verbose(6) << " with end points ";
   for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
       != OrderMap.end(); runner++)
+    {
       endpoints[Counter] = runner->second;
       cout << " " << *endpoints[Counter];
       Counter++;
+    }
   if (Counter < 3)
+    {
       cerr << "ERROR! We have a triangle with only two distinct endpoints!"
           << endl;
       //exit(1);
+    }
   cout << "." << endl;
+                int number)
+{
+        // set number
+        Nr = number;
+        // set lines
+        cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+        for (int i = 0; i < 3; i++)
+                {
+                        lines[i] = line[i];
+                        lines[i]->AddTriangle(this);
+                }
+        // get ascending order of endpoints
+        map<int, class BoundaryPointSet *> OrderMap;
+        for (int i = 0; i < 3; i++)
+                // for all three lines
+                for (int j = 0; j < 2; j++)
+                        { // for both endpoints
+                                OrderMap.insert(pair<int, class BoundaryPointSet *> (
+                                                line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+                                // and we don't care whether insertion fails
+                        }
+        // set endpoints
+        int Counter = 0;
+        cout << Verbose(6) << " with end points ";
+        for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+                        != OrderMap.end(); runner++)
+                {
+                        endpoints[Counter] = runner->second;
+                        cout << " " << *endpoints[Counter];
+                        Counter++;
+                }
+        if (Counter < 3)
+                {
+                        cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+                                        << endl;
+                        //exit(1);
+                }
+        cout << "." << endl;
+}
+;
 …
 BoundaryTriangleSet::~BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++)
+    {
+      cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+      lines[i]->triangles.erase(Nr);
+      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;
+    }
+        for (int i = 0; i < 3; i++) {
+                cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+                lines[i]->triangles.erase(Nr);
+                if (lines[i]->triangles.empty()) {
+                        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+                        if (lines[i] != NULL) {
+                                delete (lines[i]);
+                                lines[i] = NULL;
+                        } else
+                                cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+                } else
+                        cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+        }
+}
+;
 …
 BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+{
   // get normal vector
   NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
       &endpoints[2]->node->x);
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
   if (endpoints[0]->node->x.Projection(&OtherVector) > 0)
     NormalVector.Scale(-1.);
+        // get normal vector
+        NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
+                        &endpoints[2]->node->x);
+        // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+        if (NormalVector.Projection(&OtherVector) > 0)
+                NormalVector.Scale(-1.);
+}
+;
 …
 operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
       << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
   return ost;
+        ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+                        << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+        return ost;
+}
+;
 …
 GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
+{
   class BoundaryLineSet * lines[2] =
     { line1, line2 };
   class BoundaryPointSet *node = NULL;
   map<int, class BoundaryPointSet *> OrderMap;
   pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
     for (int j = 0; j < 2; j++)
       { // for both endpoints
         OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
             lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
         if (!OrderTest.second)
           { // if insertion fails, we have common endpoint
             node = OrderTest.first->second;
             cout << Verbose(5) << "Common endpoint of lines " << *line1
                 << " and " << *line2 << " is: " << *node << "." << endl;
             j = 2;
             i = 2;
             break;
+          }
+      }
   return node;
+        class BoundaryLineSet * lines[2] =
+                { line1, line2 };
+        class BoundaryPointSet *node = NULL;
+        map<int, class BoundaryPointSet *> OrderMap;
+        pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+        for (int i = 0; i < 2; i++)
+                // for both lines
+                for (int j = 0; j < 2; j++)
+                        { // for both endpoints
+                                OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
+                                                lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
+                                if (!OrderTest.second)
+                                        { // if insertion fails, we have common endpoint
+                                                node = OrderTest.first->second;
+                                                cout << Verbose(5) << "Common endpoint of lines " << *line1
+                                                                << " and " << *line2 << " is: " << *node << "." << endl;
+                                                j = 2;
+                                                i = 2;
+                                                break;
+                                        }
+                        }
+        return node;
+}
+;
 …
 GetBoundaryPoints(ofstream *out, molecule *mol)
+{
   atom *Walker = NULL;
   PointMap PointsOnBoundary;
   LineMap LinesOnBoundary;
   TriangleMap TrianglesOnBoundary;
   *out << Verbose(1) << "Finding all boundary points." << endl;
   Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
   BoundariesTestPair BoundaryTestPair;
   Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
   double radius, angle;
   // 3a. Go through every axis
   for (int axis = 0; axis < NDIM; axis++)
+    {
       AxisVector.Zero();
       AngleReferenceVector.Zero();
       AngleReferenceNormalVector.Zero();
       AxisVector.x[axis] = 1.;
       AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
       AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
       //    *out << Verbose(1) << "Axisvector is ";
       //    AxisVector.Output(out);
       //    *out << " and AngleReferenceVector is ";
       //    AngleReferenceVector.Output(out);
       //    *out << "." << endl;
       //    *out << " and AngleReferenceNormalVector is ";
       //    AngleReferenceNormalVector.Output(out);
       //    *out << "." << endl;
       // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
       Walker = mol->start;
       while (Walker->next != mol->end)
+        {
           Walker = Walker->next;
           Vector ProjectedVector;
           ProjectedVector.CopyVector(&Walker->x);
           ProjectedVector.ProjectOntoPlane(&AxisVector);
           // correct for negative side
           //if (Projection(y) < 0)
           //angle = 2.*M_PI - angle;
           radius = ProjectedVector.Norm();
           if (fabs(radius) > MYEPSILON)
             angle = ProjectedVector.Angle(&AngleReferenceVector);
           else
             angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
           //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
           if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+            {
               angle = 2. * M_PI - angle;
+            }
           //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
           //ProjectedVector.Output(out);
           //*out << endl;
           BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
               DistancePair (radius, Walker)));
           if (BoundaryTestPair.second)
             { // successfully inserted
+            }
           else
             { // same point exists, check first r, then distance of original vectors to center of gravity
               *out << Verbose(2)
                   << "Encountered two vectors whose projection onto axis "
                   << axis << " is equal: " << endl;
               *out << Verbose(2) << "Present vector: ";
               BoundaryTestPair.first->second.second->x.Output(out);
               *out << endl;
               *out << Verbose(2) << "New vector: ";
               Walker->x.Output(out);
               *out << endl;
               double tmp = ProjectedVector.Norm();
               if (tmp > BoundaryTestPair.first->second.first)
+                {
                   BoundaryTestPair.first->second.first = tmp;
                   BoundaryTestPair.first->second.second = Walker;
                   *out << Verbose(2) << "Keeping new vector." << endl;
+                }
               else if (tmp == BoundaryTestPair.first->second.first)
+                {
                   if (BoundaryTestPair.first->second.second->x.ScalarProduct(
                       &BoundaryTestPair.first->second.second->x)
                       < Walker->x.ScalarProduct(&Walker->x))
                     { // Norm() does a sqrt, which makes it a lot slower
                       BoundaryTestPair.first->second.second = Walker;
                       *out << Verbose(2) << "Keeping new vector." << endl;
+                    }
                   else
+                    {
                       *out << Verbose(2) << "Keeping present vector." << endl;
+                    }
+                }
               else
+                {
                   *out << Verbose(2) << "Keeping present vector." << endl;
+                }
+            }
+        }
       // printing all inserted for debugging
       //    {
       //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
       //      int i=0;
       //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
       //        if (runner != BoundaryPoints[axis].begin())
       //          *out << ", " << i << ": " << *runner->second.second;
       //        else
       //          *out << i << ": " << *runner->second.second;
       //        i++;
       //      }
       //      *out << endl;
       //    }
       // 3c. throw out points whose distance is less than the mean of left and right neighbours
       bool flag = false;
       do
         { // do as long as we still throw one out per round
           *out << Verbose(1)
               << "Looking for candidates to kick out by convex condition ... "
               << endl;
           flag = false;
           Boundaries::iterator left = BoundaryPoints[axis].end();
           Boundaries::iterator right = BoundaryPoints[axis].end();
           for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
               != BoundaryPoints[axis].end(); runner++)
+            {
               // set neighbours correctly
               if (runner == BoundaryPoints[axis].begin())
+                {
                   left = BoundaryPoints[axis].end();
+                }
               else
+                {
                   left = runner;
+                }
               left--;
               right = runner;
               right++;
               if (right == BoundaryPoints[axis].end())
+                {
                   right = BoundaryPoints[axis].begin();
+                }
               // check distance
               // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                {
                   Vector SideA, SideB, SideC, SideH;
                   SideA.CopyVector(&left->second.second->x);
                   SideA.ProjectOntoPlane(&AxisVector);
                   //          *out << "SideA: ";
                   //          SideA.Output(out);
                   //          *out << endl;
                   SideB.CopyVector(&right->second.second->x);
                   SideB.ProjectOntoPlane(&AxisVector);
                   //          *out << "SideB: ";
                   //          SideB.Output(out);
                   //          *out << endl;
                   SideC.CopyVector(&left->second.second->x);
                   SideC.SubtractVector(&right->second.second->x);
                   SideC.ProjectOntoPlane(&AxisVector);
                   //          *out << "SideC: ";
                   //          SideC.Output(out);
                   //          *out << endl;
                   SideH.CopyVector(&runner->second.second->x);
                   SideH.ProjectOntoPlane(&AxisVector);
                   //          *out << "SideH: ";
                   //          SideH.Output(out);
                   //          *out << endl;
                   // calculate each length
                   double a = SideA.Norm();
                   //double b = SideB.Norm();
                   //double c = SideC.Norm();
                   double h = SideH.Norm();
                   // calculate the angles
                   double alpha = SideA.Angle(&SideH);
                   double beta = SideA.Angle(&SideC);
                   double gamma = SideB.Angle(&SideH);
                   double delta = SideC.Angle(&SideH);
                   double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
                       < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
                   //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
                   //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
                   if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
                       < MYEPSILON) && (h < MinDistance))
+                    {
                       // throw out point
                       //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
                       BoundaryPoints[axis].erase(runner);
                       flag = true;
+                    }
+                }
+            }
+        }
       while (flag);
+    }
   return BoundaryPoints;
+        atom *Walker = NULL;
+        PointMap PointsOnBoundary;
+        LineMap LinesOnBoundary;
+        TriangleMap TrianglesOnBoundary;
+        *out << Verbose(1) << "Finding all boundary points." << endl;
+        Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+        BoundariesTestPair BoundaryTestPair;
+        Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+        double radius, angle;
+        // 3a. Go through every axis
+        for (int axis = 0; axis < NDIM; axis++)
+                {
+                        AxisVector.Zero();
+                        AngleReferenceVector.Zero();
+                        AngleReferenceNormalVector.Zero();
+                        AxisVector.x[axis] = 1.;
+                        AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+                        AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+                        //              *out << Verbose(1) << "Axisvector is ";
+                        //              AxisVector.Output(out);
+                        //              *out << " and AngleReferenceVector is ";
+                        //              AngleReferenceVector.Output(out);
+                        //              *out << "." << endl;
+                        //              *out << " and AngleReferenceNormalVector is ";
+                        //              AngleReferenceNormalVector.Output(out);
+                        //              *out << "." << endl;
+                        // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+                        Walker = mol->start;
+                        while (Walker->next != mol->end)
+                                {
+                                        Walker = Walker->next;
+                                        Vector ProjectedVector;
+                                        ProjectedVector.CopyVector(&Walker->x);
+                                        ProjectedVector.ProjectOntoPlane(&AxisVector);
+                                        // correct for negative side
+                                        //if (Projection(y) < 0)
+                                        //angle = 2.*M_PI - angle;
+                                        radius = ProjectedVector.Norm();
+                                        if (fabs(radius) > MYEPSILON)
+                                                angle = ProjectedVector.Angle(&AngleReferenceVector);
+                                        else
+                                                angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+                                        //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+                                        if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
+                                                {
+                                                        angle = 2. * M_PI - angle;
+                                                }
+                                        //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
+                                        //ProjectedVector.Output(out);
+                                        //*out << endl;
+                                        BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
+                                                        DistancePair (radius, Walker)));
+                                        if (BoundaryTestPair.second)
+                                                { // successfully inserted
+                                                }
+                                        else
+                                                { // same point exists, check first r, then distance of original vectors to center of gravity
+                                                        *out << Verbose(2)
+                                                                        << "Encountered two vectors whose projection onto axis "
+                                                                        << axis << " is equal: " << endl;
+                                                        *out << Verbose(2) << "Present vector: ";
+                                                        BoundaryTestPair.first->second.second->x.Output(out);
+                                                        *out << endl;
+                                                        *out << Verbose(2) << "New vector: ";
+                                                        Walker->x.Output(out);
+                                                        *out << endl;
+                                                        double tmp = ProjectedVector.Norm();
+                                                        if (tmp > BoundaryTestPair.first->second.first)
+                                                                {
+                                                                        BoundaryTestPair.first->second.first = tmp;
+                                                                        BoundaryTestPair.first->second.second = Walker;
+                                                                        *out << Verbose(2) << "Keeping new vector." << endl;
+                                                                }
+                                                        else if (tmp == BoundaryTestPair.first->second.first)
+                                                                {
+                                                                        if (BoundaryTestPair.first->second.second->x.ScalarProduct(
+                                                                                        &BoundaryTestPair.first->second.second->x)
+                                                                                        < Walker->x.ScalarProduct(&Walker->x))
+                                                                                { // Norm() does a sqrt, which makes it a lot slower
+                                                                                        BoundaryTestPair.first->second.second = Walker;
+                                                                                        *out << Verbose(2) << "Keeping new vector." << endl;
+                                                                                }
+                                                                        else
+                                                                                {
+                                                                                        *out << Verbose(2) << "Keeping present vector." << endl;
+                                                                                }
+                                                                }
+                                                        else
+                                                                {
+                                                                        *out << Verbose(2) << "Keeping present vector." << endl;
+                                                                }
+                                                }
+                                }
+                        // printing all inserted for debugging
+                        //              {
+                        //                      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+                        //                      int i=0;
+                        //                      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+                        //                              if (runner != BoundaryPoints[axis].begin())
+                        //                                      *out << ", " << i << ": " << *runner->second.second;
+                        //                              else
+                        //                                      *out << i << ": " << *runner->second.second;
+                        //                              i++;
+                        //                      }
+                        //                      *out << endl;
+                        //              }
+                        // 3c. throw out points whose distance is less than the mean of left and right neighbours
+                        bool flag = false;
+                        do
+                                { // do as long as we still throw one out per round
+                                        *out << Verbose(1)
+                                                        << "Looking for candidates to kick out by convex condition ... "
+                                                        << endl;
+                                        flag = false;
+                                        Boundaries::iterator left = BoundaryPoints[axis].end();
+                                        Boundaries::iterator right = BoundaryPoints[axis].end();
+                                        for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+                                                        != BoundaryPoints[axis].end(); runner++)
+                                                {
+                                                        // set neighbours correctly
+                                                        if (runner == BoundaryPoints[axis].begin())
+                                                                {
+                                                                        left = BoundaryPoints[axis].end();
+                                                                }
+                                                        else
+                                                                {
+                                                                        left = runner;
+                                                                }
+                                                        left--;
+                                                        right = runner;
+                                                        right++;
+                                                        if (right == BoundaryPoints[axis].end())
+                                                                {
+                                                                        right = BoundaryPoints[axis].begin();
+                                                                }
+                                                        // check distance
+                                                        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+                                                                {
+                                                                        Vector SideA, SideB, SideC, SideH;
+                                                                        SideA.CopyVector(&left->second.second->x);
+                                                                        SideA.ProjectOntoPlane(&AxisVector);
+                                                                        //                                      *out << "SideA: ";
+                                                                        //                                      SideA.Output(out);
+                                                                        //                                      *out << endl;
+                                                                        SideB.CopyVector(&right->second.second->x);
+                                                                        SideB.ProjectOntoPlane(&AxisVector);
+                                                                        //                                      *out << "SideB: ";
+                                                                        //                                      SideB.Output(out);
+                                                                        //                                      *out << endl;
+                                                                        SideC.CopyVector(&left->second.second->x);
+                                                                        SideC.SubtractVector(&right->second.second->x);
+                                                                        SideC.ProjectOntoPlane(&AxisVector);
+                                                                        //                                      *out << "SideC: ";
+                                                                        //                                      SideC.Output(out);
+                                                                        //                                      *out << endl;
+                                                                        SideH.CopyVector(&runner->second.second->x);
+                                                                        SideH.ProjectOntoPlane(&AxisVector);
+                                                                        //                                      *out << "SideH: ";
+                                                                        //                                      SideH.Output(out);
+                                                                        //                                      *out << endl;
+                                                                        // calculate each length
+                                                                        double a = SideA.Norm();
+                                                                        //double b = SideB.Norm();
+                                                                        //double c = SideC.Norm();
+                                                                        double h = SideH.Norm();
+                                                                        // calculate the angles
+                                                                        double alpha = SideA.Angle(&SideH);
+                                                                        double beta = SideA.Angle(&SideC);
+                                                                        double gamma = SideB.Angle(&SideH);
+                                                                        double delta = SideC.Angle(&SideH);
+                                                                        double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
+                                                                                        < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+                                                                        //                                      *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+                                                                        //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+                                                                        if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
+                                                                                        < MYEPSILON) && (h < MinDistance))
+                                                                                {
+                                                                                        // throw out point
+                                                                                        //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+                                                                                        BoundaryPoints[axis].erase(runner);
+                                                                                        flag = true;
+                                                                                }
+                                                                }
+                                                }
+                                }
+                        while (flag);
+                }
+        return BoundaryPoints;
+}
+;
 …
 double *
 GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
     bool IsAngstroem)
+{
   // get points on boundary of NULL was given as parameter
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
   if (BoundaryPoints == NULL)
+    {
       BoundaryFreeFlag = true;
       BoundaryPoints = GetBoundaryPoints(out, mol);
+    }
   else
+    {
       *out << Verbose(1) << "Using given boundary points set." << endl;
+    }
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
   double *GreatestDiameter = new double[NDIM];
   for (int i = 0; i < NDIM; i++)
     GreatestDiameter[i] = 0.;
   double OldComponent, tmp, w1, w2;
   Vector DistanceVector, OtherVector;
   int component, Othercomponent;
   for (int axis = 0; axis < NDIM; axis++)
     { // regard each projected plane
       //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
       for (int j = 0; j < 2; j++)
         { // and for both axis on the current plane
           component = (axis + j + 1) % NDIM;
           Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
           //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
           for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
               != BoundaryPoints[axis].end(); runner++)
+            {
               //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
               // seek for the neighbours pair where the Othercomponent sign flips
               Neighbour = runner;
               Neighbour++;
               if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                 Neighbour = BoundaryPoints[axis].begin();
               DistanceVector.CopyVector(&runner->second.second->x);
               DistanceVector.SubtractVector(&Neighbour->second.second->x);
               do
                 { // seek for neighbour pair where it flips
                   OldComponent = DistanceVector.x[Othercomponent];
                   Neighbour++;
                   if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                     Neighbour = BoundaryPoints[axis].begin();
                   DistanceVector.CopyVector(&runner->second.second->x);
                   DistanceVector.SubtractVector(&Neighbour->second.second->x);
                   //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+                }
               while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                   OldComponent) - DistanceVector.x[Othercomponent] / fabs(
                   DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
               if (runner != Neighbour)
+                {
                   OtherNeighbour = Neighbour;
                   if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
                     OtherNeighbour = BoundaryPoints[axis].end();
                   OtherNeighbour--;
                   //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                   // now we have found the pair: Neighbour and OtherNeighbour
                   OtherVector.CopyVector(&runner->second.second->x);
                   OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
                   //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                   //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                   // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                   w1 = fabs(OtherVector.x[Othercomponent]);
                   w2 = fabs(DistanceVector.x[Othercomponent]);
                   tmp = fabs((w1 * DistanceVector.x[component] + w2
                       * OtherVector.x[component]) / (w1 + w2));
                   // mark if it has greater diameter
                   //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                   GreatestDiameter[component] = (GreatestDiameter[component]
                       > tmp) ? GreatestDiameter[component] : tmp;
                 } //else
               //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+            }
+        }
+    }
   *out << Verbose(0) << "RESULT: The biggest diameters are "
       << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
       << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
       : "atomiclength") << "." << endl;
   // free reference lists
   if (BoundaryFreeFlag)
     delete[] (BoundaryPoints);
   return GreatestDiameter;
+}
+;
 /** Creates the objects in a raster3d file (renderable with a header.r3d)
+                bool IsAngstroem)
+{
+        // get points on boundary of NULL was given as parameter
+        bool BoundaryFreeFlag = false;
+        Boundaries *BoundaryPoints = BoundaryPtr;
+        if (BoundaryPoints == NULL)
+                {
+                        BoundaryFreeFlag = true;
+                        BoundaryPoints = GetBoundaryPoints(out, mol);
+                }
+        else
+                {
+                        *out << Verbose(1) << "Using given boundary points set." << endl;
+                }
+        // determine biggest "diameter" of cluster for each axis
+        Boundaries::iterator Neighbour, OtherNeighbour;
+        double *GreatestDiameter = new double[NDIM];
+        for (int i = 0; i < NDIM; i++)
+                GreatestDiameter[i] = 0.;
+        double OldComponent, tmp, w1, w2;
+        Vector DistanceVector, OtherVector;
+        int component, Othercomponent;
+        for (int axis = 0; axis < NDIM; axis++)
+                { // regard each projected plane
+                        //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
+                        for (int j = 0; j < 2; j++)
+                                { // and for both axis on the current plane
+                                        component = (axis + j + 1) % NDIM;
+                                        Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
+                                        //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
+                                        for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+                                                        != BoundaryPoints[axis].end(); runner++)
+                                                {
+                                                        //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
+                                                        // seek for the neighbours pair where the Othercomponent sign flips
+                                                        Neighbour = runner;
+                                                        Neighbour++;
+                                                        if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
+                                                                Neighbour = BoundaryPoints[axis].begin();
+                                                        DistanceVector.CopyVector(&runner->second.second->x);
+                                                        DistanceVector.SubtractVector(&Neighbour->second.second->x);
+                                                        do
+                                                                { // seek for neighbour pair where it flips
+                                                                        OldComponent = DistanceVector.x[Othercomponent];
+                                                                        Neighbour++;
+                                                                        if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
+                                                                                Neighbour = BoundaryPoints[axis].begin();
+                                                                        DistanceVector.CopyVector(&runner->second.second->x);
+                                                                        DistanceVector.SubtractVector(&Neighbour->second.second->x);
+                                                                        //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
+                                                                }
+                                                        while ((runner != Neighbour) && (fabs(OldComponent / fabs(
+                                                                        OldComponent) - DistanceVector.x[Othercomponent] / fabs(
+                                                                        DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
+                                                        if (runner != Neighbour)
+                                                                {
+                                                                        OtherNeighbour = Neighbour;
+                                                                        if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
+                                                                                OtherNeighbour = BoundaryPoints[axis].end();
+                                                                        OtherNeighbour--;
+                                                                        //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
+                                                                        // now we have found the pair: Neighbour and OtherNeighbour
+                                                                        OtherVector.CopyVector(&runner->second.second->x);
+                                                                        OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
+                                                                        //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
+                                                                        //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
+                                                                        // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
+                                                                        w1 = fabs(OtherVector.x[Othercomponent]);
+                                                                        w2 = fabs(DistanceVector.x[Othercomponent]);
+                                                                        tmp = fabs((w1 * DistanceVector.x[component] + w2
+                                                                                        * OtherVector.x[component]) / (w1 + w2));
+                                                                        // mark if it has greater diameter
+                                                                        //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
+                                                                        GreatestDiameter[component] = (GreatestDiameter[component]
+                                                                                        > tmp) ? GreatestDiameter[component] : tmp;
+                                                                } //else
+                                                        //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
+                                                }
+                                }
+                }
+        *out << Verbose(0) << "RESULT: The biggest diameters are "
+                        << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
+                        << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
+                        : "atomiclength") << "." << endl;
+        // free reference lists
+        if (BoundaryFreeFlag)
+                delete[] (BoundaryPoints);
+        return GreatestDiameter;
+}
+;
+/** Creates the objects in a VRML file.
  * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ * \param *vrmlfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
+{
+        atom *Walker = mol->start;
+        bond *Binder = mol->first;
+        int i;
+        Vector *center = mol->DetermineCenterOfAll(out);
+        if (vrmlfile != NULL) {
+                //cout << Verbose(1) << "Writing Raster3D file ... ";
+                *vrmlfile << "#VRML V2.0 utf8" << endl;
+                *vrmlfile << "#Created by molecuilder" << endl;
+                *vrmlfile << "#All atoms as spheres" << endl;
+                while (Walker->next != mol->end) {
+                        Walker = Walker->next;
+                        *vrmlfile << "Sphere {" << endl << "    "; // 2 is sphere type
+                        for (i=0;i<NDIM;i++)
+                                *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
+                        *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+                }
+                *vrmlfile << "# All bonds as vertices" << endl;
+                while (Binder->next != mol->last) {
+                        Binder = Binder->next;
+                        *vrmlfile << "3" << endl << "   "; // 2 is round-ended cylinder type
+                        for (i=0;i<NDIM;i++)
+                                *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+                        *vrmlfile << "\t0.03\t";
+                        for (i=0;i<NDIM;i++)
+                                *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+                        *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+                }
+                *vrmlfile << "# All tesselation triangles" << endl;
+                for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+                        *vrmlfile << "1" << endl << "   "; // 1 is triangle type
+                        for (i=0;i<3;i++) { // print each node
+                                for (int j=0;j<NDIM;j++)        // and for each node all NDIM coordinates
+                                        *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+                                *vrmlfile << "\t";
+                        }
+                        *vrmlfile << "1. 0. 0." << endl;        // red as colour
+                        *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+                }
+        } else {
+                cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+        }
+        delete(center);
+};
+/** Creates the objects in a raster3d file (renderable with a header.r3d).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
  * \param *Tess Tesselation structure with constructed triangles
  * \param *mol molecule structure with atom positions
 …
 void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
+{
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (rasterfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+    *rasterfile << "@header.r3d" << endl;
+    *rasterfile << "# All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *rasterfile << "2" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+    *rasterfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *rasterfile << "3" << endl << "  "; // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+    *rasterfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *rasterfile << "1" << endl << "  "; // 1 is triangle type
+      for (i=0;i<3;i++) { // print each node
+        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *rasterfile << "\t";
+      }
+      *rasterfile << "1. 0. 0." << endl;  // red as colour
+      *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+  }
+  delete(center);
+        atom *Walker = mol->start;
+        bond *Binder = mol->first;
+        int i;
+        Vector *center = mol->DetermineCenterOfAll(out);
+        if (rasterfile != NULL) {
+                //cout << Verbose(1) << "Writing Raster3D file ... ";
+                *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+                *rasterfile << "@header.r3d" << endl;
+                *rasterfile << "# All atoms as spheres" << endl;
+                while (Walker->next != mol->end) {
+                        Walker = Walker->next;
+                        *rasterfile << "2" << endl << " ";      // 2 is sphere type
+                        for (i=0;i<NDIM;i++)
+                                *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+                        *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+                }
+                *rasterfile << "# All bonds as vertices" << endl;
+                while (Binder->next != mol->last) {
+                        Binder = Binder->next;
+                        *rasterfile << "3" << endl << " ";      // 2 is round-ended cylinder type
+                        for (i=0;i<NDIM;i++)
+                                *rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+                        *rasterfile << "\t0.03\t";
+                        for (i=0;i<NDIM;i++)
+                                *rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+                        *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+                }
+                *rasterfile << "# All tesselation triangles" << endl;
+                *rasterfile << "8\n     25. -1.  1. 1. 1.        0.0            0 0 0 2\n       SOLID            1.0 0.0 0.0\n  BACKFACE        0.3 0.3 1.0      0 0\n";
+                for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+                        *rasterfile << "1" << endl << " ";      // 1 is triangle type
+                        for (i=0;i<3;i++) {     // print each node
+                                for (int j=0;j<NDIM;j++)        // and for each node all NDIM coordinates
+                                        *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+                                *rasterfile << "\t";
+                        }
+                        *rasterfile << "1. 0. 0." << endl;      // red as colour
+                        //*rasterfile << "18" << endl << "      0.5 0.5 0.5" << endl;   // 18 is transparency type for previous object
+                }
+                *rasterfile << "9\n     terminating special property\n";
+        } else {
+                cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+        }
+        delete(center);
 };
+/*
+ * This function creates the tecplot file, displaying the tesselation of the hull.
+/** This function creates the tecplot file, displaying the tesselation of the hull.
  * \param *out output stream for debugging
  * \param *tecplot output stream for tecplot data
 …
 void
 write_tecplot_file(ofstream *out, ofstream *tecplot,
     class Tesselation *TesselStruct, class molecule *mol, int N)
+{
   if (tecplot != NULL)
+    {
       *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
       *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
       *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N="
           << TesselStruct->PointsOnBoundaryCount << ", E="
           << TesselStruct->TrianglesOnBoundaryCount
           << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
       int *LookupList = new int[mol->AtomCount];
       for (int i = 0; i < mol->AtomCount; i++)
         LookupList[i] = -1;
       // print atom coordinates
       *out << Verbose(2) << "The following triangles were created:";
       int Counter = 1;
       atom *Walker = NULL;
       for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target
           != TesselStruct->PointsOnBoundary.end(); target++)
+        {
           Walker = target->second->node;
           LookupList[Walker->nr] = Counter++;
           *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " "
               << Walker->x.x[2] << " " << endl;
+        }
       *tecplot << endl;
       // print connectivity
       for (TriangleMap::iterator runner =
           TesselStruct->TrianglesOnBoundary.begin(); runner
           != TesselStruct->TrianglesOnBoundary.end(); runner++)
+        {
           *out << " " << runner->second->endpoints[0]->node->Name << "<->"
               << runner->second->endpoints[1]->node->Name << "<->"
               << runner->second->endpoints[2]->node->Name;
           *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " "
               << LookupList[runner->second->endpoints[1]->node->nr] << " "
               << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+        }
       delete[] (LookupList);
       *out << endl;
+    }
+                class Tesselation *TesselStruct, class molecule *mol, int N)
+{
+        if (tecplot != NULL)
+                {
+                        *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
+                        *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
+                        *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N="
+                                        << TesselStruct->PointsOnBoundaryCount << ", E="
+                                        << TesselStruct->TrianglesOnBoundaryCount
+                                        << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
+                        int *LookupList = new int[mol->AtomCount];
+                        for (int i = 0; i < mol->AtomCount; i++)
+                                LookupList[i] = -1;
+                        // print atom coordinates
+                        *out << Verbose(2) << "The following triangles were created:";
+                        int Counter = 1;
+                        atom *Walker = NULL;
+                        for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target
+                                        != TesselStruct->PointsOnBoundary.end(); target++)
+                                {
+                                        Walker = target->second->node;
+                                        LookupList[Walker->nr] = Counter++;
+                                        *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " "
+                                                        << Walker->x.x[2] << " " << endl;
+                                }
+                        *tecplot << endl;
+                        // print connectivity
+                        for (TriangleMap::iterator runner =
+                                        TesselStruct->TrianglesOnBoundary.begin(); runner
+                                        != TesselStruct->TrianglesOnBoundary.end(); runner++)
+                                {
+                                        *out << " " << runner->second->endpoints[0]->node->Name << "<->"
+                                                        << runner->second->endpoints[1]->node->Name << "<->"
+                                                        << runner->second->endpoints[2]->node->Name;
+                                        *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " "
+                                                        << LookupList[runner->second->endpoints[1]->node->nr] << " "
+                                                        << LookupList[runner->second->endpoints[2]->node->nr] << endl;
+                                }
+                        delete[] (LookupList);
+                        *out << endl;
+                }
+}
 …
  * Determines first the convex envelope, then tesselates it and calculates its volume.
  * \param *out output stream for debugging
  * \param *tecplot output stream for tecplot data
+ * \param *filename filename prefix for output of vertex data
  * \param *configuration needed for path to store convex envelope file
  * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
 …
  */
 double
+VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration,
+    Boundaries *BoundaryPtr, molecule *mol)
+{
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  atom *Walker = NULL;
+  struct Tesselation *TesselStruct = new Tesselation;
+  bool BoundaryFreeFlag = false;
+  Boundaries *BoundaryPoints = BoundaryPtr;
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G, h;
+  Vector x, y;
+  double a, b, c;
+  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
+  // 1. calculate center of gravity
+  *out << endl;
+  Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
+  // 2. translate all points into CoG
+  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
+  // 3. Find all points on the boundary
+  if (BoundaryPoints == NULL)
+    {
+      BoundaryFreeFlag = true;
+      BoundaryPoints = GetBoundaryPoints(out, mol);
+    }
+  else
+    {
+      *out << Verbose(1) << "Using given boundary points set." << endl;
+    }
+  // 4. fill the boundary point list
+  for (int axis = 0; axis < NDIM; axis++)
+    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+        != BoundaryPoints[axis].end(); runner++)
+      {
+        TesselStruct->AddPoint(runner->second.second);
+      }
+  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
+      << " points on the convex boundary." << endl;
+  // now we have the whole set of edge points in the BoundaryList
+  // listing for debugging
+  //  *out << Verbose(1) << "Listing PointsOnBoundary:";
+  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+  //    *out << " " << *runner->second;
+  //  }
+  //  *out << endl;
+  // 5a. guess starting triangle
+  TesselStruct->GuessStartingTriangle(out);
+  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
+  TesselStruct->TesselateOnBoundary(out, configuration, mol);
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
+      << " triangles with " << TesselStruct->LinesOnBoundaryCount
+      << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
+      << endl;
+  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+  *out << Verbose(1)
+      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
+      << endl;
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
+      != TesselStruct->TrianglesOnBoundary.end(); runner++)
+    { // go through every triangle, calculate volume of its pyramid with CoG as peak
+      x.CopyVector(&runner->second->endpoints[0]->node->x);
+      x.SubtractVector(&runner->second->endpoints[1]->node->x);
+      y.CopyVector(&runner->second->endpoints[0]->node->x);
+      y.SubtractVector(&runner->second->endpoints[2]->node->x);
+      a = sqrt(runner->second->endpoints[0]->node->x.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;
+  CenterOfGravity->Scale(-1);
+  Walker = mol->start;
+  while (Walker->next != mol->end)
+    {
+      Walker = Walker->next;
+      Walker->x.Translate(CenterOfGravity);
+    }
+  // 8. Store triangles in tecplot file
+  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[] (BoundaryPoints);
+  return volume;
+VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
+                Boundaries *BoundaryPtr, molecule *mol)
+{
+        bool IsAngstroem = configuration->GetIsAngstroem();
+        atom *Walker = NULL;
+        struct Tesselation *TesselStruct = new Tesselation;
+        bool BoundaryFreeFlag = false;
+        Boundaries *BoundaryPoints = BoundaryPtr;
+        double volume = 0.;
+        double PyramidVolume = 0.;
+        double G, h;
+        Vector x, y;
+        double a, b, c;
+        //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
+        // 1. calculate center of gravity
+        *out << endl;
+        Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
+        // 2. translate all points into CoG
+        *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
+        Walker = mol->start;
+        while (Walker->next != mol->end)
+                {
+                        Walker = Walker->next;
+                        Walker->x.Translate(CenterOfGravity);
+                }
+        // 3. Find all points on the boundary
+        if (BoundaryPoints == NULL)
+                {
+                        BoundaryFreeFlag = true;
+                        BoundaryPoints = GetBoundaryPoints(out, mol);
+                }
+        else
+                {
+                        *out << Verbose(1) << "Using given boundary points set." << endl;
+                }
+        // 4. fill the boundary point list
+        for (int axis = 0; axis < NDIM; axis++)
+                for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
+                                != BoundaryPoints[axis].end(); runner++)
+                        {
+                                TesselStruct->AddPoint(runner->second.second);
+                        }
+        *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount
+                        << " points on the convex boundary." << endl;
+        // now we have the whole set of edge points in the BoundaryList
+        // listing for debugging
+        //      *out << Verbose(1) << "Listing PointsOnBoundary:";
+        //      for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+        //              *out << " " << *runner->second;
+        //      }
+        //      *out << endl;
+        // 5a. guess starting triangle
+        TesselStruct->GuessStartingTriangle(out);
+        // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
+        TesselStruct->TesselateOnBoundary(out, configuration, mol);
+        *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
+                        << " triangles with " << TesselStruct->LinesOnBoundaryCount
+                        << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
+                        << endl;
+        // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
+        *out << Verbose(1)
+                        << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
+                        << endl;
+        for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
+                        != TesselStruct->TrianglesOnBoundary.end(); runner++)
+                { // go through every triangle, calculate volume of its pyramid with CoG as peak
+                        x.CopyVector(&runner->second->endpoints[0]->node->x);
+                        x.SubtractVector(&runner->second->endpoints[1]->node->x);
+                        y.CopyVector(&runner->second->endpoints[0]->node->x);
+                        y.SubtractVector(&runner->second->endpoints[2]->node->x);
+                        a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+                                        &runner->second->endpoints[1]->node->x));
+                        b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
+                                        &runner->second->endpoints[2]->node->x));
+                        c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
+                                        &runner->second->endpoints[1]->node->x));
+                        G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
+                        x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
+                                        &runner->second->endpoints[1]->node->x,
+                                        &runner->second->endpoints[2]->node->x);
+                        x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
+                        h = x.Norm(); // distance of CoG to triangle
+                        PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
+                        *out << Verbose(2) << "Area of triangle is " << G << " "
+                                        << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
+                                        << h << " and the volume is " << PyramidVolume << " "
+                                        << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+                        volume += PyramidVolume;
+                }
+        *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
+                        << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
+                        << endl;
+        // 7. translate all points back from CoG
+        *out << Verbose(1) << "Translating system back from Center of Gravity."
+                        << endl;
+        CenterOfGravity->Scale(-1);
+        Walker = mol->start;
+        while (Walker->next != mol->end)
+                {
+                        Walker = Walker->next;
+                        Walker->x.Translate(CenterOfGravity);
+                }
+        // 8. Store triangles in tecplot file
+        string OutputName(filename);
+        OutputName.append(TecplotSuffix);
+        ofstream *tecplot = new ofstream(OutputName.c_str());
+        write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+        tecplot->close();
+        delete(tecplot);
+        // free reference lists
+        if (BoundaryFreeFlag)
+                delete[] (BoundaryPoints);
+        return volume;
+}
+;
 …
 void
 PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
     double ClusterVolume, double celldensity)
+{
   // transform to PAS
   mol->PrincipalAxisSystem(out, true);
   // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
   Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
   double clustervolume;
   if (ClusterVolume == 0)
     clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
         BoundaryPoints, mol);
   else
     clustervolume = ClusterVolume;
   double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
       IsAngstroem);
   Vector BoxLengths;
   int repetition[NDIM] =
     { 1, 1, 1 };
   int TotalNoClusters = 1;
   for (int i = 0; i < NDIM; i++)
     TotalNoClusters *= repetition[i];
   // sum up the atomic masses
   double totalmass = 0.;
   atom *Walker = mol->start;
   while (Walker->next != mol->end)
+    {
       Walker = Walker->next;
       totalmass += Walker->type->mass;
+    }
   *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
       << totalmass << " atomicmassunit." << endl;
   *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
       << totalmass / clustervolume << " atomicmassunit/"
       << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
   // solve cubic polynomial
   *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
       << endl;
   double cellvolume;
   if (IsAngstroem)
     cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
         / clustervolume)) / (celldensity - 1);
   else
     cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
         / clustervolume)) / (celldensity - 1);
   *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
       << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
       : "atomiclength") << "^3." << endl;
   double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
       * GreatestDiameter[1] * GreatestDiameter[2]);
   *out << Verbose(1)
       << "Minimum volume of the convex envelope contained in a rectangular box is "
       << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
       : "atomiclength") << "^3." << endl;
   if (minimumvolume > cellvolume)
+    {
       cerr << Verbose(0)
           << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
           << endl;
       cout << Verbose(0)
           << "Setting Box dimensions to minimum possible, the greatest diameters."
           << endl;
       for (int i = 0; i < NDIM; i++)
         BoxLengths.x[i] = GreatestDiameter[i];
       mol->CenterEdge(out, &BoxLengths);
+    }
   else
+    {
       BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
           * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
       BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
           * GreatestDiameter[1] + repetition[0] * repetition[2]
           * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
           * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
       BoxLengths.x[2] = minimumvolume - cellvolume;
       double x0 = 0., x1 = 0., x2 = 0.;
       if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
           BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
         *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
             << " ." << endl;
       else
+        {
           *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
               << " and " << x1 << " and " << x2 << " ." << endl;
           x0 = x2; // sorted in ascending order
+        }
       cellvolume = 1;
       for (int i = 0; i < NDIM; i++)
+        {
           BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
           cellvolume *= BoxLengths.x[i];
+        }
       // set new box dimensions
       *out << Verbose(0) << "Translating to box with these boundaries." << endl;
       mol->CenterInBox((ofstream *) &cout, &BoxLengths);
+    }
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
   *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
       << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
       << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
       << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+                double ClusterVolume, double celldensity)
+{
+        // transform to PAS
+        mol->PrincipalAxisSystem(out, true);
+        // some preparations beforehand
+        bool IsAngstroem = configuration->GetIsAngstroem();
+        Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
+        double clustervolume;
+        if (ClusterVolume == 0)
+                clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration,
+                                BoundaryPoints, mol);
+        else
+                clustervolume = ClusterVolume;
+        double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
+                        IsAngstroem);
+        Vector BoxLengths;
+        int repetition[NDIM] =
+                { 1, 1, 1 };
+        int TotalNoClusters = 1;
+        for (int i = 0; i < NDIM; i++)
+                TotalNoClusters *= repetition[i];
+        // sum up the atomic masses
+        double totalmass = 0.;
+        atom *Walker = mol->start;
+        while (Walker->next != mol->end)
+                {
+                        Walker = Walker->next;
+                        totalmass += Walker->type->mass;
+                }
+        *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
+                        << totalmass << " atomicmassunit." << endl;
+        *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
+                        << totalmass / clustervolume << " atomicmassunit/"
+                        << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+        // solve cubic polynomial
+        *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
+                        << endl;
+        double cellvolume;
+        if (IsAngstroem)
+                cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
+                                / clustervolume)) / (celldensity - 1);
+        else
+                cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
+                                / clustervolume)) / (celldensity - 1);
+        *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
+                        << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
+                        : "atomiclength") << "^3." << endl;
+        double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
+                        * GreatestDiameter[1] * GreatestDiameter[2]);
+        *out << Verbose(1)
+                        << "Minimum volume of the convex envelope contained in a rectangular box is "
+                        << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
+                        : "atomiclength") << "^3." << endl;
+        if (minimumvolume > cellvolume)
+                {
+                        cerr << Verbose(0)
+                                        << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
+                                        << endl;
+                        cout << Verbose(0)
+                                        << "Setting Box dimensions to minimum possible, the greatest diameters."
+                                        << endl;
+                        for (int i = 0; i < NDIM; i++)
+                                BoxLengths.x[i] = GreatestDiameter[i];
+                        mol->CenterEdge(out, &BoxLengths);
+                }
+        else
+                {
+                        BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
+                                        * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
+                        BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
+                                        * GreatestDiameter[1] + repetition[0] * repetition[2]
+                                        * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
+                                        * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
+                        BoxLengths.x[2] = minimumvolume - cellvolume;
+                        double x0 = 0., x1 = 0., x2 = 0.;
+                        if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
+                                        BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
+                                *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
+                                                << " ." << endl;
+                        else
+                                {
+                                        *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
+                                                        << " and " << x1 << " and " << x2 << " ." << endl;
+                                        x0 = x2; // sorted in ascending order
+                                }
+                        cellvolume = 1;
+                        for (int i = 0; i < NDIM; i++)
+                                {
+                                        BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+                                        cellvolume *= BoxLengths.x[i];
+                                }
+                        // set new box dimensions
+                        *out << Verbose(0) << "Translating to box with these boundaries." << endl;
+                        mol->CenterInBox((ofstream *) &cout, &BoxLengths);
+                }
+        // update Box of atoms by boundary
+        mol->SetBoxDimension(&BoxLengths);
+        *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
+                        << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
+                        << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
+                        << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+}
+;
 …
 Tesselation::Tesselation()
+{
   PointsOnBoundaryCount = 0;
   LinesOnBoundaryCount = 0;
   TrianglesOnBoundaryCount = 0;
   TriangleFilesWritten = 0;
+        PointsOnBoundaryCount = 0;
+        LinesOnBoundaryCount = 0;
+        TrianglesOnBoundaryCount = 0;
+        TriangleFilesWritten = 0;
+}
+;
 …
 Tesselation::~Tesselation()
+{
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner
+      != TrianglesOnBoundary.end(); runner++)
+    {
+      delete (runner->second);
+    }
+        cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+        for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+                if (runner->second != NULL) {
+                        delete (runner->second);
+                        runner->second = NULL;
+                } else
+                        cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+        }
+}
+;
 …
 Tesselation::GuessStartingTriangle(ofstream *out)
+{
   // 4b. create a starting triangle
   // 4b1. create all distances
   DistanceMultiMap DistanceMMap;
   double distance, tmp;
   Vector PlaneVector, TrialVector;
   PointMap::iterator A, B, C; // three nodes of the first triangle
   A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
   // with A chosen, take each pair B,C and sort
   if (A != PointsOnBoundary.end())
+    {
       B = A;
       B++;
       for (; B != PointsOnBoundary.end(); B++)
+        {
           C = B;
           C++;
           for (; C != PointsOnBoundary.end(); C++)
+            {
               tmp = A->second->node->x.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);
+    }
+        // 4b. create a starting triangle
+        // 4b1. create all distances
+        DistanceMultiMap DistanceMMap;
+        double distance, tmp;
+        Vector PlaneVector, TrialVector;
+        PointMap::iterator A, B, C; // three nodes of the first triangle
+        A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
+        // with A chosen, take each pair B,C and sort
+        if (A != PointsOnBoundary.end())
+                {
+                        B = A;
+                        B++;
+                        for (; B != PointsOnBoundary.end(); B++)
+                                {
+                                        C = B;
+                                        C++;
+                                        for (; C != PointsOnBoundary.end(); C++)
+                                                {
+                                                        tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
+                                                        distance = tmp * tmp;
+                                                        tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
+                                                        distance += tmp * tmp;
+                                                        tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
+                                                        distance += tmp * tmp;
+                                                        DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
+                                                                        PointMap::iterator, PointMap::iterator> (B, C)));
+                                                }
+                                }
+                }
+        //              // listing distances
+        //              *out << Verbose(1) << "Listing DistanceMMap:";
+        //              for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
+        //                      *out << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
+        //              }
+        //              *out << endl;
+        // 4b2. pick three baselines forming a triangle
+        // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+        DistanceMultiMap::iterator baseline = DistanceMMap.begin();
+        for (; baseline != DistanceMMap.end(); baseline++)
+                {
+                        // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
+                        // 2. next, we have to check whether all points reside on only one side of the triangle
+                        // 3. construct plane vector
+                        PlaneVector.MakeNormalVector(&A->second->node->x,
+                                        &baseline->second.first->second->node->x,
+                                        &baseline->second.second->second->node->x);
+                        *out << Verbose(2) << "Plane vector of candidate triangle is ";
+                        PlaneVector.Output(out);
+                        *out << endl;
+                        // 4. loop over all points
+                        double sign = 0.;
+                        PointMap::iterator checker = PointsOnBoundary.begin();
+                        for (; checker != PointsOnBoundary.end(); checker++)
+                                {
+                                        // (neglecting A,B,C)
+                                        if ((checker == A) || (checker == baseline->second.first) || (checker
+                                                        == baseline->second.second))
+                                                continue;
+                                        // 4a. project onto plane vector
+                                        TrialVector.CopyVector(&checker->second->node->x);
+                                        TrialVector.SubtractVector(&A->second->node->x);
+                                        distance = TrialVector.Projection(&PlaneVector);
+                                        if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
+                                                continue;
+                                        *out << Verbose(3) << "Projection of " << checker->second->node->Name
+                                                        << " yields distance of " << distance << "." << endl;
+                                        tmp = distance / fabs(distance);
+                                        // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
+                                        if ((sign != 0) && (tmp != sign))
+                                                {
+                                                        // 4c. If so, break 4. loop and continue with next candidate in 1. loop
+                                                        *out << Verbose(2) << "Current candidates: "
+                                                                        << A->second->node->Name << ","
+                                                                        << baseline->second.first->second->node->Name << ","
+                                                                        << baseline->second.second->second->node->Name << " leave "
+                                                                        << checker->second->node->Name << " outside the convex hull."
+                                                                        << endl;
+                                                        break;
+                                                }
+                                        else
+                                                { // note the sign for later
+                                                        *out << Verbose(2) << "Current candidates: "
+                                                                        << A->second->node->Name << ","
+                                                                        << baseline->second.first->second->node->Name << ","
+                                                                        << baseline->second.second->second->node->Name << " leave "
+                                                                        << checker->second->node->Name << " inside the convex hull."
+                                                                        << endl;
+                                                        sign = tmp;
+                                                }
+                                        // 4d. Check whether the point is inside the triangle (check distance to each node
+                                        tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
+                                        int innerpoint = 0;
+                                        if ((tmp < A->second->node->x.DistanceSquared(
+                                                        &baseline->second.first->second->node->x)) && (tmp
+                                                        < A->second->node->x.DistanceSquared(
+                                                                        &baseline->second.second->second->node->x)))
+                                                innerpoint++;
+                                        tmp = checker->second->node->x.DistanceSquared(
+                                                        &baseline->second.first->second->node->x);
+                                        if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
+                                                        &A->second->node->x)) && (tmp
+                                                        < baseline->second.first->second->node->x.DistanceSquared(
+                                                                        &baseline->second.second->second->node->x)))
+                                                innerpoint++;
+                                        tmp = checker->second->node->x.DistanceSquared(
+                                                        &baseline->second.second->second->node->x);
+                                        if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
+                                                        &baseline->second.first->second->node->x)) && (tmp
+                                                        < baseline->second.second->second->node->x.DistanceSquared(
+                                                                        &A->second->node->x)))
+                                                innerpoint++;
+                                        // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+                                        if (innerpoint == 3)
+                                                break;
+                                }
+                        // 5. come this far, all on same side? Then break 1. loop and construct triangle
+                        if (checker == PointsOnBoundary.end())
+                                {
+                                        *out << "Looks like we have a candidate!" << endl;
+                                        break;
+                                }
+                }
+        if (baseline != DistanceMMap.end())
+                {
+                        BPS[0] = baseline->second.first->second;
+                        BPS[1] = baseline->second.second->second;
+                        BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+                        BPS[0] = A->second;
+                        BPS[1] = baseline->second.second->second;
+                        BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+                        BPS[0] = baseline->second.first->second;
+                        BPS[1] = A->second;
+                        BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+                        // 4b3. insert created triangle
+                        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+                        TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+                        TrianglesOnBoundaryCount++;
+                        for (int i = 0; i < NDIM; i++)
+                                {
+                                        LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+                                        LinesOnBoundaryCount++;
+                                }
+                        *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+                }
+        else
+                {
+                        *out << Verbose(1) << "No starting triangle found." << endl;
+                        exit(255);
+                }
+}
+;
 …
  * -# if the lines contains to only one triangle
  * -# We search all points in the boundary
  *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
  *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
  *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
  *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
+ *              -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors
+ *              -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
+ *                       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
+ *              -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
  * \param *out output stream for debugging
  * \param *configuration for IsAngstroem
 …
 void
 Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
     molecule *mol)
+{
   bool flag;
   PointMap::iterator winner;
   class BoundaryPointSet *peak = NULL;
   double SmallestAngle, TempAngle;
   Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
       PropagationVector;
   LineMap::iterator LineChecker[2];
   do
+    {
       flag = false;
       for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
           != LinesOnBoundary.end(); baseline++)
         if (baseline->second->TrianglesCount == 1)
+          {
             *out << Verbose(2) << "Current baseline is between "
                 << *(baseline->second) << "." << endl;
             // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
             SmallestAngle = M_PI;
             BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
             // get peak point with respect to this base line's only triangle
             for (int i = 0; i < 3; i++)
               if ((BTS->endpoints[i] != baseline->second->endpoints[0])
                   && (BTS->endpoints[i] != baseline->second->endpoints[1]))
                 peak = BTS->endpoints[i];
             *out << Verbose(3) << " and has peak " << *peak << "." << endl;
             // normal vector of triangle
             BTS->GetNormalVector(NormalVector);
             *out << Verbose(4) << "NormalVector of base triangle is ";
             NormalVector.Output(out);
             *out << endl;
             // offset to center of triangle
             CenterVector.Zero();
             for (int i = 0; i < 3; i++)
               CenterVector.AddVector(&BTS->endpoints[i]->node->x);
             CenterVector.Scale(1. / 3.);
             *out << Verbose(4) << "CenterVector of base triangle is ";
             CenterVector.Output(out);
             *out << endl;
             // vector in propagation direction (out of triangle)
             // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
             TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
             TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
             PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
             TempVector.CopyVector(&CenterVector);
             TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
             //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
             if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
               PropagationVector.Scale(-1.);
             *out << Verbose(4) << "PropagationVector of base triangle is ";
             PropagationVector.Output(out);
             *out << endl;
             winner = PointsOnBoundary.end();
             for (PointMap::iterator target = PointsOnBoundary.begin(); target
                 != PointsOnBoundary.end(); target++)
               if ((target->second != baseline->second->endpoints[0])
                   && (target->second != baseline->second->endpoints[1]))
                 { // don't take the same endpoints
                   *out << Verbose(3) << "Target point is " << *(target->second)
                       << ":";
                   bool continueflag = true;
                   VirtualNormalVector.CopyVector(
                       &baseline->second->endpoints[0]->node->x);
                   VirtualNormalVector.AddVector(
                       &baseline->second->endpoints[0]->node->x);
                   VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
                   VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
                   TempAngle = VirtualNormalVector.Angle(&PropagationVector);
                   continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
                   if (!continueflag)
+                    {
                       *out << Verbose(4)
                           << "Angle between propagation direction and base line to "
                           << *(target->second) << " is " << TempAngle
                           << ", bad direction!" << endl;
                       continue;
+                    }
                   else
                     *out << Verbose(4)
                         << "Angle between propagation direction and base line to "
                         << *(target->second) << " is " << TempAngle
                         << ", good direction!" << endl;
                   LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                       target->first);
                   LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                       target->first);
                   //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
                   //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
                   //            else
                   //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                   //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
                   //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
                   //            else
                   //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                   // check first endpoint (if any connecting line goes to target or at least not more than 1)
                   continueflag = continueflag && (((LineChecker[0]
                       == baseline->second->endpoints[0]->lines.end())
                       || (LineChecker[0]->second->TrianglesCount == 1)));
                   if (!continueflag)
+                    {
                       *out << Verbose(4) << *(baseline->second->endpoints[0])
                           << " has line " << *(LineChecker[0]->second)
                           << " to " << *(target->second)
                           << " as endpoint with "
                           << LineChecker[0]->second->TrianglesCount
                           << " triangles." << endl;
                       continue;
+                    }
                   // check second endpoint (if any connecting line goes to target or at least not more than 1)
                   continueflag = continueflag && (((LineChecker[1]
                       == baseline->second->endpoints[1]->lines.end())
                       || (LineChecker[1]->second->TrianglesCount == 1)));
                   if (!continueflag)
+                    {
                       *out << Verbose(4) << *(baseline->second->endpoints[1])
                           << " has line " << *(LineChecker[1]->second)
                           << " to " << *(target->second)
                           << " as endpoint with "
                           << LineChecker[1]->second->TrianglesCount
                           << " triangles." << endl;
                       continue;
+                    }
                   // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
                   continueflag = continueflag && (!(((LineChecker[0]
                       != baseline->second->endpoints[0]->lines.end())
                       && (LineChecker[1]
                           != baseline->second->endpoints[1]->lines.end())
                       && (GetCommonEndpoint(LineChecker[0]->second,
                           LineChecker[1]->second) == peak))));
                   if (!continueflag)
+                    {
                       *out << Verbose(4) << "Current target is peak!" << endl;
                       continue;
+                    }
                   // in case NOT both were found
                   if (continueflag)
                     { // create virtually this triangle, get its normal vector, calculate angle
                       flag = true;
                       VirtualNormalVector.MakeNormalVector(
                           &baseline->second->endpoints[0]->node->x,
                           &baseline->second->endpoints[1]->node->x,
                           &target->second->node->x);
                       // make it always point inward
                       if (baseline->second->endpoints[0]->node->x.Projection(
                           &VirtualNormalVector) > 0)
                         VirtualNormalVector.Scale(-1.);
                       // calculate angle
                       TempAngle = NormalVector.Angle(&VirtualNormalVector);
                       *out << Verbose(4) << "NormalVector is ";
                       VirtualNormalVector.Output(out);
                       *out << " and the angle is " << TempAngle << "." << endl;
                       if (SmallestAngle > TempAngle)
                         { // set to new possible winner
                           SmallestAngle = TempAngle;
                           winner = target;
+                        }
+                    }
+                }
             // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
             if (winner != PointsOnBoundary.end())
+              {
                 *out << Verbose(2) << "Winning target point is "
                     << *(winner->second) << " with angle " << SmallestAngle
                     << "." << endl;
                 // create the lins of not yet present
                 BLS[0] = baseline->second;
                 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
                 LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                     winner->first);
                 LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                     winner->first);
                 if (LineChecker[0]
                     == baseline->second->endpoints[0]->lines.end())
                   { // create
                     BPS[0] = baseline->second->endpoints[0];
                     BPS[1] = winner->second;
                     BLS[1] = new class BoundaryLineSet(BPS,
                         LinesOnBoundaryCount);
                     LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                         BLS[1]));
                     LinesOnBoundaryCount++;
+                  }
                 else
                   BLS[1] = LineChecker[0]->second;
                 if (LineChecker[1]
                     == baseline->second->endpoints[1]->lines.end())
                   { // create
                     BPS[0] = baseline->second->endpoints[1];
                     BPS[1] = winner->second;
                     BLS[2] = new class BoundaryLineSet(BPS,
                         LinesOnBoundaryCount);
                     LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                         BLS[2]));
                     LinesOnBoundaryCount++;
+                  }
                 else
                   BLS[2] = LineChecker[1]->second;
                 BTS = new class BoundaryTriangleSet(BLS,
                     TrianglesOnBoundaryCount);
                 TrianglesOnBoundary.insert(TrianglePair(
                     TrianglesOnBoundaryCount, BTS));
                 TrianglesOnBoundaryCount++;
+              }
             else
+              {
                 *out << Verbose(1)
                     << "I could not determine a winner for this baseline "
                     << *(baseline->second) << "." << endl;
+              }
             // 5d. If the set of lines is not yet empty, go to 5. and continue
+          }
         else
           *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
               << " has a triangle count of "
               << baseline->second->TrianglesCount << "." << endl;
+    }
   while (flag);
+                molecule *mol)
+{
+        bool flag;
+        PointMap::iterator winner;
+        class BoundaryPointSet *peak = NULL;
+        double SmallestAngle, TempAngle;
+        Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
+                        PropagationVector;
+        LineMap::iterator LineChecker[2];
+        do
+                {
+                        flag = false;
+                        for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
+                                        != LinesOnBoundary.end(); baseline++)
+                                if (baseline->second->TrianglesCount == 1)
+                                        {
+                                                *out << Verbose(2) << "Current baseline is between "
+                                                                << *(baseline->second) << "." << endl;
+                                                // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
+                                                SmallestAngle = M_PI;
+                                                BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+                                                // get peak point with respect to this base line's only triangle
+                                                for (int i = 0; i < 3; i++)
+                                                        if ((BTS->endpoints[i] != baseline->second->endpoints[0])
+                                                                        && (BTS->endpoints[i] != baseline->second->endpoints[1]))
+                                                                peak = BTS->endpoints[i];
+                                                *out << Verbose(3) << " and has peak " << *peak << "." << endl;
+                                                // normal vector of triangle
+                                                BTS->GetNormalVector(NormalVector);
+                                                *out << Verbose(4) << "NormalVector of base triangle is ";
+                                                NormalVector.Output(out);
+                                                *out << endl;
+                                                // offset to center of triangle
+                                                CenterVector.Zero();
+                                                for (int i = 0; i < 3; i++)
+                                                        CenterVector.AddVector(&BTS->endpoints[i]->node->x);
+                                                CenterVector.Scale(1. / 3.);
+                                                *out << Verbose(4) << "CenterVector of base triangle is ";
+                                                CenterVector.Output(out);
+                                                *out << endl;
+                                                // vector in propagation direction (out of triangle)
+                                                // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+                                                TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+                                                TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
+                                                PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
+                                                TempVector.CopyVector(&CenterVector);
+                                                TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+                                                //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
+                                                if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+                                                        PropagationVector.Scale(-1.);
+                                                *out << Verbose(4) << "PropagationVector of base triangle is ";
+                                                PropagationVector.Output(out);
+                                                *out << endl;
+                                                winner = PointsOnBoundary.end();
+                                                for (PointMap::iterator target = PointsOnBoundary.begin(); target
+                                                                != PointsOnBoundary.end(); target++)
+                                                        if ((target->second != baseline->second->endpoints[0])
+                                                                        && (target->second != baseline->second->endpoints[1]))
+                                                                { // don't take the same endpoints
+                                                                        *out << Verbose(3) << "Target point is " << *(target->second)
+                                                                                        << ":";
+                                                                        bool continueflag = true;
+                                                                        VirtualNormalVector.CopyVector(
+                                                                                        &baseline->second->endpoints[0]->node->x);
+                                                                        VirtualNormalVector.AddVector(
+                                                                                        &baseline->second->endpoints[0]->node->x);
+                                                                        VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
+                                                                        VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
+                                                                        TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+                                                                        continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
+                                                                        if (!continueflag)
+                                                                                {
+                                                                                        *out << Verbose(4)
+                                                                                                        << "Angle between propagation direction and base line to "
+                                                                                                        << *(target->second) << " is " << TempAngle
+                                                                                                        << ", bad direction!" << endl;
+                                                                                        continue;
+                                                                                }
+                                                                        else
+                                                                                *out << Verbose(4)
+                                                                                                << "Angle between propagation direction and base line to "
+                                                                                                << *(target->second) << " is " << TempAngle
+                                                                                                << ", good direction!" << endl;
+                                                                        LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                                                                                        target->first);
+                                                                        LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                                                                                        target->first);
+                                                                        //                                              if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
+                                                                        //                                                      *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+                                                                        //                                              else
+                                                                        //                                                      *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                                                                        //                                              if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
+                                                                        //                                                      *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+                                                                        //                                              else
+                                                                        //                                                      *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
+                                                                        // check first endpoint (if any connecting line goes to target or at least not more than 1)
+                                                                        continueflag = continueflag && (((LineChecker[0]
+                                                                                        == baseline->second->endpoints[0]->lines.end())
+                                                                                        || (LineChecker[0]->second->TrianglesCount == 1)));
+                                                                        if (!continueflag)
+                                                                                {
+                                                                                        *out << Verbose(4) << *(baseline->second->endpoints[0])
+                                                                                                        << " has line " << *(LineChecker[0]->second)
+                                                                                                        << " to " << *(target->second)
+                                                                                                        << " as endpoint with "
+                                                                                                        << LineChecker[0]->second->TrianglesCount
+                                                                                                        << " triangles." << endl;
+                                                                                        continue;
+                                                                                }
+                                                                        // check second endpoint (if any connecting line goes to target or at least not more than 1)
+                                                                        continueflag = continueflag && (((LineChecker[1]
+                                                                                        == baseline->second->endpoints[1]->lines.end())
+                                                                                        || (LineChecker[1]->second->TrianglesCount == 1)));
+                                                                        if (!continueflag)
+                                                                                {
+                                                                                        *out << Verbose(4) << *(baseline->second->endpoints[1])
+                                                                                                        << " has line " << *(LineChecker[1]->second)
+                                                                                                        << " to " << *(target->second)
+                                                                                                        << " as endpoint with "
+                                                                                                        << LineChecker[1]->second->TrianglesCount
+                                                                                                        << " triangles." << endl;
+                                                                                        continue;
+                                                                                }
+                                                                        // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
+                                                                        continueflag = continueflag && (!(((LineChecker[0]
+                                                                                        != baseline->second->endpoints[0]->lines.end())
+                                                                                        && (LineChecker[1]
+                                                                                                        != baseline->second->endpoints[1]->lines.end())
+                                                                                        && (GetCommonEndpoint(LineChecker[0]->second,
+                                                                                                        LineChecker[1]->second) == peak))));
+                                                                        if (!continueflag)
+                                                                                {
+                                                                                        *out << Verbose(4) << "Current target is peak!" << endl;
+                                                                                        continue;
+                                                                                }
+                                                                        // in case NOT both were found
+                                                                        if (continueflag)
+                                                                                { // create virtually this triangle, get its normal vector, calculate angle
+                                                                                        flag = true;
+                                                                                        VirtualNormalVector.MakeNormalVector(
+                                                                                                        &baseline->second->endpoints[0]->node->x,
+                                                                                                        &baseline->second->endpoints[1]->node->x,
+                                                                                                        &target->second->node->x);
+                                                                                        // make it always point inward
+                                                                                        if (baseline->second->endpoints[0]->node->x.Projection(
+                                                                                                        &VirtualNormalVector) > 0)
+                                                                                                VirtualNormalVector.Scale(-1.);
+                                                                                        // calculate angle
+                                                                                        TempAngle = NormalVector.Angle(&VirtualNormalVector);
+                                                                                        *out << Verbose(4) << "NormalVector is ";
+                                                                                        VirtualNormalVector.Output(out);
+                                                                                        *out << " and the angle is " << TempAngle << "." << endl;
+                                                                                        if (SmallestAngle > TempAngle)
+                                                                                                { // set to new possible winner
+                                                                                                        SmallestAngle = TempAngle;
+                                                                                                        winner = target;
+                                                                                                }
+                                                                                }
+                                                                }
+                                                // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
+                                                if (winner != PointsOnBoundary.end())
+                                                        {
+                                                                *out << Verbose(2) << "Winning target point is "
+                                                                                << *(winner->second) << " with angle " << SmallestAngle
+                                                                                << "." << endl;
+                                                                // create the lins of not yet present
+                                                                BLS[0] = baseline->second;
+                                                                // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
+                                                                LineChecker[0] = baseline->second->endpoints[0]->lines.find(
+                                                                                winner->first);
+                                                                LineChecker[1] = baseline->second->endpoints[1]->lines.find(
+                                                                                winner->first);
+                                                                if (LineChecker[0]
+                                                                                == baseline->second->endpoints[0]->lines.end())
+                                                                        { // create
+                                                                                BPS[0] = baseline->second->endpoints[0];
+                                                                                BPS[1] = winner->second;
+                                                                                BLS[1] = new class BoundaryLineSet(BPS,
+                                                                                                LinesOnBoundaryCount);
+                                                                                LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                                                                                                BLS[1]));
+                                                                                LinesOnBoundaryCount++;
+                                                                        }
+                                                                else
+                                                                        BLS[1] = LineChecker[0]->second;
+                                                                if (LineChecker[1]
+                                                                                == baseline->second->endpoints[1]->lines.end())
+                                                                        { // create
+                                                                                BPS[0] = baseline->second->endpoints[1];
+                                                                                BPS[1] = winner->second;
+                                                                                BLS[2] = new class BoundaryLineSet(BPS,
+                                                                                                LinesOnBoundaryCount);
+                                                                                LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
+                                                                                                BLS[2]));
+                                                                                LinesOnBoundaryCount++;
+                                                                        }
+                                                                else
+                                                                        BLS[2] = LineChecker[1]->second;
+                                                                BTS = new class BoundaryTriangleSet(BLS,
+                                                                                TrianglesOnBoundaryCount);
+                                                                TrianglesOnBoundary.insert(TrianglePair(
+                                                                                TrianglesOnBoundaryCount, BTS));
+                                                                TrianglesOnBoundaryCount++;
+                                                        }
+                                                else
+                                                        {
+                                                                *out << Verbose(1)
+                                                                                << "I could not determine a winner for this baseline "
+                                                                                << *(baseline->second) << "." << endl;
+                                                        }
+                                                // 5d. If the set of lines is not yet empty, go to 5. and continue
+                                        }
+                                else
+                                        *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
+                                                        << " has a triangle count of "
+                                                        << baseline->second->TrianglesCount << "." << endl;
+                }
+        while (flag);
+}
 …
 Tesselation::AddPoint(atom *Walker)
+{
   PointTestPair InsertUnique;
   BPS[0] = new class BoundaryPointSet(Walker);
   InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
   if (InsertUnique.second) // if new point was not present before, increase counter
     PointsOnBoundaryCount++;
+        PointTestPair InsertUnique;
+        BPS[0] = new class BoundaryPointSet(Walker);
+        InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
+        if (InsertUnique.second) // if new point was not present before, increase counter
+                PointsOnBoundaryCount++;
+}
+;
 …
 Tesselation::AddTrianglePoint(atom* Candidate, int n)
+{
   PointTestPair InsertUnique;
   TPS[n] = new class BoundaryPointSet(Candidate);
   InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
   if (InsertUnique.second) // if new point was not present before, increase counter
+    {
       PointsOnBoundaryCount++;
+    }
   else
+    {
       delete TPS[n];
       cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)
           << " gibt's schon in der PointMap." << endl;
       TPS[n] = (InsertUnique.first)->second;
+    }
+        PointTestPair InsertUnique;
+        TPS[n] = new class BoundaryPointSet(Candidate);
+        InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
+        if (InsertUnique.second) // if new point was not present before, increase counter
+                {
+                        PointsOnBoundaryCount++;
+                }
+        else
+                {
+                        delete TPS[n];
+                        cout << Verbose(2) << "Atom " << *((InsertUnique.first)->second->node)
+                                        << " gibt's schon in der PointMap." << endl;
+                        TPS[n] = (InsertUnique.first)->second;
+                }
+}
+;
 …
 void
 Tesselation::AddTriangleLine(class BoundaryPointSet *a,
     class BoundaryPointSet *b, int n)
+{
   LineMap::iterator LineWalker;
   //cout << "Manually checking endpoints for line." << endl;
   if ((a->lines.find(b->node->nr))->first == b->node->nr)
   //If a line is there, how do I recognize that beyond a shadow of a doubt?
+    {
       //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl;
       LineWalker = LinesOnBoundary.end();
       LineWalker--;
       while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr,
           b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max(
           a->node->nr, b->node->nr))
+        {
           //cout << Verbose(1) << "Looking for line which already exists"<< endl;
           LineWalker--;
+        }
       BPS[0] = LineWalker->second->endpoints[0];
       BPS[1] = LineWalker->second->endpoints[1];
       BLS[n] = LineWalker->second;
+    }
   else
+    {
       cout << Verbose(2)
           << "Adding line which has not been used before between "
           << *(a->node) << " and " << *(b->node) << "." << endl;
       BPS[0] = a;
       BPS[1] = b;
       BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
       LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
       LinesOnBoundaryCount++;
+    }
+                class BoundaryPointSet *b, int n)
+{
+        LineMap::iterator LineWalker;
+        //cout << "Manually checking endpoints for line." << endl;
+        if ((a->lines.find(b->node->nr)) != a->lines.end()) // ->first == b->node->nr)
+        //If a line is there, how do I recognize that beyond a shadow of a doubt?
+                {
+                        //cout << Verbose(2) << "Line exists already, retrieving it from LinesOnBoundarySet" << endl;
+                        LineWalker = LinesOnBoundary.end();
+                        LineWalker--;
+                        while (LineWalker->second->endpoints[0]->node->nr != min(a->node->nr,
+                                        b->node->nr) or LineWalker->second->endpoints[1]->node->nr != max(
+                                        a->node->nr, b->node->nr))
+                                {
+                                        //cout << Verbose(1) << "Looking for line which already exists"<< endl;
+                                        LineWalker--;
+                                }
+                        BPS[0] = LineWalker->second->endpoints[0];
+                        BPS[1] = LineWalker->second->endpoints[1];
+                        BLS[n] = LineWalker->second;
+                }
+        else
+                {
+                        cout << Verbose(2)
+                                        << "Adding line which has not been used before between "
+                                        << *(a->node) << " and " << *(b->node) << "." << endl;
+                        BPS[0] = a;
+                        BPS[1] = b;
+                        BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+                        LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
+                        LinesOnBoundaryCount++;
+                }
+}
+;
 …
+{
+  cout << Verbose(1) << "Adding triangle to its lines" << endl;
+  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);
+   }
+   */
+}
+;
+        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);
+         }
+         */
+}
+;
+double det_get(gsl_matrix *A, int inPlace) {
+        /*
+        inPlace = 1 => A is replaced with the LU decomposed copy.
+        inPlace = 0 => A is retained, and a copy is used for LU.
+        */
+        double det;
+        int signum;
+        gsl_permutation *p = gsl_permutation_alloc(A->size1);
+        gsl_matrix *tmpA;
+        if (inPlace)
+        tmpA = A;
+        else {
+        gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+        gsl_matrix_memcpy(tmpA , A);
+        }
+        gsl_linalg_LU_decomp(tmpA , p , &signum);
+        det = gsl_linalg_LU_det(tmpA , signum);
+        gsl_permutation_free(p);
+        if (! inPlace)
+        gsl_matrix_free(tmpA);
+        return det;
+};
+void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
+        gsl_matrix *A = gsl_matrix_calloc(3,3);
+        double m11, m12, m13, m14;
+        for(int i=0;i<3;i++) {
+                gsl_matrix_set(A, i, 0, a.x[i]);
+                gsl_matrix_set(A, i, 1, b.x[i]);
+                gsl_matrix_set(A, i, 2, c.x[i]);
+        }
+        m11 = det_get(A, 1);
+        for(int i=0;i<3;i++) {
+                gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+                gsl_matrix_set(A, i, 1, b.x[i]);
+                gsl_matrix_set(A, i, 2, c.x[i]);
+        }
+        m12 = det_get(A, 1);
+        for(int i=0;i<3;i++) {
+                gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+                gsl_matrix_set(A, i, 1, a.x[i]);
+                gsl_matrix_set(A, i, 2, c.x[i]);
+        }
+        m13 = det_get(A, 1);
+        for(int i=0;i<3;i++) {
+                gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+                gsl_matrix_set(A, i, 1, a.x[i]);
+                gsl_matrix_set(A, i, 2, b.x[i]);
+        }
+        m14 = det_get(A, 1);
+        if (fabs(m11) < MYEPSILON)
+                cerr << "ERROR: three points are colinear." << endl;
+        center->x[0] =  0.5 * m12/ m11;
+        center->x[1] = -0.5 * m13/ m11;
+        center->x[2] =  0.5 * m14/ m11;
+        if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
+                cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+        gsl_matrix_free(A);
+};
 /**
 …
  * @param Umkreisradius double radius of circumscribing circle
  */
+  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector* Direction, Vector* AlternativeDirection,
+      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+  {
+    Vector TempNormal, helper;
+    double Restradius;
+    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+    Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+    // Here we calculated center of circumscribing circle, using barycentric coordinates
+    TempNormal.CopyVector(&a);
+    TempNormal.SubtractVector(&b);
+    helper.CopyVector(&a);
+    helper.SubtractVector(&c);
+    TempNormal.VectorProduct(&helper);
+    if (fabs(HalfplaneIndicator) < MYEPSILON)
+      {
+        if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+          {
+            TempNormal.Scale(-1);
+          }
+      }
+    else
+      {
+        if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+          {
+            TempNormal.Scale(-1);
+          }
+      }
+    TempNormal.Normalize();
+    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+    TempNormal.Scale(Restradius);
+    Center->AddVector(&TempNormal);
+  }
+  ;
+void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
+                double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
+        Vector TempNormal, helper;
+        double Restradius;
+        Vector OtherCenter;
+        cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+        Center->Zero();
+        helper.CopyVector(&a);
+        helper.Scale(sin(2.*alpha));
+        Center->AddVector(&helper);
+        helper.CopyVector(&b);
+        helper.Scale(sin(2.*beta));
+        Center->AddVector(&helper);
+        helper.CopyVector(&c);
+        helper.Scale(sin(2.*gamma));
+        Center->AddVector(&helper);
+        //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+        Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+        NewUmkreismittelpunkt->CopyVector(Center);
+        cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
+        // Here we calculated center of circumscribing circle, using barycentric coordinates
+        cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
+        TempNormal.CopyVector(&a);
+        TempNormal.SubtractVector(&b);
+        helper.CopyVector(&a);
+        helper.SubtractVector(&c);
+        TempNormal.VectorProduct(&helper);
+        if (fabs(HalfplaneIndicator) < MYEPSILON)
+                {
+                        if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+                                {
+                                        TempNormal.Scale(-1);
+                                }
+                }
+        else
+                {
+                        if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+                                {
+                                        TempNormal.Scale(-1);
+                                }
+                }
+        TempNormal.Normalize();
+        Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+        cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
+        TempNormal.Scale(Restradius);
+        cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
+        Center->AddVector(&TempNormal);
+        cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
+        get_sphere(&OtherCenter, a, b, c, RADIUS);
+        cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
+        cout << Verbose(3) << "End of Get_center_of_sphere.\n";
+};
 /** This recursive function finds a third point, to form a triangle with two given ones.
  * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
  *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
  *  upon which we operate.
  *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
  *  direction and angle into Storage.
  *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
  *  with all neighbours of the candidate.
+ *      supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
+ *      upon which we operate.
+ *      If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
+ *      direction and angle into Storage.
+ *      We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
+ *      with all neighbours of the candidate.
  * @param a first point
  * @param b second point
 …
  * @param mol molecule structure with atoms and bonds
  */
+void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
+    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
+{
+  //cout << "ReferencePoint is " << ReferencePoint.x[0] << " "<< ReferencePoint.x[1] << " "<< ReferencePoint.x[2] << " "<< endl;
+  /* OldNormal is normal vector on the old triangle
+   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck;
+  Vector TempNormal, Umkreismittelpunkt;
+  Vector Mittelpunkt;
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+  double BallAngle, AlternativeSign;
+  atom *Walker; // variable atom point
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  AngleCheck.CopyVector(&(Candidate->x));
+  AngleCheck.SubtractVector(&(a->x));
+  AngleCheck.ProjectOntoPlane(Chord);
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+  if (a != Candidate and b != Candidate and c != Candidate)
+    {
+      Umkreisradius = SideA / 2.0 / sin(alpha);
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
+      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
+        {
+          cout << Verbose(1) << "Candidate is "<< *Candidate << endl;
+          sign = AngleCheck.ScalarProduct(direction1);
+          if (fabs(sign)<MYEPSILON)
+            {
+              if (AngleCheck.ScalarProduct(OldNormal)<0)
+                {
+                  sign =0;
+                  AlternativeSign=1;
+                }
+              else
+                {
+                  sign =0;
+                  AlternativeSign=-1;
+                }
+            }
+          else
+            {
+              sign /= fabs(sign);
+            }
+          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+          AngleCheck.CopyVector(&ReferencePoint);
+          AngleCheck.Scale(-1);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          AngleCheck.AddVector(&Mittelpunkt);
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          BallAngle = AngleCheck.Angle(OldNormal);
+          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+          cout << Verbose(1) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+          if (AngleCheck.ScalarProduct(direction1) >=0)
+            {
+              if (Storage[0]< -1.5) // first Candidate at all
+                {
+                  cout << "Next better candidate is " << *Candidate << " with ";
+                  Opt_Candidate = Candidate;
+                  Storage[0] = sign;
+                  Storage[1] = AlternativeSign;
+                  Storage[2] = BallAngle;
+                  cout << "Angle is " << Storage[2] << ", Halbraum ist "
+                    << Storage[0] << endl;
+                }
+              else
+                {
+                  if ( Storage[2] > BallAngle)
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = AlternativeSign;
+                      Storage[2] = BallAngle;
+                      cout << "Angle is " << Storage[2] << ", Halbraum ist "
+                        << Storage[0] << endl;
+                    }
+                  else
+                    {
+                      //if (DEBUG)
+                        cout << "Looses to better candidate" << endl;
+                    }
+                }
+              }
+            else
+              {
+                //if (DEBUG)
+                  cout << "Refused due to bad direction of ball centre." << endl;
+              }
+          }
+        else
+          {
+            //if (DEBUG)
+              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+          }
+      }
+    else
+      {
+        //if (DEBUG)
+          cout << "identisch mit Ursprungslinie" << endl;
+      }
+  if (RecursionLevel < 9) // Seven is the recursion level threshold.
+    {
+      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
+        {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
+              Candidate);
+          if (Walker == Parent)
+            { // don't go back the same bond
+              continue;
+            }
+          else
+            {
+              Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel
+                  + 1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS,
+                  mol); //call function again
+            }
+        }
+    }
+}
+;
+  /** This recursive function finds a third point, to form a triangle with two given ones.
+   * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
+   *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
+   *  upon which we operate.
+   *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
+   *  direction and angle into Storage.
+   *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
+   *  with all neighbours of the candidate.
+   * @param a first point
+   * @param b second point
+   * @param Candidate base point along whose bonds to start looking from
+   * @param Parent point to avoid during search as its wrong direction
+   * @param RecursionLevel contains current recursion depth
+   * @param Chord baseline vector of first and second point
+   * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
+   * @param OldNormal normal of the triangle which the baseline belongs to
+   * @param Opt_Candidate candidate reference to return
+   * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
+   * @param Storage array containing two angles of current Opt_Candidate
+   * @param RADIUS radius of ball
+   * @param mol molecule structure with atoms and bonds
+   */
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
+    atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
+{
+  /* OldNormal is normal vector on the old triangle
+   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
+   * Chord points from b to a!!!
+   */
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
+  Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
+  double CurrentEpsilon = 0.1;
+  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+  double BallAngle;
+  atom *Walker; // variable atom point
+  dif_a.CopyVector(&(a->x));
+  dif_a.SubtractVector(&(Candidate->x));
+  dif_b.CopyVector(&(b->x));
+  dif_b.SubtractVector(&(Candidate->x));
+  DirectionCheckPoint.CopyVector(&dif_a);
+  DirectionCheckPoint.Scale(-1);
+  DirectionCheckPoint.ProjectOntoPlane(Chord);
+  SideA = dif_b.Norm();
+  SideB = dif_a.Norm();
+  SideC = Chord->Norm();
+  //Chord->Scale(-1);
+  alpha = Chord->Angle(&dif_a);
+  beta = M_PI - Chord->Angle(&dif_b);
+  gamma = dif_a.Angle(&dif_b);
+  if (DEBUG)
+    {
+      cout << "Atom number" << Candidate->nr << endl;
+      Candidate->x.Output((ofstream *) &cout);
+      cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr]
+          << endl;
+    }
+  if (a != Candidate and b != Candidate)
+    {
+      //      alpha = dif_a.Angle(&dif_b) / 2.;
+      //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+      //      SideB = dif_a.Norm();
+      //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+      //          alpha); // note this is squared of center line length
+      //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+      // Those are remains from Freddie. Needed?
+      Umkreisradius = SideA / 2.0 / sin(alpha);
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
+      if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) //Checking whether ball will at least rest o points.
+        {
+          // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
+          Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+          Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+          TempNormal.CopyVector(&dif_a);
+          TempNormal.VectorProduct(&dif_b);
+          if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0)
+            {
+              TempNormal.Scale(-1);
+            }
+          TempNormal.Normalize();
+          Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+          TempNormal.Scale(Restradius);
+          Mittelpunkt.CopyVector(&Umkreismittelpunkt);
+          Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
+          AngleCheck.CopyVector(Chord);
+          AngleCheck.Scale(-0.5);
+          AngleCheck.SubtractVector(&(b->x));
+          AngleCheckReference.CopyVector(&AngleCheck);
+          AngleCheckReference.AddVector(Opt_Mittelpunkt);
+          AngleCheck.AddVector(&Mittelpunkt);
+          BallAngle = AngleCheck.Angle(&AngleCheckReference);
+          d1->ProjectOntoPlane(&AngleCheckReference);
+          sign = AngleCheck.ScalarProduct(d1);
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+          if (Storage[0]< -1.5) // first Candidate at all
+            {
+              cout << "Next better candidate is " << *Candidate << " with ";
+              Opt_Candidate = Candidate;
+              Storage[0] = sign;
+              Storage[1] = BallAngle;
+              Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+              cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                  << Storage[0] << endl;
+            }
+          else
+            {
+              /*
+               * removed due to change in criterium, now checking angle of ball to old normal.
+              //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
+              //within the ball.
+              Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
+              //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+              if (Distance >RADIUS) // We have no interference and may now check whether the new point is better.
+                */
+                {
+                  //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+                  if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon)))  //This will give absolute preference to those in "right-hand" quadrants
+                      //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = BallAngle;
+                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                          << Storage[0] << endl;
+                    }
+                  else
+                    {
+                      if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0
+                          && Storage[1] > BallAngle) ||
+                          (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0
+                          && Storage[1] < BallAngle))
+                      //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+                        {
+                          cout << "Next better candidate is " << *Candidate << " with ";
+                          Opt_Candidate = Candidate;
+                          Storage[0] = sign;
+                          Storage[1] = BallAngle;
+                          Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                          cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                              << Storage[0] << endl;
+                        }
+                    }
+                }
+              /*
+               * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
+               *
+                else
+                {
+                  if (sign>0 && BallAngle>0 && Storage[0]<0)
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = BallAngle;
+                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                      << Storage[0] << endl;
+//Debugging purposes only
+                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+                      cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                      cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                      cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                      cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                      cout << "Umkreisradius ist " << Umkreisradius << endl;
+                      cout << "Restradius ist " << Restradius << endl;
+                      cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                      cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                      cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                      cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                      cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                      cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+                    }
+                  else
+                    {
+                      if (DEBUG)
+                        cout << "Looses to better candidate" << endl;
+                    }
+                }
+                */
+            }
+        }
+      else
+        {
+          if (DEBUG)
+            {
+              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+            }
+        }
+    }
+  else
+    {
+      if (DEBUG)
+        cout << "identisch mit Ursprungslinie" << endl;
+    }
+  if (RecursionLevel < 9) // Five is the recursion level threshold.
+    {
+      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
+        {
+          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
+              Candidate);
+          if (Walker == Parent)
+            { // don't go back the same bond
+              continue;
+            }
+          else
+            {
+              Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel
+                  + 1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS,
+                  mol); //call function again
+            }
+        }
+    }
+}
+;
+void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
+                int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
+                atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
+{
+        cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+        cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
+        cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
+        cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
+        cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
+        cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
+        /* OldNormal is normal vector on the old triangle
+         * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+         * Chord points from b to a!!!
+         */
+        Vector dif_a; //Vector from a to candidate
+        Vector dif_b; //Vector from b to candidate
+        Vector AngleCheck;
+        Vector TempNormal, Umkreismittelpunkt;
+        Vector Mittelpunkt;
+        double CurrentEpsilon = 0.1;
+        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+        double BallAngle, AlternativeSign;
+        atom *Walker; // variable atom point
+        Vector NewUmkreismittelpunkt;
+        if (a != Candidate and b != Candidate and c != Candidate) {
+                cout << Verbose(3) << "We have a unique candidate!" << endl;
+                dif_a.CopyVector(&(a->x));
+                dif_a.SubtractVector(&(Candidate->x));
+                dif_b.CopyVector(&(b->x));
+                dif_b.SubtractVector(&(Candidate->x));
+                AngleCheck.CopyVector(&(Candidate->x));
+                AngleCheck.SubtractVector(&(a->x));
+                AngleCheck.ProjectOntoPlane(Chord);
+                SideA = dif_b.Norm();
+                SideB = dif_a.Norm();
+                SideC = Chord->Norm();
+                //Chord->Scale(-1);
+                alpha = Chord->Angle(&dif_a);
+                beta = M_PI - Chord->Angle(&dif_b);
+                gamma = dif_a.Angle(&dif_b);
+                cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+                if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
+                        cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+                        cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+                }
+                if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
+                        Umkreisradius = SideA / 2.0 / sin(alpha);
+                        //cout << Umkreisradius << endl;
+                        //cout << SideB / 2.0 / sin(beta) << endl;
+                        //cout << SideC / 2.0 / sin(gamma) << endl;
+                        if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
+                                cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+                                cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
+                                sign = AngleCheck.ScalarProduct(direction1);
+                                if (fabs(sign)<MYEPSILON) {
+                                        if (AngleCheck.ScalarProduct(OldNormal)<0) {
+                                                sign =0;
+                                                AlternativeSign=1;
+                                        } else {
+                                                sign =0;
+                                                AlternativeSign=-1;
+                                        }
+                                } else {
+                                        sign /= fabs(sign);
+                                }
+                                if (sign >= 0) {
+                                        cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+                                        Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+                                        Mittelpunkt.SubtractVector(&ReferencePoint);
+                                        cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
+                                        BallAngle = Mittelpunkt.Angle(OldNormal);
+                                        cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+                                        //cout << "direction1 is " << *direction1 << "." << endl;
+                                        //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
+                                        //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+                                        NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+                                        if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
+                                                if (Storage[0]< -1.5) { // first Candidate at all
+                                                        if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+                                                                cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
+                                                                Opt_Candidate = Candidate;
+                                                                Storage[0] = sign;
+                                                                Storage[1] = AlternativeSign;
+                                                                Storage[2] = BallAngle;
+                                                                cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
+                                                        } else
+                                                                cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+                                                } else {
+                                                        if ( Storage[2] > BallAngle) {
+                                                                if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+                                                                        cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
+                                                                        Opt_Candidate = Candidate;
+                                                                        Storage[0] = sign;
+                                                                        Storage[1] = AlternativeSign;
+                                                                        Storage[2] = BallAngle;
+                                                                        cout << " angle " << Storage[2] << " and Up/Down " << Storage[0] << endl;
+                                                                } else
+                                                                        cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+                                                        } else {
+                                                                if (DEBUG) {
+                                                                        cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
+                                                                }
+                                                        }
+                                                }
+                                        } else {
+                                                if (DEBUG) {
+                                                        cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
+                                                }
+                                        }
+                                } else {
+                                        if (DEBUG) {
+                                                cout << Verbose(3) << *Candidate << " is not in search direction." << endl;
+                                        }
+                                }
+                        } else {
+                                if (DEBUG) {
+                                        cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
+                                }
+                        }
+                } else {
+                        if (DEBUG) {
+                                cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
+                        }
+                }
+        } else {
+                if (DEBUG) {
+                        cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
+                }
+        }
+        if (RecursionLevel < 5) { // Seven is the recursion level threshold.
+                for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+                        Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+                        if (Walker == Parent) { // don't go back the same bond
+                                continue;
+                        } else {
+                                Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
+                        }
+                }
+        }
+        cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+}
+;
+/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
+ * \param *Center new center on return
+ * \param *a first point
+ * \param *b second point
+ * \param *c third point
+ */
+void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
+        Vector helper;
+        double alpha, beta, gamma;
+        Vector SideA, SideB, SideC;
+        SideA.CopyVector(b);
+        SideA.SubtractVector(c);
+        SideB.CopyVector(c);
+        SideB.SubtractVector(a);
+        SideC.CopyVector(a);
+        SideC.SubtractVector(b);
+        alpha = M_PI - SideB.Angle(&SideC);
+        beta = M_PI - SideC.Angle(&SideA);
+        gamma = M_PI - SideA.Angle(&SideB);
+        cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
+        if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
+                cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+        Center->Zero();
+        helper.CopyVector(a);
+        helper.Scale(sin(2.*alpha));
+        Center->AddVector(&helper);
+        helper.CopyVector(b);
+        helper.Scale(sin(2.*beta));
+        Center->AddVector(&helper);
+        helper.CopyVector(c);
+        helper.Scale(sin(2.*gamma));
+        Center->AddVector(&helper);
+        Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+};
+/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
+ * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
+ * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
+ * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
+ * \param CircleCenter Center of the parameter circle
+ * \param CirclePlaneNormal normal vector to plane of the parameter circle
+ * \param CircleRadius radius of the parameter circle
+ * \param NewSphereCenter new center of a circumcircle
+ * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
+ * \param NormalVector normal vector
+ * \param SearchDirection search direction to make angle unique on return.
+ * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
+ */
+double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
+        Vector helper;
+        double radius, alpha;
+        helper.CopyVector(&NewSphereCenter);
+        // test whether new center is on the parameter circle's plane
+        if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+                cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))        << "!" << endl;
+                helper.ProjectOntoPlane(&CirclePlaneNormal);
+        }
+        radius = helper.ScalarProduct(&helper);
+        // test whether the new center vector has length of CircleRadius
+        if (fabs(radius - CircleRadius) > HULLEPSILON)
+                cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+        alpha = helper.Angle(&OldSphereCenter);
+        // make the angle unique by checking the halfplanes/search direction
+        if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)      // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+                alpha = 2.*M_PI - alpha;
+        cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+        radius = helper.Distance(&OldSphereCenter);
+        helper.ProjectOntoPlane(&NormalVector);
+        // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
+        if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
+                cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+                return alpha;
+        } else {
+                cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+                return 2.*M_PI;
+        }
+};
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points
+ * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line
+ * @param OptCandidate candidate reference on return
+ * @param OptCandidateCenter candidate's sphere center on return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+// void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+// {
+//       atom *Walker = NULL;
+//       Vector CircleCenter;   // center of the circle, i.e. of the band of sphere's centers
+//       Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+//       Vector OldSphereCenter;         // center of the sphere defined by the three points of BaseTriangle
+//       Vector NewSphereCenter;         // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+//       Vector OtherNewSphereCenter;    // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+//       Vector NewNormalVector;         // normal vector of the Candidate's triangle
+//       Vector SearchDirection;         // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
+//       Vector helper;
+//       LinkedAtoms *List = NULL;
+//       double CircleRadius; // radius of this circle
+//       double radius;
+//       double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+//       double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
+//       int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+//       atom *Candidate = NULL;
+//
+//       cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
+//
+//       cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
+//
+//       // construct center of circle
+//       CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+//       CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+//       CircleCenter.Scale(0.5);
+//
+//       // construct normal vector of circle
+//       CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+//       CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+//
+//       // calculate squared radius of circle
+//       radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+//       if (radius/4. < RADIUS*RADIUS) {
+//               CircleRadius = RADIUS*RADIUS - radius/4.;
+//               CirclePlaneNormal.Normalize();
+//               cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+//
+//               // construct old center
+//               GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
+//               helper.CopyVector(&BaseTriangle->NormalVector);        // normal vector ensures that this is correct center of the two possible ones
+//               radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+//               helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//               OldSphereCenter.AddVector(&helper);
+//               OldSphereCenter.SubtractVector(&CircleCenter);
+//               cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+//
+//               // test whether old center is on the band's plane
+//               if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+//                       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+//                       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+//               }
+//               radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+//               if (fabs(radius - CircleRadius) < HULLEPSILON) {
+//
+//                       // construct SearchDirection
+//                       SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
+//                       helper.CopyVector(&BaseLine->endpoints[0]->node->x);
+//                       for(int i=0;i<3;i++)   // just take next different endpoint
+//                               if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
+//                                       helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
+//                               }
+//                       if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)     // ohoh, SearchDirection points inwards!
+//                               SearchDirection.Scale(-1.);
+//                       SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+//                       SearchDirection.Normalize();
+//                       cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+//                       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {     // rotated the wrong way!
+//                               cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+//                       }
+//
+//                       if (LC->SetIndexToVector(&CircleCenter)) {     // get cell for the starting atom
+//                               for(int i=0;i<NDIM;i++) // store indices of this cell
+//                                       N[i] = LC->n[i];
+//                               cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+//                       } else {
+//                               cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+//                               return;
+//                       }
+//                       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+//                       cout << Verbose(2) << "LC Intervals:";
+//                       for (int i=0;i<NDIM;i++) {
+//                               Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+//                               Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+//                               cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+//                       }
+//                       cout << endl;
+//                       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+//                               for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+//                                       for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+//                                               List = LC->GetCurrentCell();
+//                                               cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+//                                               if (List != NULL) {
+//                                                       for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+//                                                               Candidate = (*Runner);
+//
+//                                                               // check for three unique points
+//                                                               if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
+//                                                                       cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+//
+//                                                                       // construct both new centers
+//                                                                       GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+//                                                                       OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+//
+//                                                                       if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+//                                                                               helper.CopyVector(&NewNormalVector);
+//                                                                               cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+//                                                                               radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+//                                                                               if (radius < RADIUS*RADIUS) {
+//                                                                                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//                                                                                       cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+//                                                                                       NewSphereCenter.AddVector(&helper);
+//                                                                                       NewSphereCenter.SubtractVector(&CircleCenter);
+//                                                                                       cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+//
+//                                                                                       helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+//                                                                                       OtherNewSphereCenter.AddVector(&helper);
+//                                                                                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
+//                                                                                       cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+//
+//                                                                                       // check both possible centers
+//                                                                                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                                                                                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                                                                                       alpha = min(alpha, Otheralpha);
+//                                                                                       if (*ShortestAngle > alpha) {
+//                                                                                                       OptCandidate = Candidate;
+//                                                                                                       *ShortestAngle = alpha;
+//                                                                                                       if (alpha != Otheralpha)
+//                                                                                                               OptCandidateCenter->CopyVector(&NewSphereCenter);
+//                                                                                                       else
+//                                                                                                               OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+//                                                                                                       cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+//                                                                                       } else {
+//                                                                                               if (OptCandidate != NULL)
+//                                                                                                       cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+//                                                                                               else
+//                                                                                                       cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+//                                                                                       }
+//
+//                                                                               } else {
+//                                                                                       cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+//                                                                               }
+//                                                                       } else {
+//                                                                               cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+//                                                                       }
+//                                                               } else {
+//                                                                       cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
+//                                                               }
+//                                                       }
+//                                               }
+//                                       }
+//               } else {
+//                       cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+//               }
+//       } else {
+//               cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
+//       }
+//
+//       cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
+// };
+/** Checks whether the triangle consisting of the three atoms is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return false - triangle invalid due to edge criteria, true - triangle may be added.
+ */
+bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
+        LineMap::iterator FindLine;
+        PointMap::iterator FindPoint;
+        bool Present[3];
+        *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+        for (int i=0;i<3;i++) { // check through all endpoints
+                FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+                if (FindPoint != PointsOnBoundary.end())
+                        TPS[i] = FindPoint->second;
+                else
+                        TPS[i] = NULL;
+        }
+        // check lines
+        for (int i=0;i<3;i++)
+                if (TPS[i] != NULL)
+                        for (int j=i;j<3;j++)
+                                if (TPS[j] != NULL) {
+                                        FindLine = TPS[i]->lines.find(TPS[j]->node->nr);
+                                        if ((FindLine != TPS[i]->lines.end()) && (FindLine->second->TrianglesCount > 1)) {
+                                                *out << "WARNING: Line " << *FindLine->second << " already present with " << FindLine->second->TrianglesCount << " triangles attached." << endl;
+                                                *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+                                                return false;
+                                        }
+                                }
+        *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+        return true;
+};
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * Note that this function is for the starting triangle.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
+ * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
+ * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
+ * @param BaseLine BoundaryLineSet with the current base line
+ * @param ThirdNode third atom to avoid in search
+ * @param OptCandidate candidate reference on return
+ * @param OptCandidateCenter candidate's sphere center on return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+{
+        atom *Walker = NULL;
+        Vector CircleCenter;    // center of the circle, i.e. of the band of sphere's centers
+        Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+        Vector NewSphereCenter;  // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+        Vector OtherNewSphereCenter;     // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+        Vector NewNormalVector;  // normal vector of the Candidate's triangle
+        Vector helper;
+        LinkedAtoms *List = NULL;
+        double CircleRadius; // radius of this circle
+        double radius;
+        double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+        double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
+        int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+        atom *Candidate = NULL;
+        cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
+        cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+        // construct center of circle
+        CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+        CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+        CircleCenter.Scale(0.5);
+        // construct normal vector of circle
+        CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+        CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+        // calculate squared radius of circle
+        radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+        if (radius/4. < RADIUS*RADIUS) {
+                CircleRadius = RADIUS*RADIUS - radius/4.;
+                CirclePlaneNormal.Normalize();
+                cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                // test whether old center is on the band's plane
+                if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+                        cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+                        OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+                }
+                radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+                if (fabs(radius - CircleRadius) < HULLEPSILON) {
+                        // check SearchDirection
+                        cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                        if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {      // rotated the wrong way!
+                                cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+                        }
+                        // get cell for the starting atom
+                        if (LC->SetIndexToVector(&CircleCenter)) {
+                                        for(int i=0;i<NDIM;i++) // store indices of this cell
+                                        N[i] = LC->n[i];
+                                cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+                        } else {
+                                cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+                                return;
+                        }
+                        // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+                        cout << Verbose(2) << "LC Intervals:";
+                        for (int i=0;i<NDIM;i++) {
+                                Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+                                Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+                                cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+                        }
+                        cout << endl;
+                        for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+                                for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+                                        for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+                                                List = LC->GetCurrentCell();
+                                                //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+                                                if (List != NULL) {
+                                                        for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                                                                Candidate = (*Runner);
+                                                                // check for three unique points
+                                                                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
+                                                                        cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+                                                                        // construct both new centers
+                                                                        GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+                                                                        OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+                                                                        if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+                                                                                helper.CopyVector(&NewNormalVector);
+                                                                                cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                                                                                radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+                                                                                if (radius < RADIUS*RADIUS) {
+                                                                                        helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                                                                                        cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+                                                                                        NewSphereCenter.AddVector(&helper);
+                                                                                        NewSphereCenter.SubtractVector(&CircleCenter);
+                                                                                        cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+                                                                                        helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+                                                                                        OtherNewSphereCenter.AddVector(&helper);
+                                                                                        OtherNewSphereCenter.SubtractVector(&CircleCenter);
+                                                                                        cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+                                                                                        // check both possible centers
+                                                                                        alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                                                                                        Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                                                                                        alpha = min(alpha, Otheralpha);
+                                                                                        if (*ShortestAngle > alpha) {
+                                                                                                        OptCandidate = Candidate;
+                                                                                                        *ShortestAngle = alpha;
+                                                                                                        if (alpha != Otheralpha)
+                                                                                                                OptCandidateCenter->CopyVector(&NewSphereCenter);
+                                                                                                        else
+                                                                                                                OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+                                                                                                        cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+                                                                                        } else {
+                                                                                                if (OptCandidate != NULL)
+                                                                                                        cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+                                                                                                else
+                                                                                                        cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                                                                                        }
+                                                                                } else {
+                                                                                        cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                                                                                }
+                                                                        } else {
+                                                                                cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                                                                        }
+                                                                } else {
+                                                                        if (ThirdNode != NULL)
+                                                                                cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                                                                        else
+                                                                                cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                                                                }
+                                                        }
+                                                }
+                                        }
+                } else {
+                        cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+                }
+        } else {
+                if (ThirdNode != NULL)
+                        cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+                else
+                        cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+        }
+        cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
+};
+/** Finds the second point of starting triangle.
+ * \param *a first atom
+ * \param *Candidate pointer to candidate atom on return
+ * \param Oben vector indicating the outside
+ * \param Opt_Candidate reference to recommended candidate on return
+ * \param Storage[3] array storing angles and other candidate information
+ * \param RADIUS radius of virtual sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
+        cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
+        int i;
+        Vector AngleCheck;
+        atom* Walker;
+        double norm = -1., angle;
+        LinkedAtoms *List = NULL;
+        int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+        if (LC->SetIndexToAtom(a)) {    // get cell for the starting atom
+                for(int i=0;i<NDIM;i++) // store indices of this cell
+                        N[i] = LC->n[i];
+        } else {
+                cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+                return;
+        }
+        // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+        cout << Verbose(2) << "LC Intervals:";
+        for (int i=0;i<NDIM;i++) {
+                Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+                Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+                cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+        }
+        cout << endl;
+        for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+                for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+                        for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+                                List = LC->GetCurrentCell();
+                                //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+                                if (List != NULL) {
+                                        for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                                                Candidate = (*Runner);
+                                                                // check if we only have one unique point yet ...
+                                                                if (a != Candidate) {
+                        cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
+                        AngleCheck.CopyVector(&(Candidate->x));
+                        AngleCheck.SubtractVector(&(a->x));
+                        norm = AngleCheck.Norm();
+                                                        // second point shall have smallest angle with respect to Oben vector
+                                                        if (norm < RADIUS) {
+                                                                angle = AngleCheck.Angle(&Oben);
+                                                                if (angle < Storage[0]) {
+                                                                        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                                                                        cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
+                                                                        Opt_Candidate = Candidate;
+                                                                        Storage[0] = AngleCheck.Angle(&Oben);
+                                                                        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+                                                                } else {
+                                                                        cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+                                                                }
+                                                        } else {
+                                                                cout << "Refused due to Radius " << norm << endl;
+                                                        }
+                                                }
+                                        }
+                                }
+                        }
+        cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
+};
+/** Finds the starting triangle for find_non_convex_border().
+ * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
+ * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
+ * point are called.
+ * \param RADIUS radius of virtual rolling sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
+{
+        cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+        int i = 0;
+        LinkedAtoms *List = NULL;
+        atom* Walker;
+        atom* FirstPoint;
+        atom* SecondPoint;
+        atom* MaxAtom[NDIM];
+        double max_coordinate[NDIM];
+        Vector Oben;
+        Vector helper;
+        Vector Chord;
+        Vector SearchDirection;
+        Vector OptCandidateCenter;
+        Oben.Zero();
+        for (i = 0; i < 3; i++) {
+                MaxAtom[i] = NULL;
+                max_coordinate[i] = -1;
+        }
+        // 1. searching topmost atom with respect to each axis
+        for (int i=0;i<NDIM;i++) { // each axis
+                LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+                for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+                        for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+                                List = LC->GetCurrentCell();
+                                //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+                                if (List != NULL) {
+                                        for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+                                                cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
+                                                if ((*Runner)->x.x[i] > max_coordinate[i]) {
+                                                        max_coordinate[i] = (*Runner)->x.x[i];
+                                                        MaxAtom[i] = (*Runner);
+                                                }
+                                        }
+                                } else {
+                                        cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+                                }
+                        }
+        }
+        cout << Verbose(2) << "Found maximum coordinates: ";
+        for (int i=0;i<NDIM;i++)
+                cout << i << ": " << *MaxAtom[i] << "\t";
+        cout << endl;
+        const int k = 1;                // arbitrary choice
+        Oben.x[k] = 1.;
+        FirstPoint = MaxAtom[k];
+        cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
+        // add first point
+        AddTrianglePoint(FirstPoint, 0);
+        double ShortestAngle;
+        atom* Opt_Candidate = NULL;
+        ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+        Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+        SecondPoint = Opt_Candidate;
+        cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+        // add second point and first baseline
+        AddTrianglePoint(SecondPoint, 1);
+        AddTriangleLine(TPS[0], TPS[1], 0);
+        helper.CopyVector(&(FirstPoint->x));
+        helper.SubtractVector(&(SecondPoint->x));
+        helper.Normalize();
+        Oben.ProjectOntoPlane(&helper);
+        Oben.Normalize();
+        helper.VectorProduct(&Oben);
+        ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+        Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+        Chord.SubtractVector(&(SecondPoint->x));
+        double radius = Chord.ScalarProduct(&Chord);
+        double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+        helper.CopyVector(&Oben);
+        helper.Scale(CircleRadius);
+        // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+        cout << Verbose(2) << "Looking for third point candidates \n";
+        // look in one direction of baseline for initial candidate
+        Opt_Candidate = NULL;
+        SearchDirection.MakeNormalVector(&Chord, &Oben);        // whether we look "left" first or "right" first is not important ...
+        cout << Verbose(1) << "Looking for third point candidates ...\n";
+        Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
+        cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+        // add third point
+        AddTrianglePoint(Opt_Candidate, 2);
+        // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+        // Finally, we only have to add the found further lines
+        AddTriangleLine(TPS[1], TPS[2], 1);
+        AddTriangleLine(TPS[0], TPS[2], 2);
+        // ... and triangles to the Maps of the Tesselation class
+        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+        AddTriangleToLines();
+        // ... and calculate its normal vector (with correct orientation)
+        OptCandidateCenter.Scale(-1.);
+        cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
+        BTS->GetNormalVector(OptCandidateCenter);
+        cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+        cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+        cout << Verbose(1) << "End of Find_starting_triangle\n";
+};
 /** This function finds a triangle to a line, adjacent to an existing one.
  * @param out   output stream for debugging
  * @param mol molecule structure with all atoms and bonds
+ * @param out output stream for debugging
+ * @param *mol molecule with Atom's and Bond's
  * @param Line current baseline to search from
  * @param T current triangle which \a Line is edge of
  * @param RADIUS radius of the rolling ball
  * @param N number of found triangles
+ * @param *filename filename base for intermediate envelopes
+ * @param *LC LinkedCell structure with neighbouring atoms
  */
+void Tesselation::Find_next_suitable_triangle(ofstream *out,
+    molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+    const double& RADIUS, int N, const char *tempbasename)
+{
+  cout << Verbose(1) << "Looking for next suitable triangle \n";
+  Vector direction1;
+  Vector helper;
+  Vector Chord;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  double tmp;
+  //atom* Walker;
+  atom* OldThirdPoint;
+  double Storage[3];
+  Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
+  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+  Storage[2] = 9999999.;
+  atom* Opt_Candidate = NULL;
+  Vector Opt_Mittelpunkt;
+  cout << Verbose(1) << "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)
+        {
+          OldThirdPoint = T.endpoints[i]->node;
+          helper.SubtractVector(&T.endpoints[i]->node->x);
+          break;
+        }
+    }
+  direction1.CopyVector(&Line.endpoints[0]->node->x);
+  direction1.SubtractVector(&Line.endpoints[1]->node->x);
+  direction1.VectorProduct(&(T.NormalVector));
+  if (direction1.ScalarProduct(&helper) < 0)
+    {
+      direction1.Scale(-1);
+    }
+  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+    Vector Umkreismittelpunkt, a, b, c;
+    double alpha, beta, gamma;
+    a.CopyVector(&(T.endpoints[0]->node->x));
+    b.CopyVector(&(T.endpoints[1]->node->x));
+    c.CopyVector(&(T.endpoints[2]->node->x));
+    a.SubtractVector(&(T.endpoints[1]->node->x));
+    b.SubtractVector(&(T.endpoints[2]->node->x));
+    c.SubtractVector(&(T.endpoints[0]->node->x));
+    alpha = M_PI - a.Angle(&c);
+    beta = M_PI - b.Angle(&a);
+    gamma = M_PI - c.Angle(&b);
+    Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+    //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+    cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
+      Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
+      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
+      Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
+      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+      cout << "Letzter Winkel bei " << TrianglesOnBoundaryCount << " Winkel ist " << Storage[2] << endl;
+  if ((TrianglesOnBoundaryCount % 10) == 0) {
+    sprintf(NumberName, "-%d", TriangleFilesWritten);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(TecplotSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoRaster3DOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      NameofTempFile.append(Raster3DSuffix);
+      cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+      tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+      write_raster3d_file(out, tempstream, this, mol);
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
+  }
+      // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
+  cout << " Optimal candidate is " << *Opt_Candidate << endl;
+  AddTrianglePoint(Opt_Candidate, 0);
+  AddTrianglePoint(Line.endpoints[0]->node, 1);
+  AddTrianglePoint(Line.endpoints[1]->node, 2);
+  AddTriangleLine(TPS[0], TPS[1], 0);
+  AddTriangleLine(TPS[0], TPS[2], 1);
+  AddTriangleLine(TPS[1], TPS[2], 2);
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  AddTriangleToLines();
+  cout << "New triangle with " << *BTS << endl;
+  cout << "We have "<< TrianglesOnBoundaryCount << 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) ||
+      (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
+    {
+      BTS->NormalVector.Scale(-1);
+    };
+}
+;
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
+    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[3],
+    molecule* mol, double RADIUS)
+{
+  cout << Verbose(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[2]);
+              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[2]);
+            }
+          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;
+  atom* Walker;
+  atom* FirstPoint;
+  atom* SecondPoint;
+  atom* max_index[3];
+  double max_coordinate[3];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector CenterOfFirstLine;
+  Oben.Zero();
+  for (i = 0; i < 3; i++)
+    {
+      max_index[i] = NULL;
+      max_coordinate[i] = -1;
+    }
+  cout << Verbose(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;
+            }
+        }
+    }
+  cout << Verbose(1) << "Found maximum coordinates. " << endl;
+  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
+  const int k = 1;
+  Oben.x[k] = 1.;
+  FirstPoint = max_index[k];
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << ": "
+      << FirstPoint->x.x[0] << endl;
+  double Storage[3];
+  atom* Opt_Candidate = NULL;
+  Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+  Storage[1] = 999999.; // This will be an angle looking for the third point.
+  Storage[2] = 999999.;
+  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;
+  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[2] = 9999999.;
+  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+  Chord.SubtractVector(&(SecondPoint->x));
+  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+  cout << Verbose(1) << "Looking for third point candidates \n";
+  cout << Verbose(1) << " In direction " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << " " << endl;
+  // look in one direction of baseline for initial candidate
+  Opt_Candidate = NULL;
+  CenterOfFirstLine.CopyVector(&Chord);
+  CenterOfFirstLine.Scale(0.5);
+  CenterOfFirstLine.AddVector(&(SecondPoint->x));
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
+      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
+  // look in other direction of baseline for possible better candidate
+  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
+      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
+  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+  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);
+  // ... 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, const char *filename, 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;
+  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
+  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
+  if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
+    mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
+  Tess->Find_starting_triangle(mol, RADIUS);
+  baseline = Tess->LinesOnBoundary.begin();
+  while (baseline != Tess->LinesOnBoundary.end())
+    {
+      if (baseline->second->TrianglesCount == 1)
+        {
+          cout << Verbose(1) << "Begin of Tesselation ... " << endl;
+          Tess->Find_next_suitable_triangle(out, mol,
+              *(baseline->second),
+              *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
+          flag = true;
+          cout << Verbose(1) << "End of Tesselation ... " << endl;
+        }
+      else
+        {
+          cout << Verbose(1) << "There is a line with "
+              << baseline->second->TrianglesCount << " triangles adjacent"
+              << endl;
+        }
+      N++;
+      baseline++;
+    }
+  cout << Verbose(1) << "Writing final tecplot file\n";
+  if (DoTecplotOutput) {
+    string Name(filename);
+    Name.append(TecplotSuffix);
+    ofstream tecplot(Name.c_str(), ios::trunc);
+    write_tecplot_file(out, &tecplot, Tess, mol, -1);
+    tecplot.close();
+  }
+  if (DoRaster3DOutput) {
+    string Name(filename);
+    Name.append(Raster3DSuffix);
+    ofstream raster(Name.c_str(), ios::trunc);
+    write_raster3d_file(out, &raster, Tess, mol);
+    raster.close();
+  }
+}
+;
+bool Tesselation::Find_next_suitable_triangle(ofstream *out,
+                molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+                const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
+        cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+        ofstream *tempstream = NULL;
+        char NumberName[255];
+        double tmp;
+        atom* Opt_Candidate = NULL;
+        Vector OptCandidateCenter;
+        Vector CircleCenter;
+        Vector CirclePlaneNormal;
+        Vector OldSphereCenter;
+        Vector SearchDirection;
+        Vector helper;
+        atom *ThirdNode = NULL;
+        double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+        double radius, CircleRadius;
+        cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+        for (int i=0;i<3;i++)
+                if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+                        ThirdNode = T.endpoints[i]->node;
+        // construct center of circle
+        CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
+        CircleCenter.AddVector(&Line.endpoints[1]->node->x);
+        CircleCenter.Scale(0.5);
+        // construct normal vector of circle
+        CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
+        CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
+        // calculate squared radius of circle
+        radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+        if (radius/4. < RADIUS*RADIUS) {
+                CircleRadius = RADIUS*RADIUS - radius/4.;
+                CirclePlaneNormal.Normalize();
+                cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                // construct old center
+                GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
+                helper.CopyVector(&T.NormalVector);     // normal vector ensures that this is correct center of the two possible ones
+                radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+                helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                OldSphereCenter.AddVector(&helper);
+                OldSphereCenter.SubtractVector(&CircleCenter);
+                cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+                // construct SearchDirection
+                SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
+                helper.CopyVector(&Line.endpoints[0]->node->x);
+                helper.SubtractVector(&ThirdNode->x);
+                if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)      // ohoh, SearchDirection points inwards!
+                        SearchDirection.Scale(-1.);
+                SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+                SearchDirection.Normalize();
+                cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {      // rotated the wrong way!
+                        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+                }
+                // add third point
+                cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+                Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
+        } else {
+                cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+        }
+        if (Opt_Candidate == NULL) {
+                cerr << "WARNING: Could not find a suitable candidate." << endl;
+                return false;
+        }
+        cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << endl;
+        // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+        atom *AtomCandidates[3];
+        AtomCandidates[0] = Opt_Candidate;
+        AtomCandidates[1] = Line.endpoints[0]->node;
+        AtomCandidates[2] = Line.endpoints[1]->node;
+        bool flag = CheckPresenceOfTriangle(out, AtomCandidates);
+        if (flag) { // if so, add
+                AddTrianglePoint(Opt_Candidate, 0);
+                AddTrianglePoint(Line.endpoints[0]->node, 1);
+                AddTrianglePoint(Line.endpoints[1]->node, 2);
+                AddTriangleLine(TPS[0], TPS[1], 0);
+                AddTriangleLine(TPS[0], TPS[2], 1);
+                AddTriangleLine(TPS[1], TPS[2], 2);
+                BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+                AddTriangleToLines();
+                OptCandidateCenter.Scale(-1.);
+                BTS->GetNormalVector(OptCandidateCenter);
+                OptCandidateCenter.Scale(-1.);
+                cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
+                cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+        } else { // else, yell and do nothing
+                cout << Verbose(1) << "This triangle consisting of ";
+                cout << *Opt_Candidate << ", ";
+                cout << *Line.endpoints[0]->node << " and ";
+                cout << *Line.endpoints[1]->node << " ";
+                cout << "is invalid!" << endl;
+                return false;
+        }
+        if (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+                sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
+                if (DoTecplotOutput) {
+                        string NameofTempFile(tempbasename);
+                        NameofTempFile.append(NumberName);
+                        for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos))
+                                NameofTempFile.erase(npos, 1);
+                        NameofTempFile.append(TecplotSuffix);
+                        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+                        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+                        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+                        tempstream->close();
+                        tempstream->flush();
+                        delete(tempstream);
+                }
+                if (DoRaster3DOutput) {
+                        string NameofTempFile(tempbasename);
+                        NameofTempFile.append(NumberName);
+                        for(size_t npos = NameofTempFile.find_first_of(' '); npos != -1; npos = NameofTempFile.find(' ', npos))
+                                NameofTempFile.erase(npos, 1);
+                        NameofTempFile.append(Raster3DSuffix);
+                        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+                        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+                        write_raster3d_file(out, tempstream, this, mol);
+                        // include the current position of the virtual sphere in the temporary raster3d file
+                        // make the circumsphere's center absolute again
+                        helper.CopyVector(&Line.endpoints[0]->node->x);
+                        helper.AddVector(&Line.endpoints[1]->node->x);
+                        helper.Scale(0.5);
+                        OptCandidateCenter.AddVector(&helper);
+                        Vector *center = mol->DetermineCenterOfAll(out);
+                        OptCandidateCenter.AddVector(center);
+                        delete(center);
+                        // and add to file plus translucency object
+                        *tempstream << "# current virtual sphere\n";
+                        *tempstream << "8\n     25.0            0.6              -1.0 -1.0 -1.0          0.2                            0 0 0 0\n";
+                        *tempstream << "2\n     " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
+                        *tempstream << "9\n     terminating special property\n";
+                        tempstream->close();
+                        tempstream->flush();
+                        delete(tempstream);
+                }
+                if (DoTecplotOutput || DoRaster3DOutput)
+                        TriangleFilesWritten++;
+        }
+        cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+        return true;
+};
+/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure with Atom's and Bond's
+ * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
+ * \param *LCList linked cell list of all atoms
+ * \param *filename filename prefix for output of vertex data
+ * \para RADIUS radius of the virtual sphere
+ */
+void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
+{
+        int N = 0;
+        bool freeTess = false;
+        *out << Verbose(1) << "Entering search for non convex hull. " << endl;
+        if (Tess == NULL) {
+                *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
+                Tess = new Tesselation;
+                freeTess = true;
+        }
+        bool freeLC = false;
+        LineMap::iterator baseline;
+        *out << Verbose(0) << "Begin of Find_non_convex_border\n";
+        bool flag = false;      // marks whether we went once through all baselines without finding any without two triangles
+        bool failflag = false;
+        if (LCList == NULL) {
+                LCList = new LinkedCell(mol, 2.*RADIUS);
+                freeLC = true;
+        }
+        Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
+        baseline = Tess->LinesOnBoundary.begin();
+        while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+                if (baseline->second->TrianglesCount == 1) {
+                        failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
+                        flag = flag || failflag;
+                        if (!failflag)
+                                cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+                } else {
+                        cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
+                }
+                N++;
+                baseline++;
+                if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
+                        baseline = Tess->LinesOnBoundary.begin();        // restart if we reach end due to newly inserted lines
+                        flag = false;
+                }
+        }
+        if (1) { //failflag) {
+                *out << Verbose(1) << "Writing final tecplot file\n";
+                if (DoTecplotOutput) {
+                        string OutputName(filename);
+                        OutputName.append(TecplotSuffix);
+                        ofstream *tecplot = new ofstream(OutputName.c_str());
+                        write_tecplot_file(out, tecplot, Tess, mol, -1);
+                        tecplot->close();
+                        delete(tecplot);
+                }
+                if (DoRaster3DOutput) {
+                        string OutputName(filename);
+                        OutputName.append(Raster3DSuffix);
+                        ofstream *raster = new ofstream(OutputName.c_str());
+                        write_raster3d_file(out, raster, Tess, mol);
+                        raster->close();
+                        delete(raster);
+                }
+        } else {
+                cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+        }
+        if (freeTess)
+                delete(Tess);
+        if (freeLC)
+                delete(LCList);
+        *out << Verbose(0) << "End of Find_non_convex_border\n";
+};

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

Context Navigation

Changeset 6ac7ee for src/boundary.cpp

Legend:

src/boundary.cpp

Download in other formats: