/** \file boundary.cpp
 *
 * Implementations and super-function for envelopes
 */

#include "atom.hpp"
#include "bond.hpp"
#include "boundary.hpp"
#include "config.hpp"
#include "element.hpp"
#include "helpers.hpp"
#include "info.hpp"
#include "linkedcell.hpp"
#include "log.hpp"
#include "memoryallocator.hpp"
#include "molecule.hpp"
#include "tesselation.hpp"
#include "tesselationhelpers.hpp"
#include "World.hpp"

#include<gsl/gsl_poly.h>
#include<time.h>

// ========================================== F U N C T I O N S =================================


/** 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 *&TesselStruct Tesselation structure with triangles
 * \param IsAngstroem whether we have angstroem or atomic units
 * \return NDIM array of the diameters
 */
double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem)
{
	Info FunctionInfo(__func__);
  // get points on boundary of NULL was given as parameter
  bool BoundaryFreeFlag = false;
  double OldComponent = 0.;
  double tmp = 0.;
  double w1 = 0.;
  double w2 = 0.;
  Vector DistanceVector;
  Vector OtherVector;
  int component = 0;
  int Othercomponent = 0;
  Boundaries::const_iterator Neighbour;
  Boundaries::const_iterator OtherNeighbour;
  double *GreatestDiameter = new double[NDIM];

  const Boundaries *BoundaryPoints;
  if (BoundaryPtr == NULL) {
    BoundaryFreeFlag = true;
    BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  } else {
    BoundaryPoints = BoundaryPtr;
    Log() << Verbose(0) << "Using given boundary points set." << endl;
  }
  // determine biggest "diameter" of cluster for each axis
  for (int i = 0; i < NDIM; i++)
    GreatestDiameter[i] = 0.;
  for (int axis = 0; axis < NDIM; axis++)
    { // regard each projected plane
      //Log() << 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;
          //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl;
          for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
              //Log() << Verbose(1) << "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);
                  //Log() << Verbose(2) << "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--;
                  //Log() << Verbose(1) << "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);
                  //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                  //Log() << Verbose(1) << "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
                  //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
                  GreatestDiameter[component] = (GreatestDiameter[component]
                      > tmp) ? GreatestDiameter[component] : tmp;
                } //else
              //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl;
            }
        }
    }
  Log() << 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;
}
;


