Changes in src/boundary.cpp [8c54a3:2319ed]

File:

• src/boundary.cpp (modified) (44 diffs)

src/boundary.cpp

@@ -10 +10 @@
 #define DoSingleStepOutput 0
 #define DoTecplotOutput 1
-#define DoRaster3DOutput 0
+#define DoRaster3DOutput 1
 #define DoVRMLOutput 1
 #define TecplotSuffix ".dat"
@@ -45 +45 @@
 void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
-  cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
+  cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
+      << endl;
   if (line->endpoints[0] == this)
     {
@@ -131 +132 @@
 ;
 
+/** Checks whether we have a common endpoint with given \a *line.
+ * \param *line other line to test
+ * \return true - common endpoint present, false - not connected
+ */
+bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+{
+  if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1]))
+    return true;
+  else
+    return false;
+};
+
+/** Checks whether the adjacent triangles of a baseline are convex or not.
+ * We sum the two angles of each normal vector with a ficticious normnal vector from this baselinbe pointing outwards.
+ * If greater/equal M_PI than we are convex.
+ * \param *out output stream for debugging
+ * \return true - triangles are convex, false - concave or less than two triangles connected
+ */
+bool BoundaryLineSet::CheckConvexityCriterion(ofstream *out)
+{
+  Vector BaseLineNormal;
+  double angle = 0;
+  // get the two triangles
+  if (TrianglesCount != 2) {
+    *out << Verbose(1) << "ERROR: Baseline " << this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
+    return false;
+  }
+  // have a normal vector on the base line pointing outwards
+  BaseLineNormal.Zero();
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    BaseLineNormal.AddVector(&runner->second->NormalVector);
+  BaseLineNormal.Normalize();
+  // and calculate the sum of the angles with this normal vector and each of the triangle ones'
+  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++)
+    angle += BaseLineNormal.Angle(&runner->second->NormalVector);
+
+  if ((angle - M_PI) > -MYEPSILON)
+    return true;
+  else
+    return false;
+}
+
+/** Checks whether point is any of the two endpoints this line contains.
+ * \param *point point to test
+ * \return true - point is of the line, false - is not
+ */
+bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<2;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
+
 ostream &
 operator <<(ostream &ost, BoundaryLineSet &a)
@@ -213 +268 @@
 ;
 
-void
-BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+/** Calculates the normal vector for this triangle.
+ * Is made unique by comparison with \a OtherVector to point in the other direction.
+ * \param &OtherVector direction vector to make normal vector unique.
+ */
+void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
 {
   // get normal vector
@@ -223 +281 @@
   if (NormalVector.Projection(&OtherVector) > 0)
     NormalVector.Scale(-1.);
-}
-;
+};
+
+/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: if it's really outside, then for any endpoint of the triangle and it's opposite
+ * base line, the intersection between the line from endpoint to intersection and the base line will have a Vector::NormSquared()
+ * smaller than the first line.
+ * \param *out output stream for debugging
+ * \param *MolCenter offset vector of line
+ * \param *x second endpoint of line, minus \a *MolCenter is directional vector of line
+ * \param *Intersection intersection on plane on return
+ * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
+ */
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(ofstream *out, Vector *MolCenter, Vector *x, Vector *Intersection)
+{
+  Vector CrossPoint;
+  Vector helper;
+  int i=0;
+
+  if (Intersection->GetIntersectionWithPlane(out, &NormalVector, &endpoints[0]->node->x, MolCenter, x)) {
+    *out << Verbose(1) << "Alas! [Bronstein] failed - at least numerically - the intersection is not on the plane!" << endl;
+    return false;
+  }
+
+  // Calculate cross point between one baseline and the line from the third endpoint to intersection
+  do {
+    CrossPoint.GetIntersectionOfTwoLinesOnPlane(out, &endpoints[i%3]->node->x, &endpoints[(i+1)%3]->node->x, &endpoints[(i+2)%3]->node->x, Intersection);
+    helper.CopyVector(&endpoints[(i+1)%3]->node->x);
+    helper.SubtractVector(&endpoints[i%3]->node->x);
+    i++;
+    if (i>3)
+      break;
+  } while (CrossPoint.NormSquared() < MYEPSILON);
+  if (i>3) {
+    *out << Verbose(1) << "ERROR: Could not find any cross points, something's utterly wrong here!" << endl;
+    exit(255);
+  }
+  CrossPoint.SubtractVector(&endpoints[i%3]->node->x);
+
+  // check whether intersection is inside or not by comparing length of intersection and length of cross point
+  if ((CrossPoint.NormSquared() - helper.NormSquared()) > -MYEPSILON) { // inside
+    return true;
+  } else { // outside!
+    Intersection->Zero();
+    return false;
+  }
+};
+
+/** Checks whether lines is any of the three boundary lines this triangle contains.
+ * \param *line line to test
+ * \return true - line is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+{
+  for(int i=0;i<3;i++)
+    if (line == lines[i])
+      return true;
+  return false;
+};
+
+/** Checks whether point is any of the three endpoints this triangle contains.
+ * \param *point point to test
+ * \return true - point is of the triangle, false - is not
+ */
+bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+{
+  for(int i=0;i<3;i++)
+    if (point == endpoints[i])
+      return true;
+  return false;
+};
 
 ostream &
@@ -239 +367 @@
 
 CandidateForTesselation::CandidateForTesselation(
-        atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
+  atom *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter
 ) {
-        point = candidate;
-        BaseLine = line;
-        OptCenter.CopyVector(&OptCandidateCenter);
-        OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
+  point = candidate;
+  BaseLine = line;
+  OptCenter.CopyVector(&OptCandidateCenter);
+  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
 }
 
 CandidateForTesselation::~CandidateForTesselation() {
-        point = NULL;
-        BaseLine = NULL;
+  point = NULL;
+  BaseLine = NULL;
 }
 
@@ -301 +429 @@
   LineMap LinesOnBoundary;
   TriangleMap TrianglesOnBoundary;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  Vector helper;
 
   *out << Verbose(1) << "Finding all boundary points." << endl;
@@ -308 +438 @@
   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)
