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

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

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.);
}
;

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;

  *out << Verbose(1) << "Finding all boundary points." << endl;
  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
  BoundariesTestPair BoundaryTestPair;
  Vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
  double radius, angle;
  // 3a. Go through every axis
  for (int axis = 0; axis < NDIM; axis++)
    {
      AxisVector.Zero();
      AngleReferenceVector.Zero();
      AngleReferenceNormalVector.Zero();
      AxisVector.x[axis] = 1.;
      AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
      AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
      //    *out << Verbose(1) << "Axisvector is ";
      //    AxisVector.Output(out);
      //    *out << " and AngleReferenceVector is ";
      //    AngleReferenceVector.Output(out);
      //    *out << "." << endl;
      //    *out << " and AngleReferenceNormalVector is ";
      //    AngleReferenceNormalVector.Output(out);
      //    *out << "." << endl;
      // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
      Walker = mol->start;
      while (Walker->next != mol->end)
        {
          Walker = Walker->next;
          Vector ProjectedVector;
          ProjectedVector.CopyVector(&Walker->x);
          ProjectedVector.ProjectOntoPlane(&AxisVector);
          // correct for negative side
          //if (Projection(y) < 0)
          //angle = 2.*M_PI - angle;
          radius = ProjectedVector.Norm();
          if (fabs(radius) > MYEPSILON)
            angle = ProjectedVector.Angle(&AngleReferenceVector);
          else
            angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues

          //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
          if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0)
            {
              angle = 2. * M_PI - angle;
            }
          //*out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
          //ProjectedVector.Output(out);
          //*out << endl;
          BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle,
              DistancePair (radius, Walker)));
          if (BoundaryTestPair.second)
            { // successfully inserted
            }
          else
            { // same point exists, check first r, then distance of original vectors to center of gravity
              *out << Verbose(2)
                  << "Encountered two vectors whose projection onto axis "
                  << axis << " is equal: " << endl;
              *out << Verbose(2) << "Present vector: ";
              BoundaryTestPair.first->second.second->x.Output(out);
              *out << endl;
              *out << Verbose(2) << "New vector: ";
              Walker->x.Output(out);
              *out << endl;
              double tmp = ProjectedVector.Norm();
              if (tmp > BoundaryTestPair.first->second.first)
                {
                  BoundaryTestPair.first->second.first = tmp;
                  BoundaryTestPair.first->second.second = Walker;
                  *out << Verbose(2) << "Keeping new vector." << endl;
                }
              else if (tmp == BoundaryTestPair.first->second.first)
                {
                  if (BoundaryTestPair.first->second.second->x.ScalarProduct(
                      &BoundaryTestPair.first->second.second->x)
                      < Walker->x.ScalarProduct(&Walker->x))
                    { // Norm() does a sqrt, which makes it a lot slower
                      BoundaryTestPair.first->second.second = Walker;
                      *out << Verbose(2) << "Keeping new vector." << endl;
                    }
                  else
                    {
                      *out << Verbose(2) << "Keeping present vector." << endl;
                    }
                }
              else
                {
                  *out << Verbose(2) << "Keeping present vector." << endl;
                }
            }
        }
      // printing all inserted for debugging
      //    {
      //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
      //      int i=0;
      //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
      //        if (runner != BoundaryPoints[axis].begin())
      //          *out << ", " << i << ": " << *runner->second.second;
      //        else
      //          *out << i << ": " << *runner->second.second;
      //        i++;
      //      }
      //      *out << endl;
      //    }
      // 3c. throw out points whose distance is less than the mean of left and right neighbours
      bool flag = false;
      do
        { // do as long as we still throw one out per round
          *out << Verbose(1)
              << "Looking for candidates to kick out by convex condition ... "
              << endl;
          flag = false;
          Boundaries::iterator left = BoundaryPoints[axis].end();
          Boundaries::iterator right = BoundaryPoints[axis].end();
          for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
              != BoundaryPoints[axis].end(); runner++)
            {
              // set neighbours correctly
              if (runner == BoundaryPoints[axis].begin())
                {
                  left = BoundaryPoints[axis].end();
                }
              else
                {
                  left = runner;
                }
              left--;
              right = runner;
              right++;
              if (right == BoundaryPoints[axis].end())
                {
                  right = BoundaryPoints[axis].begin();
                }
              // check distance

              // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
                {
                  Vector SideA, SideB, SideC, SideH;
                  SideA.CopyVector(&left->second.second->x);
                  SideA.ProjectOntoPlane(&AxisVector);
                  //          *out << "SideA: ";
                  //          SideA.Output(out);
                  //          *out << endl;

                  SideB.CopyVector(&right->second.second->x);
                  SideB.ProjectOntoPlane(&AxisVector);
                  //          *out << "SideB: ";
                  //          SideB.Output(out);
                  //          *out << endl;

                  SideC.CopyVector(&left->second.second->x);
                  SideC.SubtractVector(&right->second.second->x);
                  SideC.ProjectOntoPlane(&AxisVector);
                  //          *out << "SideC: ";
                  //          SideC.Output(out);
                  //          *out << endl;

                  SideH.CopyVector(&runner->second.second->x);
                  SideH.ProjectOntoPlane(&AxisVector);
                  //          *out << "SideH: ";
                  //          SideH.Output(out);
                  //          *out << endl;

                  // calculate each length
                  double a = SideA.Norm();
                  //double b = SideB.Norm();
                  //double c = SideC.Norm();
                  double h = SideH.Norm();
                  // calculate the angles
                  double alpha = SideA.Angle(&SideH);
                  double beta = SideA.Angle(&SideC);
                  double gamma = SideB.Angle(&SideH);
                  double delta = SideC.Angle(&SideH);
                  double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha
                      < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
                  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
                  //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
                  if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance))
                      < MYEPSILON) && (h < MinDistance))
                    {
                      // throw out point
                      //*out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
                      BoundaryPoints[axis].erase(runner);
                      flag = true;
                    }
                }
            }
        }
      while (flag);
    }
  return BoundaryPoints;
}
;