/** 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
 * \param *&TesselStruct pointer to Tesselation structure
 */
Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct)
{
	Info FunctionInfo(__func__);
  atom *Walker = NULL;
  PointMap PointsOnBoundary;
  LineMap LinesOnBoundary;
  TriangleMap TrianglesOnBoundary;
  Vector *MolCenter = mol->DetermineCenterOfAll();
  Vector helper;
  BoundariesTestPair BoundaryTestPair;
  Vector AxisVector;
  Vector AngleReferenceVector;
  Vector AngleReferenceNormalVector;
  Vector ProjectedVector;
  Boundaries *BoundaryPoints = new Boundaries[NDIM]; // first is alpha, second is (r, nr)
  double angle = 0.;

  // 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.;

    Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;

    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    Walker = mol->start;
    while (Walker->next != mol->end) {
      Walker = Walker->next;
      ProjectedVector.CopyVector(&Walker->x);
      ProjectedVector.SubtractVector(MolCenter);
      ProjectedVector.ProjectOntoPlane(&AxisVector);

      // correct for negative side
      const double radius = ProjectedVector.NormSquared();
      if (fabs(radius) > MYEPSILON)
        angle = ProjectedVector.Angle(&AngleReferenceVector);
      else
        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues

      //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
      if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
        angle = 2. * M_PI - angle;
      }
      Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
        Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
        Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
        Log() << Verbose(2) << "New vector: " << *Walker << endl;
        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
          BoundaryTestPair.first->second.second = Walker;
          Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
          helper.CopyVector(&Walker->x);
          helper.SubtractVector(MolCenter);
          const double oldhelperNorm = helper.NormSquared();
          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
          helper.SubtractVector(MolCenter);
          if (helper.NormSquared() < oldhelperNorm) {
            BoundaryTestPair.first->second.second = Walker;
            Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
          } else {
            Log() << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << oldhelperNorm << "." << endl;
          }
        } else {
          Log() << Verbose(2) << "Keeping present vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
        }
      }
    }
    // printing all inserted for debugging
    //    {
    //      Log() << Verbose(1) << "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())
    //          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
    //        else
    //          Log() << Verbose(0) << i << ": " << *runner->second.second;
    //        i++;
    //      }
    //      Log() << Verbose(0) << endl;
    //    }
    // 3c. throw out points whose distance is less than the mean of left and right neighbours
    bool flag = false;
    Log() << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
    do { // do as long as we still throw one out per round
      flag = false;
      Boundaries::iterator left = BoundaryPoints[axis].end();
      Boundaries::iterator right = BoundaryPoints[axis].end();
      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
        // set neighbours correctly
        if (runner == BoundaryPoints[axis].begin()) {
          left = BoundaryPoints[axis].end();
        } else {
          left = runner;
        }
        left--;
        right = runner;
        right++;
        if (right == BoundaryPoints[axis].end()) {
          right = BoundaryPoints[axis].begin();
        }
        // check distance

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

          SideB.CopyVector(&right->second.second->x);
          SideB.SubtractVector(MolCenter);
          SideB.ProjectOntoPlane(&AxisVector);
          //          Log() << Verbose(1) << "SideB: " << SideB << endl;

          SideC.CopyVector(&left->second.second->x);
          SideC.SubtractVector(&right->second.second->x);
          SideC.ProjectOntoPlane(&AxisVector);
          //          Log() << Verbose(1) << "SideC: " << SideC << endl;

          SideH.CopyVector(&runner->second.second->x);
          SideH.SubtractVector(MolCenter);
          SideH.ProjectOntoPlane(&AxisVector);
          //          Log() << Verbose(1) << "SideH: " << SideH << endl;

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

/** Tesselates the convex boundary by finding all boundary points.
 * \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 FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename)
{
	Info FunctionInfo(__func__);
  bool BoundaryFreeFlag = false;
  Boundaries *BoundaryPoints = NULL;

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

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

// printing all inserted for debugging
  for (int axis=0; axis < NDIM; axis++)
    {
      Log() << Verbose(1) << "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())
          Log() << Verbose(0) << ", " << i << ": " << *runner->second.second;
        else
          Log() << Verbose(0) << i << ": " << *runner->second.second;
        i++;
      }
      Log() << Verbose(0) << endl;
    }

  // 2. fill the boundary point list
  for (int axis = 0; axis < NDIM; axis++)
    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++)
        if (!TesselStruct->AddBoundaryPoint(runner->second.second, 0))
          DoeLog(2) && (eLog()<< Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl);

  Log() << Verbose(0) << "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
  //  Log() << Verbose(1) << "Listing PointsOnBoundary:";
  //  for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
  //    Log() << Verbose(0) << " " << *runner->second;
  //  }
  //  Log() << Verbose(0) << endl;

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

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

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

  Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;

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

  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
  bool AllConvex = true;
  class BoundaryLineSet *line = NULL;
  do {
    AllConvex = true;
    for (LineMap::iterator LineRunner = TesselStruct->LinesOnBoundary.begin(); LineRunner != TesselStruct->LinesOnBoundary.end(); LineRunner++) {
      line = LineRunner->second;
      Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
      if (!line->CheckConvexityCriterion()) {
        Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;

        // flip the line
        if (TesselStruct->PickFarthestofTwoBaselines(line) == 0.)
          DoeLog(1) && (eLog()<< Verbose(1) << "Correction of concave baselines failed!" << endl);
        else {
          TesselStruct->FlipBaseline(line);
          Log() << Verbose(1) << "INFO: Correction of concave baselines worked." << endl;
        }
      }
    }
  } while (!AllConvex);

  // 3e. we need another correction here, for TesselPoints that are below the surface (i.e. have an odd number of concave triangles surrounding it)
//  if (!TesselStruct->CorrectConcaveTesselPoints(out))
//    Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl;

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

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


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

/** For testing removes one boundary point after another to check for leaks.
 * \param *out output stream for debugging
 * \param *TesselStruct Tesselation containing envelope with boundary points
 * \param *mol molecule
 * \param *filename name of file
 * \return true - all removed, false - something went wrong
 */
bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
{
	Info FunctionInfo(__func__);
  int i=0;
  char number[MAXSTRINGSIZE];

  if ((TesselStruct == NULL) || (TesselStruct->PointsOnBoundary.empty())) {
    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty." << endl);
    return false;
  }

  PointMap::iterator PointRunner;
  while (!TesselStruct->PointsOnBoundary.empty()) {
    Log() << Verbose(1) << "Remaining points are: ";
    for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++)
      Log() << Verbose(0) << *(PointSprinter->second) << "\t";
    Log() << Verbose(0) << endl;

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    // remove point
    TesselStruct->RemovePointFromTesselatedSurface(PointRunner->second);

    // store envelope
    sprintf(number, "-%04d", i++);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, number);
  }

  return true;
};