+  for (int axis = 0; axis < NDIM; axis++) {
+    AxisVector.Zero();
+    AngleReferenceVector.Zero();
+    AngleReferenceNormalVector.Zero();
+    AxisVector.x[axis] = 1.;
+    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
+    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+
+    *out << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
+
+    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
+    Walker = mol->start;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      Vector ProjectedVector;
+      ProjectedVector.CopyVector(&Walker->x);
+      ProjectedVector.SubtractVector(MolCenter);
+      ProjectedVector.ProjectOntoPlane(&AxisVector);
+
+      // correct for negative side
+      radius = ProjectedVector.NormSquared();
+      if (fabs(radius) > MYEPSILON)
+        angle = ProjectedVector.Angle(&AngleReferenceVector);
+      else
+        angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
+
+      //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
+      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+        angle = 2. * M_PI - angle;
+      }
+      *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;
+      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
+      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
+        *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
+        *out << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl;
+        *out << Verbose(2) << "New vector: " << *Walker << endl;
+        double tmp = ProjectedVector.NormSquared();
+        if ((tmp - BoundaryTestPair.first->second.first) > MYEPSILON) {
+          BoundaryTestPair.first->second.first = tmp;
+          BoundaryTestPair.first->second.second = Walker;
+          *out << Verbose(2) << "Keeping new vector due to larger projected distance " << tmp << "." << endl;
+        } else if (fabs(tmp - BoundaryTestPair.first->second.first) < MYEPSILON) {
+          helper.CopyVector(&Walker->x);
+          helper.SubtractVector(MolCenter);
+          tmp = helper.NormSquared();
+          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
+          helper.SubtractVector(MolCenter);
+          if (helper.NormSquared() < tmp) {
+            BoundaryTestPair.first->second.second = Walker;
+            *out << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl;
+          } else {
+            *out << Verbose(2) << "Keeping present vector due to larger distance to molecule center " << tmp << "." << endl;
+          }
+        } else {
+          *out << Verbose(2) << "Keeping present vector due to larger projected distance " << tmp << "." << endl;
+        }
+      }
+    }
+    // printing all inserted for debugging
+    //    {
+    //      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+    //      int i=0;
+    //      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+    //        if (runner != BoundaryPoints[axis].begin())
+    //          *out << ", " << i << ": " << *runner->second.second;
+    //        else
+    //          *out << i << ": " << *runner->second.second;
+    //        i++;
+    //      }
+    //      *out << endl;
+    //    }
+    // 3c. throw out points whose distance is less than the mean of left and right neighbours
+    bool flag = false;
+    *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
+    do { // do as long as we still throw one out per round
+      flag = false;
+      Boundaries::iterator left = BoundaryPoints[axis].end();
+      Boundaries::iterator right = BoundaryPoints[axis].end();
+      for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        // set neighbours correctly
+        if (runner == BoundaryPoints[axis].begin()) {
+          left = BoundaryPoints[axis].end();
+        } else {
+          left = runner;
+        }
+        left--;
+        right = runner;
+        right++;
+        if (right == BoundaryPoints[axis].end()) {
+          right = BoundaryPoints[axis].begin();
+        }
+        // check distance
+
+        // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
         {
-          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;
-                }
-            }
+          Vector SideA, SideB, SideC, SideH;
+          SideA.CopyVector(&left->second.second->x);
+          SideA.SubtractVector(MolCenter);
+          SideA.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideA: ";
+          //          SideA.Output(out);
+          //          *out << endl;
+
+          SideB.CopyVector(&right->second.second->x);
+          SideB.SubtractVector(MolCenter);
+          SideB.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideB: ";
+          //          SideB.Output(out);
+          //          *out << endl;
+
+          SideC.CopyVector(&left->second.second->x);
+          SideC.SubtractVector(&right->second.second->x);
+          SideC.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideC: ";
+          //          SideC.Output(out);
+          //          *out << endl;
+
+          SideH.CopyVector(&runner->second.second->x);
+          SideH.SubtractVector(MolCenter);
+          SideH.ProjectOntoPlane(&AxisVector);
+          //          *out << "SideH: ";
+          //          SideH.Output(out);
+          //          *out << endl;
+
+          // calculate each length
+          double a = SideA.Norm();
+          //double b = SideB.Norm();
+          //double c = SideC.Norm();
+          double h = SideH.Norm();
+          // calculate the angles
+          double alpha = SideA.Angle(&SideH);
+          double beta = SideA.Angle(&SideC);
+          double gamma = SideB.Angle(&SideH);
+          double delta = SideC.Angle(&SideH);
+          double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
+          //*out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
+          *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
+          if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) {
+            // throw out point
+            *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
+            BoundaryPoints[axis].erase(runner);
+            flag = true;
+          }
         }
-      // 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);
-    }
+      }
+    } while (flag);
+  }
+  delete(MolCenter);
   return BoundaryPoints;
 }
@@ -619 +714 @@
       *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
+        *vrmlfile << Walker->x.x[i]-center->x[i] << " ";
       *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
     }
@@ -628 +723 @@
       *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 << 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 << Binder->rightatom->x.x[i]-center->x[i] << " ";
       *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
     }
@@ -640 +735 @@
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
         *vrmlfile << "\t";
       }
@@ -673 +768 @@
       *rasterfile << "2" << endl << "  ";  // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *rasterfile << Walker->x.x[i]+center->x[i] << " ";
+        *rasterfile << Walker->x.x[i]-center->x[i] << " ";
       *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
     }
@@ -682 +777 @@
       *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 << 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 << Binder->rightatom->x.x[i]-center->x[i] << " ";
       *rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
     }
@@ -695 +790 @@
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]-center->x[j] << " ";
         *rasterfile << "\t";
       }
@@ -713 +808 @@
  * \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;
-    }
+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;
+  }
 }
 
