Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/Makefile.am	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,4 +1,4 @@
-SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp vector.cpp verbose.cpp
-HEADER = boundary.hpp defs.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp vector.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
+HEADER = atom.hpp bond.hpp boundary.hpp defs.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp
 
 bin_PROGRAMS = molecuilder joiner analyzer
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/atom.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -5,5 +5,5 @@
  */
 
-#include "molecules.hpp"
+#include "atom.hpp"
 
 /************************************* Functions for class atom *************************************/
@@ -13,13 +13,11 @@
 atom::atom()
 {
-  Name = NULL;
+  father = this;  // generally, father is itself
   previous = NULL;
   next = NULL;
-  father = this;  // generally, father is itself
   Ancestor = NULL;
   type = NULL;
   sort = NULL;
   FixedIon = 0;
-  nr = -1;
   GraphNr = -1;
   ComponentNr = NULL;
@@ -29,8 +27,9 @@
   AdaptiveOrder = 0;
   MaxOrder = false;
+  // set LCNode::Vector to our Vector
+  node = &x;
 };
 
 /** Constructor of class atom.
- *
  */
 atom::atom(atom *pointer)
@@ -41,4 +40,11 @@
   father = this;  // generally, father is itself
   Ancestor = NULL;
+  GraphNr = -1;
+  ComponentNr = NULL;
+  IsCyclic = false;
+  SeparationVertex = false;
+  LowpointNr = -1;
+  AdaptiveOrder = 0;
+  MaxOrder = false;
   type = pointer->type;  // copy element of atom
   x.CopyVector(&pointer->x); // copy coordination
@@ -54,5 +60,4 @@
 atom::~atom()
 {
-  Free((void **)&Name, "atom::~atom: *Name");
   Free((void **)&ComponentNr, "atom::~atom: *ComponentNr");
 };
Index: src/atom.hpp
===================================================================
--- src/atom.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/atom.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,55 @@
+/*
+ * atom.hpp
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#ifndef ATOM_HPP_
+#define ATOM_HPP_
+
+using namespace std;
+
+#include "periodentafel.hpp"
+#include "tesselation.hpp"
+#include "vector.hpp"
+
+/** Single atom.
+ * Class incorporates position, type
+ */
+class atom : public TesselPoint {
+  public:
+    Vector x;       //!< coordinate array of atom, giving position within cell
+    Vector v;       //!< velocity array of atom
+    element *type;  //!< pointing to element
+    atom *previous; //!< previous atom in molecule list
+    atom *next;     //!< next atom in molecule list
+    atom *father;   //!< In many-body bond order fragmentations points to originating atom
+    atom *Ancestor; //!< "Father" in Depth-First-Search
+    //char *Name;      //!< unique name used during many-body bond-order fragmentation, comes from TesselPoint
+    int FixedIon;   //!< config variable that states whether forces act on the ion or not
+    int *sort;      //!< sort criteria
+    //int nr;         //!< continuous, unique number, comes from TesselPoint
+    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
+    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
+    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
+    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
+    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
+    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
+    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
+
+  atom();
+  atom(class atom *pointer);
+  ~atom();
+
+  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
+  bool OutputXYZLine(ofstream *out) const;
+  atom *GetTrueFather();
+  bool Compare(atom &ptr);
+
+  private:
+};
+
+ostream & operator << (ostream &ost, const atom &a);
+
+#endif /* ATOM_HPP_ */
Index: src/bond.cpp
===================================================================
--- src/bond.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/bond.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -5,5 +5,5 @@
  */
 
-#include "molecules.hpp"
+#include "bond.hpp"
 
 
Index: src/bond.hpp
===================================================================
--- src/bond.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/bond.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,54 @@
+/*
+ * bond.hpp
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#ifndef BOND_HPP_
+#define BOND_HPP_
+
+using namespace std;
+
+#include "atom.hpp"
+
+/** Bonds between atoms.
+ * Class incorporates bonds between atoms in a molecule,
+ * used to derive tge fragments in many-body bond order
+ * calculations.
+ */
+class bond {
+  public:
+    atom *leftatom;    //!< first bond partner
+    atom *rightatom;  //!< second bond partner
+    bond *previous; //!< previous atom in molecule list
+    bond *next;     //!< next atom in molecule list
+    int HydrogenBond;  //!< Number of hydrogen atoms in the bond
+    int BondDegree;    //!< single, double, triple, ... bond
+    int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
+    bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
+    enum EdgeType Type;//!< whether this is a tree or back edge
+
+  atom * GetOtherAtom(atom *Atom) const;
+  bond * GetFirstBond();
+  bond * GetLastBond();
+
+  bool MarkUsed(enum Shading color);
+  enum Shading IsUsed();
+  void ResetUsed();
+  bool Contains(const atom *ptr);
+  bool Contains(const int nr);
+
+  bond();
+  bond(atom *left, atom *right);
+  bond(atom *left, atom *right, int degree);
+  bond(atom *left, atom *right, int degree, int number);
+  ~bond();
+
+  private:
+    enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
+};
+
+ostream & operator << (ostream &ost, const bond &b);
+
+#endif /* BOND_HPP_ */
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/boundary.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,592 +1,8 @@
 #include "boundary.hpp"
-#include "linkedcell.hpp"
-#include "molecules.hpp"
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_permutation.h>
-
-#define DEBUG 1
-#define DoSingleStepOutput 0
-#define DoTecplotOutput 1
-#define DoRaster3DOutput 1
-#define DoVRMLOutput 1
-#define TecplotSuffix ".dat"
-#define Raster3DSuffix ".r3d"
-#define VRMLSUffix ".wrl"
-#define HULLEPSILON 1e-7
-
-// ======================================== Points on Boundary =================================
-
-BoundaryPointSet::BoundaryPointSet()
-{
-  LinesCount = 0;
-  Nr = -1;
-}
-;
-
-BoundaryPointSet::BoundaryPointSet(atom *Walker)
-{
-  node = Walker;
-  LinesCount = 0;
-  Nr = Walker->nr;
-}
-;
-
-BoundaryPointSet::~BoundaryPointSet()
-{
-  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
-  if (!lines.empty())
-    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
-  node = NULL;
-}
-;
-
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
-{
-  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
-      << endl;
-  if (line->endpoints[0] == this)
-    {
-      lines.insert(LinePair(line->endpoints[1]->Nr, line));
-    }
-  else
-    {
-      lines.insert(LinePair(line->endpoints[0]->Nr, line));
-    }
-  LinesCount++;
-}
-;
-
-ostream &
-operator <<(ostream &ost, BoundaryPointSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.node->Name << "]";
-  return ost;
-}
-;
-
-// ======================================== Lines on Boundary =================================
-
-BoundaryLineSet::BoundaryLineSet()
-{
-  for (int i = 0; i < 2; i++)
-    endpoints[i] = NULL;
-  TrianglesCount = 0;
-  Nr = -1;
-}
-;
-
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
-{
-  // set number
-  Nr = number;
-  // set endpoints in ascending order
-  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
-  // add this line to the hash maps of both endpoints
-  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
-  Point[1]->AddLine(this); //
-  // clear triangles list
-  TrianglesCount = 0;
-  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
-}
-;
-
-BoundaryLineSet::~BoundaryLineSet()
-{
-  int Numbers[2];
-  Numbers[0] = endpoints[1]->Nr;
-  Numbers[1] = endpoints[0]->Nr;
-  for (int i = 0; i < 2; i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
-    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
-    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
-    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
-      if ((*Runner).second == this) {
-        endpoints[i]->lines.erase(Runner);
-        break;
-      }
-    if (endpoints[i]->lines.empty()) {
-      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
-      if (endpoints[i] != NULL) {
-        delete(endpoints[i]);
-        endpoints[i] = NULL;
-      } else
-        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
-  }
-  if (!triangles.empty())
-    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
-}
-;
-
-void
-BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
-{
-  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
-      << endl;
-  triangles.insert(TrianglePair(triangle->Nr, triangle));
-  TrianglesCount++;
-}
-;
-
-/** Checks whether we have a common endpoint with given \a *line.
- * \param *line other line to test
- * \return true - common endpoint present, false - not connected
- */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
-{
-  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
-    return true;
-  else
-    return false;
-};
-
-/** Checks whether the adjacent triangles of a baseline are convex or not.
- * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
- * If greater/equal M_PI than we are convex.
- * \param *out output stream for debugging
- * \return true - triangles are convex, false - concave or less than two triangles connected
- */
-bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
-{
-  Vector BaseLineNormal;
-  double angle = 0;
-  // get the two triangles
-  if (TrianglesCount != 2) {
-    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
-    return false;
-  }
-  // have a normal vector on the base line pointing outwards
-  BaseLineNormal.Zero();
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    BaseLineNormal.AddVector(&runner->second->NormalVector);
-  BaseLineNormal.Normalize();
-  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
-    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
-
-  if ((angle - M_PI) > -MYEPSILON)
-    return true;
-  else
-    return false;
-}
-
-/** Checks whether point is any of the two endpoints this line contains.
- * \param *point point to test
- * \return true - point is of the line, false - is not
- */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
-  for(int i=0;i<2;i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
-};
-
-ostream &
-operator <<(ostream &ost, BoundaryLineSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
-      << a.endpoints[1]->node->Name << "]";
-  return ost;
-}
-;
-
-// ======================================== Triangles on Boundary =================================
-
-
-BoundaryTriangleSet::BoundaryTriangleSet()
-{
-  for (int i = 0; i < 3; i++)
-    {
-      endpoints[i] = NULL;
-      lines[i] = NULL;
-    }
-  Nr = -1;
-}
-;
-
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
-{
-  // set number
-  Nr = number;
-  // set lines
-  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
-  for (int i = 0; i < 3; i++)
-    {
-      lines[i] = line[i];
-      lines[i]->AddTriangle(this);
-    }
-  // get ascending order of endpoints
-  map<int, class BoundaryPointSet *> OrderMap;
-  for (int i = 0; i < 3; i++)
-    // for all three lines
-    for (int j = 0; j < 2; j++)
-      { // for both endpoints
-        OrderMap.insert(pair<int, class BoundaryPointSet *> (
-            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
-        // and we don't care whether insertion fails
-      }
-  // set endpoints
-  int Counter = 0;
-  cout << Verbose(6) << " with end points ";
-  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
-      != OrderMap.end(); runner++)
-    {
-      endpoints[Counter] = runner->second;
-      cout << " " << *endpoints[Counter];
-      Counter++;
-    }
-  if (Counter < 3)
-    {
-      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
-          << endl;
-      //exit(1);
-    }
-  cout << "." << endl;
-}
-;
-
-BoundaryTriangleSet::~BoundaryTriangleSet()
-{
-  for (int i = 0; i < 3; i++) {
-    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-    lines[i]->triangles.erase(Nr);
-    if (lines[i]->triangles.empty()) {
-      if (lines[i] != NULL) {
-        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
-        delete (lines[i]);
-        lines[i] = NULL;
-      } else
-        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
-    } else
-      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
-  }
-}
-;
-
-/** Calculates the normal vector for this triangle.
- * Is made unique by comparison with \a OtherVector to point in the other direction.
- * \param &OtherVector direction vector to make normal vector unique.
- */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
-{
-  // get normal vector
-  NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x,
-      &endpoints[2]->node->x);
-
-  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (NormalVector.Projection(&OtherVector) > 0)
-    NormalVector.Scale(-1.);
-};
-
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
- * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
- * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
- * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
- * smaller than the first line.
- * \param *out output stream for debugging
- * \param *MolCenter offset vector of line
- * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
- * \param *Intersection intersection on plane on return
- * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
- */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
-{
-  Vector CrossPoint;
-  Vector helper;
-  int i=0;
-
-  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
-    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
-    return false;
-  }
-
-  // Calculate cross point between one baseline and the line from the third endpoint to intersection
-  do {
-    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
-    helper.CopyVector(&endpoints[(i+1)%3]->node->x);
-    helper.SubtractVector(&endpoints[i%3]->node->x);
-    i++;
-    if (i>3)
-      break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
-  if (i>3) {
-    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
-    exit(255);
-  }
-  CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
-    return true;
-  } else { // outside!
-    Intersection->Zero();
-    return false;
-  }
-};
-
-/** Checks whether lines is any of the three boundary lines this triangle contains.
- * \param *line line to test
- * \return true - line is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
-{
-  for(int i=0;i<3;i++)
-    if (line == lines[i])
-      return true;
-  return false;
-};
-
-/** Checks whether point is any of the three endpoints this triangle contains.
- * \param *point point to test
- * \return true - point is of the triangle, false - is not
- */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
-{
-  for(int i=0;i<3;i++)
-    if (point == endpoints[i])
-      return true;
-  return false;
-};
-
-ostream &
-operator <<(ostream &ost, BoundaryTriangleSet &a)
-{
-  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
-      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
-  return ost;
-}
-;
-
-
-// ============================ CandidateForTesselation =============================
-
-CandidateForTesselation::CandidateForTesselation(
-  atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
-) {
-  point = candidate;
-  BaseLine = line;
-  OptCenter.CopyVector(&OptCandidateCenter);
-  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
-}
-
-CandidateForTesselation::~CandidateForTesselation() {
-  point = NULL;
-  BaseLine = NULL;
-}
+
+#include<gsl/gsl_poly.h>
 
 // ========================================== F U N C T I O N S =================================
 