/** 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)
{
  atom *Walker = NULL;
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

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

  // 1. calculate center of gravity
  *out << endl;
  Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);

  // 2. translate all points into CoG
  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
  Walker = mol->start;
  while (Walker->next != mol->end) {
    Walker = Walker->next;
    Walker->x.Translate(CenterOfGravity);
  }

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

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

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

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

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

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

  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
      << " triangles with " << TesselStruct->LinesOnBoundaryCount
      << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
      << endl;

  // 6. translate all points back from CoG
  *out << Verbose(1) << "Translating system back from Center of Gravity."
      << endl;
  CenterOfGravity->Scale(-1);
  Walker = mol->start;
  while (Walker->next != mol->end)
    {
      Walker = Walker->next;
      Walker->x.Translate(CenterOfGravity);
    }

  // 7. Store triangles in tecplot file
  if (filename != NULL) {
    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);
    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);

};


/** 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;
}
;

// =========================================================== 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 << " leave "
                  << checker->second->node->Name << " outside the convex hull."
                  << endl;
              break;
            }
          else
            { // note the sign for later
              *out << Verbose(2) << "Current candidates: "
                  << A->second->node->Name << ","
                  << baseline->second.first->second->node->Name << ","
                  << baseline->second.second->second->node->Name << " leave "
                  << checker->second->node->Name << " inside the convex hull."
                  << endl;
              sign = tmp;
            }
          // 4d. Check whether the point is inside the triangle (check distance to each node
          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
          int innerpoint = 0;
          if ((tmp < A->second->node->x.DistanceSquared(
              &baseline->second.first->second->node->x)) && (tmp
              < A->second->node->x.DistanceSquared(
                  &baseline->second.second->second->node->x)))
            innerpoint++;
          tmp = checker->second->node->x.DistanceSquared(
              &baseline->second.first->second->node->x);
          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
              &A->second->node->x)) && (tmp
              < baseline->second.first->second->node->x.DistanceSquared(
                  &baseline->second.second->second->node->x)))
            innerpoint++;
          tmp = checker->second->node->x.DistanceSquared(
              &baseline->second.second->second->node->x);
          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
              &baseline->second.first->second->node->x)) && (tmp
              < baseline->second.second->second->node->x.DistanceSquared(
                  &A->second->node->x)))
            innerpoint++;
          // 4e. If so, break 4. loop and continue with next candidate in 1. loop
          if (innerpoint == 3)
            break;
        }
      // 5. come this far, all on same side? Then break 1. loop and construct triangle
      if (checker == PointsOnBoundary.end())
        {
          *out << "Looks like we have a candidate!" << endl;
          break;
        }
    }
  if (baseline != DistanceMMap.end())
    {
      BPS[0] = baseline->second.first->second;
      BPS[1] = baseline->second.second->second;
      BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
      BPS[0] = A->second;
      BPS[1] = baseline->second.second->second;
      BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
      BPS[0] = baseline->second.first->second;
      BPS[1] = A->second;
      BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);

      // 4b3. insert created triangle
      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
      TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
      TrianglesOnBoundaryCount++;
      for (int i = 0; i < NDIM; i++)
        {
          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
          LinesOnBoundaryCount++;
        }

      *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
    }
  else
    {
      *out << Verbose(1) << "No starting triangle found." << endl;
      exit(255);
    }
}
;

/** 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 with the baseline and the current point has the smallest of angles (comparison between normal vectors
 *    -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
 *       baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
 *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
 * \param *out output stream for debugging
 * \param *configuration for IsAngstroem
 * \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,
      PropagationVector;
  LineMap::iterator LineChecker[2];
  do
    {
      flag = false;
      for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline
          != LinesOnBoundary.end(); baseline++)
        if (baseline->second->TrianglesCount == 1)
          {
            *out << Verbose(2) << "Current baseline is between "
                << *(baseline->second) << "." << endl;
            // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
            SmallestAngle = M_PI;
            BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
            // get peak point with respect to this base line's only triangle
            for (int i = 0; i < 3; i++)
              if ((BTS->endpoints[i] != baseline->second->endpoints[0])
                  && (BTS->endpoints[i] != baseline->second->endpoints[1]))
                peak = BTS->endpoints[i];
            *out << Verbose(3) << " and has peak " << *peak << "." << endl;
            // normal vector of triangle
            BTS->GetNormalVector(NormalVector);
            *out << Verbose(4) << "NormalVector of base triangle is ";
            NormalVector.Output(out);
            *out << endl;
            // offset to center of triangle
            CenterVector.Zero();
            for (int i = 0; i < 3; i++)
              CenterVector.AddVector(&BTS->endpoints[i]->node->x);
            CenterVector.Scale(1. / 3.);
            *out << Verbose(4) << "CenterVector of base triangle is ";
            CenterVector.Output(out);
            *out << endl;
            // vector in propagation direction (out of triangle)
            // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
            TempVector.SubtractVector(&baseline->second->endpoints[1]->node->x);
            PropagationVector.MakeNormalVector(&TempVector, &NormalVector);
            TempVector.CopyVector(&CenterVector);
            TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
            //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
            if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
              PropagationVector.Scale(-1.);
            *out << Verbose(4) << "PropagationVector of base triangle is ";
            PropagationVector.Output(out);
            *out << endl;
            winner = PointsOnBoundary.end();
            for (PointMap::iterator target = PointsOnBoundary.begin(); target
                != PointsOnBoundary.end(); target++)
              if ((target->second != baseline->second->endpoints[0])
                  && (target->second != baseline->second->endpoints[1]))
                { // don't take the same endpoints
                  *out << Verbose(3) << "Target point is " << *(target->second)
                      << ":";
                  bool continueflag = true;

                  VirtualNormalVector.CopyVector(
                      &baseline->second->endpoints[0]->node->x);
                  VirtualNormalVector.AddVector(
                      &baseline->second->endpoints[0]->node->x);
                  VirtualNormalVector.Scale(-1. / 2.); // points now to center of base line
                  VirtualNormalVector.AddVector(&target->second->node->x); // points from center of base line to target
                  TempAngle = VirtualNormalVector.Angle(&PropagationVector);
                  continueflag = continueflag && (TempAngle < (M_PI/2.)); // no bends bigger than Pi/2 (90 degrees)
                  if (!continueflag)
                    {
                      *out << Verbose(4)
                          << "Angle between propagation direction and base line to "
                          << *(target->second) << " is " << TempAngle
                          << ", bad direction!" << endl;
                      continue;
                    }
                  else
                    *out << Verbose(4)
                        << "Angle between propagation direction and base line to "
                        << *(target->second) << " is " << TempAngle
                        << ", good direction!" << endl;
                  LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                      target->first);
                  LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                      target->first);
                  //            if (LineChecker[0] != baseline->second->endpoints[0]->lines.end())
                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
                  //            else
                  //              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                  //            if (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
                  //            else
                  //              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has no line to " << *(target->second) << " as endpoint." << endl;
                  // check first endpoint (if any connecting line goes to target or at least not more than 1)
                  continueflag = continueflag && (((LineChecker[0]
                      == baseline->second->endpoints[0]->lines.end())
                      || (LineChecker[0]->second->TrianglesCount == 1)));
                  if (!continueflag)
                    {
                      *out << Verbose(4) << *(baseline->second->endpoints[0])
                          << " has line " << *(LineChecker[0]->second)
                          << " to " << *(target->second)
                          << " as endpoint with "
                          << LineChecker[0]->second->TrianglesCount
                          << " triangles." << endl;
                      continue;
                    }
                  // check second endpoint (if any connecting line goes to target or at least not more than 1)
                  continueflag = continueflag && (((LineChecker[1]
                      == baseline->second->endpoints[1]->lines.end())
                      || (LineChecker[1]->second->TrianglesCount == 1)));
                  if (!continueflag)
                    {
                      *out << Verbose(4) << *(baseline->second->endpoints[1])
                          << " has line " << *(LineChecker[1]->second)
                          << " to " << *(target->second)
                          << " as endpoint with "
                          << LineChecker[1]->second->TrianglesCount
                          << " triangles." << endl;
                      continue;
                    }
                  // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
                  continueflag = continueflag && (!(((LineChecker[0]
                      != baseline->second->endpoints[0]->lines.end())
                      && (LineChecker[1]
                          != baseline->second->endpoints[1]->lines.end())
                      && (GetCommonEndpoint(LineChecker[0]->second,
                          LineChecker[1]->second) == peak))));
                  if (!continueflag)
                    {
                      *out << Verbose(4) << "Current target is peak!" << endl;
                      continue;
                    }
                  // in case NOT both were found
                  if (continueflag)
                    { // create virtually this triangle, get its normal vector, calculate angle
                      flag = true;
                      VirtualNormalVector.MakeNormalVector(
                          &baseline->second->endpoints[0]->node->x,
                          &baseline->second->endpoints[1]->node->x,
                          &target->second->node->x);
                      // make it always point inward
                      if (baseline->second->endpoints[0]->node->x.Projection(
                          &VirtualNormalVector) > 0)
                        VirtualNormalVector.Scale(-1.);
                      // calculate angle
                      TempAngle = NormalVector.Angle(&VirtualNormalVector);
                      *out << Verbose(4) << "NormalVector is ";
                      VirtualNormalVector.Output(out);
                      *out << " and the angle is " << TempAngle << "." << endl;
                      if (SmallestAngle > TempAngle)
                        { // set to new possible winner
                          SmallestAngle = TempAngle;
                          winner = target;
                        }
                    }
                }
            // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
            if (winner != PointsOnBoundary.end())
              {
                *out << Verbose(2) << "Winning target point is "
                    << *(winner->second) << " with angle " << SmallestAngle
                    << "." << endl;
                // create the lins of not yet present
                BLS[0] = baseline->second;
                // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
                LineChecker[0] = baseline->second->endpoints[0]->lines.find(
                    winner->first);
                LineChecker[1] = baseline->second->endpoints[1]->lines.find(
                    winner->first);
                if (LineChecker[0]
                    == baseline->second->endpoints[0]->lines.end())
                  { // create
                    BPS[0] = baseline->second->endpoints[0];
                    BPS[1] = winner->second;
                    BLS[1] = new class BoundaryLineSet(BPS,
                        LinesOnBoundaryCount);
                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                        BLS[1]));
                    LinesOnBoundaryCount++;
                  }
                else
                  BLS[1] = LineChecker[0]->second;
                if (LineChecker[1]
                    == baseline->second->endpoints[1]->lines.end())
                  { // create
                    BPS[0] = baseline->second->endpoints[1];
                    BPS[1] = winner->second;
                    BLS[2] = new class BoundaryLineSet(BPS,
                        LinesOnBoundaryCount);
                    LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount,
                        BLS[2]));
                    LinesOnBoundaryCount++;
                  }
                else
                  BLS[2] = LineChecker[1]->second;
                BTS = new class BoundaryTriangleSet(BLS,
                    TrianglesOnBoundaryCount);
                TrianglesOnBoundary.insert(TrianglePair(
                    TrianglesOnBoundaryCount, BTS));
                TrianglesOnBoundaryCount++;
              }
            else
              {
                *out << Verbose(1)
                    << "I could not determine a winner for this baseline "
                    << *(baseline->second) << "." << endl;
              }

            // 5d. If the set of lines is not yet empty, go to 5. and continue
          }
        else
          *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
              << " has a triangle count of "
              << baseline->second->TrianglesCount << "." << endl;
    }
  while (flag);

}
;

/** 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++;
}
;

/** Adds point to Tesselation::PointsOnBoundary if not yet present.
 * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
 * @param Candidate point to add
 * @param n index for this point in Tesselation::TPS array
 */
