source: src/boundary.cpp@ 0dbddc

#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;
}

// ========================================== 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.
 * Looks at lines parallel to one axis and where they intersect on the projected planes
 * \param *out output stream for debugging
 * \param *BoundaryPoints NDIM set of boundary points defining the convex envelope on each projected plane
 * \param *mol molecule structure representing the cluster
 * \param IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *
GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol,
    bool IsAngstroem)
{
  // get points on boundary of NULL was given as parameter
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = BoundaryPtr;
  if (BoundaryPoints == NULL)
    {
      BoundaryFreeFlag = true;
      BoundaryPoints = GetBoundaryPoints(out, mol);
    }
  else
    {
      *out << Verbose(1) << "Using given boundary points set." << endl;
    }
  // determine biggest "diameter" of cluster for each axis
  Boundaries::iterator Neighbour, OtherNeighbour;
  double *GreatestDiameter = new double[NDIM];
  for (int i = 0; i < NDIM; i++)
    GreatestDiameter[i] = 0.;
  double OldComponent, tmp, w1, w2;
  Vector DistanceVector, OtherVector;
  int component, Othercomponent;
  for (int axis = 0; axis < NDIM; axis++)
    { // regard each projected plane
      //*out << Verbose(1) << "Current axis is " << axis << "." << endl;
      for (int j = 0; j < 2; j++)
        { // and for both axis on the current plane
          component = (axis + j + 1) % NDIM;
          Othercomponent = (axis + 1 + ((j + 1) & 1)) % NDIM;
          //*out << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
              != BoundaryPoints[axis].end(); runner++)
            {
              //*out << Verbose(2) << "Current runner is " << *(runner->second.second) << "." << endl;
              // seek for the neighbours pair where the Othercomponent sign flips
              Neighbour = runner;
              Neighbour++;
              if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                Neighbour = BoundaryPoints[axis].begin();
              DistanceVector.CopyVector(&runner->second.second->x);
              DistanceVector.SubtractVector(&Neighbour->second.second->x);
              do
                { // seek for neighbour pair where it flips
                  OldComponent = DistanceVector.x[Othercomponent];
                  Neighbour++;
                  if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                    Neighbour = BoundaryPoints[axis].begin();
                  DistanceVector.CopyVector(&runner->second.second->x);
                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
                  //*out << Verbose(3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                }
              while ((runner != Neighbour) && (fabs(OldComponent / fabs(
                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
                  DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
              if (runner != Neighbour)
                {
                  OtherNeighbour = Neighbour;
                  if (OtherNeighbour == BoundaryPoints[axis].begin()) // make it wrap around
                    OtherNeighbour = BoundaryPoints[axis].end();
                  OtherNeighbour--;
                  //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                  // now we have found the pair: Neighbour and OtherNeighbour
                  OtherVector.CopyVector(&runner->second.second->x);
                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
                  //*out << Verbose(2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                  //*out << Verbose(2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                  // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
                  w1 = fabs(OtherVector.x[Othercomponent]);
                  w2 = fabs(DistanceVector.x[Othercomponent]);
                  tmp = fabs((w1 * DistanceVector.x[component] + w2
                      * OtherVector.x[component]) / (w1 + w2));
                  // mark if it has greater diameter
                  //*out << Verbose(2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                  GreatestDiameter[component] = (GreatestDiameter[component]
                      > tmp) ? GreatestDiameter[component] : tmp;
                } //else
              //*out << Verbose(2) << "Saw no sign flip, probably top or bottom node." << endl;
            }
        }
    }
  *out << Verbose(0) << "RESULT: The biggest diameters are "
      << GreatestDiameter[0] << " and " << GreatestDiameter[1] << " and "
      << GreatestDiameter[2] << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << "." << endl;

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);

  return GreatestDiameter;
}
;

/** Creates the objects in a VRML file.
 * \param *out output stream for debugging
 * \param *vrmlfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
{
  atom *Walker = mol->start;
  bond *Binder = mol->first;
  int i;
  Vector *center = mol->DetermineCenterOfAll(out);
  if (vrmlfile != NULL) {
    //cout << Verbose(1) << "Writing Raster3D file ... ";
    *vrmlfile << "#VRML V2.0 utf8" << endl;
    *vrmlfile << "#Created by molecuilder" << endl;
    *vrmlfile << "#All atoms as spheres" << endl;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
      for (i=0;i<NDIM;i++)
        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    }

    *vrmlfile << "# All bonds as vertices" << endl;
    while (Binder->next != mol->last) {
      Binder = Binder->next;
      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
      for (i=0;i<NDIM;i++)
        *vrmlfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.03\t";
      for (i=0;i<NDIM;i++)
        *vrmlfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    }

    *vrmlfile << "# All tesselation triangles" << endl;
    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
      *vrmlfile << "1" << endl << "  "; // 1 is triangle type
      for (i=0;i<3;i++) { // print each node
        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
        *vrmlfile << "\t";
      }
      *vrmlfile << "1. 0. 0." << endl;  // red as colour
      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
    }
  } else {
    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
  }
  delete(center);
};

/** Creates the objects in a raster3d file (renderable with a header.r3d).
 * \param *out output stream for debugging
 * \param *rasterfile output stream for tecplot data
 * \param *Tess Tesselation structure with constructed triangles
 * \param *mol molecule structure with atom positions
 */
void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
{
  atom *Walker = mol->start;
  bond *Binder = mol->first;
  int i;
  Vector *center = mol->DetermineCenterOfAll(out);
  if (rasterfile != NULL) {
    //cout << Verbose(1) << "Writing Raster3D file ... ";
    *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
    *rasterfile << "@header.r3d" << endl;
    *rasterfile << "# All atoms as spheres" << endl;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      *rasterfile << "2" << endl << "  ";  // 2 is sphere type
      for (i=0;i<NDIM;i++)
        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
    }

    *rasterfile << "# All bonds as vertices" << endl;
    while (Binder->next != mol->last) {
      Binder = Binder->next;
      *rasterfile << "3" << endl << "  ";  // 2 is round-ended cylinder type
      for (i=0;i<NDIM;i++)
        *rasterfile << Binder->leftatom->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.03\t";
      for (i=0;i<NDIM;i++)
        *rasterfile << Binder->rightatom->x.x[i]-center->x[i] << " ";
      *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
    }

    *rasterfile << "# All tesselation triangles" << endl;
    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
      *rasterfile << "1" << endl << "  ";  // 1 is triangle type
      for (i=0;i<3;i++) {  // print each node
        for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
        *rasterfile << "\t";
      }
      *rasterfile << "1. 0. 0." << endl;  // red as colour
      //*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
    }
    *rasterfile << "9\n  terminating special property\n";
  } else {
    cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
  }
  delete(center);
};