/** Creates a convex envelope from a given non-convex one.
 * -# First step, remove concave spots, i.e. singular "dents"
 *   -# We go through all PointsOnBoundary.
 *   -# We CheckConvexityCriterion() for all its lines.
 *   -# If all its lines are concave, it cannot be on the convex envelope.
 *   -# Hence, we remove it and re-create all its triangles from its getCircleOfConnectedPoints()
 *   -# We calculate the additional volume.
 *   -# We go over all lines until none yields a concavity anymore.
 * -# Second step, remove concave lines, i.e. line-shape "dents"
 *   -# We go through all LinesOnBoundary
 *   -# We CheckConvexityCriterion()
 *   -# If it returns concave, we flip the line in this quadrupel of points (abusing the degeneracy of the tesselation)
 *   -# We CheckConvexityCriterion(),
 *   -# if it's concave, we continue
 *   -# if not, we mark an error and stop
 * Note: This routine - for free - calculates the difference in volume between convex and
 * non-convex envelope, as the former is easy to calculate - VolumeOfConvexEnvelope() - it
 * can be used to compute volumes of arbitrary shapes.
 * \param *out output stream for debugging
 * \param *TesselStruct non-convex envelope, is changed in return!
 * \param *mol molecule
 * \param *filename name of file
 * \return volume difference between the non- and the created convex envelope
 */
double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
{
	Info FunctionInfo(__func__);
  double volume = 0;
  class BoundaryPointSet *point = NULL;
  class BoundaryLineSet *line = NULL;
  bool Concavity = false;
  char dummy[MAXSTRINGSIZE];
  PointMap::iterator PointRunner;
  PointMap::iterator PointAdvance;
  LineMap::iterator LineRunner;
  LineMap::iterator LineAdvance;
  TriangleMap::iterator TriangleRunner;
  TriangleMap::iterator TriangleAdvance;
  int run = 0;

  // check whether there is something to work on
  if (TesselStruct == NULL) {
    DoeLog(1) && (eLog()<< Verbose(1) << "TesselStruct is empty!" << endl);
    return volume;
  }

  // First step: RemovePointFromTesselatedSurface
  do {
    Concavity = false;
    sprintf(dummy, "-first-%d", run);
    //CalculateConcavityPerBoundaryPoint(TesselStruct);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);

    PointRunner = TesselStruct->PointsOnBoundary.begin();
    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
      PointAdvance++;
      point = PointRunner->second;
      Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
        line = LineRunner->second;
        Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
        if (!line->CheckConvexityCriterion()) {
          // remove the point if needed
          Log() << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
          sprintf(dummy, "-first-%d", ++run);
          StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
          Concavity = true;
          break;
        }
      }
      PointRunner = PointAdvance;
    }

    sprintf(dummy, "-second-%d", run);
    //CalculateConcavityPerBoundaryPoint(TesselStruct);
    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);

    // second step: PickFarthestofTwoBaselines
    LineRunner = TesselStruct->LinesOnBoundary.begin();
    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
      LineAdvance++;
      line = LineRunner->second;
      Log() << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
      // take highest of both lines
      if (TesselStruct->IsConvexRectangle(line) == NULL) {
        const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
        volume += tmp;
        if (tmp != 0.) {
          TesselStruct->FlipBaseline(line);
          Concavity = true;
        }
      }
      LineRunner = LineAdvance;
    }
    run++;
  } while (Concavity);
  //CalculateConcavityPerBoundaryPoint(TesselStruct);
  //StoreTrianglesinFile(mol, filename, "-third");

  // third step: IsConvexRectangle