void
Tesselation::AddTrianglePoint(atom* Candidate, int n)
{
  PointTestPair InsertUnique;
  TPS[n] = new class BoundaryPointSet(Candidate);
  InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
  if (InsertUnique.second) { // if new point was not present before, increase counter
    PointsOnBoundaryCount++;
  } else {
    delete TPS[n];
    cout << Verbose(3) << "Atom " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;
    TPS[n] = (InsertUnique.first)->second;
  }
}
;

/** Function tries to add line from current Points in BPS to BoundaryLineSet.
 * If successful it raises the line count and inserts the new line into the BLS,
 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n) {
  bool insertNewLine = true;

  if (a->lines.find(b->node->nr) != a->lines.end()) {
    LineMap::iterator FindLine;
    pair<LineMap::iterator,LineMap::iterator> FindPair;
    FindPair = a->lines.equal_range(b->node->nr);

    for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
      // If there is a line with less than two attached triangles, we don't need a new line.
      if (FindLine->second->TrianglesCount < 2) {
        insertNewLine = false;
        cout << Verbose(3) << "Using existing line " << *FindLine->second << endl;

        BPS[0] = FindLine->second->endpoints[0];
        BPS[1] = FindLine->second->endpoints[1];
        BLS[n] = FindLine->second;

        break;
      }
    }
  }

  if (insertNewLine) {
    AlwaysAddTriangleLine(a, b, n);
  }
}
;

/**
 * Adds lines from each of the current points in the BPS to BoundaryLineSet.
 * Raises the line count and inserts the new line into the BLS.
 *
 * @param *a first endpoint
 * @param *b second endpoint
 * @param n index of Tesselation::BLS giving the line with both endpoints
 */
void Tesselation::AlwaysAddTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, int n)
{
  cout << Verbose(3) << "Adding line between " << *(a->node) << " and " << *(b->node) << "." << endl;
  BPS[0] = a;
  BPS[1] = b;
  BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);  // this also adds the line to the local maps
  // add line to global map
  LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
  // increase counter
  LinesOnBoundaryCount++;
};

/** Function tries to add Triangle just created to Triangle and remarks if already existent (Failure of algorithm).
 * Furthermore it adds the triangle to all of its lines, in order to recognize those which are saturated later.
 */
void
Tesselation::AddTriangle()
{
  cout << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl;

  // add triangle to global map
  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
  TrianglesOnBoundaryCount++;

  // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
}
;


double det_get(gsl_matrix *A, int inPlace) {
  /*
  inPlace = 1 => A is replaced with the LU decomposed copy.
  inPlace = 0 => A is retained, and a copy is used for LU.
  */

  double det;
  int signum;
  gsl_permutation *p = gsl_permutation_alloc(A->size1);
  gsl_matrix *tmpA;

  if (inPlace)
  tmpA = A;
  else {
  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
  gsl_matrix_memcpy(tmpA , A);
  }


  gsl_linalg_LU_decomp(tmpA , p , &signum);
  det = gsl_linalg_LU_det(tmpA , signum);
  gsl_permutation_free(p);
  if (! inPlace)
  gsl_matrix_free(tmpA);

  return det;
};