/** This function creates the tecplot file, displaying the tesselation of the hull.
 * \param *out output stream for debugging
 * \param *tecplot output stream for tecplot data
 * \param N arbitrary number to differentiate various zones in the tecplot format
 */
void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
{
  if (tecplot != NULL) {
    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\"" << endl;
    *tecplot << "ZONE T=\"TRIANGLES" << N << "\", N=" << TesselStruct->PointsOnBoundaryCount << ", E=" << TesselStruct->TrianglesOnBoundaryCount << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    int *LookupList = new int[mol->AtomCount];
    for (int i = 0; i < mol->AtomCount; i++)
      LookupList[i] = -1;

    // print atom coordinates
    *out << Verbose(2) << "The following triangles were created:";
    int Counter = 1;
    atom *Walker = NULL;
    for (PointMap::iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
      Walker = target->second->node;
      LookupList[Walker->nr] = Counter++;
      *tecplot << Walker->x.x[0] << " " << Walker->x.x[1] << " " << Walker->x.x[2] << " " << endl;
    }
    *tecplot << endl;
    // print connectivity
    for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
      *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    }
    delete[] (LookupList);
    *out << endl;
  }
}

/** Tesselates the convex boundary by finding all boundary points.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell list
 * \param *filename filename prefix for output of vertex data
 * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
 */