-/** Determines the volume of a cluster.
- * Determines first the convex envelope, then tesselates it and calculates its volume.
+/** 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
- * \param *configuration needed for path to store convex envelope file
- * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
- * \param *mol molecule structure representing the cluster
- * \return determined volume of the cluster in cubed config:GetIsAngstroem()
- */
-double
-VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
-    Boundaries *BoundaryPtr, molecule *mol)
-{
-  bool IsAngstroem = configuration->GetIsAngstroem();
-  atom *Walker = NULL;
-  struct Tesselation *TesselStruct = new Tesselation;
+ * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
+ */
+void Find_convex_border(ofstream *out, molecule* mol, class Tesselation *&TesselStruct, class LinkedCell *LCList, const char *filename)
+{
   bool BoundaryFreeFlag = false;
-  Boundaries *BoundaryPoints = BoundaryPtr;
-  double volume = 0.;
-  double PyramidVolume = 0.;
-  double G, h;
-  Vector x, y;
-  double a, b, c;
-
-  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
-
-  // 1. calculate center of gravity
-  *out << endl;
-  Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
-
-  // 2. translate all points into CoG
-  *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
-  Walker = mol->start;
-  while (Walker->next != mol->end)
-    {
-      Walker = Walker->next;
-      Walker->x.Translate(CenterOfGravity);
-    }
-
-  // 3. Find all points on the boundary
-  if (BoundaryPoints == NULL)
-    {
+  Boundaries *BoundaryPoints = NULL;
+
+  cout << Verbose(1) << "Begin of find_convex_border" << endl;
+
+  if (TesselStruct != NULL) // free if allocated
+    delete(TesselStruct);
+  TesselStruct = new class Tesselation;
+
+  // 1. Find all points on the boundary
+  if (BoundaryPoints == NULL) {
       BoundaryFreeFlag = true;
       BoundaryPoints = GetBoundaryPoints(out, mol);
-    }
-  else
+  } else {
+      *out << Verbose(1) << "Using given boundary points set." << endl;
+  }
+
+// printing all inserted for debugging
+  for (int axis=0; axis < NDIM; axis++)
     {
-      *out << Verbose(1) << "Using given boundary points set." << endl;
-    }
-
-  // 4. fill the boundary point list
+      *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
+      int i=0;
+      for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
+        if (runner != BoundaryPoints[axis].begin())
+          *out << ", " << i << ": " << *runner->second.second;
+        else
+          *out << i << ": " << *runner->second.second;
+        i++;
+      }
+      *out << endl;
+    }
+
+  // 2. fill the boundary point list
   for (int axis = 0; axis < NDIM; axis++)
-    for (Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner
-        != BoundaryPoints[axis].end(); runner++)
-      {
+    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;
+
+  *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
 
@@ -828 +895 @@
   //  *out << endl;
 
-  // 5a. guess starting triangle
+  // 3a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
 
-  // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
-  TesselStruct->TesselateOnBoundary(out, configuration, mol);
-
-  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount
-      << " triangles with " << TesselStruct->LinesOnBoundaryCount
-      << " lines and " << TesselStruct->PointsOnBoundaryCount << " points."
-      << endl;
+  // 3b. go through all lines, that are not yet part of two triangles (only of one so far)
+  TesselStruct->TesselateOnBoundary(out, mol);
+
+  // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point
+  if (!TesselStruct->InsertStraddlingPoints(out, mol))
+    *out << Verbose(1) << "Insertion of straddling points failed!" << endl;
+
+  // 3d. check all baselines whether the peaks of the two adjacent triangles with respect to center of baseline are convex, if not, make the baseline between the two peaks and baseline endpoints become the new peaks
+  if (!TesselStruct->CorrectConcaveBaselines(out))
+    *out << Verbose(1) << "Correction of concave baselines failed!" << endl;
+
+  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;
+
+  // 4. Store triangles in tecplot file
+  if (filename != NULL) {
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
+      write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(Raster3DSuffix);
+      ofstream *rasterplot = new ofstream(OutputName.c_str());
+      write_raster3d_file(out, rasterplot, TesselStruct, mol);
+      rasterplot->close();
+      delete(rasterplot);
+    }
+  }
+
+  // free reference lists
+  if (BoundaryFreeFlag)
+    delete[] (BoundaryPoints);
+
+  cout << Verbose(1) << "End of find_convex_border" << endl;
+};
+
+
+/** Determines the volume of a cluster.
+ * Determines first the convex envelope, then tesselates it and calculates its volume.
+ * \param *out output stream for debugging
+ * \param *TesselStruct Tesselation filled with points, lines and triangles on boundary on return
+ * \param *configuration needed for path to store convex envelope file
+ * \return determined volume of the cluster in cubed config:GetIsAngstroem()
+ */
+double VolumeOfConvexEnvelope(ofstream *out, class Tesselation *TesselStruct, class config *configuration)
+{
+  bool IsAngstroem = configuration->GetIsAngstroem();
+  double volume = 0.;
+  double PyramidVolume = 0.;
+  double G, h;
+  Vector x, y;
+  double a, b, c;
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
@@ -873 +989 @@
       << endl;
 
-  // 7. translate all points back from CoG
-  *out << Verbose(1) << "Translating system back from Center of Gravity."
-      << endl;
-  CenterOfGravity->Scale(-1);
-  Walker = mol->start;
-  while (Walker->next != mol->end)
-    {
-      Walker = Walker->next;
-      Walker->x.Translate(CenterOfGravity);
-    }
-
-  // 8. Store triangles in tecplot file
-  string OutputName(filename);
-  OutputName.append(TecplotSuffix);
-  ofstream *tecplot = new ofstream(OutputName.c_str());
-  write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
-  tecplot->close();
-  delete(tecplot);
-
-  // free reference lists
-  if (BoundaryFreeFlag)
-    delete[] (BoundaryPoints);
-
   return volume;
 }
@@ -918 +1011 @@
   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, NULL, configuration,
-        BoundaryPoints, mol);
+    clustervolume = VolumeOfConvexEnvelope(out, TesselStruct, configuration);
   else
     clustervolume = ClusterVolume;
-  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol,
-      IsAngstroem);
+  double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
   Vector BoxLengths;
   int repetition[NDIM] =
@@ -1010 +1104 @@
       // set new box dimensions
       *out << Verbose(0) << "Translating to box with these boundaries." << endl;
-      mol->CenterInBox((ofstream *) &cout, &BoxLengths);
+      mol->SetBoxDimension(&BoxLengths);
+      mol->CenterInBox((ofstream *) &cout);
     }
   // update Box of atoms by boundary
@@ -1020 +1115 @@
 }
 ;