void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
{
  gsl_matrix *A = gsl_matrix_calloc(3,3);
  double m11, m12, m13, m14;

  for(int i=0;i<3;i++) {
    gsl_matrix_set(A, i, 0, a.x[i]);
    gsl_matrix_set(A, i, 1, b.x[i]);
    gsl_matrix_set(A, i, 2, c.x[i]);
  }
  m11 = det_get(A, 1);

  for(int i=0;i<3;i++) {
    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    gsl_matrix_set(A, i, 1, b.x[i]);
    gsl_matrix_set(A, i, 2, c.x[i]);
  }
  m12 = det_get(A, 1);

  for(int i=0;i<3;i++) {
    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    gsl_matrix_set(A, i, 1, a.x[i]);
    gsl_matrix_set(A, i, 2, c.x[i]);
  }
  m13 = det_get(A, 1);

  for(int i=0;i<3;i++) {
    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
    gsl_matrix_set(A, i, 1, a.x[i]);
    gsl_matrix_set(A, i, 2, b.x[i]);
  }
  m14 = det_get(A, 1);

  if (fabs(m11) < MYEPSILON)
    cerr << "ERROR: three points are colinear." << endl;

  center->x[0] =  0.5 * m12/ m11;
  center->x[1] = -0.5 * m13/ m11;
  center->x[2] =  0.5 * m14/ m11;

  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;

  gsl_matrix_free(A);
};



/**
 * Function returns center of sphere with RADIUS, which rests on points a, b, c
 * @param Center this vector will be used for return
 * @param a vector first point of triangle
 * @param b vector second point of triangle
 * @param c vector third point of triangle
 * @param *Umkreismittelpunkt new cneter point of circumference
 * @param Direction vector indicates up/down
 * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
 * @param Halfplaneindicator double indicates whether Direction is up or down
 * @param AlternativeIndicator doube indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
 * @param alpha double angle at a
 * @param beta double, angle at b
 * @param gamma, double, angle at c
 * @param Radius, double
 * @param Umkreisradius double radius of circumscribing circle
 */
void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
{
  Vector TempNormal, helper;
  double Restradius;
  Vector OtherCenter;
  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
  Center->Zero();
  helper.CopyVector(&a);
  helper.Scale(sin(2.*alpha));
  Center->AddVector(&helper);
  helper.CopyVector(&b);
  helper.Scale(sin(2.*beta));
  Center->AddVector(&helper);
  helper.CopyVector(&c);
  helper.Scale(sin(2.*gamma));
  Center->AddVector(&helper);
  //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
  NewUmkreismittelpunkt->CopyVector(Center);
  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
  // Here we calculated center of circumscribing circle, using barycentric coordinates
  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";

  TempNormal.CopyVector(&a);
  TempNormal.SubtractVector(&b);
  helper.CopyVector(&a);
  helper.SubtractVector(&c);
  TempNormal.VectorProduct(&helper);
  if (fabs(HalfplaneIndicator) < MYEPSILON)
    {
      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
        {
          TempNormal.Scale(-1);
        }
    }
  else
    {
      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
        {
          TempNormal.Scale(-1);
        }
    }

  TempNormal.Normalize();
  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
  TempNormal.Scale(Restradius);
  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";

  Center->AddVector(&TempNormal);
  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
  get_sphere(&OtherCenter, a, b, c, RADIUS);
  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
};


/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
 * \param *Center new center on return
 * \param *a first point
 * \param *b second point
 * \param *c third point
 */
void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
{
  Vector helper;
  double alpha, beta, gamma;
  Vector SideA, SideB, SideC;
  SideA.CopyVector(b);
  SideA.SubtractVector(c);
  SideB.CopyVector(c);
  SideB.SubtractVector(a);
  SideC.CopyVector(a);
  SideC.SubtractVector(b);
  alpha = M_PI - SideB.Angle(&SideC);
  beta = M_PI - SideC.Angle(&SideA);
  gamma = M_PI - SideA.Angle(&SideB);
  //cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
    cerr << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;

  Center->Zero();
  helper.CopyVector(a);
  helper.Scale(sin(2.*alpha));
  Center->AddVector(&helper);
  helper.CopyVector(b);
  helper.Scale(sin(2.*beta));
  Center->AddVector(&helper);
  helper.CopyVector(c);
  helper.Scale(sin(2.*gamma));
  Center->AddVector(&helper);
  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
};

/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
 * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
 * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
 * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
 * \param CircleCenter Center of the parameter circle
 * \param CirclePlaneNormal normal vector to plane of the parameter circle
 * \param CircleRadius radius of the parameter circle
 * \param NewSphereCenter new center of a circumcircle
 * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
 * \param NormalVector normal vector
 * \param SearchDirection search direction to make angle unique on return.
 * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
 */
double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
{
  Vector helper;
  double radius, alpha;

  helper.CopyVector(&NewSphereCenter);
  // test whether new center is on the parameter circle's plane
  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
    helper.ProjectOntoPlane(&CirclePlaneNormal);
  }
  radius = helper.ScalarProduct(&helper);
  // test whether the new center vector has length of CircleRadius
  if (fabs(radius - CircleRadius) > HULLEPSILON)
    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
  alpha = helper.Angle(&OldSphereCenter);
  // make the angle unique by checking the halfplanes/search direction
  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    alpha = 2.*M_PI - alpha;
  //cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
  radius = helper.Distance(&OldSphereCenter);
  helper.ProjectOntoPlane(&NormalVector);
  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    //cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
    return alpha;
  } else {
    //cout << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
    return 2.*M_PI;
  }
};


/** Checks whether the triangle consisting of the three atoms is already present.
 * Searches for the points in Tesselation::PointsOnBoundary and checks their
 * lines. If any of the three edges already has two triangles attached, false is
 * returned.
 * \param *out output stream for debugging
 * \param *Candidates endpoints of the triangle candidate
 * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
 *                 triangles exist which is the maximum for three points
 */