void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
{
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

  cout << Verbose(1) << "Begin of find_convex_border" << endl;

  if (TesselStruct != NULL) // free if allocated
    delete(TesselStruct);
  TesselStruct = new class Tesselation;

  // 1. Find all points on the boundary
  if (BoundaryPoints == NULL) {
      BoundaryFreeFlag = true;
      BoundaryPoints = GetBoundaryPoints(out, mol);
  } else {
      *out << Verbose(1) << "Using given boundary points set." << endl;
  }

// printing all inserted for debugging
  for (int axis=0; axis < NDIM; axis++)
    {
      *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;
    }

  // 2. fill the boundary point list
  for (int axis = 0; axis < NDIM; axis++)
    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
        TesselStruct->AddPoint(runner->second.second);

  *out << Verbose(2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;
  // now we have the whole set of edge points in the BoundaryList

  // listing for debugging
  //  *out << Verbose(1) << "Listing PointsOnBoundary:";
  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
  //    *out << " " << *runner->second;
  //  }
  //  *out << endl;

  // 3a. guess starting triangle
  TesselStruct->GuessStartingTriangle(out);

  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
  TesselStruct->TesselateOnBoundary(out, mol);

  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
  if (!TesselStruct->InsertStraddlingPoints(out, mol))
    *out << Verbose(1) << "Insertion of straddling points failed!" << endl;

  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
  if (!TesselStruct->CorrectConcaveBaselines(out))
    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;

  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;

  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      write_raster3d_file(out, rasterplot, TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }

  // free reference lists
  if (BoundaryFreeFlag)
    delete[] (BoundaryPoints);

  cout << Verbose(1) << "End of find_convex_border" << endl;
};


/** Determines the volume of a cluster.
 * Determines first the convex envelope, then tesselates it and calculates its volume.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *configuration needed for path to store convex envelope file
 * \return determined volume of the cluster in cubed config:GetIsAngstroem()
 */
double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
{
  bool IsAngstroem = configuration->GetIsAngstroem();
  double volume = 0.;
  double PyramidVolume = 0.;
  double G, h;
  Vector x, y;
  double a, b, c;

  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
  *out << Verbose(1)
      << "Calculating the volume of the pyramids formed out of triangles and center of gravity."
      << endl;
  for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner
      != TesselStruct->TrianglesOnBoundary.end(); runner++)
    { // go through every triangle, calculate volume of its pyramid with CoG as peak
      x.CopyVector(&runner->second->endpoints[0]->node->x);
      x.SubtractVector(&runner->second->endpoints[1]->node->x);
      y.CopyVector(&runner->second->endpoints[0]->node->x);
      y.SubtractVector(&runner->second->endpoints[2]->node->x);
      a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
          &runner->second->endpoints[1]->node->x));
      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
          &runner->second->endpoints[2]->node->x));
      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
          &runner->second->endpoints[1]->node->x));
      G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
      x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
          &runner->second->endpoints[1]->node->x,
          &runner->second->endpoints[2]->node->x);
      x.Scale(runner->second->endpoints[1]->node->x.Projection(&x));
      h = x.Norm(); // distance of CoG to triangle
      PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
      *out << Verbose(2) << "Area of triangle is " << G << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
          << h << " and the volume is " << PyramidVolume << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
      volume += PyramidVolume;
    }
  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10)
      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
      << endl;

  return volume;
}
;

/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
 * We get cluster volume by VolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
 * \param *out output stream for debugging
 * \param *configuration needed for path to store convex envelope file
 * \param *mol molecule structure representing the cluster
 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
 * \param celldensity desired average density in final cell
 */