+
+
+/** Fills the empty space of the simulation box with water/
+ * \param *out output stream for debugging
+ * \param *List list of molecules already present in box
+ * \param *TesselStruct contains tesselated surface
+ * \param *filler molecule which the box is to be filled with
+ * \param configuration contains box dimensions
+ * \param distance[NDIM] distance between filling molecules in each direction
+ * \param RandAtomDisplacement maximum distance for random displacement per atom
+ * \param RandMolDisplacement maximum distance for random displacement per filler molecule
+ * \param DoRandomRotation true - do random rotiations, false - don't
+ * \return *mol pointer to new molecule with filled atoms
+ */
+molecule * FillBoxWithMolecule(ofstream *out, MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
+{
+  molecule *Filling = new molecule(filler->elemente);
+  Vector CurrentPosition;
+  int N[NDIM];
+  int n[NDIM];
+  double *M =  filler->ReturnFullMatrixforSymmetric(filler->cell_size);
+  double Rotations[NDIM*NDIM];
+  Vector AtomTranslations;
+  Vector FillerTranslations;
+  Vector FillerDistance;
+  double FillIt = false;
+  atom *Walker = NULL, *Runner = NULL;
+  bond *Binder = NULL;
+
+  // Center filler at origin
+  filler->CenterOrigin(out);
+  filler->Center.Zero();
+
+  // calculate filler grid in [0,1]^3
+  FillerDistance.Init(distance[0], distance[1], distance[2]);
+  FillerDistance.InverseMatrixMultiplication(M);
+  for(int i=0;i<NDIM;i++)
+    N[i] = ceil(1./FillerDistance.x[i]);
+
+  // go over [0,1]^3 filler grid
+  for (n[0] = 0; n[0] < N[0]; n[0]++)
+    for (n[1] = 0; n[1] < N[1]; n[1]++)
+      for (n[2] = 0; n[2] < N[2]; n[2]++) {
+        // calculate position of current grid vector in untransformed box
+        CurrentPosition.Init(n[0], n[1], n[2]);
+        CurrentPosition.MatrixMultiplication(M);
+        // Check whether point is in- or outside
+        FillIt = true;
+        for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
+          FillIt = FillIt && (!(*ListRunner)->TesselStruct->IsInside(&CurrentPosition));
+        }
+
+        if (FillIt) {
+          // fill in Filler
+          *out << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
+
+          // create molecule random translation vector ...
+          for (int i=0;i<NDIM;i++)
+            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+
+          // go through all atoms
+          Walker = filler->start;
+          while (Walker != filler->end) {
+            Walker = Walker->next;
+            // copy atom ...
+            *Runner = new atom(Walker);
+
+            // create atomic random translation vector ...
+            for (int i=0;i<NDIM;i++)
+              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+
+            // ... and rotation matrix
+            if (DoRandomRotation) {
+              double phi[NDIM];
+              for (int i=0;i<NDIM;i++) {
+                phi[i] = rand()/(RAND_MAX/(2.*M_PI));
+              }
+
+              Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
+              Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
+              Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
+              Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
+              Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
+              Rotations[7] =               sin(phi[1])                                                ;
+              Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
+              Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
+              Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            }
+
+            // ... and put at new position
+            if (DoRandomRotation)
+              Runner->x.MatrixMultiplication(Rotations);
+            Runner->x.AddVector(&AtomTranslations);
+            Runner->x.AddVector(&FillerTranslations);
+            Runner->x.AddVector(&CurrentPosition);
+            // insert into Filling
+            Filling->AddAtom(Runner);
+          }
+
+          // go through all bonds and add as well
+          Binder = filler->first;
+          while(Binder != filler->last) {
+            Binder = Binder->next;
+          }
+        } else {
+          // leave empty
+          *out << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
+        }
+      }
+  return Filling;
+};
+
 
 // =========================================================== class TESSELATION ===========================================
@@ -1133 +1340 @@
                   << A->second->node->Name << ","
                   << baseline->second.first->second->node->Name << ","
-                  << baseline->second.second->second->node->Name << " leave "
+                  << baseline->second.second->second->node->Name << " leaves "
                   << checker->second->node->Name << " outside the convex hull."
                   << endl;
@@ -1218 +1425 @@
  * -# 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)
+ *    -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
  *    -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
  * \param *out output stream for debugging
@@ -1226 +1433 @@
  * \param *mol the cluster as a molecule structure
  */
-void
-Tesselation::TesselateOnBoundary(ofstream *out, config *configuration,
-    molecule *mol)
+void Tesselation::TesselateOnBoundary(ofstream *out, molecule *mol)
 {
   bool flag;
@@ -1234 +1439 @@
   class BoundaryPointSet *peak = NULL;
   double SmallestAngle, TempAngle;
-  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector,
-      PropagationVector;
+  Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *MolCenter = NULL;
   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)
+
+  MolCenter = mol->DetermineCenterOfAll(out);
+  // create a first tesselation with the given BoundaryPoints
+  do {
+    flag = false;
+    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
+      if (baseline->second->TrianglesCount == 1) {
+        // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
+        SmallestAngle = M_PI;
+
+        // get peak point with respect to this base line's only triangle
+        BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
+        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
+        for (int i = 0; i < 3; i++)
+          if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
+            peak = BTS->endpoints[i];
+        *out << Verbose(3) << " and has peak " << *peak << "." << endl;
+
+        // prepare some auxiliary vectors
+        Vector BaseLineCenter, BaseLine;
+        BaseLineCenter.CopyVector(&baseline->second->endpoints[0]->node->x);
+        BaseLineCenter.AddVector(&baseline->second->endpoints[1]->node->x);
+        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
+        BaseLine.CopyVector(&baseline->second->endpoints[0]->node->x);
+        BaseLine.SubtractVector(&baseline->second->endpoints[1]->node->x);
+
+        // offset to center of triangle
+        CenterVector.Zero();
+        for (int i = 0; i < 3; i++)
+          CenterVector.AddVector(&BTS->endpoints[i]->node->x);
+        CenterVector.Scale(1. / 3.);
+        *out << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl;
+
+        // normal vector of triangle
+        NormalVector.CopyVector(MolCenter);
+        NormalVector.SubtractVector(&CenterVector);
+        BTS->GetNormalVector(NormalVector);
+        NormalVector.CopyVector(&BTS->NormalVector);
+        *out << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl;
+
+        // vector in propagation direction (out of triangle)
+        // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
+        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
+        TempVector.CopyVector(&CenterVector);
+        TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
+        if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+          PropagationVector.Scale(-1.);
+        *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
+        winner = PointsOnBoundary.end();
+
+        // loop over all points and calculate angle between normal vector of new and present triangle
+        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
+          if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
+            *out << Verbose(3) << "Target point is " << *(target->second) << ":" << endl;
+
+            // first check direction, so that triangles don't intersect
+            VirtualNormalVector.CopyVector(&target->second->node->x);
+            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
+            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
+            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+            *out << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
+            if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
+              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;
+              continue;
+            } else
+              *out << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;
+
+            // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
+            LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
+            LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
+            if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->TrianglesCount == 2))) {
+              *out << Verbose(4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->TrianglesCount << " triangles." << endl;
+              continue;
+            }
+            if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->TrianglesCount == 2))) {
+              *out << Verbose(4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->TrianglesCount << " triangles." << endl;
+              continue;
+            }
+
+            // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
+            if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
+              *out << Verbose(4) << "Current target is peak!" << endl;
+              continue;
+            }
+
+            // check for linear dependence
             TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