//  LineRunner = TesselStruct->LinesOnBoundary.begin();
//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
//    LineAdvance++;
//    line = LineRunner->second;
//    Log() << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
//      //Log() << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
//    if (!line->CheckConvexityCriterion(out)) {
//      Log() << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
//
//      // take highest of both lines
//      point = TesselStruct->IsConvexRectangle(line);
//      if (point != NULL)
//        volume += TesselStruct->RemovePointFromTesselatedSurface(point);
//    }
//    LineRunner = LineAdvance;
//  }

  CalculateConcavityPerBoundaryPoint(TesselStruct);
  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");

  // end
  Log() << Verbose(0) << "Volume is " << volume << "." << endl;
  return volume;
};


/** 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(class Tesselation *TesselStruct, class config *configuration)
{
	Info FunctionInfo(__func__);
  bool IsAngstroem = configuration->GetIsAngstroem();
  double volume = 0.;
  Vector x;
  Vector y;

  // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
  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->node);
      x.SubtractVector(runner->second->endpoints[1]->node->node);
      y.CopyVector(runner->second->endpoints[0]->node->node);
      y.SubtractVector(runner->second->endpoints[2]->node->node);
      const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
      const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
      const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
      const double 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->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
      const double h = x.Norm(); // distance of CoG to triangle
      const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
      Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is "
          << h << " and the volume is " << PyramidVolume << " "
          << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
      volume += PyramidVolume;
    }
  Log() << Verbose(0) << "RESULT: The summed volume is " << setprecision(6)
      << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."
      << endl;

  return volume;
};

/** Stores triangles to file.
 * \param *out output stream for debugging
 * \param *mol molecule with atoms and bonds
 * \param *TesselStruct Tesselation with boundary triangles
 * \param *filename prefix of filename
 * \param *extraSuffix intermediate suffix
 */
void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
{
	Info FunctionInfo(__func__);
  // 4. Store triangles in tecplot file
  if (filename != NULL) {
    if (DoTecplotOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(TecplotSuffix);
      ofstream *tecplot = new ofstream(OutputName.c_str());
      WriteTecplotFile(tecplot, TesselStruct, mol, -1);
      tecplot->close();
      delete(tecplot);
    }
    if (DoRaster3DOutput) {
      string OutputName(filename);
      OutputName.append(extraSuffix);
      OutputName.append(Raster3DSuffix);
      ofstream *rasterplot = new ofstream(OutputName.c_str());
      WriteRaster3dFile(rasterplot, TesselStruct, mol);
      rasterplot->close();
      delete(rasterplot);
    }
  }
};

/** 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 *&TesselStruct Tesselation structure with triangles on return
 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used VolumeOfConvexEnvelope() instead.
 * \param celldensity desired average density in final cell
 */