void
PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol,
    double ClusterVolume, double celldensity)
{
  // transform to PAS
  mol->PrincipalAxisSystem(out, true);

  // some preparations beforehand
  bool IsAngstroem = configuration->GetIsAngstroem();
  Boundaries *BoundaryPoints = GetBoundaryPoints(out, mol);
  class Tesselation *TesselStruct = NULL;
  LinkedCell LCList(mol, 10.);
  Find_convex_border(out, mol, TesselStruct, &LCList, NULL);
  double clustervolume;
  if (ClusterVolume == 0)
    clustervolume = VolumeOfConvexEnvelope(out, TesselStruct, configuration);
  else
    clustervolume = ClusterVolume;
  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
  Vector BoxLengths;
  int repetition[NDIM] =
    { 1, 1, 1 };
  int TotalNoClusters = 1;
  for (int i = 0; i < NDIM; i++)
    TotalNoClusters *= repetition[i];

  // sum up the atomic masses
  double totalmass = 0.;
  atom *Walker = mol->start;
  while (Walker->next != mol->end)
    {
      Walker = Walker->next;
      totalmass += Walker->type->mass;
    }
  *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10)
      << totalmass << " atomicmassunit." << endl;

  *out << Verbose(0) << "RESULT: The average density is " << setprecision(10)
      << totalmass / clustervolume << " atomicmassunit/"
      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;

  // solve cubic polynomial
  *out << Verbose(1) << "Solving equidistant suspension in water problem ..."
      << endl;
  double cellvolume;
  if (IsAngstroem)
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass
        / clustervolume)) / (celldensity - 1);
  else
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass
        / clustervolume)) / (celldensity - 1);
  *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity
      << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom"
      : "atomiclength") << "^3." << endl;

  double minimumvolume = TotalNoClusters * (GreatestDiameter[0]
      * GreatestDiameter[1] * GreatestDiameter[2]);
  *out << Verbose(1)
      << "Minimum volume of the convex envelope contained in a rectangular box is "
      << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom"
      : "atomiclength") << "^3." << endl;
  if (minimumvolume > cellvolume)
    {
      cerr << Verbose(0)
          << "ERROR: the containing box already has a greater volume than the envisaged cell volume!"
          << endl;
      cout << Verbose(0)
          << "Setting Box dimensions to minimum possible, the greatest diameters."
          << endl;
      for (int i = 0; i < NDIM; i++)
        BoxLengths.x[i] = GreatestDiameter[i];
      mol->CenterEdge(out, &BoxLengths);
    }
  else
    {
      BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1]
          * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
      BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0]
          * GreatestDiameter[1] + repetition[0] * repetition[2]
          * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1]
          * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
      BoxLengths.x[2] = minimumvolume - cellvolume;
      double x0 = 0., x1 = 0., x2 = 0.;
      if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1],
          BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
        *out << Verbose(0) << "RESULT: The resulting spacing is: " << x0
            << " ." << endl;
      else
        {
          *out << Verbose(0) << "RESULT: The resulting spacings are: " << x0
              << " and " << x1 << " and " << x2 << " ." << endl;
          x0 = x2; // sorted in ascending order
        }

      cellvolume = 1;
      for (int i = 0; i < NDIM; i++)
        {
          BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
          cellvolume *= BoxLengths.x[i];
        }

      // set new box dimensions
      *out << Verbose(0) << "Translating to box with these boundaries." << endl;
      mol->SetBoxDimension(&BoxLengths);
      mol->CenterInBox((ofstream *) &cout);
    }
  // update Box of atoms by boundary
  mol->SetBoxDimension(&BoxLengths);
  *out << Verbose(0) << "RESULT: The resulting cell dimensions are: "
      << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and "
      << BoxLengths.x[2] << " with total volume of " << cellvolume << " "
      << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
}
;


/** Fills the empty space of the simulation box with water/
 * \param *out output stream for debugging
 * \param *List list of molecules already present in box
 * \param *TesselStruct contains tesselated surface
 * \param *filler molecule which the box is to be filled with
 * \param configuration contains box dimensions
 * \param distance[NDIM] distance between filling molecules in each direction
 * \param RandAtomDisplacement maximum distance for random displacement per atom
 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 * \param DoRandomRotation true - do random rotiations, false - don't
 * \return *mol pointer to new molecule with filled atoms
 */
molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
{
  molecule *Filling = new molecule(filler->elemente);
  Vector CurrentPosition;
  int N[NDIM];
  int n[NDIM];
  double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
  double Rotations[NDIM*NDIM];
  Vector AtomTranslations;
  Vector FillerTranslations;
  Vector FillerDistance;
  double FillIt = false;
  atom *Walker = NULL, *Runner = NULL;
  bond *Binder = NULL;

  // Center filler at origin
  filler->CenterOrigin(out);
  filler->Center.Zero();

  // calculate filler grid in [0,1]^3
  FillerDistance.Init(distance[0], distance[1], distance[2]);
  FillerDistance.InverseMatrixMultiplication(M);
  for(int i=0;i<NDIM;i++)
    N[i] = (int) ceil(1./FillerDistance.x[i]);

  // go over [0,1]^3 filler grid
  for (n[0] = 0; n[0] < N[0]; n[0]++)
    for (n[1] = 0; n[1] < N[1]; n[1]++)
      for (n[2] = 0; n[2] < N[2]; n[2]++) {
        // calculate position of current grid vector in untransformed box
        CurrentPosition.Init(n[0], n[1], n[2]);
        CurrentPosition.MatrixMultiplication(M);
        // Check whether point is in- or outside
        FillIt = true;
        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
          //FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
        }

        if (FillIt) {
          // fill in Filler
          *out << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;

          // create molecule random translation vector ...
          for (int i=0;i<NDIM;i++)
            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);

          // go through all atoms
          Walker = filler->start;
          while (Walker != filler->end) {
            Walker = Walker->next;
            // copy atom ...
            *Runner = new atom(Walker);

            // create atomic random translation vector ...
            for (int i=0;i<NDIM;i++)
              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);

            // ... and rotation matrix
            if (DoRandomRotation) {
              double phi[NDIM];
              for (int i=0;i<NDIM;i++) {
                phi[i] = rand()/(RAND_MAX/(2.*M_PI));
              }

              Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
              Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
              Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
              Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
              Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
              Rotations[7] =               sin(phi[1])                                                ;
              Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
              Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
              Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
            }

            // ... and put at new position
            if (DoRandomRotation)
              Runner->x.MatrixMultiplication(Rotations);
            Runner->x.AddVector(&AtomTranslations);
            Runner->x.AddVector(&FillerTranslations);
            Runner->x.AddVector(&CurrentPosition);
            // insert into Filling
            Filling->AddAtom(Runner);
          }

          // go through all bonds and add as well
          Binder = filler->first;
          while(Binder != filler->last) {
            Binder = Binder->next;
          }
        } else {
          // leave empty
          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
        }
      }
  return Filling;
};