-/** Finds the endpoint two lines are sharing.
- * \param *line1 first line
- * \param *line2 second line
- * \return point which is shared or NULL if none
- */
-class BoundaryPointSet *
-GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2)
-{
-  class BoundaryLineSet * lines[2] =
-    { line1, line2 };
-  class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
-  for (int i = 0; i < 2; i++)
-    // for both lines
-    for (int j = 0; j < 2; j++)
-      { // for both endpoints
-        OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (
-            lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
-        if (!OrderTest.second)
-          { // if insertion fails, we have common endpoint
-            node = OrderTest.first->second;
-            cout << Verbose(5) << "Common endpoint of lines " << *line1
-                << " and " << *line2 << " is: " << *node << "." << endl;
-            j = 2;
-            i = 2;
-            break;
-          }
-      }
-  return node;
-}
-;
-
-/** Determines the boundary points of a cluster.
- * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
- * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
- * center and first and last point in the triple, it is thrown out.
- * \param *out output stream for debugging
- * \param *mol molecule structure representing the cluster
- */
-Boundaries *
-GetBoundaryPoints(ofstream *out, molecule *mol)
-{
-  atom *Walker = NULL;
-  PointMap PointsOnBoundary;
-  LineMap LinesOnBoundary;
-  TriangleMap TrianglesOnBoundary;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  Vector helper;
-
-  *out << Verbose(1) << "Finding all boundary points." << endl;
-  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
-  BoundariesTestPair BoundaryTestPair;
-  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
-  double radius, angle;
-  // 3a. Go through every axis
-  for (int axis = 0; axis < NDIM; axis++) {
-    AxisVector.Zero();
-    AngleReferenceVector.Zero();
-    AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
-
-    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
-
-    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
-    Walker = mol->start;
-    while (Walker->next != mol->end) {
-      Walker = Walker->next;
-      Vector ProjectedVector;
-      ProjectedVector.CopyVector(&Walker->x);
-      ProjectedVector.SubtractVector(MolCenter);
-      ProjectedVector.ProjectOntoPlane(&AxisVector);
-
-      // correct for negative side
-      radius = ProjectedVector.NormSquared();
-      if (fabs(radius) > MYEPSILON)
-        angle = ProjectedVector.Angle(&AngleReferenceVector);
-      else
-        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
-
-      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
-        angle = 2. * M_PI - angle;
-      }
-      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
-      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
-      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
-        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
-        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
-        *out << Verbose(2) << "New vector: " << *Walker << endl;
-        double tmp = ProjectedVector.NormSquared();
-        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
-          BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker;
-          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
-        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
-          helper.SubtractVector(MolCenter);
-          tmp = helper.NormSquared();
-          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
-          helper.SubtractVector(MolCenter);
-          if (helper.NormSquared() < tmp) {
-            BoundaryTestPair.first->second.second = Walker;
-            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
-          } else {
-            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
-          }
-        } else {
-          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
-        }
-      }
-    }
-    // printing all inserted for debugging
-    //    {
-    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
-    //      int i=0;
-    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-    //        if (runner != BoundaryPoints[axis].begin())
-    //          *out << ", " << i << ": " << *runner->second.second;
-    //        else
-    //          *out << i << ": " << *runner->second.second;
-    //        i++;
-    //      }
-    //      *out << endl;
-    //    }
-    // 3c. throw out points whose distance is less than the mean of left and right neighbours
-    bool flag = false;
-    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
-    do { // do as long as we still throw one out per round
-      flag = false;
-      Boundaries::iterator left = BoundaryPoints[axis].end();
-      Boundaries::iterator right = BoundaryPoints[axis].end();
-      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
-        // set neighbours correctly
-        if (runner == BoundaryPoints[axis].begin()) {
-          left = BoundaryPoints[axis].end();
-        } else {
-          left = runner;
-        }
-        left--;
-        right = runner;
-        right++;
-        if (right == BoundaryPoints[axis].end()) {
-          right = BoundaryPoints[axis].begin();
-        }
-        // check distance
-
-        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
-        {
-          Vector SideA, SideB, SideC, SideH;
-          SideA.CopyVector(&left->second.second->x);
-          SideA.SubtractVector(MolCenter);
-          SideA.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideA: ";
-          //          SideA.Output(out);
-          //          *out << endl;
-
-          SideB.CopyVector(&right->second.second->x);
-          SideB.SubtractVector(MolCenter);
-          SideB.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideB: ";
-          //          SideB.Output(out);
-          //          *out << endl;
-
-          SideC.CopyVector(&left->second.second->x);
-          SideC.SubtractVector(&right->second.second->x);
-          SideC.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideC: ";
-          //          SideC.Output(out);
-          //          *out << endl;
-
-          SideH.CopyVector(&runner->second.second->x);
-          SideH.SubtractVector(MolCenter);
-          SideH.ProjectOntoPlane(&AxisVector);
-          //          *out << "SideH: ";
-          //          SideH.Output(out);
-          //          *out << endl;
-
-          // calculate each length
-          double a = SideA.Norm();
-          //double b = SideB.Norm();
-          //double c = SideC.Norm();
-          double h = SideH.Norm();
-          // calculate the angles
-          double alpha = SideA.Angle(&SideH);
-          double beta = SideA.Angle(&SideC);
-          double gamma = SideB.Angle(&SideH);
-          double delta = SideC.Angle(&SideH);
-          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
-          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
-          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
-          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
-            // throw out point
-            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
-            BoundaryPoints[axis].erase(runner);
-            flag = true;
-          }
-        }
-      }
-    } while (flag);
-  }
-  delete(MolCenter);
-  return BoundaryPoints;
-}
-;
 
 /** Determines greatest diameters of a cluster defined by its convex envelope.
@@ -605,13 +21,10 @@
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = BoundaryPtr;
-  if (BoundaryPoints == NULL)
-    {
-      BoundaryFreeFlag = true;
-      BoundaryPoints = GetBoundaryPoints(out, mol);
-    }
-  else
-    {
-      *out << Verbose(1) << "Using given boundary points set." << endl;
-    }
+  if (BoundaryPoints == NULL) {
+    BoundaryFreeFlag = true;
+    BoundaryPoints = GetBoundaryPoints(out, mol);
+  } else {
+    *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -735,5 +148,5 @@
       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 << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *vrmlfile << "\t";
       }
@@ -790,5 +203,5 @@
       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 << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
         *rasterfile << "\t";
       }
@@ -821,9 +234,9 @@
     *out << Verbose(2) << "The following triangles were created:";
     int Counter = 1;
-    atom *Walker = NULL;
+    TesselPoint *Walker = NULL;
     for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
+      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << endl;
     }
     *tecplot << endl;
@@ -837,4 +250,178 @@
   }
 }
+
+
+/** Determines the boundary points of a cluster.
+ * Does a projection per axis onto the orthogonal plane, transforms into spherical coordinates, sorts them by the angle
+ * and looks at triples: if the middle has less a distance than the allowed maximum height of the triangle formed by the plane's
+ * center and first and last point in the triple, it is thrown out.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure representing the cluster
+ */
+Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol)
+{
+  atom *Walker = NULL;
+  PointMap PointsOnBoundary;
+  LineMap LinesOnBoundary;
+  TriangleMap TrianglesOnBoundary;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  Vector helper;
+
+  *out << Verbose(1) << "Finding all boundary points." << endl;
+  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
+  BoundariesTestPair BoundaryTestPair;
+  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
+  double radius, angle;
+  // 3a. Go through every axis
+  for (int axis = 0; axis < NDIM; axis++) {
+    AxisVector.Zero();
+    AngleReferenceVector.Zero();
+    AngleReferenceNormalVector.Zero();
+    AxisVector.x[axis] = 1.;
+    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+
+    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
+
+    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+    Walker = mol->start;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      Vector ProjectedVector;
+      ProjectedVector.CopyVector(&Walker->x);
+      ProjectedVector.SubtractVector(MolCenter);
+      ProjectedVector.ProjectOntoPlane(&AxisVector);
+
+      // correct for negative side
+      radius = ProjectedVector.NormSquared();
+      if (fabs(radius) > MYEPSILON)
+        angle = ProjectedVector.Angle(&AngleReferenceVector);
+      else
+        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+
+      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+        angle = 2. * M_PI - angle;
+      }
+      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
+        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
+        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
+        *out << Verbose(2) << "New vector: " << *Walker << endl;
+        double tmp = ProjectedVector.NormSquared();
+        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
+          BoundaryTestPair.first->second.first = tmp;
+          BoundaryTestPair.first->second.second = Walker;
+          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
+        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          helper.CopyVector(&Walker->x);
+          helper.SubtractVector(MolCenter);
+          tmp = helper.NormSquared();
+          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
+          helper.SubtractVector(MolCenter);
+          if (helper.NormSquared() < tmp) {
+            BoundaryTestPair.first->second.second = Walker;
+            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
+          } else {
+            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+          }
+        } else {
+          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
+        }
+      }
+    }
+    // printing all inserted for debugging
+    //    {
+    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+    //      int i=0;
+    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+    //        if (runner != BoundaryPoints[axis].begin())
+    //          *out << ", " << i << ": " << *runner->second.second;
+    //        else
+    //          *out << i << ": " << *runner->second.second;
+    //        i++;
+    //      }
+    //      *out << endl;
+    //    }
+    // 3c. throw out points whose distance is less than the mean of left and right neighbours
+    bool flag = false;
+    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
+    do { // do as long as we still throw one out per round
+      flag = false;
+      Boundaries::iterator left = BoundaryPoints[axis].end();
+      Boundaries::iterator right = BoundaryPoints[axis].end();
+      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        // set neighbours correctly
+        if (runner == BoundaryPoints[axis].begin()) {
+          left = BoundaryPoints[axis].end();
+        } else {
+          left = runner;
+        }
+        left--;
+        right = runner;
+        right++;
+        if (right == BoundaryPoints[axis].end()) {
+          right = BoundaryPoints[axis].begin();
+        }
+        // check distance
+
+        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
+        {
+          Vector SideA, SideB, SideC, SideH;
+          SideA.CopyVector(&left->second.second->x);
+          SideA.SubtractVector(MolCenter);
+          SideA.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideA: ";
+          //          SideA.Output(out);
+          //          *out << endl;
+
+          SideB.CopyVector(&right->second.second->x);
+          SideB.SubtractVector(MolCenter);
+          SideB.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideB: ";
+          //          SideB.Output(out);
+          //          *out << endl;
+
+          SideC.CopyVector(&left->second.second->x);
+          SideC.SubtractVector(&right->second.second->x);
+          SideC.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideC: ";
+          //          SideC.Output(out);
+          //          *out << endl;
+
+          SideH.CopyVector(&runner->second.second->x);
+          SideH.SubtractVector(MolCenter);
+          SideH.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideH: ";
+          //          SideH.Output(out);
+          //          *out << endl;
+
+          // calculate each length
+          double a = SideA.Norm();
+          //double b = SideB.Norm();
+          //double c = SideC.Norm();
+          double h = SideH.Norm();
+          // calculate the angles
+          double alpha = SideA.Angle(&SideH);
+          double beta = SideA.Angle(&SideC);
+          double gamma = SideB.Angle(&SideH);
+          double delta = SideC.Angle(&SideH);
+          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
+            // throw out point
+            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+            BoundaryPoints[axis].erase(runner);
+            flag = true;
+          }
+        }
+      }
+    } while (flag);
+  }
+  delete(MolCenter);
+  return BoundaryPoints;
+};
 
 /** Tesselates the convex boundary by finding all boundary points.
@@ -959,22 +546,16 @@
       << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
       << endl;
-  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
-      != TesselStruct->TrianglesOnBoundary.end(); runner++)
+  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
-      x.CopyVector(&runner->second->endpoints[0]->node->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));
+      x.CopyVector(runner->second->endpoints[0]->node->node);
+      x.SubtractVector(runner->second->endpoints[1]->node->node);
+      y.CopyVector(runner->second->endpoints[0]->node->node);
+      y.SubtractVector(runner->second->endpoints[2]->node->node);
+      a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
+      b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
+      c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
       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));
+      x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
+      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
       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)
@@ -1227,1940 +808,4 @@
 };
 
-
-// =========================================================== class TESSELATION ===========================================
-
-/** Constructor of class Tesselation.
- */
-Tesselation::Tesselation()
-{
-  PointsOnBoundaryCount = 0;
-  LinesOnBoundaryCount = 0;
-  TrianglesOnBoundaryCount = 0;
-  TriangleFilesWritten = 0;
-}
-;
-
-/** Constructor of class Tesselation.
- * We have to free all points, lines and triangles.
- */
-Tesselation::~Tesselation()
-{
-  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
-    if (runner->second != NULL) {
-      delete (runner->second);
-      runner->second = NULL;
-    } else
-      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
-  }
-}
-;
-
-/** Gueses first starting triangle of the convex envelope.
- * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
- * \param *out output stream for debugging
- * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
- */
-void
-Tesselation::GuessStartingTriangle(ofstream *out)
-{
-  // 4b. create a starting triangle
-  // 4b1. create all distances
-  DistanceMultiMap DistanceMMap;
-  double distance, tmp;
-  Vector PlaneVector, TrialVector;
-  PointMap::iterator A, B, C; // three nodes of the first triangle
-  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
-
-  // with A chosen, take each pair B,C and sort
-  if (A != PointsOnBoundary.end())
-    {
-      B = A;
-      B++;
-      for (; B != PointsOnBoundary.end(); B++)
-        {
-          C = B;
-          C++;
-          for (; C != PointsOnBoundary.end(); C++)
-            {
-              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
-              distance = tmp * tmp;
-              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
-              distance += tmp * tmp;
-              DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
-                  PointMap::iterator, PointMap::iterator> (B, C)));
-            }
-        }
-    }
-  //    // 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 << " leaves "
-                  << checker->second->node->Name << " outside the convex hull."
-                  << endl;
-              break;
-            }
-          else
-            { // note the sign for later
-              *out << Verbose(2) << "Current candidates: "
-                  << A->second->node->Name << ","
-                  << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leave "
-                  << checker->second->node->Name << " inside the convex hull."
-                  << endl;
-              sign = tmp;
-            }
-          // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
-          int innerpoint = 0;
-          if ((tmp < A->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < A->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x);
-          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
-              &A->second->node->x)) && (tmp
-              < baseline->second.first->second->node->x.DistanceSquared(
-                  &baseline->second.second->second->node->x)))
-            innerpoint++;
-          tmp = checker->second->node->x.DistanceSquared(
-              &baseline->second.second->second->node->x);
-          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
-              &baseline->second.first->second->node->x)) && (tmp
-              < baseline->second.second->second->node->x.DistanceSquared(
-                  &A->second->node->x)))
-            innerpoint++;
-          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
-          if (innerpoint == 3)
-            break;
-        }
-      // 5. come this far, all on same side? Then break 1. loop and construct triangle
-      if (checker == PointsOnBoundary.end())
-        {
-          *out << "Looks like we have a candidate!" << endl;
-          break;
-        }
-    }
-  if (baseline != DistanceMMap.end())
-    {
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = A->second;
-      BPS[1] = baseline->second.second->second;
-      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-      BPS[0] = baseline->second.first->second;
-      BPS[1] = A->second;
-      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-
-      // 4b3. insert created triangle
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-      TrianglesOnBoundaryCount++;
-      for (int i = 0; i < NDIM; i++)
-        {
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
-          LinesOnBoundaryCount++;
-        }
-
-      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
-    }
-  else
-    {
-      *out << Verbose(1) << "No starting triangle found." << endl;
-      exit(255);
-    }
-}
-;
-
-/** Tesselates the convex envelope of a cluster from a single starting triangle.
- * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
- * 2 triangles. Hence, we go through all current lines:
- * -# if the lines contains to only one triangle
- * -# We search all points in the boundary
- *    -# if the triangle 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)
- *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
- *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
- * \param *out output stream for debugging
- * \param *configuration for IsAngstroem
- * \param *mol the cluster as a molecule structure
- */
-void Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
-{
-  bool flag;
-  PointMap::iterator winner;
-  class BoundaryPointSet *peak = NULL;
-  double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *MolCenter = NULL;
-  LineMap::iterator LineChecker[2];
-
-  MolCenter = mol->DetermineCenterOfAll(out);
-  // create a first tesselation with the given BoundaryPoints
-  do {
-    flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
-      if (baseline->second->TrianglesCount == 1) {
-        // 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;
-
-        // get peak point with respect to this base line's only triangle
-        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
-        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;
-
-        // prepare some auxiliary vectors
-        Vector BaseLineCenter, BaseLine;
-        BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLineCenter.AddVector(&baseline->second->endpoints[1]->node->x);
-        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
-        BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x);
-        BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x);
-
-        // 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 << endl;
-
-        // normal vector of triangle
-        NormalVector.CopyVector(MolCenter);
-        NormalVector.SubtractVector(&CenterVector);
-        BTS->GetNormalVector(NormalVector);
-        NormalVector.CopyVector(&BTS->NormalVector);
-        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
-
-        // vector in propagation direction (out of triangle)
-        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector.MakeNormalVector(&BaseLine, &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 << endl;
-        winner = PointsOnBoundary.end();
-
-        // loop over all points and calculate angle between normal vector of new and present triangle
-        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) << ":" << endl;
-
-            // first check direction, so that triangles don't intersect
-            VirtualNormalVector.CopyVector(&target->second->node->x);
-            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
-            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
-            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
-            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
-            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
-              continue;
-            } else
-              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
-
-            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
-            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()) && (LineChecker[0]->second->TrianglesCount == 2))) {
-              *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;
-            }
-            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
-              *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)
-            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
-              *out << Verbose(4) << "Current target is peak!" << endl;
-              continue;
-            }
-
-            // check for linear dependence
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.SubtractVector(&target->second->node->x);
-            helper.CopyVector(&baseline->second->endpoints[1]->node->x);
-            helper.SubtractVector(&target->second->node->x);
-            helper.ProjectOntoPlane(&TempVector);
-            if (fabs(helper.NormSquared()) < MYEPSILON) {
-              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
-              continue;
-            }
-
-            // in case NOT both were found, 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);
-            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.AddVector(&baseline->second->endpoints[1]->node->x);
-            TempVector.AddVector(&target->second->node->x);
-            TempVector.Scale(1./3.);
-            TempVector.SubtractVector(MolCenter);
-            // make it always point outward
-            if (VirtualNormalVector.Projection(&TempVector) < 0)
-              VirtualNormalVector.Scale(-1.);
-            // calculate angle
-            TempAngle = NormalVector.Angle(&VirtualNormalVector);
-            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
-            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
-              SmallestAngle = TempAngle;
-              winner = target;
-              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
-            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
-              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
-              helper.CopyVector(&target->second->node->x);
-              helper.SubtractVector(&BaseLineCenter);
-              helper.ProjectOntoPlane(&BaseLine);
-              // ...the one with the smaller angle is the better candidate
-              TempVector.CopyVector(&target->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              TempAngle = TempVector.Angle(&helper);
-              TempVector.CopyVector(&winner->second->node->x);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              if (TempAngle < TempVector.Angle(&helper)) {
-                TempAngle = NormalVector.Angle(&VirtualNormalVector);
-                SmallestAngle = TempAngle;
-                winner = target;
-                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
-              } else
-                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
-            } else
-              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
-          }
-        } // end of loop over all boundary points
-
-        // 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);
-
-  // exit
-  delete(MolCenter);
-};
-
-/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
- * \param *out output stream for debugging
- * \param *mol molecule structure with atoms
- * \return true - all straddling points insert, false - something went wrong
- */
-bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
-{
-  Vector Intersection;
-  atom *Walker = mol->start;
-  Vector *MolCenter = mol->DetermineCenterOfAll(out);
-  while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
-    // get the next triangle
-    BTS = FindClosestTriangleToPoint(out, &Walker->x);
-    if (BTS == NULL) {
-      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
-      return false;
-    }
-    // get the intersection point
-    if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
-      // we have the intersection, check whether in- or outside of boundary
-      if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
-        // inside, next!
-        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
-      } else {
-        // outside!
-        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
-        class BoundaryLineSet *OldLines[3], *NewLines[3];
-        class BoundaryPointSet *OldPoints[3], *NewPoint;
-        // store the three old lines and old points
-        for (int i=0;i<3;i++) {
-          OldLines[i] = BTS->lines[i];
-          OldPoints[i] = BTS->endpoints[i];
-        }
-        // add Walker to boundary points
-        AddPoint(Walker);
-        if (BPS[0] == NULL)
-          NewPoint = BPS[0];
-        else
-          continue;
-        // remove triangle
-        TrianglesOnBoundary.erase(BTS->Nr);
-        // create three new boundary lines
-        for (int i=0;i<3;i++) {
-          BPS[0] = NewPoint;
-          BPS[1] = OldPoints[i];
-          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
-          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
-          LinesOnBoundaryCount++;
-        }
-        // create three new triangle with new point
-        for (int i=0;i<3;i++) { // find all baselines
-          BLS[0] = OldLines[i];
-          int n = 1;
-          for (int j=0;j<3;j++) {
-            if (NewLines[j]->IsConnectedTo(BLS[0])) {
-              if (n>2) {
-                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
-                return false;
-              } else
-                BLS[n++] = NewLines[j];
-            }
-          }
-          // create the triangle
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-          TrianglesOnBoundaryCount++;
-        }
-      }
-    } else { // something is wrong with FindClosestTriangleToPoint!
-      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
-      return false;
-    }
-  }
-
-  // exit
-  delete(MolCenter);
-  return true;
-};
-
-/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
- * \param *out output stream for debugging
- * \return true - all baselines were corrected, false - there are still concave pieces
- */
-bool Tesselation::CorrectConcaveBaselines(ofstream *out)
-{
-  class BoundaryTriangleSet *triangle[2];
-  class BoundaryLineSet *OldLines[4], *NewLine;
-  class BoundaryPointSet *OldPoints[2];
-  Vector BaseLineNormal;
-  class BoundaryLineSet *Base = NULL;
-  int OldTriangles[2], OldBaseLine;
-  int i;
-  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
-    Base = baseline->second;
-
-    // check convexity
-    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
-      *out << Verbose(3) << Base << " has two convex triangles." << endl;
-    } else { // not convex!
-      // get the two triangles
-      i=0;
-      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
-        triangle[i++] = runner->second;
-      // gather four endpoints and four lines
-      for (int j=0;j<4;j++)
-        OldLines[j] = NULL;
-      for (int j=0;j<2;j++)
-        OldPoints[j] = NULL;
-      i=0;
-      for (int m=0;m<2;m++) { // go over both triangles
-        for (int j=0;j<3;j++) { // all of their endpoints and baselines
-          if (triangle[m]->lines[j] != Base) // pick not the central baseline
-            OldLines[i++] = triangle[m]->lines[j];
-          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
-            OldPoints[m] = triangle[m]->endpoints[j];
-        }
-      }
-      if (i<4) {
-        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-        return false;
-      }
-      for (int j=0;j<4;j++)
-        if (OldLines[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
-          return false;
-        }
-      for (int j=0;j<2;j++)
-        if (OldPoints[j] == NULL) {
-          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
-          return false;
-        }
-
-      // remove triangles
-      for (int j=0;j<2;j++) {
-        OldTriangles[j] = triangle[j]->Nr;
-        TrianglesOnBoundary.erase(OldTriangles[j]);
-        triangle[j] = NULL;
-      }
-
-      // remove baseline
-      OldBaseLine = Base->Nr;
-      LinesOnBoundary.erase(OldBaseLine);
-      Base = NULL;
-
-      // construct new baseline (with same number as old one)
-      BPS[0] = OldPoints[0];
-      BPS[1] = OldPoints[1];
-      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
-      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
-
-      // construct new triangles with flipped baseline
-      i=-1;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=2;
-      if (BLS[0]->IsConnectedTo(OldLines[2]))
-        i=3;
-      if (i!=-1) {
-        BLS[0] = OldLines[0];
-        BLS[1] = OldLines[i];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
-
-        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
-        BLS[1] = OldLines[1];
-        BLS[2] = NewLine;
-        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
-        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
-      } else {
-        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
-        return false;
-      }
-    }
-  }
-  return true;
-};
-
-
-/** States whether point is in- or outside of a tesselated surface.
- * \param *pointer point to be checked
- * \return true - is inside, false - is outside
- */
-bool Tesselation::IsInside(Vector *pointer)
-{
-
-  // hier kommt dann Saskias Routine hin...
-
-  return true;
-};
-
-/** Finds the closest triangle to a given point.
- * \param *out output stream for debugging
- * \param *x second endpoint of line
- * \return pointer triangle that is closest, NULL if none was found
- */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
-{
-  class BoundaryTriangleSet *triangle = NULL;
-
-  // hier kommt dann Saskias Routine hin...
-
-  return triangle;
-};
-
-/** Adds an atom to the tesselation::PointsOnBoundary list.
- * \param *Walker atom to add
- */
-void
-Tesselation::AddPoint(atom *Walker)
-{
-  PointTestPair InsertUnique;
-  BPS[0] = new class BoundaryPointSet(Walker);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[0]));
-  if (InsertUnique.second) // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  else {
-    delete(BPS[0]);
-    BPS[0] = NULL;
-  }
-}
-;
-
-/** Adds point to Tesselation::PointsOnBoundary if not yet present.
- * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
- * @param Candidate point to add
- * @param n index for this point in Tesselation::TPS array
- */
-void
-Tesselation::AddTrianglePoint(atom* Candidate, int n)
-{
-  PointTestPair InsertUnique;
-  TPS[n] = new class BoundaryPointSet(Candidate);
-  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
-  if (InsertUnique.second) { // if new point was not present before, increase counter
-    PointsOnBoundaryCount++;
-  } else {
-    delete TPS[n];
-    cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
-    TPS[n] = (InsertUnique.first)->second;
-  }
-}
-;
-
-/** Function tries to add line from current Points in BPS to BoundaryLineSet.
- * If successful it raises the line count and inserts the new line into the BLS,
- * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
-  bool insertNewLine = true;
-
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine;
-    pair<LineMap::iterator,LineMap::iterator> FindPair;
-    FindPair = a->lines.equal_range(b->node->nr);
-
-    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-      // If there is a line with less than two attached triangles, we don't need a new line.
-      if (FindLine->second->TrianglesCount < 2) {
-        insertNewLine = false;
-        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
-
-        BPS[0] = FindLine->second->endpoints[0];
-        BPS[1] = FindLine->second->endpoints[1];
-        BLS[n] = FindLine->second;
-
-        break;
-      }
-    }
-  }
-
-  if (insertNewLine) {
-    AlwaysAddTriangleLine(a, b, n);
-  }
-}
-;
-
-/**
- * Adds lines from each of the current points in the BPS to BoundaryLineSet.
- * Raises the line count and inserts the new line into the BLS.
- *
- * @param *a first endpoint
- * @param *b second endpoint
- * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
-{
-  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
-  BPS[0] = a;
-  BPS[1] = b;
-  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
-  // add line to global map
-  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
-  // increase counter
-  LinesOnBoundaryCount++;
-};
-
-/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
- * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
- */
-void
-Tesselation::AddTriangle()
-{
-  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
-
-  // add triangle to global map
-  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
-  TrianglesOnBoundaryCount++;
-
-  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
-}
-;
-
-
-double det_get(gsl_matrix *A, int inPlace) {
-  /*
-  inPlace = 1 => A is replaced with the LU decomposed copy.
-  inPlace = 0 => A is retained, and a copy is used for LU.
-  */
-
-  double det;
-  int signum;
-  gsl_permutation *p = gsl_permutation_alloc(A->size1);
-  gsl_matrix *tmpA;
-
-  if (inPlace)
-  tmpA = A;
-  else {
-  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
-  gsl_matrix_memcpy(tmpA , A);
-  }
-
-
-  gsl_linalg_LU_decomp(tmpA , p , &signum);
-  det = gsl_linalg_LU_det(tmpA , signum);
-  gsl_permutation_free(p);
-  if (! inPlace)
-  gsl_matrix_free(tmpA);
-
-  return det;
-};
-
-void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
-{
-  gsl_matrix *A = gsl_matrix_calloc(3,3);
-  double m11, m12, m13, m14;
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m11 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m12 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
-  }
-  m13 = det_get(A, 1);
-
-  for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, b.x[i]);
-  }
-  m14 = det_get(A, 1);
-
-  if (fabs(m11) < MYEPSILON)
-    cerr << "ERROR: three points are colinear." << endl;
-
-  center->x[0] =  0.5 * m12/ m11;
-  center->x[1] = -0.5 * m13/ m11;
-  center->x[2] =  0.5 * m14/ m11;
-
-  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
-
-  gsl_matrix_free(A);
-};
-
-
-
-/**
- * Function returns center of sphere with RADIUS, which rests on points a, b, c
- * @param Center this vector will be used for return
- * @param a vector first point of triangle
- * @param b vector second point of triangle
- * @param c vector third point of triangle
- * @param *Umkreismittelpunkt new cneter point of circumference
- * @param Direction vector indicates up/down
- * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
- * @param Halfplaneindicator double indicates whether Direction is up or down
- * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
- * @param alpha double angle at a
- * @param beta double, angle at b
- * @param gamma, double, angle at c
- * @param Radius, double
- * @param Umkreisradius double radius of circumscribing circle
- */
-void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
-    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
-{
-  Vector TempNormal, helper;
-  double Restradius;
-  Vector OtherCenter;
-  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
-  Center->Zero();
-  helper.CopyVector(&a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(&b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(&c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-  NewUmkreismittelpunkt->CopyVector(Center);
-  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
-  // Here we calculated center of circumscribing circle, using barycentric coordinates
-  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
-
-  TempNormal.CopyVector(&a);
-  TempNormal.SubtractVector(&b);
-  helper.CopyVector(&a);
-  helper.SubtractVector(&c);
-  TempNormal.VectorProduct(&helper);
-  if (fabs(HalfplaneIndicator) < MYEPSILON)
-    {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
-        {
-          TempNormal.Scale(-1);
-        }
-    }
-  else
-    {
-      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
-        {
-          TempNormal.Scale(-1);
-        }
-    }
-
-  TempNormal.Normalize();
-  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
-  TempNormal.Scale(Restradius);
-  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
-
-  Center->AddVector(&TempNormal);
-  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
-  get_sphere(&OtherCenter, a, b, c, RADIUS);
-  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
-  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
-};
-
-
-/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
- * \param *Center new center on return
- * \param *a first point
- * \param *b second point
- * \param *c third point
- */
-void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
-{
-  Vector helper;
-  double alpha, beta, gamma;
-  Vector SideA, SideB, SideC;
-  SideA.CopyVector(b);
-  SideA.SubtractVector(c);
-  SideB.CopyVector(c);
-  SideB.SubtractVector(a);
-  SideC.CopyVector(a);
-  SideC.SubtractVector(b);
-  alpha = M_PI - SideB.Angle(&SideC);
-  beta = M_PI - SideC.Angle(&SideA);
-  gamma = M_PI - SideA.Angle(&SideB);
-  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
-  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
-    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
-
-  Center->Zero();
-  helper.CopyVector(a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
-  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-};
-
-/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
- * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
- * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
- * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
- * \param CircleCenter Center of the parameter circle
- * \param CirclePlaneNormal normal vector to plane of the parameter circle
- * \param CircleRadius radius of the parameter circle
- * \param NewSphereCenter new center of a circumcircle
- * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
- * \param NormalVector normal vector
- * \param SearchDirection search direction to make angle unique on return.
- * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
- */
-double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
-{
-  Vector helper;
-  double radius, alpha;
-
-  helper.CopyVector(&NewSphereCenter);
-  // test whether new center is on the parameter circle's plane
-  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
-    helper.ProjectOntoPlane(&CirclePlaneNormal);
-  }
-  radius = helper.ScalarProduct(&helper);
-  // test whether the new center vector has length of CircleRadius
-  if (fabs(radius - CircleRadius) > HULLEPSILON)
-    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&OldSphereCenter);
-  // make the angle unique by checking the halfplanes/search direction
-  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
-    alpha = 2.*M_PI - alpha;
-  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&OldSphereCenter);
-  helper.ProjectOntoPlane(&NormalVector);
-  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
-  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
-    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
-    return alpha;
-  } else {
-    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
-    return 2.*M_PI;
-  }
-};
-
-
-/** Checks whether the triangle consisting of the three atoms is already present.
- * Searches for the points in Tesselation::PointsOnBoundary and checks their
- * lines. If any of the three edges already has two triangles attached, false is
- * returned.
- * \param *out output stream for debugging
- * \param *Candidates endpoints of the triangle candidate
- * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
- *                 triangles exist which is the maximum for three points
- */
-int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
-  int adjacentTriangleCount = 0;
-  class BoundaryPointSet *Points[3];
-
-  //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
-  // builds a triangle point set (Points) of the end points
-  for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      Points[i] = FindPoint->second;
-    } else {
-      Points[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (Points[i] != NULL) {
-      for (int j = i; j < 3; j++) {
-        if (Points[j] != NULL) {
-          FindLine = Points[i]->lines.find(Points[j]->node->nr);
-          if (FindLine != Points[i]->lines.end()) {
-            for (; FindLine->first == Points[j]->node->nr; FindLine++) {
-              FindTriangle = FindLine->second->triangles.begin();
-              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == Points[0])
-                    || (FindTriangle->second->endpoints[0] == Points[1])
-                    || (FindTriangle->second->endpoints[0] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == Points[0])
-                    || (FindTriangle->second->endpoints[1] == Points[1])
-                    || (FindTriangle->second->endpoints[1] == Points[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == Points[0])
-                    || (FindTriangle->second->endpoints[2] == Points[1])
-                    || (FindTriangle->second->endpoints[2] == Points[2])
-                  )
-                ) {
-                  adjacentTriangleCount++;
-                }
-              }
-            }
-            // Only one of the triangle lines must be considered for the triangle count.
-            *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-            return adjacentTriangleCount;
-
-          }
-        }
-      }
-    }
-  }
-
-  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
-  return adjacentTriangleCount;
-};
-
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * Note that this function is for the starting triangle.
- * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
- * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
- * the center of the sphere is still fixed up to a single parameter. The band of possible values
- * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
- * us the "null" on this circle, the new center of the candidate point will be some way along this
- * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
- * by the normal vector of the base triangle that always points outwards by construction.
- * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
- * We construct the normal vector that defines the plane this circle lies in, it is just in the
- * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
- * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
- * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
- * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
- * both.
- * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
- * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
- * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
- * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
- * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
- * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
- * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
- * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
- * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
- * @param BaseLine BoundaryLineSet with the current base line
- * @param ThirdNode third atom to avoid in search
- * @param candidates list of equally good candidates to return
- * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
- * @param RADIUS radius of sphere
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_third_point_for_Tesselation(
-  Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
-  class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
-  double *ShortestAngle, const double RADIUS, LinkedCell *LC
-) {
-  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
-  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
-  Vector SphereCenter;
-  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-  Vector NewNormalVector;   // normal vector of the Candidate's triangle
-  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
-  LinkedAtoms *List = NULL;
-  double CircleRadius; // radius of this circle
-  double radius;
-  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-  atom *Candidate = NULL;
-  CandidateForTesselation *optCandidate = NULL;
-
-  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
-
-  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
-
-  // construct center of circle
-  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
-  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
-
-  // calculate squared radius atom *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-
-    // test whether old center is on the band's plane
-    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-    }
-    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
-    if (fabs(radius - CircleRadius) < HULLEPSILON) {
-
-      // check SearchDirection
-      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
-      }
-
-      // get cell for the starting atom
-      if (LC->SetIndexToVector(&CircleCenter)) {
-        for(int i=0;i<NDIM;i++) // store indices of this cell
-        N[i] = LC->n[i];
-        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-      } else {
-        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
-        return;
-      }
-      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-      //cout << Verbose(2) << "LC Intervals:";
-      for (int i=0;i<NDIM;i++) {
-        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-      }
-      //cout << endl;
-      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-            List = LC->GetCurrentCell();
-            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-            if (List != NULL) {
-              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-                Candidate = (*Runner);
-
-                // check for three unique points
-                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
-                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
-
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-
-                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
-                  ) {
-                    helper.CopyVector(&NewNormalVector);
-                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
-                    if (radius < RADIUS*RADIUS) {
-                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
-                      NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
-                      // OtherNewSphereCenter is created by the same vector just in the other direction
-                      helper.Scale(-1.);
-                      OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
-                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
-
-                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
-                      alpha = min(alpha, Otheralpha);
-                      // if there is a better candidate, drop the current list and add the new candidate
-                      // otherwise ignore the new candidate and keep the list
-                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
-                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
-                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
-                        } else {
-                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
-                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
-                        }
-                        // if there is an equal candidate, add it to the list without clearing the list
-                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
-                        } else {
-                          // remove all candidates from the list and then the list itself
-                          class CandidateForTesselation *remover = NULL;
-                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
-                            remover = *it;
-                            delete(remover);
-                          }
-                          candidates->clear();
-                          candidates->push_back(optCandidate);
-                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
-                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
-                        }
-                        *ShortestAngle = alpha;
-                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
-                      } else {
-                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
-                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
-                        } else {
-                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
-                        }
-                      }
-
-                    } else {
-                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
-                    }
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
-                  }
-                } else {
-                  if (ThirdNode != NULL) {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
-                  } else {
-                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
-                  }
-                }
-              }
-            }
-          }
-    } else {
-      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-    }
-  } else {
-    if (ThirdNode != NULL)
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
-    else
-      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
-  }
-
-  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
-  if (candidates->size() > 1) {
-    candidates->unique();
-    candidates->sort(sortCandidates);
-  }
-
-  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
-};
-
-struct Intersection {
-  Vector x1;
-  Vector x2;
-  Vector x3;
-  Vector x4;
-};
-
-/**
- * Intersection calculation function.
- *
- * @param x to find the result for
- * @param function parameter
- */
-double MinIntersectDistance(const gsl_vector * x, void *params) {
-  double retval = 0;
-  struct Intersection *I = (struct Intersection *)params;
-  Vector intersection;
-  Vector SideA,SideB,HeightA, HeightB;
-  for (int i=0;i<NDIM;i++)
-    intersection.x[i] = gsl_vector_get(x, i);
-
-  SideA.CopyVector(&(I->x1));
-  SideA.SubtractVector(&I->x2);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&I->x1);
-  HeightA.ProjectOntoPlane(&SideA);
-
-  SideB.CopyVector(&I->x3);
-  SideB.SubtractVector(&I->x4);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&I->x3);
-  HeightB.ProjectOntoPlane(&SideB);
-
-  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
-  //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
-
-  return retval;
-};
-
-
-/**
- * Calculates whether there is an intersection between two lines. The first line
- * always goes through point 1 and point 2 and the second line is given by the
- * connection between point 4 and point 5.
- *
- * @param point 1 of line 1
- * @param point 2 of line 1
- * @param point 1 of line 2
- * @param point 2 of line 2
- *
- * @return true if there is an intersection between the given lines, false otherwise
- */
-bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
-  bool result;
-
-  struct Intersection par;
-    par.x1.CopyVector(&point1);
-    par.x2.CopyVector(&point2);
-    par.x3.CopyVector(&point3);
-    par.x4.CopyVector(&point4);
-
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minex_func;
-
-    size_t iter = 0;
-    int status;
-    double size;
-
-    /* Starting point */
-    x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1.x[0]);
-    gsl_vector_set(x, 1, point1.x[1]);
-    gsl_vector_set(x, 2, point1.x[2]);
-
-    /* Set initial step sizes to 1 */
-    ss = gsl_vector_alloc(NDIM);
-    gsl_vector_set_all(ss, 1.0);
-
-    /* Initialize method and iterate */
-    minex_func.n = NDIM;
-    minex_func.f = &MinIntersectDistance;
-    minex_func.params = (void *)&par;
-
-    s = gsl_multimin_fminimizer_alloc(T, NDIM);
-    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
-
-    do {
-        iter++;
-        status = gsl_multimin_fminimizer_iterate(s);
-
-        if (status) {
-          break;
-        }
-
-        size = gsl_multimin_fminimizer_size(s);
-        status = gsl_multimin_test_size(size, 1e-2);
-
-        if (status == GSL_SUCCESS) {
-          cout << Verbose(2) << "converged to minimum" <<  endl;
-        }
-    } while (status == GSL_CONTINUE && iter < 100);
-
-    // check whether intersection is in between or not
-  Vector intersection, SideA, SideB, HeightA, HeightB;
-  double t1, t2;
-  for (int i = 0; i < NDIM; i++) {
-    intersection.x[i] = gsl_vector_get(s->x, i);
-  }
-
-  SideA.CopyVector(&par.x2);
-  SideA.SubtractVector(&par.x1);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&par.x1);
-
-  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
-
-  SideB.CopyVector(&par.x4);
-  SideB.SubtractVector(&par.x3);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&par.x3);
-
-  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
-
-  cout << Verbose(2) << "Intersection " << intersection << " is at "
-    << t1 << " for (" << point1 << "," << point2 << ") and at "
-    << t2 << " for (" << point3 << "," << point4 << "): ";
-
-  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
-    cout << "true intersection." << endl;
-    result = true;
-  } else {
-    cout << "intersection out of region of interest." << endl;
-    result = false;
-  }
-
-  // free minimizer stuff
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
-    gsl_multimin_fminimizer_free(s);
-
-  return result;
-}
-
-/** Finds the second point of starting triangle.
- * \param *a first atom
- * \param *Candidate pointer to candidate atom on return
- * \param Oben vector indicating the outside
- * \param Opt_Candidate reference to recommended candidate on return
- * \param Storage[3] array storing angles and other candidate information
- * \param RADIUS radius of virtual sphere
- * \param *LC LinkedCell structure with neighbouring atoms
- */
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
-{
-  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
-  Vector AngleCheck;
-  double norm = -1., angle;
-  LinkedAtoms *List = NULL;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
-    for(int i=0;i<NDIM;i++) // store indices of this cell
-      N[i] = LC->n[i];
-  } else {
-    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
-    return;
-  }
-  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-  cout << Verbose(3) << "LC Intervals from [";
-  for (int i=0;i<NDIM;i++) {
-  cout << " " << N[i] << "<->" << LC->N[i];
-  }
-  cout << "] :";
-  for (int i=0;i<NDIM;i++) {
-    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-  }
-  cout << endl;
-
-
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            Candidate = (*Runner);
-            // check if we only have one unique point yet ...
-            if (a != Candidate) {
-              // Calculate center of the circle with radius RADIUS through points a and Candidate
-              Vector OrthogonalizedOben, a_Candidate, Center;
-              double distance, scaleFactor;
-
-              OrthogonalizedOben.CopyVector(&Oben);
-              a_Candidate.CopyVector(&(a->x));
-              a_Candidate.SubtractVector(&(Candidate->x));
-              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
-              OrthogonalizedOben.Normalize();
-              distance = 0.5 * a_Candidate.Norm();
-              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
-              OrthogonalizedOben.Scale(scaleFactor);
-
-              Center.CopyVector(&(Candidate->x));
-              Center.AddVector(&(a->x));
-              Center.Scale(0.5);
-              Center.AddVector(&OrthogonalizedOben);
-
-              AngleCheck.CopyVector(&Center);
-              AngleCheck.SubtractVector(&(a->x));
-              norm = a_Candidate.Norm();
-              // second point shall have smallest angle with respect to Oben vector
-              if (norm < RADIUS*2.) {
-                angle = AngleCheck.Angle(&Oben);
-                if (angle < Storage[0]) {
-                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
-                  Opt_Candidate = Candidate;
-                  Storage[0] = angle;
-                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-                } else {
-                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
-                }
-              } else {
-                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
-              }
-            } else {
-              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
-            }
-          }
-        } else {
-          cout << Verbose(3) << "Linked cell list is empty." << 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* FirstPoint = NULL;
-  atom* SecondPoint = NULL;
-  atom* MaxAtom[NDIM];
-  double max_coordinate[NDIM];
-  Vector Oben;
-  Vector helper;
-  Vector Chord;
-  Vector SearchDirection;
-
-  Oben.Zero();
-
-  for (i = 0; i < 3; i++) {
-    MaxAtom[i] = NULL;
-    max_coordinate[i] = -1;
-  }
-
-  // 1. searching topmost atom with respect to each axis
-  for (int i=0;i<NDIM;i++) { // each axis
-    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
-    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
-      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
-        List = LC->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if ((*Runner)->x.x[i] > max_coordinate[i]) {
-              cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
-              max_coordinate[i] = (*Runner)->x.x[i];
-              MaxAtom[i] = (*Runner);
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-  }
-
-  cout << Verbose(2) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *MaxAtom[i] << "\t";
-  cout << endl;
-
-  BTS = NULL;
-  CandidateList *Opt_Candidates = new CandidateList();
-  for (int k=0;k<NDIM;k++) {
-    Oben.x[k] = 1.;
-    FirstPoint = MaxAtom[k];
-    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
-
-    double ShortestAngle;
-    atom* Opt_Candidate = NULL;
-    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-
-    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
-    SecondPoint = Opt_Candidate;
-    if (SecondPoint == NULL)  // have we found a second point?
-      continue;
-    else
-      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-    helper.CopyVector(&(FirstPoint->x));
-    helper.SubtractVector(&(SecondPoint->x));
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
-    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-
-    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-    Chord.SubtractVector(&(SecondPoint->x));
-    double radius = Chord.ScalarProduct(&Chord);
-    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
-    helper.CopyVector(&Oben);
-    helper.Scale(CircleRadius);
-    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-    // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
-
-    // adding point 1 and point 2 and the line between them
-    AddTrianglePoint(FirstPoint, 0);
-    AddTrianglePoint(SecondPoint, 1);
-    AddTriangleLine(TPS[0], TPS[1], 0);
-
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
-    Find_third_point_for_Tesselation(
-      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
-    );
-    cout << Verbose(1) << "List of third Points is ";
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        cout << " " << *(*it)->point;
-    }
-    cout << endl;
-
-    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-      // add third triangle point
-      AddTrianglePoint((*it)->point, 2);
-      // add the second and third line
-      AddTriangleLine(TPS[1], TPS[2], 1);
-      AddTriangleLine(TPS[0], TPS[2], 2);
-      // ... and triangles to the Maps of the Tesselation class
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      // ... and calculate its normal vector (with correct orientation)
-      (*it)->OptCenter.Scale(-1.);
-      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
-      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
-      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
-      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
-
-      // if we do not reach the end with the next step of iteration, we need to setup a new first line
-      if (it != Opt_Candidates->end()--) {
-        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
-        SecondPoint = (*it)->point;
-        // adding point 1 and point 2 and the line between them
-        AddTrianglePoint(FirstPoint, 0);
-        AddTrianglePoint(SecondPoint, 1);
-        AddTriangleLine(TPS[0], TPS[1], 0);
-      }
-      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
-    }
-    if (BTS != NULL) // we have created one starting triangle
-      break;
-    else {
-      // remove all candidates from the list and then the list itself
-      class CandidateForTesselation *remover = NULL;
-      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-        remover = *it;
-        delete(remover);
-      }
-      Opt_Candidates->clear();
-    }
-  }
-
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
-  }
-  delete(Opt_Candidates);
-  cout << Verbose(1) << "End of Find_starting_triangle\n";
-};
-
-/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
- * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
- * make it bigger (i.e. closing one (the baseline) and opening two new ones).
- * \param TPS[3] nodes of the triangle
- * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
- */
-bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
-{
-  bool result = false;
-  int counter = 0;
-
-  // check all three points
-  for (int i=0;i<3;i++)
-    for (int j=i+1; j<3; j++) {
-      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
-        LineMap::iterator FindLine;
-        pair<LineMap::iterator,LineMap::iterator> FindPair;
-        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
-        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
-          // If there is a line with less than two attached triangles, we don't need a new line.
-          if (FindLine->second->TrianglesCount < 2) {
-            counter++;
-            break;  // increase counter only once per edge
-          }
-        }
-      } else { // no line
-        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
-        result = true;
-      }
-    }
-  if (counter > 1) {
-    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
-    result = true;
-  }
-  return result;
-};
-
-
-/** This function finds a triangle to a line, adjacent to an existing one.
- * @param out output stream for debugging
- * @param *mol molecule with Atom's and Bond's
- * @param Line current baseline to search from
- * @param T current triangle which \a Line is edge of
- * @param RADIUS radius of the rolling ball
- * @param N number of found triangles
- * @param *filename filename base for intermediate envelopes
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-bool Tesselation::Find_next_suitable_triangle(ofstream *out,
-    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
-{
-  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
-  ofstream *tempstream = NULL;
-  char NumberName[255];
-  bool result = true;
-  CandidateList *Opt_Candidates = new CandidateList();
-
-  Vector CircleCenter;
-  Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
-  Vector SearchDirection;
-  Vector helper;
-  atom *ThirdNode = NULL;
-  LineMap::iterator testline;
-  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
-  double radius, CircleRadius;
-
-  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-  for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
-      ThirdNode = T.endpoints[i]->node;
-
-  // construct center of circle
-  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
-  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
-  CircleCenter.Scale(0.5);
-
-  // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
-  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
-
-  // calculate squared radius of circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
-    CirclePlaneNormal.Normalize();
-    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(&Line.endpoints[0]->node->x);
-    helper.SubtractVector(&ThirdNode->x);
-    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
-      SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
-    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
-      // rotated the wrong way!
-      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
-    }
-
-    // add third point
-    Find_third_point_for_Tesselation(
-      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
-      &ShortestAngle, RADIUS, LC
-    );
-
-  } else {
-    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
-  }
-
-  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
-    cerr << "WARNING: Could not find a suitable candidate." << endl;
-    return false;
-  }
-  cout << Verbose(1) << "Third Points are ";
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << " " << *(*it)->point;
-  }
-  cout << endl;
-
-  BoundaryLineSet *BaseRay = &Line;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
-    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
-    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
-
-    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-    atom *AtomCandidates[3];
-    AtomCandidates[0] = (*it)->point;
-    AtomCandidates[1] = BaseRay->endpoints[0]->node;
-    AtomCandidates[2] = BaseRay->endpoints[1]->node;
-    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
-
-    BTS = NULL;
-    // If there is no triangle, add it regularly.
-    if (existentTrianglesCount == 0) {
-      AddTrianglePoint((*it)->point, 0);
-      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-
-      AddTriangleLine(TPS[0], TPS[1], 0);
-      AddTriangleLine(TPS[0], TPS[2], 1);
-      AddTriangleLine(TPS[1], TPS[2], 2);
-
-      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-      AddTriangle();
-      (*it)->OptCenter.Scale(-1.);
-      BTS->GetNormalVector((*it)->OptCenter);
-      (*it)->OptCenter.Scale(-1.);
-
-      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-        << " for this triangle ... " << endl;
-      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
-    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
-        AddTrianglePoint((*it)->point, 0);
-        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
-        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
-
-        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
-        // i.e. at least one of the three lines must be present with TriangleCount <= 1
-        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
-          AddTriangleLine(TPS[0], TPS[1], 0);
-          AddTriangleLine(TPS[0], TPS[2], 1);
-          AddTriangleLine(TPS[1], TPS[2], 2);
-
-          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-          AddTriangle();
-
-          (*it)->OtherOptCenter.Scale(-1.);
-          BTS->GetNormalVector((*it)->OtherOptCenter);
-          (*it)->OtherOptCenter.Scale(-1.);
-
-          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
-          << " for this triangle ... " << endl;
-          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
-        } else {
-          cout << Verbose(1) << "WARNING: This triangle consisting of ";
-          cout << *(*it)->point << ", ";
-          cout << *BaseRay->endpoints[0]->node << " and ";
-          cout << *BaseRay->endpoints[1]->node << " ";
-          cout << "exists and is not added, as it does not seem helpful!" << endl;
-          result = false;
-        }
-    } else {
-      cout << Verbose(1) << "This triangle consisting of ";
-      cout << *(*it)->point << ", ";
-      cout << *BaseRay->endpoints[0]->node << " and ";
-      cout << *BaseRay->endpoints[1]->node << " ";
-      cout << "is invalid!" << endl;
-      result = false;
-    }
-
-    if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
-      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
-      if (DoTecplotOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(TecplotSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-
-      if (DoRaster3DOutput) {
-        string NameofTempFile(tempbasename);
-        NameofTempFile.append(NumberName);
-        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
-        NameofTempFile.erase(npos, 1);
-        NameofTempFile.append(Raster3DSuffix);
-        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
-        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-        write_raster3d_file(out, tempstream, this, mol);
-        // include the current position of the virtual sphere in the temporary raster3d file
-        // make the circumsphere's center absolute again
-        helper.CopyVector(&BaseRay->endpoints[0]->node->x);
-        helper.AddVector(&BaseRay->endpoints[1]->node->x);
-        helper.Scale(0.5);
-        (*it)->OptCenter.AddVector(&helper);
-        Vector *center = mol->DetermineCenterOfAll(out);
-        (*it)->OptCenter.SubtractVector(center);
-        delete(center);
-        // and add to file plus translucency object
-        *tempstream << "# current virtual sphere\n";
-        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
-          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
-          << "\t" << RADIUS << "\t1 0 0\n";
-        *tempstream << "9\n  terminating special property\n";
-        tempstream->close();
-        tempstream->flush();
-        delete(tempstream);
-      }
-      if (DoTecplotOutput || DoRaster3DOutput)
-        TriangleFilesWritten++;
-    }
-
-    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
-    BaseRay = BLS[0];
-  }
-
-  // remove all candidates from the list and then the list itself
-  class CandidateForTesselation *remover = NULL;
-  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-    remover = *it;
-    delete(remover);
-  }
-  delete(Opt_Candidates);
-  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
-  return result;
-};
-
-/**
- * Sort function for the candidate list.
- */
-bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
-  Vector BaseLineVector, OrthogonalVector, helper;
-  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-    //return false;
-    exit(1);
-  }
-  // create baseline vector
-  BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
-  BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  BaseLineVector.Normalize();
-
-  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
-  helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  helper.SubtractVector(&(candidate1->point->x));
-  OrthogonalVector.CopyVector(&helper);
-  helper.VectorProduct(&BaseLineVector);
-  OrthogonalVector.SubtractVector(&helper);
-  OrthogonalVector.Normalize();
-
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(&(candidate1->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double phi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-  helper.CopyVector(&(candidate2->point->x));
-  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-  double psi = BaseLineVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    psi = 2.*M_PI - psi;
-  }
-
-  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-
-  // return comparison
-  return phi < psi;
-}
-
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
  * \param *out output stream for debugging