-            TempVector.SubtractVector(&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
+            TempVector.SubtractVector(&target->second->node->x);
+            helper.CopyVector(&baseline->second->endpoints[1]->node->x);
+            helper.SubtractVector(&target->second->node->x);
+            helper.ProjectOntoPlane(&TempVector);
+            if (fabs(helper.NormSquared()) < MYEPSILON) {
+              *out << Verbose(4) << "Chosen set of vectors is linear dependent." << endl;
+              continue;
+            }
+
+            // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
+            flag = true;
+            VirtualNormalVector.MakeNormalVector(&baseline->second->endpoints[0]->node->x, &baseline->second->endpoints[1]->node->x, &target->second->node->x);
+            TempVector.CopyVector(&baseline->second->endpoints[0]->node->x);
+            TempVector.AddVector(&baseline->second->endpoints[1]->node->x);
+            TempVector.AddVector(&target->second->node->x);
+            TempVector.Scale(1./3.);
+            TempVector.SubtractVector(MolCenter);
+            // make it always point outward
+            if (VirtualNormalVector.Projection(&TempVector) < 0)
+              VirtualNormalVector.Scale(-1.);
+            // calculate angle
+            TempAngle = NormalVector.Angle(&VirtualNormalVector);
+            *out << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
+            if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
+              SmallestAngle = TempAngle;
+              winner = target;
+              *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
+            } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
+              // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
+              helper.CopyVector(&target->second->node->x);
+              helper.SubtractVector(&BaseLineCenter);
+              helper.ProjectOntoPlane(&BaseLine);
+              // ...the one with the smaller angle is the better candidate
+              TempVector.CopyVector(&target->second->node->x);
+              TempVector.SubtractVector(&BaseLineCenter);
+              TempVector.ProjectOntoPlane(&VirtualNormalVector);
+              TempAngle = TempVector.Angle(&helper);
+              TempVector.CopyVector(&winner->second->node->x);
+              TempVector.SubtractVector(&BaseLineCenter);
+              TempVector.ProjectOntoPlane(&VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(&helper)) {
+                TempAngle = NormalVector.Angle(&VirtualNormalVector);
+                SmallestAngle = TempAngle;
+                winner = target;
+                *out << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;
+              } else
+                *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;
+            } else
+              *out << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;
           }