// =========================================================== 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]) {
// TODO: use findTriangles here and return result.size();
  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) {
  // TODO: use get angle and remove duplicate code
  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;
}

/**
 * Checks whether the provided atom is within the tesselation structure.
 *
 * @param atom of which to check the position
 * @param tesselation structure
 *
 * @return true if the atom is inside the tesselation structure, false otherwise
 */
bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC) {
  if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
    cout << Verbose(0) << "Error: There is no tesselation structure to compare the atom with, "
        << "please create one first.";
    exit(1);
  }

  class atom *trianglePoints[3];
  trianglePoints[0] = findClosestAtom(Atom, LC);
  // check whether closest atom is "too close" :), then it's inside
  if (trianglePoints[0]->x.DistanceSquared(&Atom->x) < MYEPSILON)
    return true;
  list<atom*> *connectedClosestAtoms = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
  trianglePoints[1] = connectedClosestAtoms->front();
  trianglePoints[2] = connectedClosestAtoms->back();
  for (int i=0;i<3;i++) {
    if (trianglePoints[i] == NULL) {
      cout << Verbose(1) << "IsInnerAtom encounters serious error, point " << i << " not found." << endl;
    }

    cout << Verbose(1) << "List of possible atoms:" << endl;
    cout << *trianglePoints[i] << endl;

//    for(list<atom*>::iterator runner = connectedClosestAtoms->begin(); runner != connectedClosestAtoms->end(); runner++)
//      cout << Verbose(2) << *(*runner) << endl;
  }
  delete(connectedClosestAtoms);

  list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);

  if (triangles->empty()) {
    cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
    exit(1);
  }

  Vector helper;
  helper.CopyVector(&Atom->x);

  // Only in case of degeneration, there will be two different scalar products.
  // If at least one scalar product is positive, the point is considered to be outside.
  bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
      && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
  delete(triangles);
  return status;
}

/**
 * Finds the atom which is closest to the provided one.
 *
 * @param atom to which to find the closest other atom
 * @param linked cell structure
 *
 * @return atom which is closest to the provided one
 */
atom* findClosestAtom(const atom* Atom, LinkedCell* LC) {
  LinkedAtoms *List = NULL;
  atom* closestAtom = NULL;
  double distance = 1e16;
  Vector helper;
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];

  LC->SetIndexToVector(&Atom->x); // ignore status as we calculate bounds below sensibly
  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;

  LC->GetNeighbourBounds(Nlower, Nupper);
  //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 << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
        if (List != NULL) {
          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            helper.CopyVector(&Atom->x);
            helper.SubtractVector(&(*Runner)->x);
            double currentNorm = helper. Norm();
            if (currentNorm < distance) {
              distance = currentNorm;
              closestAtom = (*Runner);
            }
          }
        } else {
          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
            << LC->n[2] << " is invalid!" << endl;
        }
      }

  return closestAtom;
}

/**
 * Gets all atoms connected to the provided atom by triangulation lines.
 *
 * @param atom of which get all connected atoms
 * @param atom to be checked whether it is an inner atom
 *
 * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
 */