int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
  LineMap::iterator FindLine;
  PointMap::iterator FindPoint;
  TriangleMap::iterator FindTriangle;
  int adjacentTriangleCount = 0;
  class BoundaryPointSet *Points[3];

  //*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
  // builds a triangle point set (Points) of the end points
  for (int i = 0; i < 3; i++) {
    FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
    if (FindPoint != PointsOnBoundary.end()) {
      Points[i] = FindPoint->second;
    } else {
      Points[i] = NULL;
    }
  }

  // checks lines between the points in the Points for their adjacent triangles
  for (int i = 0; i < 3; i++) {
    if (Points[i] != NULL) {
      for (int j = i; j < 3; j++) {
        if (Points[j] != NULL) {
          FindLine = Points[i]->lines.find(Points[j]->node->nr);
          if (FindLine != Points[i]->lines.end()) {
            for (; FindLine->first == Points[j]->node->nr; FindLine++) {
              FindTriangle = FindLine->second->triangles.begin();
              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
                if ((
                  (FindTriangle->second->endpoints[0] == Points[0])
                    || (FindTriangle->second->endpoints[0] == Points[1])
                    || (FindTriangle->second->endpoints[0] == Points[2])
                  ) && (
                    (FindTriangle->second->endpoints[1] == Points[0])
                    || (FindTriangle->second->endpoints[1] == Points[1])
                    || (FindTriangle->second->endpoints[1] == Points[2])
                  ) && (
                    (FindTriangle->second->endpoints[2] == Points[0])
                    || (FindTriangle->second->endpoints[2] == Points[1])
                    || (FindTriangle->second->endpoints[2] == Points[2])
                  )
                ) {
                  adjacentTriangleCount++;
                }
              }
            }
            // Only one of the triangle lines must be considered for the triangle count.
            *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
            return adjacentTriangleCount;

          }
        }
      }
    }
  }

  *out << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
  return adjacentTriangleCount;
};

/** This recursive function finds a third point, to form a triangle with two given ones.
 * Note that this function is for the starting triangle.
 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
 * the center of the sphere is still fixed up to a single parameter. The band of possible values
 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
 * us the "null" on this circle, the new center of the candidate point will be some way along this
 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
 * by the normal vector of the base triangle that always points outwards by construction.
 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
 * We construct the normal vector that defines the plane this circle lies in, it is just in the
 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
 * both.
 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
 * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
 * @param BaseLine BoundaryLineSet with the current base line
 * @param ThirdNode third atom to avoid in search
 * @param candidates list of equally good candidates to return
 * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
 * @param RADIUS radius of sphere
 * @param *LC LinkedCell structure with neighbouring atoms
 */
void Find_third_point_for_Tesselation(
  Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter,
  class BoundaryLineSet *BaseLine, atom *ThirdNode, CandidateList* &candidates,
  double *ShortestAngle, const double RADIUS, LinkedCell *LC
) {
  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
  Vector SphereCenter;
  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
  Vector NewNormalVector;   // normal vector of the Candidate's triangle
  Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
  LinkedAtoms *List = NULL;
  double CircleRadius; // radius of this circle
  double radius;
  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
  atom *Candidate = NULL;
  CandidateForTesselation *optCandidate = NULL;

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

  //cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;

  // construct center of circle
  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
  CircleCenter.Scale(0.5);

  // construct normal vector of circle
  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);

  // calculate squared radius atom *ThirdNode,f circle
  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
  if (radius/4. < RADIUS*RADIUS) {
    CircleRadius = RADIUS*RADIUS - radius/4.;
    CirclePlaneNormal.Normalize();
    //cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

    // test whether old center is on the band's plane
    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
    }
    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
    if (fabs(radius - CircleRadius) < HULLEPSILON) {

      // check SearchDirection
      //cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
      }

      // get cell for the starting atom
      if (LC->SetIndexToVector(&CircleCenter)) {
        for(int i=0;i<NDIM;i++) // store indices of this cell
        N[i] = LC->n[i];
        //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
      } else {
        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
        return;
      }
      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
      //cout << Verbose(2) << "LC Intervals:";
      for (int i=0;i<NDIM;i++) {
        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
        //cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
      }
      //cout << endl;
      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
            List = LC->GetCurrentCell();
            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
            if (List != NULL) {
              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
                Candidate = (*Runner);

                // check for three unique points
                //cout << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){

                  // construct both new centers
                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);

                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
                  ) {
                    helper.CopyVector(&NewNormalVector);
                    //cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
                    if (radius < RADIUS*RADIUS) {
                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
                      //cout << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
                      NewSphereCenter.AddVector(&helper);
                      NewSphereCenter.SubtractVector(&CircleCenter);
                      //cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;

                      // OtherNewSphereCenter is created by the same vector just in the other direction
                      helper.Scale(-1.);
                      OtherNewSphereCenter.AddVector(&helper);
                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                      //cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;

                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                      alpha = min(alpha, Otheralpha);
                      // if there is a better candidate, drop the current list and add the new candidate
                      // otherwise ignore the new candidate and keep the list
                      if (*ShortestAngle > (alpha - HULLEPSILON)) {
                        optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter);
                        if (fabs(alpha - Otheralpha) > MYEPSILON) {
                          optCandidate->OptCenter.CopyVector(&NewSphereCenter);
                          optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);
                        } else {
                          optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);
                          optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);
                        }
                        // if there is an equal candidate, add it to the list without clearing the list
                        if ((*ShortestAngle - HULLEPSILON) < alpha) {
                          candidates->push_back(optCandidate);
                          cout << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "
                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
                        } else {
                          // remove all candidates from the list and then the list itself
                          class CandidateForTesselation *remover = NULL;
                          for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) {
                            remover = *it;
                            delete(remover);
                          }
                          candidates->clear();
                          candidates->push_back(optCandidate);
                          cout << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with "
                            << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;
                        }
                        *ShortestAngle = alpha;
                        //cout << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;
                      } else {
                        if ((optCandidate != NULL) && (optCandidate->point != NULL)) {
                          //cout << Verbose(2) << "REJECT: Old candidate: " << *(optCandidate->point) << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
                        } else {
                          //cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
                        }
                      }

                    } else {
                      //cout << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
                    }
                  } else {
                    //cout << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                  }
                } else {
                  if (ThirdNode != NULL) {
                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
                  } else {
                    //cout << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
                  }
                }
              }
            }
          }
    } else {
      cerr << Verbose(2) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    }
  } else {
    if (ThirdNode != NULL)
      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
    else
      cout << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
  }

  //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
  if (candidates->size() > 1) {
    candidates->unique();
    candidates->sort(sortCandidates);
  }

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

struct Intersection {
  Vector x1;
  Vector x2;
  Vector x3;
  Vector x4;
};

/**
 * Intersection calculation function.
 *
 * @param x to find the result for
 * @param function parameter
 */
double MinIntersectDistance(const gsl_vector * x, void *params) {
  double retval = 0;
  struct Intersection *I = (struct Intersection *)params;
  Vector intersection;
  Vector SideA,SideB,HeightA, HeightB;
  for (int i=0;i<NDIM;i++)
    intersection.x[i] = gsl_vector_get(x, i);

  SideA.CopyVector(&(I->x1));
  SideA.SubtractVector(&I->x2);
  HeightA.CopyVector(&intersection);
  HeightA.SubtractVector(&I->x1);
  HeightA.ProjectOntoPlane(&SideA);

  SideB.CopyVector(&I->x3);
  SideB.SubtractVector(&I->x4);
  HeightB.CopyVector(&intersection);
  HeightB.SubtractVector(&I->x3);
  HeightB.ProjectOntoPlane(&SideB);

  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
  //cout << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;

  return retval;
};