void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
{
	Info FunctionInfo(__func__);
  bool IsAngstroem = true;
  double *GreatestDiameter = NULL;
  Boundaries *BoundaryPoints = NULL;
  class Tesselation *TesselStruct = NULL;
  Vector BoxLengths;
  int repetition[NDIM] = { 1, 1, 1 };
  int TotalNoClusters = 1;
  atom *Walker = NULL;
  double totalmass = 0.;
  double clustervolume = 0.;
  double cellvolume = 0.;

  // transform to PAS
  mol->PrincipalAxisSystem(true);

  IsAngstroem = configuration->GetIsAngstroem();
  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
  LinkedCell LCList(mol, 10.);
  FindConvexBorder(mol, TesselStruct, &LCList, NULL);

  // some preparations beforehand
  if (ClusterVolume == 0)
    clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
  else
    clustervolume = ClusterVolume;

  for (int i = 0; i < NDIM; i++)
    TotalNoClusters *= repetition[i];

  // sum up the atomic masses
  Walker = mol->start;
  while (Walker->next != mol->end) {
      Walker = Walker->next;
      totalmass += Walker->type->mass;
  }
  Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
  Log() << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;

  // solve cubic polynomial
  Log() << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
  if (IsAngstroem)
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
  else
    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
  Log() << 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]);
  Log() << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
  if (minimumvolume > cellvolume) {
    DoeLog(1) && (eLog()<< Verbose(1) << "the containing box already has a greater volume than the envisaged cell volume!" << endl);
    Log() << 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(&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.;
    double x1 = 0.;
    double 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
      Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
    else {
      Log() << 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
    Log() << Verbose(0) << "Translating to box with these boundaries." << endl;
    mol->SetBoxDimension(&BoxLengths);
    mol->CenterInBox();
  }
  // update Box of atoms by boundary
  mol->SetBoxDimension(&BoxLengths);
  Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
};


/** Fills the empty space of the simulation box with water/
 * \param *out output stream for debugging
 * \param *List list of molecules already present in box
 * \param *TesselStruct contains tesselated surface
 * \param *filler molecule which the box is to be filled with
 * \param configuration contains box dimensions
 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
 * \param distance[NDIM] distance between filling molecules in each direction
 * \param boundary length of boundary zone between molecule and filling mollecules
 * \param epsilon distance to surface which is not filled
 * \param RandAtomDisplacement maximum distance for random displacement per atom
 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
 * \param DoRandomRotation true - do random rotiations, false - don't
 * \return *mol pointer to new molecule with filled atoms
 */
molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
{
	Info FunctionInfo(__func__);
  molecule *Filling = new molecule(filler->elemente);
  Vector CurrentPosition;
  int N[NDIM];
  int n[NDIM];
  double *M =  ReturnFullMatrixforSymmetric(World::get()->cell_size);
  double Rotations[NDIM*NDIM];
  double *MInverse = InverseMatrix(M);
  Vector AtomTranslations;
  Vector FillerTranslations;
  Vector FillerDistance;
  Vector Inserter;
  double FillIt = false;
  atom *Walker = NULL;
  bond *Binder = NULL;
  double phi[NDIM];
  map<molecule *, Tesselation *> TesselStruct;
  map<molecule *, LinkedCell *> LCList;

  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
    if ((*ListRunner)->AtomCount > 0) {
      Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
      Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
      TesselStruct[(*ListRunner)] = NULL;
      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
    }

  // Center filler at origin
  filler->CenterEdge(&Inserter);
  filler->Center.Zero();
  Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl;
  Binder = filler->first;
  while(Binder->next != filler->last) {
    Binder = Binder->next;
    Log() << Verbose(2) << "  " << *Binder << endl;
  }

  filler->CountAtoms();
  atom * CopyAtoms[filler->AtomCount];

  // calculate filler grid in [0,1]^3
  FillerDistance.Init(distance[0], distance[1], distance[2]);
  FillerDistance.InverseMatrixMultiplication(M);
  for(int i=0;i<NDIM;i++)
    N[i] = (int) ceil(1./FillerDistance.x[i]);
  Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;

  // initialize seed of random number generator to current time
  srand ( time(NULL) );

  // go over [0,1]^3 filler grid
  for (n[0] = 0; n[0] < N[0]; n[0]++)
    for (n[1] = 0; n[1] < N[1]; n[1]++)
      for (n[2] = 0; n[2] < N[2]; n[2]++) {
        // calculate position of current grid vector in untransformed box
        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
        CurrentPosition.MatrixMultiplication(M);
        // create molecule random translation vector ...
        for (int i=0;i<NDIM;i++)
          FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
        Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl;

        // go through all atoms
        for (int i=0;i<filler->AtomCount;i++)
          CopyAtoms[i] = NULL;
        Walker = filler->start;
        while (Walker->next != filler->end) {
          Walker = Walker->next;

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

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

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

          // ... and put at new position
          Inserter.CopyVector(&(Walker->x));
          if (DoRandomRotation)
            Inserter.MatrixMultiplication(Rotations);
          Inserter.AddVector(&AtomTranslations);
          Inserter.AddVector(&FillerTranslations);
          Inserter.AddVector(&CurrentPosition);

          // check whether inserter is inside box
          Inserter.MatrixMultiplication(MInverse);
          FillIt = true;
          for (int i=0;i<NDIM;i++)
            FillIt = FillIt && (Inserter.x[i] >= -MYEPSILON) && ((Inserter.x[i]-1.) <= MYEPSILON);
          Inserter.MatrixMultiplication(M);

          // Check whether point is in- or outside
          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
            // get linked cell list
            if (TesselStruct[(*ListRunner)] != NULL) {
              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
            }
          }
          // insert into Filling
          if (FillIt) {
            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl;
            // copy atom ...
            CopyAtoms[Walker->nr] = new atom(Walker);
            CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
            Filling->AddAtom(CopyAtoms[Walker->nr]);
            Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
          } else {
            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl;
            CopyAtoms[Walker->nr] = NULL;
            continue;
          }
        }
        // go through all bonds and add as well
        Binder = filler->first;
        while(Binder->next != filler->last) {
          Binder = Binder->next;
          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
          }
        }
      }
  Free(&M);
  Free(&MInverse);

  return Filling;
};