@@ -3176,4 +821,8 @@
   bool freeTess = false;
   bool freeLC = false;
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+  int TriangleFilesWritten = 0;
+
   *out << Verbose(1) << "Entering search for non convex hull. " << endl;
   if (Tess == NULL) {
@@ -3193,10 +842,10 @@
   }
 
-  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
+  Tess->Find_starting_triangle(out, 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.
+      failflag = Tess->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
@@ -3214,5 +863,58 @@
       flag = false;
     }
-  }
+
+    // write temporary envelope
+    if ((DoSingleStepOutput && (Tess->TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+      class BoundaryTriangleSet *triangle = (Tess->TrianglesOnBoundary.end()--)->second;
+      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, triangle->endpoints[0]->node->Name, triangle->endpoints[1]->node->Name, triangle->endpoints[2]->node->Name);
+      if (DoTecplotOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(TecplotSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_tecplot_file(out, tempstream, Tess, mol, TriangleFilesWritten);
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+
+      if (DoRaster3DOutput) {
+        string NameofTempFile(filename);
+        NameofTempFile.append(NumberName);
+        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+        NameofTempFile.append(Raster3DSuffix);
+        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+        write_raster3d_file(out, tempstream, Tess, mol);
+//        // include the current position of the virtual sphere in the temporary raster3d file
+//        // make the circumsphere's center absolute again
+//        helper.CopyVector(BaseRay->endpoints[0]->node->node);
+//        helper.AddVector(BaseRay->endpoints[1]->node->node);
+//        helper.Scale(0.5);
+//        (*it)->OptCenter.AddVector(&helper);
+//        Vector *center = mol->DetermineCenterOfAll(out);
+//        (*it)->OptCenter.SubtractVector(center);
+//        delete(center);
+//        // and add to file plus translucency object
+//        *tempstream << "# current virtual sphere\n";
+//        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+//        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
+//          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
+//          << "\t" << RADIUS << "\t1 0 0\n";
+//        *tempstream << "9\n  terminating special property\n";
+        tempstream->close();
+        tempstream->flush();
+        delete(tempstream);
+      }
+      if (DoTecplotOutput || DoRaster3DOutput)
+        TriangleFilesWritten++;
+    }
+  }
+
+  // write final envelope
   if (filename != 0) {
     *out << Verbose(1) << "Writing final tecplot file\n";
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/boundary.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,9 +1,4 @@
 #ifndef BOUNDARY_HPP_
 #define BOUNDARY_HPP_
-
-class BoundaryPointSet;
-class BoundaryLineSet;
-class BoundaryTriangleSet;
-class CandidateForTesselation;
 
 // include config.h
@@ -17,120 +12,17 @@
 #include <deque>
 
-#include <gsl/gsl_poly.h>
-
 #include "linkedcell.hpp"
 #include "molecules.hpp"
+#include "tesselation.hpp"
+#include "tesselationhelpers.hpp"
 
-template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
-{
-  if (endpoint1->Nr < endpoint2->Nr) {
-    endpoints[0] = endpoint1;
-    endpoints[1] = endpoint2;
-  } else {
-    endpoints[0] = endpoint2;
-    endpoints[1] = endpoint1;
-  }
-};
-
-class BoundaryPointSet {
-  public:
-    BoundaryPointSet();
-    BoundaryPointSet(atom *Walker);
-    ~BoundaryPointSet();
-
-    void AddLine(class BoundaryLineSet *line);
-
-    LineMap lines;
-    int LinesCount;
-    atom *node;
-    int Nr;
-};
-
-class BoundaryLineSet {
-  public:
-    BoundaryLineSet();
-    BoundaryLineSet(class BoundaryPointSet *Point[2], int number);
-    ~BoundaryLineSet();
-
-    void AddTriangle(class BoundaryTriangleSet *triangle);
-    bool IsConnectedTo(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-    bool CheckConvexityCriterion(ofstream *out);
-
-    class BoundaryPointSet *endpoints[2];
-    TriangleMap triangles;
-    int TrianglesCount;
-    int Nr;
-};
-
-class BoundaryTriangleSet {
-  public:
-    BoundaryTriangleSet();
-    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
-    ~BoundaryTriangleSet();
-
-    void GetNormalVector(Vector &NormalVector);
-    bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
-    bool ContainsBoundaryLine(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-
-    class BoundaryPointSet *endpoints[3];
-    class BoundaryLineSet *lines[3];
-    Vector NormalVector;
-    int Nr;
-};
-
-
-class CandidateForTesselation {
-  public :
-  CandidateForTesselation(atom* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
-  ~CandidateForTesselation();
-  atom *point;
-  BoundaryLineSet *BaseLine;
-  Vector OptCenter;
-  Vector OtherOptCenter;
-};
-
-
-class Tesselation {
-  public:
-
-    Tesselation();
-    ~Tesselation();
-
-    void TesselateOnBoundary(ofstream *out, molecule *mol);
-    void GuessStartingTriangle(ofstream *out);
-    void AddPoint(atom * Walker);
-    void AddTrianglePoint(atom* Candidate, int n);
-    void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
-    void AddTriangle();
-    void Find_starting_triangle(ofstream *out, molecule* mol, const double RADIUS, LinkedCell *LC);
-    bool Find_next_suitable_triangle(ofstream *out, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename, LinkedCell *LC);
-    int CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]);
-    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);
-    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x);
-    bool IsInside(Vector *pointer);
-    bool InsertStraddlingPoints(ofstream *out, molecule *mol);
-    bool CorrectConcaveBaselines(ofstream *out);
-
-    PointMap PointsOnBoundary;
-    LineMap LinesOnBoundary;
-    TriangleMap TrianglesOnBoundary;
-    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
-    class BoundaryPointSet *BPS[2];
-    class BoundaryLineSet *BLS[3];
-    class BoundaryTriangleSet *BTS;
-    int PointsOnBoundaryCount;
-    int LinesOnBoundaryCount;
-    int TrianglesOnBoundaryCount;
-    int TriangleFilesWritten;
-};
-
-
-ostream & operator << (ostream &ost, BoundaryPointSet &a);
-ostream & operator << (ostream &ost, BoundaryLineSet &a);
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
-
+#define DEBUG 1
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
+#define DoRaster3DOutput 1
+#define DoVRMLOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
 
 double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration);
@@ -141,7 +33,5 @@
 void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *T, class LinkedCell *LC, const char *tempbasename, const double RADIUS);
 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
-bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess);
-bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
-bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2);
+Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
 
 #endif /*BOUNDARY_HPP_*/
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/builder.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1229,5 +1229,6 @@
   }
   molecules->SimpleMultiAdd(mol, src, N);
-  delete[](src);
+  delete(src);
+
   // ... and translate back
   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/config.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -276,7 +276,7 @@
 
 /** Readin of Thermostat related values from parameter file.
- * \param *source parameter file
- */
-void config::InitThermostats(ifstream *source)
+ * \param *fb file buffer containing the config file
+ */
+void config::InitThermostats(class ConfigFileBuffer *fb)
 {
   char *thermo = MallocString(12, "IonsInitRead: thermo");
@@ -284,5 +284,5 @@
 
   // read desired Thermostat from file along with needed additional parameters
-  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
+  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
       if (ThermostatImplemented[0] == 1) {
@@ -295,5 +295,5 @@
       if (ThermostatImplemented[1] == 1) {
         Thermostat = Woodcock;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -303,5 +303,5 @@
       if (ThermostatImplemented[2] == 1) {
         Thermostat = Gaussian;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -311,6 +311,6 @@
       if (ThermostatImplemented[3] == 1) {
         Thermostat = Langevin;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
-        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
+        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
           cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
         } else {
@@ -324,5 +324,5 @@
       if (ThermostatImplemented[4] == 1) {
         Thermostat = Berendsen;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
       } else {
         cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
@@ -332,5 +332,5 @@
       if (ThermostatImplemented[5] == 1) {
         Thermostat = NoseHoover;
-        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
+        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
         alpha = 0.;
       } else {
@@ -704,4 +704,6 @@
     return;
   }
+  file->close();
+  delete(file);
   RetrieveConfigPathAndName(filename);
 
@@ -722,5 +724,5 @@
   double value[3];
   
-  InitThermostats(file);
+  InitThermostats(FileBuffer);
   
   /* Namen einlesen */
Index: src/ellipsoid.cpp
===================================================================
--- src/ellipsoid.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/ellipsoid.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -6,4 +6,8 @@
  */
 
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_vector.h>
+
+#include "boundary.hpp"
 #include "ellipsoid.hpp"
 
@@ -221,6 +225,6 @@
   set<int>::iterator current;
   int index;
-  atom *Candidate = NULL;
-  LinkedAtoms *List = NULL;
+  TesselPoint *Candidate = NULL;
+  LinkedNodes *List = NULL;
   *out << Verbose(2) << "Begin of PickRandomPointSet" << endl;
 
@@ -288,9 +292,9 @@
 //            else
 //              *out << Verbose(2) << "Cell is empty ... " << endl;
-            for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
               if ((current != PickedAtomNrs.end()) && (*current == index)) {
                 Candidate = (*Runner);
                 *out << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl;
-                x[PointsPicked++].CopyVector(&(Candidate->x));    // we have one more atom picked
+                x[PointsPicked++].CopyVector(Candidate->node);    // we have one more atom picked
                 current++;    // next pre-picked atom
               }
@@ -338,5 +342,5 @@
       //*out << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
       if (value > threshold) {
-        x[PointsPicked].CopyVector(&(Runner->second->node->x));
+        x[PointsPicked].CopyVector(Runner->second->node->node);
         PointsPicked++;
         //*out << "IN." << endl;
@@ -375,5 +379,5 @@
   Center.Zero();
   for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
-    Center.AddVector(&Runner->second->node->x);
+    Center.AddVector(Runner->second->node->node);
   Center.Scale(1./T->PointsOnBoundaryCount);
   *out << Verbose(1) << "Center is at " << Center << "." << endl;
Index: src/ellipsoid.hpp
===================================================================
--- src/ellipsoid.hpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/ellipsoid.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -14,7 +14,5 @@
 #endif
 
-#include <cstdlib>
-
-#include "boundary.hpp"
+#include "linkedcell.hpp"
 #include "vector.hpp"
 
Index: src/linkedcell.cpp
===================================================================
--- src/linkedcell.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/linkedcell.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,4 +1,8 @@
 #include "linkedcell.hpp"
 #include "molecules.hpp"
+#include "tesselation.hpp"
+
+// ========================================================= class LinkedCell ===========================================
+
 
 /** Constructor for class LinkedCell.
@@ -16,10 +20,10 @@
 
 /** Puts all atoms in \a *mol into a linked cell list with cell's lengths of \a RADIUS
- * \param *mol molecule structure with all Atom's
+ * \param *set LCNodeSet class with all LCNode's
  * \param RADIUS edge length of cells
  */
-LinkedCell::LinkedCell(molecule *mol, double radius)
+LinkedCell::LinkedCell(PointCloud *set, double radius)
 {
-  atom *Walker = NULL;
+  TesselPoint *Walker = NULL;
 
   RADIUS = radius;
@@ -31,26 +35,29 @@
   min.Zero();
   cout << Verbose(1) << "Begin of LinkedCell" << endl;
-  if (mol->start->next == mol->end) {
-    cerr << "ERROR: molecule contains no atoms!" << endl;
+  if (set->IsEmpty()) {
+    cerr << "ERROR: set contains no linked cell nodes!" << endl;
     return;
   }
   // 1. find max and min per axis of atoms
-  Walker = mol->start->next;
+  set->GoToFirst();
+  Walker = set->GetPoint();
   for (int i=0;i<NDIM;i++) {
-    max.x[i] = Walker->x.x[i];
-    min.x[i] = Walker->x.x[i];
+    max.x[i] = Walker->node->x[i];
+    min.x[i] = Walker->node->x[i];
   }
-  while (Walker != mol->end) {
+  set->GoToFirst();
+  while (!set->IsLast()) {
+    Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      if (max.x[i] < Walker->x.x[i])
-        max.x[i] = Walker->x.x[i];
-      if (min.x[i] > Walker->x.x[i])
-        min.x[i] = Walker->x.x[i];
+      if (max.x[i] < Walker->node->x[i])
+        max.x[i] = Walker->node->x[i];
+      if (min.x[i] > Walker->node->x[i])
+        min.x[i] = Walker->node->x[i];
     }
-    Walker = Walker->next;
+    set->GoToNext();
   }
   cout << Verbose(2) << "Bounding box is " << min << " and " << max << "." << endl;
 
-  // 2. find then umber of cells per axis
+  // 2. find then number of cells per axis
   for (int i=0;i<NDIM;i++) {
     N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
@@ -64,5 +71,5 @@
     return;
   }
-  LC = new LinkedAtoms[N[0]*N[1]*N[2]];
+  LC = new LinkedNodes[N[0]*N[1]*N[2]];
   for (index=0;index<N[0]*N[1]*N[2];index++) {
     LC [index].clear();
@@ -72,13 +79,14 @@
   // 4. put each atom into its respective cell
   cout << Verbose(2) << "Filling cells ... ";
-  Walker = mol->start;
-  while (Walker->next != mol->end) {
-    Walker = Walker->next;
+  set->GoToFirst();
+  while (!set->IsLast()) {
+    Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
+      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
     }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
     LC[index].push_back(Walker);
     //cout << Verbose(2) << *Walker << " goes into cell " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << index << "." << endl;
+    set->GoToNext();
   }
   cout << "done."  << endl;
@@ -118,5 +126,5 @@
  * \return LinkedAtoms pointer to current cell, NULL if LinkedCell::n[] are out of bounds.
  */
-LinkedAtoms* LinkedCell::GetCurrentCell()
+LinkedNodes* LinkedCell::GetCurrentCell()
 {
   if (CheckBounds()) {
@@ -128,21 +136,21 @@
 };
 
-/** Calculates the index for a given atom *Walker.
- * \param *Walker atom to set index to
+/** Calculates the index for a given LCNode *Walker.
+ * \param *Walker LCNode to set index tos
  * \return if the atom is also found in this cell - true, else - false
  */
-bool LinkedCell::SetIndexToAtom(atom *Walker)
+bool LinkedCell::SetIndexToNode(TesselPoint *Walker)
 {
   bool status = false;
   for (int i=0;i<NDIM;i++) {
-    n[i] = (int)floor((Walker->x.x[i] - min.x[i])/RADIUS);
+    n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
   }
   index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
   if (CheckBounds()) {
-    for (LinkedAtoms::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
+    for (LinkedNodes::iterator Runner = LC[index].begin(); Runner != LC[index].end(); Runner++)
       status = status || ((*Runner) == Walker);
     return status;
   } else {
-    cerr << Verbose(1) << "ERROR: Atom "<< *Walker << " at " << Walker->x << " is out of bounds." << endl;
+    cerr << Verbose(1) << "ERROR: Node at " << *Walker << " is out of bounds." << endl;
     return false;
   }
Index: src/linkedcell.hpp
===================================================================
--- src/linkedcell.hpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/linkedcell.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,4 +1,16 @@
+/*
+ * linkedcell.hpp
+ *
+ *  If the linked cell should be usable, the class has to inherit LCNodeSet and the nodes (containing the Vectors) have to inherit LCNode. This works well
+ *  for molecule and atom classes.
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
 #ifndef LINKEDCELL_HPP_
 #define LINKEDCELL_HPP_
+
+using namespace std;
 
 // include config.h
@@ -7,13 +19,22 @@
 #endif
 
-#include "molecules.hpp"
+#include <list>
 
-#define LinkedAtoms list <atom *>
+#include "defs.hpp"
+#include "helpers.hpp"
+#include "vector.hpp"
 
-class LinkedCell{
+class TesselPoint;
+class PointCloud;
+
+#define LinkedNodes list<TesselPoint *>
+
+/** Linked Cell class for containing Vectors in real space efficiently.
+ */
+class LinkedCell {
   public:
     Vector max;       // upper boundary
     Vector min;       // lower boundary
-    LinkedAtoms *LC;  // linked cell list
+    LinkedNodes *LC;  // linked cell list
     double RADIUS;    // cell edge length
     int N[NDIM];      // number of cells per axis
@@ -22,15 +43,15 @@
 
     LinkedCell();
-    LinkedCell(molecule *mol, double RADIUS);
+    LinkedCell(PointCloud *set, double RADIUS);
     ~LinkedCell();
-    LinkedAtoms* GetCurrentCell();
-    bool SetIndexToAtom(atom *Walker);
+    LinkedNodes* GetCurrentCell();
+    bool SetIndexToNode(TesselPoint *Walker);
     bool SetIndexToVector(Vector *x);
     bool CheckBounds();
 
     // not implemented yet
-    bool AddAtom(atom *Walker);
-    bool DeleteAtom(atom *Walker);
-    bool MoveAtom(atom *Walker);
+    bool AddNode(Vector *Walker);
+    bool DeleteNode(Vector *Walker);
+    bool MoveNode(Vector *Walker);
 };
 
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/molecules.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -42,4 +42,5 @@
   end->father = NULL;
   link(start,end);
+  InternalPointer = start;
   // init bond chain list
   first = new bond(start, end, 1, -1);
@@ -82,4 +83,82 @@
   delete(end);
   delete(start);
+};
+
+
+/** Determine center of all atoms.
+ * \param *out output stream for debugging
+ * \return pointer to allocated with central coordinates
+ */
+Vector *molecule::GetCenter(ofstream *out)
+{
+  Vector *center = DetermineCenterOfAll(out);
+  return center;
+};
+
+/** Return current atom in the list.
+ * \return pointer to atom or NULL if none present
+ */
+TesselPoint *molecule::GetPoint()
+{
+  if ((InternalPointer != start) && (InternalPointer != end))
+    return InternalPointer;
+  else
+    return NULL;
+};
+
+/** Return pointer to one after last atom in the list.
+ * \return pointer to end marker
+ */
+TesselPoint *molecule::GetTerminalPoint()
+{
+  return end;
+};
+
+/** Go to next atom.
+ * Stops at last one.
+ */
+void molecule::GoToNext()
+{
+  if (InternalPointer->next != end)
+    InternalPointer = InternalPointer->next;
+};
+
+/** Go to previous atom.
+ * Stops at first one.
+ */
+void molecule::GoToPrevious()
+{
+  if (InternalPointer->previous != start)
+    InternalPointer = InternalPointer->previous;
+};
+
+/** Goes to first atom.
+ */
+void molecule::GoToFirst()
+{
+  InternalPointer = start->next;
+};
+
+/** Goes to last atom.
+ */
+void molecule::GoToLast()
+{
+  InternalPointer = end->previous;
+};
+
+/** Checks whether we have any atoms in molecule.
+ * \return true - no atoms, false - not empty
+ */
+bool molecule::IsEmpty()
+{
+  return (start->next == end);
+};
+
+/** Checks whether we are at the last atom
+ * \return true - current atom is last one, false - is not last one
+ */
+bool molecule::IsLast()
+{
+  return (InternalPointer->next == end);
 };
 
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/molecules.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -25,12 +25,14 @@
 #include <vector>
 
+#include "atom.hpp"
+#include "bond.hpp"
 #include "helpers.hpp"
+#include "linkedcell.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
+#include "tesselation.hpp"
 #include "vector.hpp"
 
-class atom;
-class bond;
 class config;
 class molecule;
@@ -38,5 +40,4 @@
 class MoleculeListClass;
 class Verbose;
-class Tesselation;
 
 /******************************** Some definitions for easier reading **********************************/
@@ -58,20 +59,4 @@
 #define BoundariesTestPair pair< Boundaries::iterator, bool>
 
-#define PointMap map < int, class BoundaryPointSet * >
-#define PointPair pair < int, class BoundaryPointSet * >
-#define PointTestPair pair < PointMap::iterator, bool >
-#define CandidateList list <class CandidateForTesselation *>
-
-#define LineMap multimap < int, class BoundaryLineSet * >
-#define LinePair pair < int, class BoundaryLineSet * >
-#define LineTestPair pair < LineMap::iterator, bool >
-
-#define TriangleMap map < int, class BoundaryTriangleSet * >
-#define TrianglePair pair < int, class BoundaryTriangleSet * >
-#define TriangleTestPair pair < TrianglePair::iterator, bool >
-
-#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
-#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
-
 #define MoleculeList list <molecule *>
 #define MoleculeListTest pair <MoleculeList::iterator, bool>
@@ -125,81 +110,4 @@
 };
 
-/** Single atom.
- * Class incoporates position, type
- */
-class atom {
-  public:
-    Vector x;       //!< coordinate array of atom, giving position within cell
-    Vector v;       //!< velocity array of atom
-    element *type;  //!< pointing to element
-    atom *previous; //!< previous atom in molecule list
-    atom *next;     //!< next atom in molecule list
-    atom *father;   //!< In many-body bond order fragmentations points to originating atom
-    atom *Ancestor; //!< "Father" in Depth-First-Search
-    char *Name;      //!< unique name used during many-body bond-order fragmentation
-    int FixedIon;   //!< config variable that states whether forces act on the ion or not
-    int *sort;      //!< sort criteria
-    int nr;         //!< continuous, unique number
-    int GraphNr;      //!< unique number, given in DepthFirstSearchAnalysis()
-    int *ComponentNr;//!< belongs to this nonseparable components, given in DepthFirstSearchAnalysis() (if more than one, then is SeparationVertex)
-    int LowpointNr; //!< needed in DepthFirstSearchAnalysis() to detect nonseparable components, is the lowest possible number of an atom to reach via tree edges only followed by at most one back edge.
-    bool SeparationVertex; //!< whether this atom separates off subsets of atoms or not, determined in DepthFirstSearchAnalysis()
-    bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
-    unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
-    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
-
-  atom();
-  atom(class atom *pointer);
-  ~atom();
-
-  bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
-  bool OutputXYZLine(ofstream *out) const;
-  atom *GetTrueFather();
-  bool Compare(atom &ptr);
-
-  private:
-};
-
-ostream & operator << (ostream &ost, const atom &a);
-
-/** Bonds between atoms.
- * Class incorporates bonds between atoms in a molecule,
- * used to derive tge fragments in many-body bond order
- * calculations.
- */
-class bond {
-  public:
-    atom *leftatom;    //!< first bond partner
-    atom *rightatom;  //!< second bond partner
-    bond *previous; //!< previous atom in molecule list
-    bond *next;     //!< next atom in molecule list
-    int HydrogenBond;  //!< Number of hydrogen atoms in the bond
-    int BondDegree;    //!< single, double, triple, ... bond
-    int nr;           //!< unique number in a molecule, updated by molecule::CreateAdjacencyList()
-    bool Cyclic;      //!< flag whether bond is part of a cycle or not, given in DepthFirstSearchAnalysis()
-    enum EdgeType Type;//!< whether this is a tree or back edge
-
-  atom * GetOtherAtom(atom *Atom) const;
-  bond * GetFirstBond();
-  bond * GetLastBond();
-
-  bool MarkUsed(enum Shading color);
-  enum Shading IsUsed();
-  void ResetUsed();
-  bool Contains(const atom *ptr);
-  bool Contains(const int nr);
-
-  bond();
-  bond(atom *left, atom *right);
-  bond(atom *left, atom *right, int degree);
-  bond(atom *left, atom *right, int degree, int number);
-  ~bond();
-
-  private:
-    enum Shading Used;        //!< marker in depth-first search, DepthFirstSearchAnalysis()
-};
-
-
-ostream & operator << (ostream &ost, const bond &b);
 
 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
@@ -210,5 +118,5 @@
  * Class incorporates number of types
  */
-class molecule {
+class molecule : public PointCloud {
   public:
     double cell_size[6];//!< cell size
@@ -234,7 +142,19 @@
     char name[MAXSTRINGSIZE];         //!< arbitrary name
     int IndexNr;        //!< index of molecule in a MoleculeListClass
+    class Tesselation *TesselStruct;
 
   molecule(periodentafel *teil);
   ~molecule();
+
+  // re-definition of virtual functions from PointCloud
+  Vector *GetCenter(ofstream *out);
+  TesselPoint *GetPoint();
+  TesselPoint *GetTerminalPoint();
+  void GoToNext();
+  void GoToPrevious();
+  void GoToFirst();
+  void GoToLast();
+  bool IsEmpty();
+  bool IsLast();
 
   /// remove atoms from molecule.
@@ -350,4 +270,5 @@
   private:
   int last_atom;      //!< number given to last atom
+  atom *InternalPointer;  //!< internal pointer for PointCloud
 };
 
@@ -524,5 +445,5 @@
   char *GetDefaultPath() const;
   void SetDefaultPath(const char *path);
-  void InitThermostats(ifstream *source);
+  void InitThermostats(class ConfigFileBuffer *fb);
 };
 
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/tesselation.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,1955 @@
+/*
+ * tesselation.cpp
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#include "tesselation.hpp"
+
+// ======================================== Points on Boundary =================================
+
+BoundaryPointSet::BoundaryPointSet()
+{
+  LinesCount = 0;
+  Nr = -1;
+}
+;
+
+BoundaryPointSet::BoundaryPointSet(TesselPoint *Walker)
+{
+  node = Walker;
+  LinesCount = 0;
+  Nr = Walker->nr;
+}
+;
+
+BoundaryPointSet::~BoundaryPointSet()
+{
+  cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  if (!lines.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
+  node = NULL;
+}
+;
+
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+{
+  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+      << endl;
+  if (line->endpoints[0] == this)
+    {
+      lines.insert(LinePair(line->endpoints[1]->Nr, line));
+    }
+  else
+    {
+      lines.insert(LinePair(line->endpoints[0]->Nr, line));
+    }
+  LinesCount++;
+}
+;
+
+ostream &
+operator <<(ostream &ost, BoundaryPointSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.node->Name << "]";
+  return ost;
+}
+;
+
+// ======================================== Lines on Boundary =================================
+
+BoundaryLineSet::BoundaryLineSet()
+{
+  for (int i = 0; i < 2; i++)
+    endpoints[i] = NULL;
+  TrianglesCount = 0;
+  Nr = -1;
+}
+;
+
+BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], int number)
+{
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point[0], Point[1]);
+  // add this line to the hash maps of both endpoints
+  Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+  Point[1]->AddLine(this); //
+  // clear triangles list
+  TrianglesCount = 0;
+  cout << Verbose(5) << "New Line with endpoints " << *this << "." << endl;
+}
+;
+
+BoundaryLineSet::~BoundaryLineSet()
+{
+  int Numbers[2];
+  Numbers[0] = endpoints[1]->Nr;
+  Numbers[1] = endpoints[0]->Nr;
+  for (int i = 0; i < 2; i++) {
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+    // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
+    pair<LineMap::iterator, LineMap::iterator> erasor = endpoints[i]->lines.equal_range(Numbers[i]);
+    for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
+      if ((*Runner).second == this) {
+        endpoints[i]->lines.erase(Runner);
+        break;
+      }
+    if (endpoints[i]->lines.empty()) {
+      cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+      if (endpoints[i] != NULL) {
+        delete(endpoints[i]);
+        endpoints[i] = NULL;
+      } else
+        cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+  }
+  if (!triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
+}
+;
+
+void
+BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+{
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
+      << endl;
+  triangles.insert(TrianglePair(triangle->Nr, triangle));
+  TrianglesCount++;
+}
+;
+
+/** Checks whether we have a common endpoint with given \a *line.
+ * \param *line other line to test
+ * \return true - common endpoint present, false - not connected
+ */
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+{
+  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
+    return true;
+  else
+    return false;
+};
+
+/** Checks whether the adjacent triangles of a baseline are convex or not.
+ * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * If greater/equal M_PI than we are convex.
+ * \param *out output stream for debugging
+ * \return true - triangles are convex, false - concave or less than two triangles connected
+ */
+bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
+{
+  Vector BaseLineNormal;
+  double angle = 0;
+  // get the two triangles
+  if (TrianglesCount != 2) {
+    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    return false;
+  }
+  // have a normal vector on the base line pointing outwards
+  BaseLineNormal.Zero();
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    BaseLineNormal.AddVector(&runner->second->NormalVector);
+  BaseLineNormal.Normalize();
+  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
+
+  if ((angle - M_PI) > -MYEPSILON)
+    return true;
+  else
+    return false;
+}
+
+/** Checks whether point is any of the two endpoints this line contains.
+ * \param *point point to test
+ * \return true - point is of the line, false - is not
+ */
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<2;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+
+ostream &
+operator <<(ostream &ost, BoundaryLineSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
+  return ost;
+}
+;
+
+// ======================================== Triangles on Boundary =================================
+
+
+BoundaryTriangleSet::BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++)
+    {
+      endpoints[i] = NULL;
+      lines[i] = NULL;
+    }
+  Nr = -1;
+}
+;
+
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number)
+{
+  // set number
+  Nr = number;
+  // set lines
+  cout << Verbose(5) << "New triangle " << Nr << ":" << endl;
+  for (int i = 0; i < 3; i++)
+    {
+      lines[i] = line[i];
+      lines[i]->AddTriangle(this);
+    }
+  // get ascending order of endpoints
+  map<int, class BoundaryPointSet *> OrderMap;
+  for (int i = 0; i < 3; i++)
+    // for all three lines
+    for (int j = 0; j < 2; j++)
+      { // for both endpoints
+        OrderMap.insert(pair<int, class BoundaryPointSet *> (
+            line[i]->endpoints[j]->Nr, line[i]->endpoints[j]));
+        // and we don't care whether insertion fails
+      }
+  // set endpoints
+  int Counter = 0;
+  cout << Verbose(6) << " with end points ";
+  for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner
+      != OrderMap.end(); runner++)
+    {
+      endpoints[Counter] = runner->second;
+      cout << " " << *endpoints[Counter];
+      Counter++;
+    }
+  if (Counter < 3)
+    {
+      cerr << "ERROR! We have a triangle with only two distinct endpoints!"
+          << endl;
+      //exit(1);
+    }
+  cout << "." << endl;
+}
+;
+
+BoundaryTriangleSet::~BoundaryTriangleSet()
+{
+  for (int i = 0; i < 3; i++) {
+    cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+    lines[i]->triangles.erase(Nr);
+    if (lines[i]->triangles.empty()) {
+      if (lines[i] != NULL) {
+        cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+        delete (lines[i]);
+        lines[i] = NULL;
+      } else
+        cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+    } else
+      cout << Verbose(5) << *lines[i] << " is still attached to another triangle." << endl;
+  }
+}
+;
+
+/** Calculates the normal vector for this triangle.
+ * Is made unique by comparison with \a OtherVector to point in the other direction.
+ * \param &OtherVector direction vector to make normal vector unique.
+ */
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+{
+  // get normal vector
+  NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
+
+  // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
+  if (NormalVector.Projection(&OtherVector) > 0)
+    NormalVector.Scale(-1.);
+};
+
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
+ * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
+ * smaller than the first line.
+ * \param *out output stream for debugging
+ * \param *MolCenter offset vector of line
+ * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
+ * \param *Intersection intersection on plane on return
+ * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
+ */
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+  Vector CrossPoint;
+  Vector helper;
+  int i=0;
+
+  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, endpoints[0]->node->node, MolCenter, x)) {
+    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
+    return false;
+  }
+
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  do {
+    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection);
+    helper.CopyVector(endpoints[(i+1)%3]->node->node);
+    helper.SubtractVector(endpoints[i%3]->node->node);
+    i++;
+    if (i>3)
+      break;
+  } while (CrossPoint.NormSquared() < MYEPSILON);
+  if (i>3) {
+    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
+    exit(255);
+  }
+  CrossPoint.SubtractVector(endpoints[i%3]->node->node);
+
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+    return true;
+  } else { // outside!
+    Intersection->Zero();
+    return false;
+  }
+};
+
+/** Checks whether lines is any of the three boundary lines this triangle contains.
+ * \param *line line to test
+ * \return true - line is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+{
+  for(int i=0;i<3;i++)
+    if (line == lines[i])
+      return true;
+  return false;
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point point to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+
+/** Checks whether three given \a *Points coincide with triangle's endpoints.
+ * \param *Points[3] pointer to BoundaryPointSet
+ * \return true - is the very triangle, false - is not
+ */
+bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
+{
+  return (((endpoints[0] == Points[0])
+            || (endpoints[0] == Points[1])
+            || (endpoints[0] == Points[2])
+          ) && (
+            (endpoints[1] == Points[0])
+            || (endpoints[1] == Points[1])
+            || (endpoints[1] == Points[2])
+          ) && (
+            (endpoints[2] == Points[0])
+            || (endpoints[2] == Points[1])
+            || (endpoints[2] == Points[2])
+          ));
+};
+
+ostream &
+operator <<(ostream &ost, BoundaryTriangleSet &a)
+{
+  ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << ","
+      << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
+  return ost;
+}
+;
+
+// =========================================================== class TESSELPOINT ===========================================
+
+/** Constructor of class TesselPoint.
+ */
+TesselPoint::TesselPoint()
+{
+  node = NULL;
+  nr = -1;
+  Name =  NULL;
+};
+
+/** Destructor for class TesselPoint.
+ */
+TesselPoint::~TesselPoint()
+{
+  Free((void **)&Name, "TesselPoint::~TesselPoint: *Name");
+};
+
+/** Prints LCNode to screen.
+ */
+ostream & operator << (ostream &ost, const TesselPoint &a)
+{
+  ost << "[" << (a.Name) << "|" << &a << "]";
+  return ost;
+};
+
+
+// =========================================================== class POINTCLOUD ============================================
+
+/** Constructor of class PointCloud.
+ */
+PointCloud::PointCloud()
+{
+
+};
+
+/** Destructor for class PointCloud.
+ */
+PointCloud::~PointCloud()
+{
+
+};
+
+// ============================ CandidateForTesselation =============================
+
+/** Constructor of class CandidateForTesselation.
+ */
+CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) {
+  point = candidate;
+  BaseLine = line;
+  OptCenter.CopyVector(&OptCandidateCenter);
+  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
+};
+
+/** Destructor for class CandidateForTesselation.
+ */
+CandidateForTesselation::~CandidateForTesselation() {
+  point = NULL;
+  BaseLine = NULL;
+};
+
+// =========================================================== class TESSELATION ===========================================
+
+/** Constructor of class Tesselation.
+ */
+Tesselation::Tesselation()
+{
+  PointsOnBoundaryCount = 0;
+  LinesOnBoundaryCount = 0;
+  TrianglesOnBoundaryCount = 0;
+}
+;
+
+/** Destructor of class Tesselation.
+ * We have to free all points, lines and triangles.
+ */
+Tesselation::~Tesselation()
+{
+  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+    if (runner->second != NULL) {
+      delete (runner->second);
+      runner->second = NULL;
+    } else
+      cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+  }
+}
+;
+
+/** Gueses first starting triangle of the convex envelope.
+ * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
+ * \param *out output stream for debugging
+ * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
+ */
+void
+Tesselation::GuessStartingTriangle(ofstream *out)
+{
+  // 4b. create a starting triangle
+  // 4b1. create all distances
+  DistanceMultiMap DistanceMMap;
+  double distance, tmp;
+  Vector PlaneVector, TrialVector;
+  PointMap::iterator A, B, C; // three nodes of the first triangle
+  A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
+
+  // with A chosen, take each pair B,C and sort
+  if (A != PointsOnBoundary.end())
+    {
+      B = A;
+      B++;
+      for (; B != PointsOnBoundary.end(); B++)
+        {
+          C = B;
+          C++;
+          for (; C != PointsOnBoundary.end(); C++)
+            {
+              tmp = A->second->node->node->DistanceSquared(B->second->node->node);
+              distance = tmp * tmp;
+              tmp = A->second->node->node->DistanceSquared(C->second->node->node);
+              distance += tmp * tmp;
+              tmp = B->second->node->node->DistanceSquared(C->second->node->node);
+              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->node,
+          baseline->second.first->second->node->node,
+          baseline->second.second->second->node->node);
+      *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->node);
+          TrialVector.SubtractVector(A->second->node->node);
+          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 << " leaves "
+                  << 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->node->DistanceSquared(A->second->node->node);
+          int innerpoint = 0;
+          if ((tmp < A->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node)) && (tmp
+              < A->second->node->node->DistanceSquared(
+                  baseline->second.second->second->node->node)))
+            innerpoint++;
+          tmp = checker->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node);
+          if ((tmp < baseline->second.first->second->node->node->DistanceSquared(
+              A->second->node->node)) && (tmp
+              < baseline->second.first->second->node->node->DistanceSquared(
+                  baseline->second.second->second->node->node)))
+            innerpoint++;
+          tmp = checker->second->node->node->DistanceSquared(
+              baseline->second.second->second->node->node);
+          if ((tmp < baseline->second.second->second->node->node->DistanceSquared(
+              baseline->second.first->second->node->node)) && (tmp
+              < baseline->second.second->second->node->node->DistanceSquared(
+                  A->second->node->node)))
+            innerpoint++;
+          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
+          if (innerpoint == 3)
+            break;
+        }
+      // 5. come this far, all on same side? Then break 1. loop and construct triangle
+      if (checker == PointsOnBoundary.end())
+        {
+          *out << "Looks like we have a candidate!" << endl;
+          break;
+        }
+    }
+  if (baseline != DistanceMMap.end())
+    {
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = A->second;
+      BPS[1] = baseline->second.second->second;
+      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+      BPS[0] = baseline->second.first->second;
+      BPS[1] = A->second;
+      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+
+      // 4b3. insert created triangle
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+      TrianglesOnBoundaryCount++;
+      for (int i = 0; i < NDIM; i++)
+        {
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
+          LinesOnBoundaryCount++;
+        }
+
+      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
+    }
+  else
+    {
+      *out << Verbose(1) << "No starting triangle found." << endl;
+      exit(255);
+    }
+}
+;
+
+/** Tesselates the convex envelope of a cluster from a single starting triangle.
+ * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
+ * 2 triangles. Hence, we go through all current lines:
+ * -# if the lines contains to only one triangle
+ * -# We search all points in the boundary
+ *    -# if the triangle 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)
+ *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
+ *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
+ * \param *out output stream for debugging
+ * \param *configuration for IsAngstroem
+ * \param *cloud cluster of points
+ */
+void Tesselation::TesselateOnBoundary(ofstream *out, PointCloud *cloud)
+{
+  bool flag;
+  PointMap::iterator winner;
+  class BoundaryPointSet *peak = NULL;
+  double SmallestAngle, TempAngle;
+  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *Center = NULL;
+  LineMap::iterator LineChecker[2];
+
+  Center = cloud->GetCenter(out);
+  // create a first tesselation with the given BoundaryPoints
+  do {
+    flag = false;
+    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
+      if (baseline->second->TrianglesCount == 1) {
+        // 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;
+
+        // get peak point with respect to this base line's only triangle
+        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
+        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;
+
+        // prepare some auxiliary vectors
+        Vector BaseLineCenter, BaseLine;
+        BaseLineCenter.CopyVector(baseline->second->endpoints[0]->node->node);
+        BaseLineCenter.AddVector(baseline->second->endpoints[1]->node->node);
+        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
+        BaseLine.CopyVector(baseline->second->endpoints[0]->node->node);
+        BaseLine.SubtractVector(baseline->second->endpoints[1]->node->node);
+
+        // offset to center of triangle
+        CenterVector.Zero();
+        for (int i = 0; i < 3; i++)
+          CenterVector.AddVector(BTS->endpoints[i]->node->node);
+        CenterVector.Scale(1. / 3.);
+        *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
+
+        // normal vector of triangle
+        NormalVector.CopyVector(Center);
+        NormalVector.SubtractVector(&CenterVector);
+        BTS->GetNormalVector(NormalVector);
+        NormalVector.CopyVector(&BTS->NormalVector);
+        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
+
+        // vector in propagation direction (out of triangle)
+        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
+        TempVector.CopyVector(&CenterVector);
+        TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // 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 << endl;
+        winner = PointsOnBoundary.end();
+
+        // loop over all points and calculate angle between normal vector of new and present triangle
+        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) << ":" << endl;
+
+            // first check direction, so that triangles don't intersect
+            VirtualNormalVector.CopyVector(target->second->node->node);
+            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
+            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
+            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
+            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
+              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
+              continue;
+            } else
+              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
+
+            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
+            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()) && (LineChecker[0]->second->TrianglesCount == 2))) {
+              *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;
+            }
+            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
+              *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)
+            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
+              *out << Verbose(4) << "Current target is peak!" << endl;
+              continue;
+            }
+
+            // check for linear dependence
+            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
+            TempVector.SubtractVector(target->second->node->node);
+            helper.CopyVector(baseline->second->endpoints[1]->node->node);
+            helper.SubtractVector(target->second->node->node);
+            helper.ProjectOntoPlane(&TempVector);
+            if (fabs(helper.NormSquared()) < MYEPSILON) {
+              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
+              continue;
+            }
+
+            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
+            flag = true;
+            VirtualNormalVector.MakeNormalVector(baseline->second->endpoints[0]->node->node, baseline->second->endpoints[1]->node->node, target->second->node->node);
+            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
+            TempVector.AddVector(baseline->second->endpoints[1]->node->node);
+            TempVector.AddVector(target->second->node->node);
+            TempVector.Scale(1./3.);
+            TempVector.SubtractVector(Center);
+            // make it always point outward
+            if (VirtualNormalVector.Projection(&TempVector) < 0)
+              VirtualNormalVector.Scale(-1.);
+            // calculate angle
+            TempAngle = NormalVector.Angle(&VirtualNormalVector);
+            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
+            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
+              SmallestAngle = TempAngle;
+              winner = target;
+              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
+            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
+              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
+              helper.CopyVector(target->second->node->node);
+              helper.SubtractVector(&BaseLineCenter);
+              helper.ProjectOntoPlane(&BaseLine);
+              // ...the one with the smaller angle is the better candidate
+              TempVector.CopyVector(target->second->node->node);
+              TempVector.SubtractVector(&BaseLineCenter);
+              TempVector.ProjectOntoPlane(&VirtualNormalVector);
+              TempAngle = TempVector.Angle(&helper);
+              TempVector.CopyVector(winner->second->node->node);
+              TempVector.SubtractVector(&BaseLineCenter);
+              TempVector.ProjectOntoPlane(&VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(&helper)) {
+                TempAngle = NormalVector.Angle(&VirtualNormalVector);
+                SmallestAngle = TempAngle;
+                winner = target;
+                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
+              } else
+                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
+            } else
+              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
+          }
+        } // end of loop over all boundary points
+
+        // 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);
+
+  // exit
+  delete(Center);
+};
+
+/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
+ * \param *out output stream for debugging
+ * \param *cloud cluster of points
+ * \return true - all straddling points insert, false - something went wrong
+ */
+bool Tesselation::InsertStraddlingPoints(ofstream *out, PointCloud *cloud)
+{
+  Vector Intersection;
+  TesselPoint *Walker = NULL;
+  Vector *Center = cloud->GetCenter(out);
+
+  cloud->GoToFirst();
+  while (!cloud->IsLast()) {  // we only have to go once through all points, as boundary can become only bigger
+    Walker = cloud->GetPoint();
+    // get the next triangle
+    BTS = FindClosestTriangleToPoint(out, Walker->node);
+    if (BTS == NULL) {
+      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
+      return false;
+    }
+    // get the intersection point
+    if (BTS->GetIntersectionInsideTriangle(out, Center, Walker->node, &Intersection)) {
+      // we have the intersection, check whether in- or outside of boundary
+      if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
+        // inside, next!
+        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
+      } else {
+        // outside!
+        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
+        class BoundaryLineSet *OldLines[3], *NewLines[3];
+        class BoundaryPointSet *OldPoints[3], *NewPoint;
+        // store the three old lines and old points
+        for (int i=0;i<3;i++) {
+          OldLines[i] = BTS->lines[i];
+          OldPoints[i] = BTS->endpoints[i];
+        }
+        // add Walker to boundary points
+        AddPoint(Walker);
+        if (BPS[0] == NULL)
+          NewPoint = BPS[0];
+        else
+          continue;
+        // remove triangle
+        TrianglesOnBoundary.erase(BTS->Nr);
+        // create three new boundary lines
+        for (int i=0;i<3;i++) {
+          BPS[0] = NewPoint;
+          BPS[1] = OldPoints[i];
+          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
+          LinesOnBoundaryCount++;
+        }
+        // create three new triangle with new point
+        for (int i=0;i<3;i++) { // find all baselines
+          BLS[0] = OldLines[i];
+          int n = 1;
+          for (int j=0;j<3;j++) {
+            if (NewLines[j]->IsConnectedTo(BLS[0])) {
+              if (n>2) {
+                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
+                return false;
+              } else
+                BLS[n++] = NewLines[j];
+            }
+          }
+          // create the triangle
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+          TrianglesOnBoundaryCount++;
+        }
+      }
+    } else { // something is wrong with FindClosestTriangleToPoint!
+      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
+      return false;
+    }
+    cloud->GoToNext();
+  }
+
+  // exit
+  delete(Center);
+  return true;
+};
+
+/** Adds an atom to the tesselation::PointsOnBoundary list.
+ * \param *Walker atom to add
+ */
+void
+Tesselation::AddPoint(TesselPoint *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++;
+  else {
+    delete(BPS[0]);
+    BPS[0] = NULL;
+  }
+}
+;
+
+/** Adds point to Tesselation::PointsOnBoundary if not yet present.
+ * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
+ * @param Candidate point to add
+ * @param n index for this point in Tesselation::TPS array
+ */
+void
+Tesselation::AddTrianglePoint(TesselPoint* Candidate, int n)
+{
+  PointTestPair InsertUnique;
+  TPS[n] = new class BoundaryPointSet(Candidate);
+  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
+  if (InsertUnique.second) { // if new point was not present before, increase counter
+    PointsOnBoundaryCount++;
+  } else {
+    delete TPS[n];
+    cout << Verbose(3) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
+    TPS[n] = (InsertUnique.first)->second;
+  }
+}
+;
+
+/** Function tries to add line from current Points in BPS to BoundaryLineSet.
+ * If successful it raises the line count and inserts the new line into the BLS,
+ * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
+ * @param *a first endpoint
+ * @param *b second endpoint
+ * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
+  bool insertNewLine = true;
+
+  if (a->lines.find(b->node->nr) != a->lines.end()) {
+    LineMap::iterator FindLine;
+    pair<LineMap::iterator,LineMap::iterator> FindPair;
+    FindPair = a->lines.equal_range(b->node->nr);
+
+    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+      // If there is a line with less than two attached triangles, we don't need a new line.
+      if (FindLine->second->TrianglesCount < 2) {
+        insertNewLine = false;
+        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;
+
+        BPS[0] = FindLine->second->endpoints[0];
+        BPS[1] = FindLine->second->endpoints[1];
+        BLS[n] = FindLine->second;
+
+        break;
+      }
+    }
+  }
+
+  if (insertNewLine) {
+    AlwaysAddTriangleLine(a, b, n);
+  }
+}
+;
+
+/**
+ * Adds lines from each of the current points in the BPS to BoundaryLineSet.
+ * Raises the line count and inserts the new line into the BLS.
+ *
+ * @param *a first endpoint
+ * @param *b second endpoint
+ * @param n index of Tesselation::BLS giving the line with both endpoints
+ */
+void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
+{
+  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
+  BPS[0] = a;
+  BPS[1] = b;
+  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
+  // add line to global map
+  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
+  // increase counter
+  LinesOnBoundaryCount++;
+};
+
+/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
+ * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
+ */
+void
+Tesselation::AddTriangle()
+{
+  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;
+
+  // add triangle to global map
+  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+  TrianglesOnBoundaryCount++;
+
+  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
+}
+;
+
+
+
+/** Checks whether the triangle consisting of the three atoms is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *Candidates endpoints of the triangle candidate
+ * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
+ *                 triangles exist which is the maximum for three points
+ */
+int Tesselation::CheckPresenceOfTriangle(ofstream *out, TesselPoint *Candidates[3]) {
+  int adjacentTriangleCount = 0;
+  class BoundaryPointSet *Points[3];
+
+  *out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+  // builds a triangle point set (Points) of the end points
+  for (int i = 0; i < 3; i++) {
+    PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
+    if (FindPoint != PointsOnBoundary.end()) {
+      Points[i] = FindPoint->second;
+    } else {
+      Points[i] = NULL;
+    }
+  }
+
+  // checks lines between the points in the Points for their adjacent triangles
+  for (int i = 0; i < 3; i++) {
+    if (Points[i] != NULL) {
+      for (int j = i; j < 3; j++) {
+        if (Points[j] != NULL) {
+          LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
+          for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
+            TriangleMap *triangles = &FindLine->second->triangles;
+            *out << Verbose(3) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;
+            for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
+              if (FindTriangle->second->IsPresentTupel(Points)) {
+                adjacentTriangleCount++;
+              }
+            }
+            *out << Verbose(3) << "end." << endl;
+          }
+          // Only one of the triangle lines must be considered for the triangle count.
+          *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+          return adjacentTriangleCount;
+        }
+      }
+    }
+  }
+
+  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
+  *out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+  return adjacentTriangleCount;
+};
+
+
+/** 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 *out output stream for debugging
+ * \param RADIUS radius of virtual rolling sphere
+ * \param *LC LinkedCell structure with neighbouring TesselPoint's
+ */
+void Tesselation::Find_starting_triangle(ofstream *out, const double RADIUS, LinkedCell *LC)
+{
+  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+  int i = 0;
+  LinkedNodes *List = NULL;
+  TesselPoint* FirstPoint = NULL;
+  TesselPoint* SecondPoint = NULL;
+  TesselPoint* MaxAtom[NDIM];
+  double max_coordinate[NDIM];
+  Vector Oben;
+  Vector helper;
+  Vector Chord;
+  Vector SearchDirection;
+
+  Oben.Zero();
+
+  for (i = 0; i < 3; i++) {
+    MaxAtom[i] = NULL;
+    max_coordinate[i] = -1;
+  }
+
+  // 1. searching topmost atom with respect to each axis
+  for (int i=0;i<NDIM;i++) { // each axis
+    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin();Runner != List->end();Runner++) {
+            if ((*Runner)->node->x[i] > max_coordinate[i]) {
+              cout << Verbose(2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
+              max_coordinate[i] = (*Runner)->node->x[i];
+              MaxAtom[i] = (*Runner);
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+  }
+
+  cout << Verbose(2) << "Found maximum coordinates: ";
+  for (int i=0;i<NDIM;i++)
+    cout << i << ": " << *MaxAtom[i] << "\t";
+  cout << endl;
+
+  BTS = NULL;
+  CandidateList *Opt_Candidates = new CandidateList();
+  for (int k=0;k<NDIM;k++) {
+    Oben.x[k] = 1.;
+    FirstPoint = MaxAtom[k];
+    cout << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;
+
+    double ShortestAngle;
+    TesselPoint* Opt_Candidate = NULL;
+    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+
+    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    SecondPoint = Opt_Candidate;
+    if (SecondPoint == NULL)  // have we found a second point?
+      continue;
+    else
+      cout << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n";
+
+    helper.CopyVector(FirstPoint->node);
+    helper.SubtractVector(SecondPoint->node);
+    helper.Normalize();
+    Oben.ProjectOntoPlane(&helper);
+    Oben.Normalize();
+    helper.VectorProduct(&Oben);
+    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+
+    Chord.CopyVector(FirstPoint->node); // bring into calling function
+    Chord.SubtractVector(SecondPoint->node);
+    double radius = Chord.ScalarProduct(&Chord);
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+    helper.CopyVector(&Oben);
+    helper.Scale(CircleRadius);
+    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+
+    // look in one direction of baseline for initial candidate
+    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+
+    // adding point 1 and point 2 and the line between them
+    AddTrianglePoint(FirstPoint, 0);
+    AddTrianglePoint(SecondPoint, 1);
+    AddTriangleLine(TPS[0], TPS[1], 0);
+
+    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
+    Find_third_point_for_Tesselation(
+      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
+    );
+    cout << Verbose(1) << "List of third Points is ";
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        cout << " " << *(*it)->point;
+    }
+    cout << endl;
+
+    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+      // add third triangle point
+      AddTrianglePoint((*it)->point, 2);
+      // add the second and third line
+      AddTriangleLine(TPS[1], TPS[2], 1);
+      AddTriangleLine(TPS[0], TPS[2], 2);
+      // ... and triangles to the Maps of the Tesselation class
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      AddTriangle();
+      // ... and calculate its normal vector (with correct orientation)
+      (*it)->OptCenter.Scale(-1.);
+      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
+      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
+      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
+      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";
+
+      // if we do not reach the end with the next step of iteration, we need to setup a new first line
+      if (it != Opt_Candidates->end()--) {
+        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
+        SecondPoint = (*it)->point;
+        // adding point 1 and point 2 and the line between them
+        AddTrianglePoint(FirstPoint, 0);
+        AddTrianglePoint(SecondPoint, 1);
+        AddTriangleLine(TPS[0], TPS[1], 0);
+      }
+      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+    }
+    if (BTS != NULL) // we have created one starting triangle
+      break;
+    else {
+      // remove all candidates from the list and then the list itself
+      class CandidateForTesselation *remover = NULL;
+      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+        remover = *it;
+        delete(remover);
+      }
+      Opt_Candidates->clear();
+    }
+  }
+
+  // remove all candidates from the list and then the list itself
+  class CandidateForTesselation *remover = NULL;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    remover = *it;
+    delete(remover);
+  }
+  delete(Opt_Candidates);
+  cout << Verbose(1) << "End of Find_starting_triangle\n";
+};
+
+
+/** This function finds a triangle to a line, adjacent to an existing one.
+ * @param out output stream for debugging
+ * @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 *LC LinkedCell structure with neighbouring atoms
+ */
+bool Tesselation::Find_next_suitable_triangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC)
+{
+  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
+  bool result = true;
+  CandidateList *Opt_Candidates = new CandidateList();
+
+  Vector CircleCenter;
+  Vector CirclePlaneNormal;
+  Vector OldSphereCenter;
+  Vector SearchDirection;
+  Vector helper;
+  TesselPoint *ThirdNode = NULL;
+  LineMap::iterator testline;
+  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+  double radius, CircleRadius;
+
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+  for (int i=0;i<3;i++)
+    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+      ThirdNode = T.endpoints[i]->node;
+
+  // construct center of circle
+  CircleCenter.CopyVector(Line.endpoints[0]->node->node);
+  CircleCenter.AddVector(Line.endpoints[1]->node->node);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(Line.endpoints[0]->node->node);
+  CirclePlaneNormal.SubtractVector(Line.endpoints[1]->node->node);
+
+  // 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->node, T.endpoints[1]->node->node, T.endpoints[2]->node->node);
+    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+    radius = Line.endpoints[0]->node->node->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->node);
+    helper.SubtractVector(ThirdNode->node);
+    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+      SearchDirection.Scale(-1.);
+    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+    SearchDirection.Normalize();
+    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+      // rotated the wrong way!
+      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+    }
+
+    // add third point
+    Find_third_point_for_Tesselation(
+      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
+      &ShortestAngle, RADIUS, LC
+    );
+
+  } else {
+    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+  }
+
+  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << "Third Points are ";
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    cout << " " << *(*it)->point;
+  }
+  cout << endl;
+
+  BoundaryLineSet *BaseRay = &Line;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
+    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
+    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;
+
+    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+    TesselPoint *AtomCandidates[3];
+    AtomCandidates[0] = (*it)->point;
+    AtomCandidates[1] = BaseRay->endpoints[0]->node;
+    AtomCandidates[2] = BaseRay->endpoints[1]->node;
+    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);
+
+    BTS = NULL;
+    // If there is no triangle, add it regularly.
+    if (existentTrianglesCount == 0) {
+      AddTrianglePoint((*it)->point, 0);
+      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
+      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+
+      AddTriangleLine(TPS[0], TPS[1], 0);
+      AddTriangleLine(TPS[0], TPS[2], 1);
+      AddTriangleLine(TPS[1], TPS[2], 2);
+
+      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+      AddTriangle();
+      (*it)->OptCenter.Scale(-1.);
+      BTS->GetNormalVector((*it)->OptCenter);
+      (*it)->OptCenter.Scale(-1.);
+
+      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
+        << " for this triangle ... " << endl;
+      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
+    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
+        AddTrianglePoint((*it)->point, 0);
+        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
+        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);
+
+        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
+        // i.e. at least one of the three lines must be present with TriangleCount <= 1
+        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
+          AddTriangleLine(TPS[0], TPS[1], 0);
+          AddTriangleLine(TPS[0], TPS[2], 1);
+          AddTriangleLine(TPS[1], TPS[2], 2);
+
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          AddTriangle();
+
+          (*it)->OtherOptCenter.Scale(-1.);
+          BTS->GetNormalVector((*it)->OtherOptCenter);
+          (*it)->OtherOptCenter.Scale(-1.);
+
+          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
+          << " for this triangle ... " << endl;
+          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
+        } else {
+          cout << Verbose(1) << "WARNING: This triangle consisting of ";
+          cout << *(*it)->point << ", ";
+          cout << *BaseRay->endpoints[0]->node << " and ";
+          cout << *BaseRay->endpoints[1]->node << " ";
+          cout << "exists and is not added, as it does not seem helpful!" << endl;
+          result = false;
+        }
+    } else {
+      cout << Verbose(1) << "This triangle consisting of ";
+      cout << *(*it)->point << ", ";
+      cout << *BaseRay->endpoints[0]->node << " and ";
+      cout << *BaseRay->endpoints[1]->node << " ";
+      cout << "is invalid!" << endl;
+      result = false;
+    }
+
+    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
+    BaseRay = BLS[0];
+  }
+
+  // remove all candidates from the list and then the list itself
+  class CandidateForTesselation *remover = NULL;
+  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
+    remover = *it;
+    delete(remover);
+  }
+  delete(Opt_Candidates);
+  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
+  return result;
+};
+
+
+/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
+ * \param *out output stream for debugging
+ * \return true - all baselines were corrected, false - there are still concave pieces
+ */
+bool Tesselation::CorrectConcaveBaselines(ofstream *out)
+{
+  class BoundaryTriangleSet *triangle[2];
+  class BoundaryLineSet *OldLines[4], *NewLine;
+  class BoundaryPointSet *OldPoints[2];
+  Vector BaseLineNormal;
+  class BoundaryLineSet *Base = NULL;
+  int OldTriangles[2], OldBaseLine;
+  int i;
+  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
+    Base = baseline->second;
+
+    // check convexity
+    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
+      *out << Verbose(3) << Base << " has two convex triangles." << endl;
+    } else { // not convex!
+      // get the two triangles
+      i=0;
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+        triangle[i++] = runner->second;
+      // gather four endpoints and four lines
+      for (int j=0;j<4;j++)
+        OldLines[j] = NULL;
+      for (int j=0;j<2;j++)
+        OldPoints[j] = NULL;
+      i=0;
+      for (int m=0;m<2;m++) { // go over both triangles
+        for (int j=0;j<3;j++) { // all of their endpoints and baselines
+          if (triangle[m]->lines[j] != Base) // pick not the central baseline
+            OldLines[i++] = triangle[m]->lines[j];
+          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
+            OldPoints[m] = triangle[m]->endpoints[j];
+        }
+      }
+      if (i<4) {
+        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+        return false;
+      }
+      for (int j=0;j<4;j++)
+        if (OldLines[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+          return false;
+        }
+      for (int j=0;j<2;j++)
+        if (OldPoints[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
+          return false;
+        }
+
+      // remove triangles
+      for (int j=0;j<2;j++) {
+        OldTriangles[j] = triangle[j]->Nr;
+        TrianglesOnBoundary.erase(OldTriangles[j]);
+        triangle[j] = NULL;
+      }
+
+      // remove baseline
+      OldBaseLine = Base->Nr;
+      LinesOnBoundary.erase(OldBaseLine);
+      Base = NULL;
+
+      // construct new baseline (with same number as old one)
+      BPS[0] = OldPoints[0];
+      BPS[1] = OldPoints[1];
+      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
+      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
+
+      // construct new triangles with flipped baseline
+      i=-1;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=2;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=3;
+      if (i!=-1) {
+        BLS[0] = OldLines[0];
+        BLS[1] = OldLines[i];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
+
+        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
+        BLS[1] = OldLines[1];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
+      } else {
+        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
+        return false;
+      }
+    }
+  }
+  return true;
+};
+
+
+/** States whether point is in- or outside of a tesselated surface.
+ * \param *pointer point to be checked
+ * \return true - is inside, false - is outside
+ */
+bool Tesselation::IsInside(Vector *pointer)
+{
+
+  // hier kommt dann Saskias Routine hin...
+
+  return true;
+};
+
+
+/** Finds the closest triangle to a given point.
+ * \param *out output stream for debugging
+ * \param *x second endpoint of line
+ * \return pointer triangle that is closest, NULL if none was found
+ */
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
+{
+  class BoundaryTriangleSet *triangle = NULL;
+
+  // hier kommt dann Saskias Routine hin...
+
+  return triangle;
+};
+
+
+/** Finds the second point of starting triangle.
+ * \param *a first node
+ * \param *Candidate pointer to candidate node 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 Tesselation::Find_second_point_for_Tesselation(TesselPoint* a, TesselPoint* Candidate, Vector Oben, TesselPoint*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
+  Vector AngleCheck;
+  double norm = -1., angle;
+  LinkedNodes *List = NULL;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LC->SetIndexToNode(a)) {  // get cell for the starting atom
+    for(int i=0;i<NDIM;i++) // store indices of this cell
+      N[i] = LC->n[i];
+  } else {
+    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+    return;
+  }
+  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+  cout << Verbose(3) << "LC Intervals from [";
+  for (int i=0;i<NDIM;i++) {
+  cout << " " << N[i] << "<->" << LC->N[i];
+  }
+  cout << "] :";
+  for (int i=0;i<NDIM;i++) {
+    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+  }
+  cout << endl;
+
+
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            Candidate = (*Runner);
+            // check if we only have one unique point yet ...
+            if (a != Candidate) {
+              // Calculate center of the circle with radius RADIUS through points a and Candidate
+              Vector OrthogonalizedOben, a_Candidate, Center;
+              double distance, scaleFactor;
+
+              OrthogonalizedOben.CopyVector(&Oben);
+              a_Candidate.CopyVector(a->node);
+              a_Candidate.SubtractVector(Candidate->node);
+              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
+              OrthogonalizedOben.Normalize();
+              distance = 0.5 * a_Candidate.Norm();
+              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
+              OrthogonalizedOben.Scale(scaleFactor);
+
+              Center.CopyVector(Candidate->node);
+              Center.AddVector(a->node);
+              Center.Scale(0.5);
+              Center.AddVector(&OrthogonalizedOben);
+
+              AngleCheck.CopyVector(&Center);
+              AngleCheck.SubtractVector(a->node);
+              norm = a_Candidate.Norm();
+              // second point shall have smallest angle with respect to Oben vector
+              if (norm < RADIUS*2.) {
+                angle = AngleCheck.Angle(&Oben);
+                if (angle < Storage[0]) {
+                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
+                  Opt_Candidate = Candidate;
+                  Storage[0] = angle;
+                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+                } else {
+                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+                }
+              } else {
+                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
+              }
+            } else {
+              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
+            }
+          }
+        } else {
+          cout << Verbose(3) << "Linked cell list is empty." << endl;
+        }
+      }
+  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
+};
+
+
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * Note that this function is for the starting triangle.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
+ * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
+ * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
+ * @param BaseLine BoundaryLineSet with the current base line
+ * @param ThirdNode third atom to avoid in search
+ * @param candidates list of equally good candidates to return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+void Tesselation::Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint  *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+{
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector SphereCenter;
+  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  Vector NewNormalVector;   // normal vector of the Candidate's triangle
+  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  LinkedNodes *List = NULL;
+  double CircleRadius; // radius of this circle
+  double radius;
+  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  TesselPoint *Candidate = NULL;
+  CandidateForTesselation *optCandidate = NULL;
+
+  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
+
+  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+
+  // construct center of circle
+  CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node);
+  CircleCenter.AddVector(BaseLine->endpoints[1]->node->node);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node);
+  CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node);
+
+  // calculate squared radius atom *ThirdNode,f circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    // test whether old center is on the band's plane
+    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(radius - CircleRadius) < HULLEPSILON) {
+
+      // check SearchDirection
+      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
+      }
+
+      // get cell for the starting atom
+      if (LC->SetIndexToVector(&CircleCenter)) {
+        for(int i=0;i<NDIM;i++) // store indices of this cell
+        N[i] = LC->n[i];
+        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+      } else {
+        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+        return;
+      }
+      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+      //cout << Verbose(2) << "LC Intervals:";
+      for (int i=0;i<NDIM;i++) {
+        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+      }
+      //cout << endl;
+      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+            List = LC->GetCurrentCell();
+            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+            if (List != NULL) {
+              for (LinkedNodes::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                Candidate = (*Runner);
+
+                // check for three unique points
+                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){
+
+                  // construct both new centers
+                  GetCenterofCircumcircle(&NewSphereCenter, BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node);
+                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+
+                  if ((NewNormalVector.MakeNormalVector(BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node))
+                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  ) {
+                    helper.CopyVector(&NewNormalVector);
+                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    if (radius < RADIUS*RADIUS) {
+                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      NewSphereCenter.AddVector(&helper);
+                      NewSphereCenter.SubtractVector(&CircleCenter);
+                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+
+                      // OtherNewSphereCenter is created by the same vector just in the other direction
+                      helper.Scale(-1.);
+                      OtherNewSphereCenter.AddVector(&helper);
+                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
+                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+
+                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      alpha = min(alpha, Otheralpha);
+                      // if there is a better candidate, drop the current list and add the new candidate
+                      // otherwise ignore the new candidate and keep the list
+                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
+                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
+                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
+                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
+                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                        } else {
+                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
+                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
+                        }
+                        // if there is an equal candidate, add it to the list without clearing the list
+                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
+                          candidates->push_back(optCandidate);
+                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
+                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        } else {
+                          // remove all candidates from the list and then the list itself
+                          class CandidateForTesselation *remover = NULL;
+                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
+                            remover = *it;
+                            delete(remover);
+                          }
+                          candidates->clear();
+                          candidates->push_back(optCandidate);
+                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
+                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
+                        }
+                        *ShortestAngle = alpha;
+                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
+                      } else {
+                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
+                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+                        } else {
+                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                        }
+                      }
+
+                    } else {
+                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                    }
+                  } else {
+                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                  }
+                } else {
+                  if (ThirdNode != NULL) {
+                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                  } else {
+                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                  }
+                }
+              }
+            }
+          }
+    } else {
+      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+    }
+  } else {
+    if (ThirdNode != NULL)
+      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+    else
+      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+  }
+
+  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
+  if (candidates->size() > 1) {
+    candidates->unique();
+    candidates->sort(sortCandidates);
+  }
+
+  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
+};
+
+/** Finds the endpoint two lines are sharing.
+ * \param *line1 first line
+ * \param *line2 second line
+ * \return point which is shared or NULL if none
+ */
+class BoundaryPointSet *Tesselation::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;
+};
+
+/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
+ * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
+ * make it bigger (i.e. closing one (the baseline) and opening two new ones).
+ * \param TPS[3] nodes of the triangle
+ * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
+ */
+bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
+{
+  bool result = false;
+  int counter = 0;
+
+  // check all three points
+  for (int i=0;i<3;i++)
+    for (int j=i+1; j<3; j++) {
+      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
+        LineMap::iterator FindLine;
+        pair<LineMap::iterator,LineMap::iterator> FindPair;
+        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
+        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
+          // If there is a line with less than two attached triangles, we don't need a new line.
+          if (FindLine->second->TrianglesCount < 2) {
+            counter++;
+            break;  // increase counter only once per edge
+          }
+        }
+      } else { // no line
+        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
+        result = true;
+      }
+    }
+  if (counter > 1) {
+    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
+    result = true;
+  }
+  return result;
+};
+
+
+/** Sort function for the candidate list.
+ */
+bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
+  Vector BaseLineVector, OrthogonalVector, helper;
+  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
+    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
+    //return false;
+    exit(1);
+  }
+  // create baseline vector
+  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
+  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  BaseLineVector.Normalize();
+
+  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
+  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
+  helper.SubtractVector(candidate1->point->node);
+  OrthogonalVector.CopyVector(&helper);
+  helper.VectorProduct(&BaseLineVector);
+  OrthogonalVector.SubtractVector(&helper);
+  OrthogonalVector.Normalize();
+
+  // calculate both angles and correct with in-plane vector
+  helper.CopyVector(candidate1->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double phi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    phi = 2.*M_PI - phi;
+  }
+  helper.CopyVector(candidate2->point->node);
+  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
+  double psi = BaseLineVector.Angle(&helper);
+  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
+    psi = 2.*M_PI - psi;
+  }
+
+  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
+  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
+
+  // return comparison
+  return phi < psi;
+};
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/tesselation.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,228 @@
+/*
+ * tesselation.hpp
+ *
+ *  The tesselation class is meant to contain the envelope (concave, convex or neither) of a set of Vectors. As we actually mean this stuff for atoms, we have to encapsulate it all a bit.
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#ifndef TESSELATION_HPP_
+#define TESSELATION_HPP_
+
+using namespace std;
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <map>
+#include <list>
+
+
+#include "helpers.hpp"
+#include "linkedcell.hpp"
+#include "tesselationhelpers.hpp"
+#include "vector.hpp"
+
+class BoundaryPointSet;
+class BoundaryLineSet;
+class BoundaryTriangleSet;
+class TesselPoint;
+class PointCloud;
+class Tesselation;
+
+// ======================================================= some template functions =========================================
+
+#define PointMap map < int, class BoundaryPointSet * >
+#define PointPair pair < int, class BoundaryPointSet * >
+#define PointTestPair pair < PointMap::iterator, bool >
+#define CandidateList list <class CandidateForTesselation *>
+
+#define LineMap multimap < int, class BoundaryLineSet * >
+#define LinePair pair < int, class BoundaryLineSet * >
+#define LineTestPair pair < LineMap::iterator, bool >
+
+#define TriangleMap map < int, class BoundaryTriangleSet * >
+#define TrianglePair pair < int, class BoundaryTriangleSet * >
+#define TriangleTestPair pair < TrianglePair::iterator, bool >
+
+#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
+#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
+
+
+
+template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
+{
+  if (endpoint1->Nr < endpoint2->Nr) {
+    endpoints[0] = endpoint1;
+    endpoints[1] = endpoint2;
+  } else {
+    endpoints[0] = endpoint2;
+    endpoints[1] = endpoint1;
+  }
+};
+
+// ======================================================== class BoundaryPointSet =========================================
+
+class BoundaryPointSet {
+  public:
+    BoundaryPointSet();
+    BoundaryPointSet(TesselPoint *Walker);
+    ~BoundaryPointSet();
+
+    void AddLine(class BoundaryLineSet *line);
+
+    LineMap lines;
+    int LinesCount;
+    TesselPoint *node;
+    int Nr;
+};
+
+ostream & operator << (ostream &ost, BoundaryPointSet &a);
+
+// ======================================================== class BoundaryLineSet ==========================================
+
+class BoundaryLineSet {
+  public:
+    BoundaryLineSet();
+    BoundaryLineSet(class BoundaryPointSet *Point[2], int number);
+    ~BoundaryLineSet();
+
+    void AddTriangle(class BoundaryTriangleSet *triangle);
+    bool IsConnectedTo(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool CheckConvexityCriterion(ofstream *out);
+
+    class BoundaryPointSet *endpoints[2];
+    TriangleMap triangles;
+    int TrianglesCount;
+    int Nr;
+};
+
+ostream & operator << (ostream &ost, BoundaryLineSet &a);
+
+// ======================================================== class BoundaryTriangleSet =======================================
+
+class BoundaryTriangleSet {
+  public:
+    BoundaryTriangleSet();
+    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
+    ~BoundaryTriangleSet();
+
+    void GetNormalVector(Vector &NormalVector);
+    bool GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection);
+    bool ContainsBoundaryLine(class BoundaryLineSet *line);
+    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
+    bool IsPresentTupel(class BoundaryPointSet *Points[3]);
+
+    class BoundaryPointSet *endpoints[3];
+    class BoundaryLineSet *lines[3];
+    Vector NormalVector;
+    int Nr;
+};
+
+ostream & operator << (ostream &ost, BoundaryTriangleSet &a);
+
+// =========================================================== class TESSELPOINT ===========================================
+
+/** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
+ */
+class TesselPoint {
+public:
+  TesselPoint();
+  ~TesselPoint();
+
+  Vector *node;   // pointer to position of the dot in space
+  int nr;       // index to easierly identify
+  char *Name;   // some name to reference to on output
+};
+
+ostream & operator << (ostream &ost, const TesselPoint &a);
+
+// =========================================================== class POINTCLOUD ============================================
+
+/** Super-class for all point clouds structures, also molecules. They have to inherit this structure and implement the virtual function to access the Vectors.
+ * This basically encapsulates a list structure.
+ */
+class PointCloud {
+public:
+  PointCloud();
+  ~PointCloud();
+
+  virtual Vector *GetCenter(ofstream *out) { return NULL; };
+  virtual TesselPoint *GetPoint() { return NULL; };
+  virtual TesselPoint *GetTerminalPoint() { return NULL; };
+  virtual void GoToNext() {};
+  virtual void GoToPrevious() {};
+  virtual void GoToFirst() {};
+  virtual void GoToLast() {};
+  virtual bool IsEmpty() { return false; };
+  virtual bool IsLast() { return false; };
+};
+
+// ======================================================== class CandidateForTesselation =========================================
+
+class CandidateForTesselation {
+  public :
+  CandidateForTesselation(TesselPoint* candidate, BoundaryLineSet* currentBaseLine, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
+  ~CandidateForTesselation();
+
+  TesselPoint *point;
+  BoundaryLineSet *BaseLine;
+  Vector OptCenter;
+  Vector OtherOptCenter;
+};
+
+// =========================================================== class TESSELATION ===========================================
+
+/** Contains the envelope to a PointCloud.
+ */
+class Tesselation {
+  public:
+
+    Tesselation();
+    ~Tesselation();
+
+    void AddPoint(TesselPoint *Walker);
+    void AddTrianglePoint(TesselPoint* Candidate, int n);
+    void AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
+    void AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n);
+    void AddTriangle();
+    bool IsInside(Vector *pointer);
+    class BoundaryPointSet *GetCommonEndpoint(class BoundaryLineSet * line1, class BoundaryLineSet * line2);
+
+    // concave envelope
+    void Find_starting_triangle(ofstream *out, const double RADIUS, class LinkedCell *LC);
+    void Find_second_point_for_Tesselation(class TesselPoint* a, class TesselPoint* Candidate, Vector Oben, class TesselPoint*& Opt_Candidate, double Storage[3], double RADIUS, class LinkedCell *LC);
+    void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, class TesselPoint *ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, class LinkedCell *LC);
+    bool Find_next_suitable_triangle(ofstream *out, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, LinkedCell *LC);
+    int CheckPresenceOfTriangle(ofstream *out, class TesselPoint *Candidates[3]);
+
+    // convex envelope
+    class BoundaryTriangleSet * FindClosestTriangleToPoint(ofstream *out, Vector *x);
+    void TesselateOnBoundary(ofstream *out, PointCloud *cloud);
+    void GuessStartingTriangle(ofstream *out);
+    bool InsertStraddlingPoints(ofstream *out, PointCloud *cloud);
+    bool CorrectConcaveBaselines(ofstream *out);
+
+    PointMap PointsOnBoundary;
+    LineMap LinesOnBoundary;
+    TriangleMap TrianglesOnBoundary;
+    int PointsOnBoundaryCount;
+    int LinesOnBoundaryCount;
+    int TrianglesOnBoundaryCount;
+
+  private:
+    class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions
+    class BoundaryPointSet *BPS[2];
+    class BoundaryLineSet *BLS[3];
+    class BoundaryTriangleSet *BTS;
+};
+
+bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3]);
+bool sortCandidates(class CandidateForTesselation* candidate1, class CandidateForTesselation* candidate2);
+
+
+#endif /* TESSELATION_HPP_ */
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/tesselationhelpers.cpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,389 @@
+/*
+ * TesselationHelpers.cpp
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#include "tesselationhelpers.hpp"
+
+double det_get(gsl_matrix *A, int inPlace) {
+  /*
+  inPlace = 1 => A is replaced with the LU decomposed copy.
+  inPlace = 0 => A is retained, and a copy is used for LU.
+  */
+
+  double det;
+  int signum;
+  gsl_permutation *p = gsl_permutation_alloc(A->size1);
+  gsl_matrix *tmpA;
+
+  if (inPlace)
+  tmpA = A;
+  else {
+  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+  gsl_matrix_memcpy(tmpA , A);
+  }
+
+
+  gsl_linalg_LU_decomp(tmpA , p , &signum);
+  det = gsl_linalg_LU_det(tmpA , signum);
+  gsl_permutation_free(p);
+  if (! inPlace)
+  gsl_matrix_free(tmpA);
+
+  return det;
+};
+
+void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
+  gsl_matrix *A = gsl_matrix_calloc(3,3);
+  double m11, m12, m13, m14;
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m11 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m12 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m13 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, b.x[i]);
+  }
+  m14 = det_get(A, 1);
+
+  if (fabs(m11) < MYEPSILON)
+    cerr << "ERROR: three points are colinear." << endl;
+
+  center->x[0] =  0.5 * m12/ m11;
+  center->x[1] = -0.5 * m13/ m11;
+  center->x[2] =  0.5 * m14/ m11;
+
+  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
+    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+
+  gsl_matrix_free(A);
+};
+
+
+
+/**
+ * Function returns center of sphere with RADIUS, which rests on points a, b, c
+ * @param Center this vector will be used for return
+ * @param a vector first point of triangle
+ * @param b vector second point of triangle
+ * @param c vector third point of triangle
+ * @param *Umkreismittelpunkt new cneter point of circumference
+ * @param Direction vector indicates up/down
+ * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
+ * @param Halfplaneindicator double indicates whether Direction is up or down
+ * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
+ * @param alpha double angle at a
+ * @param beta double, angle at b
+ * @param gamma, double, angle at c
+ * @param Radius, double
+ * @param Umkreisradius double radius of circumscribing circle
+ */
+void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
+    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
+  Vector TempNormal, helper;
+  double Restradius;
+  Vector OtherCenter;
+  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+  Center->Zero();
+  helper.CopyVector(&a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(&b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(&c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+  NewUmkreismittelpunkt->CopyVector(Center);
+  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
+  // Here we calculated center of circumscribing circle, using barycentric coordinates
+  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
+
+  TempNormal.CopyVector(&a);
+  TempNormal.SubtractVector(&b);
+  helper.CopyVector(&a);
+  helper.SubtractVector(&c);
+  TempNormal.VectorProduct(&helper);
+  if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  else
+    {
+      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+
+  TempNormal.Normalize();
+  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
+  TempNormal.Scale(Restradius);
+  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
+
+  Center->AddVector(&TempNormal);
+  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
+  get_sphere(&OtherCenter, a, b, c, RADIUS);
+  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
+  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
+};
+
+
+/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
+ * \param *Center new center on return
+ * \param *a first point
+ * \param *b second point
+ * \param *c third point
+ */
+void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
+  Vector helper;
+  double alpha, beta, gamma;
+  Vector SideA, SideB, SideC;
+  SideA.CopyVector(b);
+  SideA.SubtractVector(c);
+  SideB.CopyVector(c);
+  SideB.SubtractVector(a);
+  SideC.CopyVector(a);
+  SideC.SubtractVector(b);
+  alpha = M_PI - SideB.Angle(&SideC);
+  beta = M_PI - SideC.Angle(&SideA);
+  gamma = M_PI - SideA.Angle(&SideB);
+  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
+  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
+    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+
+  Center->Zero();
+  helper.CopyVector(a);
+  helper.Scale(sin(2.*alpha));
+  Center->AddVector(&helper);
+  helper.CopyVector(b);
+  helper.Scale(sin(2.*beta));
+  Center->AddVector(&helper);
+  helper.CopyVector(c);
+  helper.Scale(sin(2.*gamma));
+  Center->AddVector(&helper);
+  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+};
+
+/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
+ * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
+ * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
+ * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
+ * \param CircleCenter Center of the parameter circle
+ * \param CirclePlaneNormal normal vector to plane of the parameter circle
+ * \param CircleRadius radius of the parameter circle
+ * \param NewSphereCenter new center of a circumcircle
+ * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
+ * \param NormalVector normal vector
+ * \param SearchDirection search direction to make angle unique on return.
+ * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
+ */
+double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
+  Vector helper;
+  double radius, alpha;
+
+  helper.CopyVector(&NewSphereCenter);
+  // test whether new center is on the parameter circle's plane
+  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
+  radius = helper.ScalarProduct(&helper);
+  // test whether the new center vector has length of CircleRadius
+  if (fabs(radius - CircleRadius) > HULLEPSILON)
+    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+  alpha = helper.Angle(&OldSphereCenter);
+  // make the angle unique by checking the halfplanes/search direction
+  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+    alpha = 2.*M_PI - alpha;
+  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&OldSphereCenter);
+  helper.ProjectOntoPlane(&NormalVector);
+  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
+  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
+    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    return alpha;
+  } else {
+    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    return 2.*M_PI;
+  }
+};
+
+struct Intersection {
+  Vector x1;
+  Vector x2;
+  Vector x3;
+  Vector x4;
+};
+
+/**
+ * Intersection calculation function.
+ *
+ * @param x to find the result for
+ * @param function parameter
+ */
+double MinIntersectDistance(const gsl_vector * x, void *params)
+{
+  double retval = 0;
+  struct Intersection *I = (struct Intersection *)params;
+  Vector intersection;
+  Vector SideA,SideB,HeightA, HeightB;
+  for (int i=0;i<NDIM;i++)
+    intersection.x[i] = gsl_vector_get(x, i);
+
+  SideA.CopyVector(&(I->x1));
+  SideA.SubtractVector(&I->x2);
+  HeightA.CopyVector(&intersection);
+  HeightA.SubtractVector(&I->x1);
+  HeightA.ProjectOntoPlane(&SideA);
+
+  SideB.CopyVector(&I->x3);
+  SideB.SubtractVector(&I->x4);
+  HeightB.CopyVector(&intersection);
+  HeightB.SubtractVector(&I->x3);
+  HeightB.ProjectOntoPlane(&SideB);
+
+  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
+  //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
+
+  return retval;
+};
+
+
+/**
+ * Calculates whether there is an intersection between two lines. The first line
+ * always goes through point 1 and point 2 and the second line is given by the
+ * connection between point 4 and point 5.
+ *
+ * @param point 1 of line 1
+ * @param point 2 of line 1
+ * @param point 1 of line 2
+ * @param point 2 of line 2
+ *
+ * @return true if there is an intersection between the given lines, false otherwise
+ */
+bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4)
+{
+  bool result;
+
+  struct Intersection par;
+    par.x1.CopyVector(&point1);
+    par.x2.CopyVector(&point2);
+    par.x3.CopyVector(&point3);
+    par.x4.CopyVector(&point4);
+
+    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
+    gsl_multimin_fminimizer *s = NULL;
+    gsl_vector *ss, *x;
+    gsl_multimin_function minex_func;
+
+    size_t iter = 0;
+    int status;
+    double size;
+
+    /* Starting point */
+    x = gsl_vector_alloc(NDIM);
+    gsl_vector_set(x, 0, point1.x[0]);
+    gsl_vector_set(x, 1, point1.x[1]);
+    gsl_vector_set(x, 2, point1.x[2]);
+
+    /* Set initial step sizes to 1 */
+    ss = gsl_vector_alloc(NDIM);
+    gsl_vector_set_all(ss, 1.0);
+
+    /* Initialize method and iterate */
+    minex_func.n = NDIM;
+    minex_func.f = &MinIntersectDistance;
+    minex_func.params = (void *)&par;
+
+    s = gsl_multimin_fminimizer_alloc(T, NDIM);
+    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);
+
+    do {
+        iter++;
+        status = gsl_multimin_fminimizer_iterate(s);
+
+        if (status) {
+          break;
+        }
+
+        size = gsl_multimin_fminimizer_size(s);
+        status = gsl_multimin_test_size(size, 1e-2);
+
+        if (status == GSL_SUCCESS) {
+          cout << Verbose(2) << "converged to minimum" <<  endl;
+        }
+    } while (status == GSL_CONTINUE && iter < 100);
+
+    // check whether intersection is in between or not
+  Vector intersection, SideA, SideB, HeightA, HeightB;
+  double t1, t2;
+  for (int i = 0; i < NDIM; i++) {
+    intersection.x[i] = gsl_vector_get(s->x, i);
+  }
+
+  SideA.CopyVector(&par.x2);
+  SideA.SubtractVector(&par.x1);
+  HeightA.CopyVector(&intersection);
+  HeightA.SubtractVector(&par.x1);
+
+  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
+
+  SideB.CopyVector(&par.x4);
+  SideB.SubtractVector(&par.x3);
+  HeightB.CopyVector(&intersection);
+  HeightB.SubtractVector(&par.x3);
+
+  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
+
+  cout << Verbose(2) << "Intersection " << intersection << " is at "
+    << t1 << " for (" << point1 << "," << point2 << ") and at "
+    << t2 << " for (" << point3 << "," << point4 << "): ";
+
+  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
+    cout << "true intersection." << endl;
+    result = true;
+  } else {
+    cout << "intersection out of region of interest." << endl;
+    result = false;
+  }
+
+  // free minimizer stuff
+    gsl_vector_free(x);
+    gsl_vector_free(ss);
+    gsl_multimin_fminimizer_free(s);
+
+  return result;
+}
+
Index: src/tesselationhelpers.hpp
===================================================================
--- src/tesselationhelpers.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
+++ src/tesselationhelpers.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -0,0 +1,41 @@
+/*
+ * TesselationHelpers.hpp
+ *
+ *  Auxiliary functions for the tesselation.
+ *
+ *  Created on: Aug 3, 2009
+ *      Author: heber
+ */
+
+#ifndef TESSELATIONHELPERS_HPP_
+#define TESSELATIONHELPERS_HPP_
+
+using namespace std;
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#define HULLEPSILON 1e-7
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
+
+#include "vector.hpp"
+
+double det_get(gsl_matrix *A, int inPlace);
+void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS);
+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);
+void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c);
+double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection);
+double MinIntersectDistance(const gsl_vector * x, void *params);
+bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4);
+
+//bool Choose_preferable_third_point(atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine, atom *ThirdNode, Tesselation *Tess);
+
+
+#endif /* TESSELATIONHELPERS_HPP_ */
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision 2319ed7194f122bf4ff43ec6b7ea796f88df7fb3)
+++ src/vector.hpp	(revision 357fba6ccb2cf9e2d97f501d5968558a706028a5)
@@ -1,4 +1,8 @@
 #ifndef VECTOR_HPP_
 #define VECTOR_HPP_
+
+using namespace std;
+
+#include "helpers.hpp"
 
 class Vector;