/**
 * Calculates whether there is an intersection between two lines. The first line
 * always goes through point 1 and point 2 and the second line is given by the
 * connection between point 4 and point 5.
 *
 * @param point 1 of line 1
 * @param point 2 of line 1
 * @param point 1 of line 2
 * @param point 2 of line 2
 *
 * @return true if there is an intersection between the given lines, false otherwise
 */
bool existsIntersection(Vector point1, Vector point2, Vector point3, Vector point4) {
  bool result;

  struct Intersection par;
    par.x1.CopyVector(&point1);
    par.x2.CopyVector(&point2);
    par.x3.CopyVector(&point3);
    par.x4.CopyVector(&point4);

    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *s = NULL;
    gsl_vector *ss, *x;
    gsl_multimin_function minex_func;

    size_t iter = 0;
    int status;
    double size;

    /* Starting point */
    x = gsl_vector_alloc(NDIM);
    gsl_vector_set(x, 0, point1.x[0]);
    gsl_vector_set(x, 1, point1.x[1]);
    gsl_vector_set(x, 2, point1.x[2]);

    /* Set initial step sizes to 1 */
    ss = gsl_vector_alloc(NDIM);
    gsl_vector_set_all(ss, 1.0);

    /* Initialize method and iterate */
    minex_func.n = NDIM;
    minex_func.f = &MinIntersectDistance;
    minex_func.params = (void *)&par;

    s = gsl_multimin_fminimizer_alloc(T, NDIM);
    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);

    do {
        iter++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status) {
          break;
        }

        size = gsl_multimin_fminimizer_size(s);
        status = gsl_multimin_test_size(size, 1e-2);

        if (status == GSL_SUCCESS) {
          cout << Verbose(2) << "converged to minimum" <<  endl;
        }
    } while (status == GSL_CONTINUE && iter < 100);

    // check whether intersection is in between or not
  Vector intersection, SideA, SideB, HeightA, HeightB;
  double t1, t2;
  for (int i = 0; i < NDIM; i++) {
    intersection.x[i] = gsl_vector_get(s->x, i);
  }

  SideA.CopyVector(&par.x2);
  SideA.SubtractVector(&par.x1);
  HeightA.CopyVector(&intersection);
  HeightA.SubtractVector(&par.x1);

  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);

  SideB.CopyVector(&par.x4);
  SideB.SubtractVector(&par.x3);
  HeightB.CopyVector(&intersection);
  HeightB.SubtractVector(&par.x3);

  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);

  cout << Verbose(2) << "Intersection " << intersection << " is at "
    << t1 << " for (" << point1 << "," << point2 << ") and at "
    << t2 << " for (" << point3 << "," << point4 << "): ";

  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    cout << "true intersection." << endl;
    result = true;
  } else {
    cout << "intersection out of region of interest." << endl;
    result = false;
  }

  // free minimizer stuff
    gsl_vector_free(x);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free(s);

  return result;
}

/** Finds the second point of starting triangle.
 * \param *a first atom
 * \param *Candidate pointer to candidate atom on return
 * \param Oben vector indicating the outside
 * \param Opt_Candidate reference to recommended candidate on return
 * \param Storage[3] array storing angles and other candidate information
 * \param RADIUS radius of virtual sphere
 * \param *LC LinkedCell structure with neighbouring atoms
 */
void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
{
  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
  Vector AngleCheck;
  double norm = -1., angle;
  LinkedAtoms *List = NULL;
  int N[NDIM], Nlower[NDIM], Nupper[NDIM];

  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
    for(int i=0;i<NDIM;i++) // store indices of this cell
      N[i] = LC->n[i];
  } else {
    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
    return;
  }
  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
  cout << Verbose(3) << "LC Intervals from [";
  for (int i=0;i<NDIM;i++) {
  cout << " " << N[i] << "<->" << LC->N[i];
  }
  cout << "] :";
  for (int i=0;i<NDIM;i++) {
    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
  }
  cout << endl;


  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
        List = LC->GetCurrentCell();
        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            Candidate = (*Runner);
            // check if we only have one unique point yet ...
            if (a != Candidate) {
              // Calculate center of the circle with radius RADIUS through points a and Candidate
              Vector OrthogonalizedOben, a_Candidate, Center;
              double distance, scaleFactor;

              OrthogonalizedOben.CopyVector(&Oben);
              a_Candidate.CopyVector(&(a->x));
              a_Candidate.SubtractVector(&(Candidate->x));
              OrthogonalizedOben.ProjectOntoPlane(&a_Candidate);
              OrthogonalizedOben.Normalize();
              distance = 0.5 * a_Candidate.Norm();
              scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
              OrthogonalizedOben.Scale(scaleFactor);

              Center.CopyVector(&(Candidate->x));
              Center.AddVector(&(a->x));
              Center.Scale(0.5);
              Center.AddVector(&OrthogonalizedOben);

              AngleCheck.CopyVector(&Center);
              AngleCheck.SubtractVector(&(a->x));
              norm = a_Candidate.Norm();
              // second point shall have smallest angle with respect to Oben vector
              if (norm < RADIUS*2.) {
                angle = AngleCheck.Angle(&Oben);
                if (angle < Storage[0]) {
                  //cout << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
                  cout << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";
                  Opt_Candidate = Candidate;
                  Storage[0] = angle;
                  //cout << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
                } else {
                  //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
                }
              } else {
                //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
              }
            } else {
              //cout << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
            }
          }
        } else {
          cout << Verbose(3) << "Linked cell list is empty." << endl;
        }
      }
  cout << Verbose(2) << "End of Find_second_point_for_Tesselation" << endl;
};

/** Finds the starting triangle for find_non_convex_border().
 * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
 * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
 * point are called.
 * \param RADIUS radius of virtual rolling sphere
 * \param *LC LinkedCell structure with neighbouring atoms
 */