/** 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 *&TesselStruct Tesselation filled with points, lines and triangles on boundary on return
 * \param *&LCList atoms in LinkedCell list
 * \param RADIUS radius of the virtual sphere
 * \param *filename filename prefix for output of vertex data
 * \return true - tesselation successful, false - tesselation failed
 */
bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL)
{
	Info FunctionInfo(__func__);
  bool freeLC = false;
  bool status = false;
  CandidateForTesselation *baseline = NULL;
  bool OneLoopWithoutSuccessFlag = true;  // marks whether we went once through all baselines without finding any without two triangles
  bool TesselationFailFlag = false;

  if (TesselStruct == NULL) {
    Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl;
    TesselStruct= new Tesselation;
  } else {
    delete(TesselStruct);
    Log() << Verbose(1) << "Re-Allocating Tesselation struct ..." << endl;
    TesselStruct = new Tesselation;
  }

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

  // 1. get starting triangle
  TesselStruct->FindStartingTriangle(RADIUS, LCList);
  if (filename != NULL) {
    if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
      TesselStruct->Output(filename, mol);
    }
  }

  // 2. expand from there
  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
    // 2a. print OpenLines
    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      Log() << Verbose(1) << " " << *(Runner->second) << endl;

    // 2b. find best candidate for each OpenLine
    TesselationFailFlag = TesselStruct->FindCandidatesforOpenLines(RADIUS, LCList);

    // 2c. print OpenLines with candidates again
    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
      Log() << Verbose(1) << " " << *(Runner->second) << endl;

    // 2d. search for smallest ShortestAngle among all candidates
    double ShortestAngle = 4.*M_PI;
    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
      if (Runner->second->ShortestAngle < ShortestAngle) {
        baseline = Runner->second;
        ShortestAngle = baseline->ShortestAngle;
        Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *(*baseline->pointlist.begin()) << " and angle " << baseline->ShortestAngle << endl;
      }
    }
    // 2e. if we found one, add candidate
    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
      OneLoopWithoutSuccessFlag = false;
    else {
      TesselStruct->AddCandidatePolygon(*baseline, RADIUS, LCList);
    }

    // 2f. write temporary envelope
    if (filename != NULL) {
      if ((DoSingleStepOutput && ((TesselStruct->TrianglesOnBoundary.size() % SingleStepWidth == 0)))) { // if we have a new triangle and want to output each new triangle configuration
        TesselStruct->Output(filename, mol);
      }
    }
  }
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);
//
//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
//  //->InsertStraddlingPoints(mol, LCList);
//  mol->GoToFirst();
//  class TesselPoint *Runner = NULL;
//  while (!mol->IsEnd()) {
//    Runner = mol->GetPoint();
//    Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
//    if (!->IsInnerPoint(Runner, LCList)) {
//      Log() << Verbose(2) << Runner->Name << " is outside of envelope, adding via degenerated triangles." << endl;
//      ->AddBoundaryPointByDegeneratedTriangle(Runner, LCList);
//    } else {
//      Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
//    }
//    mol->GoToNext();
//  }

//  // Purges surplus triangles.
//  TesselStruct->RemoveDegeneratedTriangles();

  // check envelope for consistency
  status = CheckListOfBaselines(TesselStruct);

  // store before correction
  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");

//  // correct degenerated polygons
//  TesselStruct->CorrectAllDegeneratedPolygons();
//
//  // check envelope for consistency
//  status = CheckListOfBaselines(TesselStruct);

  // write final envelope
  CalculateConcavityPerBoundaryPoint(TesselStruct);
  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");

  if (freeLC)
    delete(LCList);

  return status;
};


/** Finds a hole of sufficient size in \a *mols to embed \a *srcmol into it.
 * \param *out output stream for debugging
 * \param *mols molecules in the domain to embed in between
 * \param *srcmol embedding molecule
 * \return *Vector new center of \a *srcmol for embedding relative to \a this
 */
Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol)
{
	Info FunctionInfo(__func__);
  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;
};