+        } // end of loop over all boundary points
+
+        // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
+        if (winner != PointsOnBoundary.end()) {
+          *out << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;
+          // create the lins of not yet present
+          BLS[0] = baseline->second;
+          // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
+          LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
+          LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
+          if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
+            BPS[0] = baseline->second->endpoints[0];
+            BPS[1] = winner->second;
+            BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
+            LinesOnBoundaryCount++;
+          } else
+            BLS[1] = LineChecker[0]->second;
+          if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
+            BPS[0] = baseline->second->endpoints[1];
+            BPS[1] = winner->second;
+            BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+            LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
+            LinesOnBoundaryCount++;
+          } else
+            BLS[2] = LineChecker[1]->second;
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+          TrianglesOnBoundaryCount++;
+        } else {
+          *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
+        }
+
+        // 5d. If the set of lines is not yet empty, go to 5. and continue
+      } else
+        *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
+  } while (flag);
+
+  // exit
+  delete(MolCenter);
+};
+
+/** Inserts all atoms outside of the tesselated surface into it by adding new triangles.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure with atoms
+ * \return true - all straddling points insert, false - something went wrong
+ */
+bool Tesselation::InsertStraddlingPoints(ofstream *out, molecule *mol)
+{
+  Vector Intersection;
+  atom *Walker = mol->start;
+  Vector *MolCenter = mol->DetermineCenterOfAll(out);
+  while (Walker->next != mol->end) {  // we only have to go once through all points, as boundary can become only bigger
+    // get the next triangle
+    BTS = FindClosestTriangleToPoint(out, &Walker->x);
+    if (BTS == NULL) {
+      *out << Verbose(1) << "ERROR: No triangle closest to " << Walker << " was found." << endl;
+      return false;
+    }
+    // get the intersection point
+    if (BTS->GetIntersectionInsideTriangle(out, MolCenter, &Walker->x, &Intersection)) {
+      // we have the intersection, check whether in- or outside of boundary
+      if ((MolCenter->DistanceSquared(&Walker->x) - MolCenter->DistanceSquared(&Intersection)) < -MYEPSILON) {
+        // inside, next!
+        *out << Verbose(4) << Walker << " is inside wrt triangle " << BTS << "." << endl;
+      } else {
+        // outside!
+        *out << Verbose(3) << Walker << " is outside wrt triangle " << BTS << "." << endl;
+        class BoundaryLineSet *OldLines[3], *NewLines[3];
+        class BoundaryPointSet *OldPoints[3], *NewPoint;
+        // store the three old lines and old points
+        for (int i=0;i<3;i++) {
+          OldLines[i] = BTS->lines[i];
+          OldPoints[i] = BTS->endpoints[i];
+        }
+        // add Walker to boundary points
+        AddPoint(Walker);
+        if (BPS[0] == NULL)
+          NewPoint = BPS[0];
         else
-          *out << Verbose(2) << "Baseline candidate " << *(baseline->second)
-              << " has a triangle count of "
-              << baseline->second->TrianglesCount << "." << endl;
-    }
-  while (flag);
-
-}
-;
+          continue;
+        // remove triangle
+        TrianglesOnBoundary.erase(BTS->Nr);
+        // create three new boundary lines
+        for (int i=0;i<3;i++) {
+          BPS[0] = NewPoint;
+          BPS[1] = OldPoints[i];
+          NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
+          LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
+          LinesOnBoundaryCount++;
+        }
+        // create three new triangle with new point
+        for (int i=0;i<3;i++) { // find all baselines
+          BLS[0] = OldLines[i];
+          int n = 1;
+          for (int j=0;j<3;j++) {
+            if (NewLines[j]->IsConnectedTo(BLS[0])) {
+              if (n>2) {
+                *out << Verbose(1) << "ERROR: " << BLS[0] << " connects to all of the new lines?!" << endl;
+                return false;
+              } else
+                BLS[n++] = NewLines[j];
+            }
+          }
+          // create the triangle
+          BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+          TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
+          TrianglesOnBoundaryCount++;
+        }
+      }
+    } else { // something is wrong with FindClosestTriangleToPoint!
+      *out << Verbose(1) << "ERROR: The closest triangle did not produce an intersection!" << endl;
+      return false;
+    }
+  }
+
+  // exit
+  delete(MolCenter);
+  return true;
+};
+
+/** Goes over all baselines and checks whether adjacent triangles and convex to each other.
+ * \param *out output stream for debugging
+ * \return true - all baselines were corrected, false - there are still concave pieces
+ */
+bool Tesselation::CorrectConcaveBaselines(ofstream *out)
+{
+  class BoundaryTriangleSet *triangle[2];
+  class BoundaryLineSet *OldLines[4], *NewLine;
+  class BoundaryPointSet *OldPoints[2];
+  Vector BaseLineNormal;
+  class BoundaryLineSet *Base = NULL;
+  int OldTriangles[2], OldBaseLine;
+  int i;
+  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) {
+    Base = baseline->second;
+
+    // check convexity
+    if (Base->CheckConvexityCriterion(out)) { // triangles are convex
+      *out << Verbose(3) << Base << " has two convex triangles." << endl;
+    } else { // not convex!
+      // get the two triangles
+      i=0;
+      for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
+        triangle[i++] = runner->second;
+      // gather four endpoints and four lines
+      for (int j=0;j<4;j++)
+        OldLines[j] = NULL;
+      for (int j=0;j<2;j++)
+        OldPoints[j] = NULL;
+      i=0;
+      for (int m=0;m<2;m++) { // go over both triangles
+        for (int j=0;j<3;j++) { // all of their endpoints and baselines
+          if (triangle[m]->lines[j] != Base) // pick not the central baseline
+            OldLines[i++] = triangle[m]->lines[j];
+          if (!Base->ContainsBoundaryPoint(triangle[m]->endpoints[j])) // and neither of its endpoints
+            OldPoints[m] = triangle[m]->endpoints[j];
+        }
+      }
+      if (i<4) {
+        *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+        return false;
+      }
+      for (int j=0;j<4;j++)
+        if (OldLines[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough baselines!" << endl;
+          return false;
+        }
+      for (int j=0;j<2;j++)
+        if (OldPoints[j] == NULL) {
+          *out << Verbose(1) << "ERROR: We have not gathered enough endpoints!" << endl;
+          return false;
+        }
+
+      // remove triangles
+      for (int j=0;j<2;j++) {
+        OldTriangles[j] = triangle[j]->Nr;
+        TrianglesOnBoundary.erase(OldTriangles[j]);
+        triangle[j] = NULL;
+      }
+
+      // remove baseline
+      OldBaseLine = Base->Nr;
+      LinesOnBoundary.erase(OldBaseLine);
+      Base = NULL;
+
+      // construct new baseline (with same number as old one)
+      BPS[0] = OldPoints[0];
+      BPS[1] = OldPoints[1];
+      NewLine = new class BoundaryLineSet(BPS, OldBaseLine);
+      LinesOnBoundary.insert(LinePair(OldBaseLine, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
+
+      // construct new triangles with flipped baseline
+      i=-1;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=2;
+      if (BLS[0]->IsConnectedTo(OldLines[2]))
+        i=3;
+      if (i!=-1) {
+        BLS[0] = OldLines[0];
+        BLS[1] = OldLines[i];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[0]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[0], BTS));
+
+        BLS[0] = (i==2 ? OldLines[3] : OldLines[2]);
+        BLS[1] = OldLines[1];
+        BLS[2] = NewLine;
+        BTS = new class BoundaryTriangleSet(BLS, OldTriangles[1]);
+        TrianglesOnBoundary.insert(TrianglePair(OldTriangles[1], BTS));
+      } else {
+        *out << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;
+        return false;
+      }
+    }
+  }
+  return true;
+};
+
+
+/** States whether point is in- or outside of a tesselated surface.
+ * \param *pointer point to be checked
+ * \return true - is inside, false - is outside
+ */
+bool Tesselation::IsInside(Vector *pointer)
+{
+
+  // hier kommt dann Saskias Routine hin...
+
+  return true;
+};
+
+/** Finds the closest triangle to a given point.
+ * \param *out output stream for debugging
+ * \param *x second endpoint of line
+ * \return pointer triangle that is closest, NULL if none was found
+ */
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(ofstream *out, Vector *x)
+{
+  class BoundaryTriangleSet *triangle = NULL;
+
+  // hier kommt dann Saskias Routine hin...
+
+  return triangle;
+};
 
 /** Adds an atom to the tesselation::PointsOnBoundary list.
@@ -1463 +1836 @@
   if (InsertUnique.second) // if new point was not present before, increase counter
     PointsOnBoundaryCount++;
+  else {
+    delete(BPS[0]);
+    BPS[0] = NULL;
+  }
 }
 ;
@@ -1804 +2181 @@
  */
 int Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *Candidates[3]) {
-// TODO: use findTriangles here and return result.size();
   LineMap::iterator FindLine;
   PointMap::iterator FindPoint;
@@ -1904 +2280 @@
   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 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;
@@ -1985 +2361 @@
 
                   if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x)))
-                        && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
                   ) {
                     helper.CopyVector(&NewNormalVector);
@@ -2072 +2448 @@
   //cout << Verbose(2) << "INFO: Sorting candidate list ..." << endl;
   if (candidates->size() > 1) {
-          candidates->unique();
-          candidates->sort(sortCandidates);
+    candidates->unique();
+    candidates->sort(sortCandidates);
   }
 
@@ -2080 +2456 @@
 
 struct Intersection {
-        Vector x1;
-        Vector x2;
-        Vector x3;
-        Vector x4;
+  Vector x1;
+  Vector x2;
+  Vector x3;
+  Vector x4;
 };
 
@@ -2093 +2469 @@
  */
 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;
+  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;
 };
 
@@ -2132 +2508 @@
  */
 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);
+  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;
@@ -2179 +2555 @@
 
         if (status == GSL_SUCCESS) {
-                cout << Verbose(2) << "converged to minimum" <<  endl;
+          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
+  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;
+  return result;
 }
 
@@ -2241 +2617 @@
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
 
-  if (LC->SetIndexToAtom(*a)) {  // get cell for the starting atom
+  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];
@@ -2603 +2979 @@
   cout << Verbose(1) << "Third Points are ";
   for (CandidateList::iterator it = Opt_Candidates->begin(); it != Opt_Candidates->end(); ++it) {
-          cout << " " << *(*it)->point;
+    cout << " " << *(*it)->point;
   }
   cout << endl;
@@ -2609 +2985 @@
   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) {
+    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);
@@ -2640 +3016 @@
         << " 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.
+    } 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);
@@ -2679 +3055 @@
     }
 
-          if ((result) && (existentTrianglesCount < 2) && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 1 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+    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) {
@@ -2711 +3087 @@
         (*it)->OptCenter.AddVector(&helper);
         Vector *center = mol->DetermineCenterOfAll(out);
-        (*it)->OptCenter.AddVector(center);
+        (*it)->OptCenter.SubtractVector(center);
         delete(center);
         // and add to file plus translucency object
@@ -2726 +3102 @@
       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];