void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
{
  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
  int i = 0;
  LinkedAtoms *List = NULL;
  atom* FirstPoint = NULL;
  atom* SecondPoint = NULL;
  atom* MaxAtom[NDIM];
  double max_coordinate[NDIM];
  Vector Oben;
  Vector helper;
  Vector Chord;
  Vector SearchDirection;

  Oben.Zero();

  for (i = 0; i < 3; i++) {
    MaxAtom[i] = NULL;
    max_coordinate[i] = -1;
  }

  // 1. searching topmost atom with respect to each axis
  for (int i=0;i<NDIM;i++) { // each axis
    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
        List = LC->GetCurrentCell();
        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
        if (List != NULL) {
          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
            if ((*Runner)->x.x[i] > max_coordinate[i]) {
              cout << Verbose(2) << "New maximal for axis " << i << " atom is " << *(*Runner) << " at " << (*Runner)->x << "." << endl;
              max_coordinate[i] = (*Runner)->x.x[i];
              MaxAtom[i] = (*Runner);
            }
          }
        } else {
          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
        }
      }
  }

  cout << Verbose(2) << "Found maximum coordinates: ";
  for (int i=0;i<NDIM;i++)
    cout << i << ": " << *MaxAtom[i] << "\t";
  cout << endl;

  BTS = NULL;
  CandidateList *Opt_Candidates = new CandidateList();
  for (int k=0;k<NDIM;k++) {
    Oben.x[k] = 1.;
    FirstPoint = MaxAtom[k];
    cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;

    double ShortestAngle;
    atom* Opt_Candidate = NULL;
    ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.

    Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
    SecondPoint = Opt_Candidate;
    if (SecondPoint == NULL)  // have we found a second point?
      continue;
    else
      cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";

    helper.CopyVector(&(FirstPoint->x));
    helper.SubtractVector(&(SecondPoint->x));
    helper.Normalize();
    Oben.ProjectOntoPlane(&helper);
    Oben.Normalize();
    helper.VectorProduct(&Oben);
    ShortestAngle = 2.*M_PI; // This will indicate the quadrant.

    Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
    Chord.SubtractVector(&(SecondPoint->x));
    double radius = Chord.ScalarProduct(&Chord);
    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
    helper.CopyVector(&Oben);
    helper.Scale(CircleRadius);
    // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)

    // look in one direction of baseline for initial candidate
    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...

    // adding point 1 and point 2 and the line between them
    AddTrianglePoint(FirstPoint, 0);
    AddTrianglePoint(SecondPoint, 1);
    AddTriangleLine(TPS[0], TPS[1], 0);

    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n";
    Find_third_point_for_Tesselation(
      Oben, SearchDirection, helper, BLS[0], NULL, *&Opt_Candidates, &ShortestAngle, RADIUS, LC
    );
    cout << Verbose(1) << "List of third Points is ";
    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
        cout << " " << *(*it)->point;
    }
    cout << endl;

    for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
      // add third triangle point
      AddTrianglePoint((*it)->point, 2);
      // add the second and third line
      AddTriangleLine(TPS[1], TPS[2], 1);
      AddTriangleLine(TPS[0], TPS[2], 2);
      // ... and triangles to the Maps of the Tesselation class
      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
      AddTriangle();
      // ... and calculate its normal vector (with correct orientation)
      (*it)->OptCenter.Scale(-1.);
      cout << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl;
      BTS->GetNormalVector((*it)->OptCenter);  // vector to compare with should point inwards
      cout << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and "
      << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n";

      // if we do not reach the end with the next step of iteration, we need to setup a new first line
      if (it != Opt_Candidates->end()--) {
        FirstPoint = (*it)->BaseLine->endpoints[0]->node;
        SecondPoint = (*it)->point;
        // adding point 1 and point 2 and the line between them
        AddTrianglePoint(FirstPoint, 0);
        AddTrianglePoint(SecondPoint, 1);
        AddTriangleLine(TPS[0], TPS[1], 0);
      }
      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
    }
    if (BTS != NULL) // we have created one starting triangle
      break;
    else {
      // remove all candidates from the list and then the list itself
      class CandidateForTesselation *remover = NULL;
      for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
        remover = *it;
        delete(remover);
      }
      Opt_Candidates->clear();
    }
  }

  // remove all candidates from the list and then the list itself
  class CandidateForTesselation *remover = NULL;
  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    remover = *it;
    delete(remover);
  }
  delete(Opt_Candidates);
  cout << Verbose(1) << "End of Find_starting_triangle\n";
};

/** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
 * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
 * make it bigger (i.e. closing one (the baseline) and opening two new ones).
 * \param TPS[3] nodes of the triangle
 * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
 */
bool CheckLineCriteriaforDegeneratedTriangle(class BoundaryPointSet *nodes[3])
{
  bool result = false;
  int counter = 0;

  // check all three points
  for (int i=0;i<3;i++)
    for (int j=i+1; j<3; j++) {
      if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) {  // there already is a line
        LineMap::iterator FindLine;
        pair<LineMap::iterator,LineMap::iterator> FindPair;
        FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
        for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
          // If there is a line with less than two attached triangles, we don't need a new line.
          if (FindLine->second->TrianglesCount < 2) {
            counter++;
            break;  // increase counter only once per edge
          }
        }
      } else { // no line
        cout << Verbose(1) << "ERROR: The line between " << nodes[i] << " and " << nodes[j] << " is not yet present, hence no need for a degenerate triangle!" << endl;
        result = true;
      }
    }
  if (counter > 1) {
    cout << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    result = true;
  }
  return result;
};


/** This function finds a triangle to a line, adjacent to an existing one.
 * @param out output stream for debugging
 * @param *mol molecule with Atom's and Bond's
 * @param Line current baseline to search from
 * @param T current triangle which \a Line is edge of
 * @param RADIUS radius of the rolling ball
 * @param N number of found triangles
 * @param *filename filename base for intermediate envelopes
 * @param *LC LinkedCell structure with neighbouring atoms
 */