list<atom*> * Tesselation::getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck) {
  list<atom*> connectedAtoms;
  map<double, atom*> anglesOfAtoms;
  map<double, atom*>::iterator runner;
  LineMap::iterator findLines = LinesOnBoundary.begin();
  list<atom*>::iterator listRunner;
  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
  atom* current;
  bool takeAtom = false;

  planeNorm.CopyVector(&Atom->x);
  planeNorm.SubtractVector(&AtomToCheck->x);
  planeNorm.Normalize();

  while (findLines != LinesOnBoundary.end()) {
    takeAtom = false;

    if (findLines->second->endpoints[0]->Nr == Atom->nr) {
      takeAtom = true;
      current = findLines->second->endpoints[1]->node;
    } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
      takeAtom = true;
      current = findLines->second->endpoints[0]->node;
    }

    if (takeAtom) {
      connectedAtoms.push_back(current);
      currentPoint.CopyVector(&current->x);
      currentPoint.ProjectOntoPlane(&planeNorm);
      center.AddVector(&currentPoint);
    }

    findLines++;
  }

  cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
    << "; scale factor " << 1.0/connectedAtoms.size();

  center.Scale(1.0/connectedAtoms.size());
  listRunner = connectedAtoms.begin();

  cout << " calculated center " << center <<  endl;

  // construct one orthogonal vector
  helper.CopyVector(&AtomToCheck->x);
  helper.ProjectOntoPlane(&planeNorm);
  OrthogonalVector.MakeNormalVector(&center, &helper, &(*listRunner)->x);
  while (listRunner != connectedAtoms.end()) {
    double angle = getAngle((*listRunner)->x, AtomToCheck->x, center, OrthogonalVector);
    cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
    anglesOfAtoms.insert(pair<double, atom*>(angle, (*listRunner)));
    listRunner++;
  }

  list<atom*> *result = new list<atom*>;
  runner = anglesOfAtoms.begin();
  cout << "First value is " << *runner->second << endl;
  result->push_back(runner->second);
  runner = anglesOfAtoms.end();
  runner--;
  cout << "Second value is " << *runner->second << endl;
  result->push_back(runner->second);

  cout << "List of closest atoms has " << result->size() << " elements, which are "
    << *(result->front()) << " and " << *(result->back()) << endl;

  return result;
}

/**
 * Finds triangles belonging to the three provided atoms.
 *
 * @param atom list, is expected to contain three atoms
 *
 * @return triangles which belong to the provided atoms, will be empty if there are none,
 *         will usually be one, in case of degeneration, there will be two
 */
list<BoundaryTriangleSet*> *Tesselation::FindTriangles(atom* Points[3]) {
  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
  LineMap::iterator FindLine;
  PointMap::iterator FindPoint;
  TriangleMap::iterator FindTriangle;
  class BoundaryPointSet *TrianglePoints[3];

  for (int i = 0; i < 3; i++) {
    FindPoint = PointsOnBoundary.find(Points[i]->nr);
    if (FindPoint != PointsOnBoundary.end()) {
      TrianglePoints[i] = FindPoint->second;
    } else {
      TrianglePoints[i] = NULL;
    }
  }

  // checks lines between the points in the Points for their adjacent triangles
  for (int i = 0; i < 3; i++) {
    if (TrianglePoints[i] != NULL) {
      for (int j = i; j < 3; j++) {
        if (TrianglePoints[j] != NULL) {
          FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
          if (FindLine != TrianglePoints[i]->lines.end()) {
            for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
              FindTriangle = FindLine->second->triangles.begin();
              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
                if ((
                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
                  ) && (
                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
                  ) && (
                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
                  )
                ) {
                  result->push_back(FindTriangle->second);
                }
              }
            }
            // Is it sufficient to consider one of the triangle lines for this.
            return result;

          }
        }
      }
    }
  }

  return result;
}

/**
 * Gets the angle between a point and a reference relative to the provided center.
 *
 * @param point to calculate the angle for
 * @param reference to which to calculate the angle
 * @param center for which to calculate the angle between the vectors
 * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
 *
 * @return angle between point and reference
 */
double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) {
  Vector ReferenceVector, helper;
  cout << Verbose(2) << reference << " is our reference " << endl;
  cout << Verbose(2) << center << " is our center " << endl;
  // create baseline vector
  ReferenceVector.CopyVector(&reference);
  ReferenceVector.SubtractVector(&center);
  if (ReferenceVector.IsNull())
    return M_PI;

  // calculate both angles and correct with in-plane vector
  helper.CopyVector(&point);
  helper.SubtractVector(&center);
  if (helper.IsNull())
    return M_PI;
  double phi = ReferenceVector.Angle(&helper);
  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    phi = 2.*M_PI - phi;
  }

  cout << Verbose(2) << point << " has angle " << phi << endl;

  return phi;
}

