Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 0936451cc73a5f07a9a4d57aa947d09abf3898a4)
+++ src/boundary.cpp	(revision 0077b5c0708da6c63d98d637cee5272400fe8565)
@@ -620,4 +620,7 @@
   }
 
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-first");
+
   // First step: RemovePointFromTesselatedSurface
   PointRunner = TesselStruct->PointsOnBoundary.begin();
@@ -641,29 +644,6 @@
   }
 
-
-//  // print all lines
-//  LineRunner = TesselStruct->LinesOnBoundary.begin();
-//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary lines for debugging." << endl;
-//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
-//    LineAdvance++;
-//    line = LineRunner->second;
-//    *out << Verbose(2) << "INFO: Current line is " << *line << "." << endl;
-//    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
-//    LineRunner = LineAdvance;
-//  }
-//
-//  // print all triangles
-//  TriangleRunner = TesselStruct->TrianglesOnBoundary.begin();
-//  TriangleAdvance = TriangleRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary triangles for debugging." << endl;
-//  while (TriangleRunner != TesselStruct->TrianglesOnBoundary.end()) {
-//    TriangleAdvance++;
-//    *out << Verbose(2) << "INFO: Current triangle is " << *(TriangleRunner->second) << "." << endl;
-//    if (TriangleAdvance != TesselStruct->TrianglesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next triangle will be " << *(TriangleAdvance->second) << "." << endl;
-//    TriangleRunner = TriangleAdvance;
-//  }
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-second");
 
   // second step: PickFarthestofTwoBaselines
@@ -674,44 +654,12 @@
     line = LineRunner->second;
     *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
-    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
     // take highest of both lines
-    TesselStruct->PickFarthestofTwoBaselines(out, line);
+    if (TesselStruct->IsConvexRectangle(out, line) == NULL)
+      TesselStruct->PickFarthestofTwoBaselines(out, line);
     LineRunner = LineAdvance;
   }
 
-  // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    point->value = 0;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
-      line = LineRunner->second;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      if (!line->CheckConvexityCriterion(out))
-        point->value += 1;
-    }
-  }
-
-  // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append("_intermed");
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append("_intermed");
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-third");
 
   // third step: IsConvexRectangle
@@ -735,6 +683,22 @@
   }
 
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+
+  // end
+  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
+  return volume;
+};
+
+/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
+ * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
+ * \param *out output stream for debugging
+ * \param *TesselStruct pointer to Tesselation structure
+ */
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+{
+  class BoundaryPointSet *point = NULL;
+  class BoundaryLineSet *line = NULL;
   // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     point = PointRunner->second;
     *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
@@ -747,9 +711,19 @@
     }
   }
-
+};
+
+/** Stores triangles to file.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and bonds
+ * \param *filename prefix of filename
+ * \param *extraSuffix intermediate suffix
+ */
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+{
   // 4. Store triangles in tecplot file
   if (filename != NULL) {
     if (DoTecplotOutput) {
       string OutputName(filename);
+      OutputName.append(extraSuffix);
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
@@ -760,4 +734,5 @@
     if (DoRaster3DOutput) {
       string OutputName(filename);
+      OutputName.append(extraSuffix);
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
@@ -767,8 +742,4 @@
     }
   }
-
-  // end
-  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
-  return volume;
 };
 
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision 0936451cc73a5f07a9a4d57aa947d09abf3898a4)
+++ src/boundary.hpp	(revision 0077b5c0708da6c63d98d637cee5272400fe8565)
@@ -42,4 +42,7 @@
 void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
 Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct);
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix);
+
 
 #endif /*BOUNDARY_HPP_*/
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision 0936451cc73a5f07a9a4d57aa947d09abf3898a4)
+++ src/tesselation.cpp	(revision 0077b5c0708da6c63d98d637cee5272400fe8565)
@@ -1751,5 +1751,5 @@
   class BoundaryPointSet *Spot = NULL;
   class BoundaryLineSet *OtherBase;
-  Vector *ClosestPoint[2];
+  Vector *ClosestPoint;
 
   int m=0;
@@ -1764,31 +1764,30 @@
 
   // get the closest point on each line to the other line
-  ClosestPoint[0] = GetClosestPointBetweenLine(out, Base, OtherBase);
-  ClosestPoint[1] = GetClosestPointBetweenLine(out, OtherBase, Base);
+  ClosestPoint = GetClosestPointBetweenLine(out, Base, OtherBase);
 
   // delete the temporary other base line
+  delete(ClosestPoint);
   delete(OtherBase);
 
   // get the distance vector from Base line to OtherBase line
-  Vector DistanceToIntersection, BaseLine;
+  Vector DistanceToIntersection[2], BaseLine;
+  double distance[2];
   BaseLine.CopyVector(Base->endpoints[1]->node->node);
   BaseLine.SubtractVector(Base->endpoints[0]->node->node);