bool Tesselation::Find_next_suitable_triangle(ofstream *out,
    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
{
  cout << Verbose(0) << "Begin of Find_next_suitable_triangle\n";
  ofstream *tempstream = NULL;
  char NumberName[255];
  bool result = true;
  CandidateList *Opt_Candidates = new CandidateList();

  Vector CircleCenter;
  Vector CirclePlaneNormal;
  Vector OldSphereCenter;
  Vector SearchDirection;
  Vector helper;
  atom *ThirdNode = NULL;
  LineMap::iterator testline;
  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
  double radius, CircleRadius;

  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
  for (int i=0;i<3;i++)
    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
      ThirdNode = T.endpoints[i]->node;

  // construct center of circle
  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
  CircleCenter.Scale(0.5);

  // construct normal vector of circle
  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);

  // calculate squared radius of circle
  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
  if (radius/4. < RADIUS*RADIUS) {
    CircleRadius = RADIUS*RADIUS - radius/4.;
    CirclePlaneNormal.Normalize();
    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;

    // construct old center
    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
    helper.Scale(sqrt(RADIUS*RADIUS - radius));
    OldSphereCenter.AddVector(&helper);
    OldSphereCenter.SubtractVector(&CircleCenter);
    //cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;

    // construct SearchDirection
    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
    helper.CopyVector(&Line.endpoints[0]->node->x);
    helper.SubtractVector(&ThirdNode->x);
    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
      SearchDirection.Scale(-1.);
    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
    SearchDirection.Normalize();
    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
      // rotated the wrong way!
      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
    }

    // add third point
    Find_third_point_for_Tesselation(
      T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidates,
      &ShortestAngle, RADIUS, LC
    );

  } else {
    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
  }

  if (Opt_Candidates->begin() == Opt_Candidates->end()) {
    cerr << "WARNING: Could not find a suitable candidate." << endl;
    return false;
  }
  cout << Verbose(1) << "Third Points are ";
  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    cout << " " << *(*it)->point;
  }
  cout << endl;

  BoundaryLineSet *BaseRay = &Line;
  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    cout << Verbose(1) << " Third point candidate is " << *(*it)->point
    << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
    cout << Verbose(1) << " Baseline is " << *BaseRay << endl;

    // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
    atom *AtomCandidates[3];
    AtomCandidates[0] = (*it)->point;
    AtomCandidates[1] = BaseRay->endpoints[0]->node;
    AtomCandidates[2] = BaseRay->endpoints[1]->node;
    int existentTrianglesCount = CheckPresenceOfTriangle(out, AtomCandidates);

    BTS = NULL;
    // If there is no triangle, add it regularly.
    if (existentTrianglesCount == 0) {
      AddTrianglePoint((*it)->point, 0);
      AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
      AddTrianglePoint(BaseRay->endpoints[1]->node, 2);

      AddTriangleLine(TPS[0], TPS[1], 0);
      AddTriangleLine(TPS[0], TPS[2], 1);
      AddTriangleLine(TPS[1], TPS[2], 2);

      BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
      AddTriangle();
      (*it)->OptCenter.Scale(-1.);
      BTS->GetNormalVector((*it)->OptCenter);
      (*it)->OptCenter.Scale(-1.);

      cout << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector
        << " for this triangle ... " << endl;
      //cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl;
    } else if (existentTrianglesCount == 1) { // If there is a planar region within the structure, we need this triangle a second time.
        AddTrianglePoint((*it)->point, 0);
        AddTrianglePoint(BaseRay->endpoints[0]->node, 1);
        AddTrianglePoint(BaseRay->endpoints[1]->node, 2);

        // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
        // i.e. at least one of the three lines must be present with TriangleCount <= 1
        if (CheckLineCriteriaforDegeneratedTriangle(TPS)) {
          AddTriangleLine(TPS[0], TPS[1], 0);
          AddTriangleLine(TPS[0], TPS[2], 1);
          AddTriangleLine(TPS[1], TPS[2], 2);

          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
          AddTriangle();

          (*it)->OtherOptCenter.Scale(-1.);
          BTS->GetNormalVector((*it)->OtherOptCenter);
          (*it)->OtherOptCenter.Scale(-1.);

          cout << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector
          << " for this triangle ... " << endl;
          cout << Verbose(1) << "We have "<< BaseRay->TrianglesCount << " for line " << BaseRay << "." << endl;
        } else {
          cout << Verbose(1) << "WARNING: This triangle consisting of ";
          cout << *(*it)->point << ", ";
          cout << *BaseRay->endpoints[0]->node << " and ";
          cout << *BaseRay->endpoints[1]->node << " ";
          cout << "exists and is not added, as it does not seem helpful!" << endl;
          result = false;
        }
    } else {
      cout << Verbose(1) << "This triangle consisting of ";
      cout << *(*it)->point << ", ";
      cout << *BaseRay->endpoints[0]->node << " and ";
      cout << *BaseRay->endpoints[1]->node << " ";
      cout << "is invalid!" << endl;
      result = false;
    }

    if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
      sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
      if (DoTecplotOutput) {
        string NameofTempFile(tempbasename);
        NameofTempFile.append(NumberName);
        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
        NameofTempFile.erase(npos, 1);
        NameofTempFile.append(TecplotSuffix);
        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
        write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
        tempstream->close();
        tempstream->flush();
        delete(tempstream);
      }

      if (DoRaster3DOutput) {
        string NameofTempFile(tempbasename);
        NameofTempFile.append(NumberName);
        for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
        NameofTempFile.erase(npos, 1);
        NameofTempFile.append(Raster3DSuffix);
        cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
        tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
        write_raster3d_file(out, tempstream, this, mol);
        // include the current position of the virtual sphere in the temporary raster3d file
        // make the circumsphere's center absolute again
        helper.CopyVector(&BaseRay->endpoints[0]->node->x);
        helper.AddVector(&BaseRay->endpoints[1]->node->x);
        helper.Scale(0.5);
        (*it)->OptCenter.AddVector(&helper);
        Vector *center = mol->DetermineCenterOfAll(out);
        (*it)->OptCenter.AddVector(center);
        delete(center);
        // and add to file plus translucency object
        *tempstream << "# current virtual sphere\n";
        *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
        *tempstream << "2\n  " << (*it)->OptCenter.x[0] << " "
          << (*it)->OptCenter.x[1] << " " << (*it)->OptCenter.x[2]
          << "\t" << RADIUS << "\t1 0 0\n";
        *tempstream << "9\n  terminating special property\n";
        tempstream->close();
        tempstream->flush();
        delete(tempstream);
      }
      if (DoTecplotOutput || DoRaster3DOutput)
        TriangleFilesWritten++;
    }

    // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point))
    BaseRay = BLS[0];
  }

  // remove all candidates from the list and then the list itself
  class CandidateForTesselation *remover = NULL;
  for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
    remover = *it;
    delete(remover);
  }
  delete(Opt_Candidates);
  cout << Verbose(0) << "End of Find_next_suitable_triangle\n";
  return result;
};

/**
 * Sort function for the candidate list.
 */
bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
  Vector BaseLineVector, OrthogonalVector, helper;
  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    //return false;
    exit(1);
  }
  // create baseline vector
  BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
  BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
  BaseLineVector.Normalize();

  // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
  helper.CopyVector(&(candidate1->BaseLine->endpoints[0]->node->x));
  helper.SubtractVector(&(candidate1->point->x));
  OrthogonalVector.CopyVector(&helper);
  helper.VectorProduct(&BaseLineVector);
  OrthogonalVector.SubtractVector(&helper);
  OrthogonalVector.Normalize();

  // calculate both angles and correct with in-plane vector
  helper.CopyVector(&(candidate1->point->x));
  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
  double phi = BaseLineVector.Angle(&helper);
  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    phi = 2.*M_PI - phi;
  }
  helper.CopyVector(&(candidate2->point->x));
  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
  double psi = BaseLineVector.Angle(&helper);
  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    psi = 2.*M_PI - psi;
  }

  cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
  cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;

  // return comparison
  return phi < psi;
}

/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
 * \param *out output stream for debugging
 * \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 (1) { //failflag) {
    *out << Verbose(1) << "Writing final tecplot file\n";
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      write_tecplot_file(out, tecplot, Tess, mol, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(Raster3DSuffix);
      ofstream *raster = new ofstream(OutputName.c_str());
      write_raster3d_file(out, raster, Tess, mol);
      raster->close();
      delete(raster);
    }
  } else {
    cerr << "ERROR: Could 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;

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