/**
 * Checks whether the provided point is within the tesselation structure.
 *
 * This is a wrapper function for IsInnerAtom, so it can be used with a simple
 * vector, too.
 *
 * @param point of which to check the position
 * @param tesselation structure
 *
 * @return true if the point is inside the tesselation structure, false otherwise
 */
bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) {
  atom *temp_atom = new atom;
  bool status = false;
  temp_atom->x.CopyVector(&Point);

  status = IsInnerAtom(temp_atom, Tess, LC);
  delete(temp_atom);

  return status;
}

/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
 * \param *out output stream for debugging
 * \param *mol molecule structure with Atom's and Bond's
 * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
 * \param *LCList atoms in LinkedCell list
 * \param *filename filename prefix for output of vertex data
 * \para RADIUS radius of the virtual sphere
 */
void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, const char *filename, const double RADIUS)
{
  int N = 0;
  bool freeTess = false;
  bool freeLC = false;
  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
  if (Tess == NULL) {
    *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    Tess = new Tesselation;
    freeTess = true;
  }
  LineMap::iterator baseline;
  LineMap::iterator testline;
  *out << Verbose(0) << "Begin of Find_non_convex_border\n";
  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
  bool failflag = false;

  if (LCList == NULL) {
    LCList = new LinkedCell(mol, 2.*RADIUS);
    freeLC = true;
  }

  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);

  baseline = Tess->LinesOnBoundary.begin();
  while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
    if (baseline->second->TrianglesCount == 1) {
      failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
      flag = flag || failflag;
      if (!failflag)
        cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
    } else {
      //cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
      if (baseline->second->TrianglesCount != 2)
        cout << Verbose(1) << "ERROR: TESSELATION FINISHED WITH INVALID TRIANGLE COUNT!" << endl;
    }

    N++;
    baseline++;
    if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
      baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
      flag = false;
    }
  }
  if (filename != 0) {
    *out << Verbose(1) << "Writing final tecplot file\n";
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      write_tecplot_file(out, tecplot, Tess, mol, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *raster = new ofstream(OutputName.c_str());
      write_raster3d_file(out, raster, Tess, mol);
      raster->close();
      delete(raster);
    }
  } else {
    cerr << "ERROR: Could definitively not find all necessary triangles!" << endl;
  }

  cout << Verbose(2) << "Check: List of Baselines with not two connected triangles:" << endl;
  int counter = 0;
  for (testline = Tess->LinesOnBoundary.begin(); testline != Tess->LinesOnBoundary.end(); testline++) {
    if (testline->second->TrianglesCount != 2) {
      cout << Verbose(2) << *testline->second << "\t" << testline->second->TrianglesCount << endl;
      counter++;
    }
  }
  if (counter == 0)
    cout << Verbose(2) << "None." << endl;

  // Tests the IsInnerAtom() function.
  Vector x (0, 0, 0);
  cout << Verbose(0) << "Point to check: " << x << endl;
  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
    << "for vector " << x << "." << endl;
  atom* a = Tess->PointsOnBoundary.begin()->second->node;
  cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
    << "for atom " << a << " (on boundary)." << endl;
  LinkedAtoms *List = NULL;
  for (int i=0;i<NDIM;i++) { // each axis
    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
        List = LCList->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 (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
              a = *Runner;
              i=3;
              for (int j=0;j<NDIM;j++)
                LCList->n[j] = LCList->N[j];
              break;
            }
          }
        }
      }
  }
  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
    << "for atom " << a << " (inside)." << endl;


  if (freeTess)
    delete(Tess);
  if (freeLC)
    delete(LCList);
  *out << Verbose(0) << "End of Find_non_convex_border\n";
};

/** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
 * \param *out output stream for debugging
 * \param *srcmol molecule to embed into
 * \return *Vector new center of \a *srcmol for embedding relative to \a this
 */
Vector* molecule::FindEmbeddingHole(ofstream *out, molecule *srcmol)
{
  Vector *Center = new Vector;
  Center->Zero();
  // calculate volume/shape of \a *srcmol

  // find embedding holes

  // if more than one, let user choose

  // return embedding center
  return Center;
};