-  DistanceToIntersection.CopyVector(ClosestPoint[0]);
-  DistanceToIntersection.SubtractVector(Base->endpoints[0]->node->node);
-  Spot = Base->endpoints[1];
-  if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
-    DistanceToIntersection.CopyVector(ClosestPoint[0]);
-    DistanceToIntersection.SubtractVector(Base->endpoints[1]->node->node);
-    Spot = Base->endpoints[0];
-    if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
-      *out << Verbose(3) << "ERROR: Something's very wrong here, both SKPs return negative?!" << endl;
-      return NULL;
-    }
-  }
-  if ((BaseLine.NormSquared() - DistanceToIntersection.NormSquared()) < -MYEPSILON) {  // distance is smaller: concave
-    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave: Base line squared length " << BaseLine.NormSquared() << " against Distance to intersection squared " << DistanceToIntersection.NormSquared() << "." << endl;
+  for (int i=0;i<2;i++) {
+    DistanceToIntersection[i].CopyVector(ClosestPoint);
+    DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node);
+    distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]);
+  }
+  if ((distance[0] * distance[1]) > MYEPSILON)  { // have same sign?
+    *out << Verbose(3) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1]  << ". " << *Base << " is concave." << endl;
+    if (distance[0] < distance[1]) {
+      Spot = Base->endpoints[0];
+    } else {
+      Spot = Base->endpoints[1];
+    }
     return Spot;
-  } else {  // base line is longer: convex
-    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave." << endl;
+  } else {  // different sign, i.e. we are in between
+    *out << Verbose(3) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
     return NULL;
   }
@@ -1796,4 +1795,27 @@
 };
 
+void Tesselation::PrintAllBoundaryPoints(ofstream *out)
+{
+  // print all lines
+  *out << Verbose(1) << "Printing all boundary points for debugging:" << endl;
+  for (PointMap::iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
+    *out << Verbose(2) << *(PointRunner->second) << endl;
+};
+
+void Tesselation::PrintAllBoundaryLines(ofstream *out)
+{
+  // print all lines
+  *out << Verbose(1) << "Printing all boundary lines for debugging:" << endl;
+  for (LineMap::iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
+    *out << Verbose(2) << *(LineRunner->second) << endl;
+};
+
+void Tesselation::PrintAllBoundaryTriangles(ofstream *out)
+{
+  // print all triangles
+  *out << Verbose(1) << "Printing all boundary triangles for debugging:" << endl;
+  for (TriangleMap::iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
+    *out << Verbose(2) << *(TriangleRunner->second) << endl;
+};
 
 /** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
@@ -1826,5 +1848,7 @@
   Distance.SubtractVector(ClosestPoint[0]);
 
-  // delete the temporary other base line
+  // delete the temporary other base line and the closest points
+  delete(ClosestPoint[0]);
+  delete(ClosestPoint[1]);
   delete(OtherBase);
 
@@ -1843,5 +1867,5 @@
       BaseLineNormal.AddVector(&(runner->second->NormalVector));
     }
-    BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+    BaseLineNormal.Scale(1./2.);
 
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
@@ -1866,23 +1890,26 @@
   // construct the plane of the two baselines (i.e. take both their directional vectors)
   Vector Normal;
-  Vector OtherBaseline;
-  Normal.CopyVector(Base->endpoints[1]->node->node);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  Vector Baseline, OtherBaseline;
+  Baseline.CopyVector(Base->endpoints[1]->node->node);
+  Baseline.SubtractVector(Base->endpoints[0]->node->node);
   OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
   OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.CopyVector(&Baseline);
   Normal.VectorProduct(&OtherBaseline);
   Normal.Normalize();
-
-  // project one offset point of OtherBase onto this plane
+  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
+
+  // project one offset point of OtherBase onto this plane (and add plane offset vector)
   Vector NewOffset;
   NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
   NewOffset.ProjectOntoPlane(&Normal);
-  Vector NewDirection;
-  NewDirection.CopyVector(OtherBase->endpoints[1]->node->node);
-  NewDirection.ProjectOntoPlane(&Normal);
+  NewOffset.AddVector(Base->endpoints[0]->node->node);
 
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &OtherBaseline, &Normal);
+  Normal.CopyVector(Intersection);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
 
   return Intersection;
@@ -2938,9 +2965,9 @@
 double getAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
 {
-  if (reference.IsNull())
+  if (reference.IsZero())
     return M_PI;
 
   // calculate both angles and correct with in-plane vector
-  if (point.IsNull())
+  if (point.IsZero())
     return M_PI;
   double phi = point.Angle(&reference);
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision 0936451cc73a5f07a9a4d57aa947d09abf3898a4)
+++ src/tesselation.hpp	(revision 0077b5c0708da6c63d98d637cee5272400fe8565)
@@ -227,4 +227,10 @@
     bool IsInnerPoint(ofstream *out, TesselPoint *Point, LinkedCell* LC);
     bool AddBoundaryPoint(TesselPoint *Walker, int n);
+
+    // print for debugging
+    void PrintAllBoundaryPoints(ofstream *out);
+    void PrintAllBoundaryLines(ofstream *out);
+    void PrintAllBoundaryTriangles(ofstream *out);
+
 
     PointMap PointsOnBoundary;