+    BaseRay = BLS[0];
   }
 
@@ -2747 +3123 @@
  */
 bool sortCandidates(CandidateForTesselation* candidate1, CandidateForTesselation* candidate2) {
-  // TODO: use get angle and remove duplicate code
   Vector BaseLineVector, OrthogonalVector, helper;
-        if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
-          cout << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
-          //return false;
-          exit(1);
-        }
-        // create baseline vector
-        BaseLineVector.CopyVector(&(candidate1->BaseLine->endpoints[1]->node->x));
-        BaseLineVector.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
-        BaseLineVector.Normalize();
+  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();
+  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);
+  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);
+  helper.SubtractVector(&(candidate1->BaseLine->endpoints[0]->node->x));
+  double psi = BaseLineVector.Angle(&helper);
   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-                psi = 2.*M_PI - psi;
-        }
-
-        cout << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
-        cout << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
-
-        // return comparison
-        return phi < psi;
-}
-
-/**
- * Checks whether the provided atom is within the tesselation structure.
- *
- * @param atom of which to check the position
- * @param tesselation structure
- *
- * @return true if the atom is inside the tesselation structure, false otherwise
- */
-bool IsInnerAtom(atom *Atom, class Tesselation *Tess, LinkedCell* LC) {
-  if (Tess->LinesOnBoundary.begin() == Tess->LinesOnBoundary.end()) {
-    cout << Verbose(0) << "Error: There is no tesselation structure to compare the atom with, "
-        << "please create one first.";
-    exit(1);
-  }
-
-  class atom *trianglePoints[3];
-  trianglePoints[0] = findClosestAtom(Atom, LC);
-  // check whether closest atom is "too close" :), then it's inside
-  if (trianglePoints[0]->x.DistanceSquared(&Atom->x) < MYEPSILON)
-    return true;
-  list<atom*> *connectedClosestAtoms = Tess->getClosestConnectedAtoms(trianglePoints[0], Atom);
-  trianglePoints[1] = connectedClosestAtoms->front();
-  trianglePoints[2] = connectedClosestAtoms->back();
-  for (int i=0;i<3;i++) {
-    if (trianglePoints[i] == NULL) {
-      cout << Verbose(1) << "IsInnerAtom encounters serious error, point " << i << " not found." << endl;
-    }
-
-    cout << Verbose(1) << "List of possible atoms:" << endl;
-    cout << *trianglePoints[i] << endl;
-
-//    for(list<atom*>::iterator runner = connectedClosestAtoms->begin(); runner != connectedClosestAtoms->end(); runner++)
-//      cout << Verbose(2) << *(*runner) << endl;
-  }
-  delete(connectedClosestAtoms);
-
-  list<BoundaryTriangleSet*> *triangles = Tess->FindTriangles(trianglePoints);
-
-  if (triangles->empty()) {
-    cout << Verbose(0) << "Error: There is no nearest triangle. Please check the tesselation structure.";
-    exit(1);
-  }
-
-  Vector helper;
-  helper.CopyVector(&Atom->x);
-
-  // Only in case of degeneration, there will be two different scalar products.
-  // If at least one scalar product is positive, the point is considered to be outside.
-  bool status = (helper.ScalarProduct(&triangles->front()->NormalVector) < 0)
-      && (helper.ScalarProduct(&triangles->back()->NormalVector) < 0);
-  delete(triangles);
-  return status;
-}
-
-/**
- * Finds the atom which is closest to the provided one.
- *
- * @param atom to which to find the closest other atom
- * @param linked cell structure
- *
- * @return atom which is closest to the provided one
- */
-atom* findClosestAtom(const atom* Atom, LinkedCell* LC) {
-  LinkedAtoms *List = NULL;
-  atom* closestAtom = NULL;
-  double distance = 1e16;
-  Vector helper;
-  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-
-  LC->SetIndexToVector(&Atom->x); // ignore status as we calculate bounds below sensibly
-  for(int i=0;i<NDIM;i++) // store indices of this cell
-    N[i] = LC->n[i];
-  //cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-
-  LC->GetNeighbourBounds(Nlower, Nupper);
-  //cout << endl;
-  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-        List = LC->GetCurrentCell();
-        //cout << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(&Atom->x);
-            helper.SubtractVector(&(*Runner)->x);
-            double currentNorm = helper. Norm();
-            if (currentNorm < distance) {
-              distance = currentNorm;
-              closestAtom = (*Runner);
-            }
-          }
-        } else {
-          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << ","
-            << LC->n[2] << " is invalid!" << endl;
-        }
-      }
-
-  return closestAtom;
-}
-
-/**
- * Gets all atoms connected to the provided atom by triangulation lines.
- *
- * @param atom of which get all connected atoms
- * @param atom to be checked whether it is an inner atom
- *
- * @return list of the two atoms linked to the provided one and closest to the atom to be checked,
- */
-list<atom*> * Tesselation::getClosestConnectedAtoms(atom* Atom, atom* AtomToCheck) {
-  list<atom*> connectedAtoms;
-  map<double, atom*> anglesOfAtoms;
-  map<double, atom*>::iterator runner;
-  LineMap::iterator findLines = LinesOnBoundary.begin();
-  list<atom*>::iterator listRunner;
-  Vector center, planeNorm, currentPoint, OrthogonalVector, helper;
-  atom* current;
-  bool takeAtom = false;
-
-  planeNorm.CopyVector(&Atom->x);
-  planeNorm.SubtractVector(&AtomToCheck->x);
-  planeNorm.Normalize();
-
-  while (findLines != LinesOnBoundary.end()) {
-    takeAtom = false;
-
-    if (findLines->second->endpoints[0]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[1]->node;
-    } else if (findLines->second->endpoints[1]->Nr == Atom->nr) {
-      takeAtom = true;
-      current = findLines->second->endpoints[0]->node;
-    }
-
-    if (takeAtom) {
-      connectedAtoms.push_back(current);
-      currentPoint.CopyVector(&current->x);
-      currentPoint.ProjectOntoPlane(&planeNorm);
-      center.AddVector(&currentPoint);
-    }
-
-    findLines++;
-  }
-
-  cout << "Summed vectors " << center << "; number of atoms " << connectedAtoms.size()
-    << "; scale factor " << 1.0/connectedAtoms.size();
-
-  center.Scale(1.0/connectedAtoms.size());
-  listRunner = connectedAtoms.begin();
-
-  cout << " calculated center " << center <<  endl;
-
-  // construct one orthogonal vector
-  helper.CopyVector(&AtomToCheck->x);
-  helper.ProjectOntoPlane(&planeNorm);
-  OrthogonalVector.MakeNormalVector(&center, &helper, &(*listRunner)->x);
-  while (listRunner != connectedAtoms.end()) {
-    double angle = getAngle((*listRunner)->x, AtomToCheck->x, center, OrthogonalVector);
-    cout << "Calculated angle " << angle << " for atom " << **listRunner << endl;
-    anglesOfAtoms.insert(pair<double, atom*>(angle, (*listRunner)));
-    listRunner++;
-  }
-
-  list<atom*> *result = new list<atom*>;
-  runner = anglesOfAtoms.begin();
-  cout << "First value is " << *runner->second << endl;
-  result->push_back(runner->second);
-  runner = anglesOfAtoms.end();
-  runner--;
-  cout << "Second value is " << *runner->second << endl;
-  result->push_back(runner->second);
-
-  cout << "List of closest atoms has " << result->size() << " elements, which are "
-    << *(result->front()) << " and " << *(result->back()) << endl;
-
-  return result;
-}
-
-/**
- * Finds triangles belonging to the three provided atoms.
- *
- * @param atom list, is expected to contain three atoms
- *
- * @return triangles which belong to the provided atoms, will be empty if there are none,
- *         will usually be one, in case of degeneration, there will be two
- */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(atom* Points[3]) {
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
-  LineMap::iterator FindLine;
-  PointMap::iterator FindPoint;
-  TriangleMap::iterator FindTriangle;
-  class BoundaryPointSet *TrianglePoints[3];
-
-  for (int i = 0; i < 3; i++) {
-    FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
-    } else {
-      TrianglePoints[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr);
-          if (FindLine != TrianglePoints[i]->lines.end()) {
-            for (; FindLine->first == TrianglePoints[j]->node->nr; FindLine++) {
-              FindTriangle = FindLine->second->triangles.begin();
-              for (; FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
-                if ((
-                  (FindTriangle->second->endpoints[0] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[0] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[1] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[1] == TrianglePoints[2])
-                  ) && (
-                    (FindTriangle->second->endpoints[2] == TrianglePoints[0])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[1])
-                    || (FindTriangle->second->endpoints[2] == TrianglePoints[2])
-                  )
-                ) {
-                  result->push_back(FindTriangle->second);
-                }
-              }
-            }
-            // Is it sufficient to consider one of the triangle lines for this.
-            return result;
-
-          }
-        }
-      }
-    }
-  }
-
-  return result;
-}
-
-/**
- * Gets the angle between a point and a reference relative to the provided center.
- *
- * @param point to calculate the angle for
- * @param reference to which to calculate the angle
- * @param center for which to calculate the angle between the vectors
- * @param OrthogonalVector helps discern between [0,pi] and [pi,2pi] in acos()
- *
- * @return angle between point and reference
- */
-double getAngle(Vector point, Vector reference, Vector center, Vector OrthogonalVector) {
-  Vector ReferenceVector, helper;
-  cout << Verbose(2) << reference << " is our reference " << endl;
-  cout << Verbose(2) << center << " is our center " << endl;
-  // create baseline vector
-  ReferenceVector.CopyVector(&reference);
-  ReferenceVector.SubtractVector(&center);
-  if (ReferenceVector.IsNull())
-    return M_PI;
-
-  // calculate both angles and correct with in-plane vector
-  helper.CopyVector(&point);
-  helper.SubtractVector(&center);
-  if (helper.IsNull())
-    return M_PI;
-  double phi = ReferenceVector.Angle(&helper);
-  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
-    phi = 2.*M_PI - phi;
-  }
-
-  cout << Verbose(2) << point << " has angle " << phi << endl;
-
-  return phi;
-}
-
-/**
- * Checks whether the provided point is within the tesselation structure.
- *
- * This is a wrapper function for IsInnerAtom, so it can be used with a simple
- * vector, too.
- *
- * @param point of which to check the position
- * @param tesselation structure
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool IsInnerPoint(Vector Point, class Tesselation *Tess, LinkedCell* LC) {
-  atom *temp_atom = new atom;
-  bool status = false;
-  temp_atom->x.CopyVector(&Point);
-
-  status = IsInnerAtom(temp_atom, Tess, LC);
-  delete(temp_atom);
-
-  return status;
+    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;
 }
 
@@ -3090 +3167 @@
  * \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
@@ -3137 +3215 @@
     }
   }
-  if (1) { //failflag) {
+  if (filename != 0) {
     *out << Verbose(1) << "Writing final tecplot file\n";
     if (DoTecplotOutput) {
@@ -3170 +3248 @@
     cout << Verbose(2) << "None." << endl;
 
-  // Tests the IsInnerAtom() function.
-  Vector x (0, 0, 0);
-  cout << Verbose(0) << "Point to check: " << x << endl;
-  cout << Verbose(0) << "Check: IsInnerPoint() returns " << IsInnerPoint(x, Tess, LCList)
-    << "for vector " << x << "." << endl;
-  atom* a = Tess->PointsOnBoundary.begin()->second->node;
-  cout << Verbose(0) << "Point to check: " << *a << " (on boundary)." << endl;
-  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
-    << "for atom " << a << " (on boundary)." << endl;
-  LinkedAtoms *List = NULL;
-  for (int i=0;i<NDIM;i++) { // each axis
-    LCList->n[i] = LCList->N[i]-1; // current axis is topmost cell
-    for (LCList->n[(i+1)%NDIM]=0;LCList->n[(i+1)%NDIM]<LCList->N[(i+1)%NDIM];LCList->n[(i+1)%NDIM]++)
-      for (LCList->n[(i+2)%NDIM]=0;LCList->n[(i+2)%NDIM]<LCList->N[(i+2)%NDIM];LCList->n[(i+2)%NDIM]++) {
-        List = LCList->GetCurrentCell();
-        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-        if (List != NULL) {
-          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if (Tess->PointsOnBoundary.find((*Runner)->nr) == Tess->PointsOnBoundary.end()) {
-              a = *Runner;
-              i=3;
-              for (int j=0;j<NDIM;j++)
-                LCList->n[j] = LCList->N[j];
-              break;
-            }
-          }
-        }
-      }
-  }
-  cout << Verbose(0) << "Check: IsInnerAtom() returns " << IsInnerAtom(a, Tess, LCList)
-    << "for atom " << a << " (inside)." << endl;
-
-
   if (freeTess)
     delete(Tess);