Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 02bfdee48bf20f43a0f4e7f6b4c779daa4163c92)
+++ src/boundary.cpp	(revision a98603750c1b5ea8639c4f89919340c98b730445)
@@ -3,4 +3,8 @@
 
 #define DEBUG 1
+#define DoTecplotOutput 0
+#define DoRaster3DOutput 1
+#define TecplotSuffix ".dat"
+#define Raster3DSuffix ".r3d"
 
 // ======================================== Points on Boundary =================================
@@ -25,4 +29,5 @@
   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
   node = NULL;
+  lines.clear();
 }
 ;
@@ -81,21 +86,19 @@
 BoundaryLineSet::~BoundaryLineSet()
 {
-  for (int i = 0; i < 2; i++)
-    {
-      cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point "
-          << *endpoints[i] << "." << endl;
-      endpoints[i]->lines.erase(Nr);
-      LineMap::iterator tester = endpoints[i]->lines.begin();
-      tester++;
-      if (tester == endpoints[i]->lines.end())
-        {
-          cout << Verbose(5) << *endpoints[i]
-              << " has no more lines it's attached to, erasing." << endl;
-          //delete(endpoints[i]);
-        }
-      else
-        cout << Verbose(5) << *endpoints[i]
-            << " has still lines it's attached to." << endl;
-    }
+	for (int i = 0; i < 2; i++) {
+		cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
+		endpoints[i]->lines.erase(Nr);
+		LineMap::iterator tester = endpoints[i]->lines.begin();
+		tester++;
+		if (tester == endpoints[i]->lines.end()) {
+			cout << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;
+			if (endpoints[i] != NULL) {
+				delete(endpoints[i]);
+				endpoints[i] = NULL;
+			} else
+				cerr << "ERROR: Endpoint " << i << " has already been free'd." << endl;
+		} else
+			cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
+	}
 }
 ;
@@ -178,20 +181,17 @@
 BoundaryTriangleSet::~BoundaryTriangleSet()
 {
-  for (int i = 0; i < 3; i++)
-    {
-      cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
-      lines[i]->triangles.erase(Nr);
-      TriangleMap::iterator tester = lines[i]->triangles.begin();
-      tester++;
-      if (tester == lines[i]->triangles.end())
-        {
-          cout << Verbose(5) << *lines[i]
-              << " is no more attached to any triangle, erasing." << endl;
-          delete (lines[i]);
-        }
-      else
-        cout << Verbose(5) << *lines[i] << " is still attached to a triangle."
-            << endl;
-    }
+	for (int i = 0; i < 3; i++) {
+		cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
+		lines[i]->triangles.erase(Nr);
+		if (lines[i]->triangles.empty()) {
+			cout << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;
+			if (lines[i] != NULL) {
+				delete (lines[i]);
+				lines[i] = NULL;
+			} else
+				cerr << "ERROR: This line " << i << " has already been free'd." << endl;
+		} else
+			cout << Verbose(5) << *lines[i] << " is still attached to a triangle." << endl;
+	}
 }
 ;
@@ -565,6 +565,59 @@
 ;
 
-/*
- * This function creates the tecplot file, displaying the tesselation of the hull.
+/** Creates the objects in a raster3d file (renderable with a header.r3d)
+ * \param *out output stream for debugging
+ * \param *tecplot output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
+{
+	atom *Walker = mol->start;
+	bond *Binder = mol->first;
+	int i;
+	Vector *center = mol->DetermineCenterOfAll(out);
+	if (rasterfile != NULL) {
+		//cout << Verbose(1) << "Writing Raster3D file ... ";
+		*rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
+		*rasterfile << "@header.r3d" << endl;
+		*rasterfile << "# All atoms as spheres" << endl;
+		while (Walker->next != mol->end) {
+			Walker = Walker->next;
+			*rasterfile << "2" << endl << "  ";	// 2 is sphere type
+			for (i=0;i<NDIM;i++)
+				*rasterfile << Walker->x.x[i]+center->x[i] << " ";
+			*rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+		}
+
+		*rasterfile << "# All bonds as vertices" << endl;
+		while (Binder->next != mol->last) {
+			Binder = Binder->next;
+			*rasterfile << "3" << endl << "  ";	// 2 is round-ended cylinder type
+			for (i=0;i<NDIM;i++)
+				*rasterfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+			*rasterfile << "\t0.03\t";
+			for (i=0;i<NDIM;i++)
+				*rasterfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+			*rasterfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+		}
+
+		*rasterfile << "# All tesselation triangles" << endl;
+		for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+			*rasterfile << "1" << endl << "  ";	// 1 is triangle type
+			for (i=0;i<3;i++) {	// print each node
+				for (int j=0;j<NDIM;j++)	// and for each node all NDIM coordinates
+					*rasterfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+				*rasterfile << "\t";
+			}
+			*rasterfile << "1. 0. 0." << endl;	// red as colour
+			*rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;	// 18 is transparency type for previous object
+		}
+	} else {
+		cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
+	}
+	delete(center);
+};
+
+/** This function creates the tecplot file, displaying the tesselation of the hull.
  * \param *out output stream for debugging
  * \param *tecplot output stream for tecplot data
@@ -893,10 +946,26 @@
 Tesselation::~Tesselation()
 {
-  cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
-  for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner
-      != TrianglesOnBoundary.end(); runner++)
-    {
-      delete (runner->second);
-    }
+	cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
+	for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
+		if (runner->second != NULL) {
+			delete (runner->second);
+			runner->second = NULL;
+		} else
+			cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+	}
+	for (LineMap::iterator runner = LinesOnBoundary.begin(); runner != LinesOnBoundary.end(); runner++) {
+		if (runner->second != NULL) {
+			delete (runner->second);
+			runner->second = NULL;
+		} else
+			cerr << "ERROR: The line " << runner->first << " has already been free'd." << endl;
+	}
+	for (PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) {
+		if (runner->second != NULL) {
+			delete (runner->second);
+			runner->second = NULL;
+		} else
+			cerr << "ERROR: The point " << runner->first << " has already been free'd." << endl;
+	}
 }
 ;
@@ -1439,5 +1508,15 @@
     double Restradius;
     cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
-    *Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
+    Center->Zero();
+    helper.CopyVector(&a);
+    helper.Scale(sin(2.*alpha));
+    Center->AddVector(&helper);
+    helper.CopyVector(&b);
+    helper.Scale(sin(2.*beta));
+    Center->AddVector(&helper);
+    helper.CopyVector(&c);
+    helper.Scale(sin(2.*gamma));
+    Center->AddVector(&helper);
+    //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
     Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
     NewUmkreismittelpunkt->CopyVector(Center);
@@ -1503,184 +1582,161 @@
  */
 
-void Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
+void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
     atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
 {
-  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
-  /* OldNormal is normal vector on the old triangle
-   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
-   * Chord points from b to a!!!
-   */
-  Vector dif_a; //Vector from a to candidate
-  Vector dif_b; //Vector from b to candidate
-  Vector AngleCheck;
-  Vector TempNormal, Umkreismittelpunkt;
-  Vector Mittelpunkt;
-
-  double CurrentEpsilon = 0.1;
-  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-  double BallAngle, AlternativeSign;
-  atom *Walker; // variable atom point
-
-  Vector NewUmkreismittelpunkt;
-
-
-  if (a != Candidate and b != Candidate and c != Candidate)
-    {
-      cout << Verbose(3) << "We have a unique candidate!" << endl;
-      dif_a.CopyVector(&(a->x));
-      dif_a.SubtractVector(&(Candidate->x));
-      dif_b.CopyVector(&(b->x));
-      dif_b.SubtractVector(&(Candidate->x));
-      AngleCheck.CopyVector(&(Candidate->x));
-      AngleCheck.SubtractVector(&(a->x));
-      AngleCheck.ProjectOntoPlane(Chord);
-
-      SideA = dif_b.Norm();
-      SideB = dif_a.Norm();
-      SideC = Chord->Norm();
-      //Chord->Scale(-1);
-
-      alpha = Chord->Angle(&dif_a);
-      beta = M_PI - Chord->Angle(&dif_b);
-      gamma = dif_a.Angle(&dif_b);
-
-      cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-      if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-        cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-      Umkreisradius = SideA / 2.0 / sin(alpha);
-      //cout << Umkreisradius << endl;
-      //cout << SideB / 2.0 / sin(beta) << endl;
-      //cout << SideC / 2.0 / sin(gamma) << endl;
-
-      if (Umkreisradius < RADIUS) //Checking whether ball will at least rest on points.
-        {
-          cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
-          cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
-          sign = AngleCheck.ScalarProduct(direction1);
-          if (fabs(sign)<MYEPSILON)
-            {
-              if (AngleCheck.ScalarProduct(OldNormal)<0)
-                {
-                  sign =0;
-                  AlternativeSign=1;
-                }
-              else
-                {
-                  sign =0;
-                  AlternativeSign=-1;
-                }
-            }
-          else
-            {
-              sign /= fabs(sign);
-              cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
-            }
-
-          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
-
-          AngleCheck.CopyVector(&ReferencePoint);
-          AngleCheck.Scale(-1);
-          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-          AngleCheck.AddVector(&Mittelpunkt);
-          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-          cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
-
-          BallAngle = AngleCheck.Angle(OldNormal);
-          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-
-          //cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
-          //cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
-
-          cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
-
-          NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
-
-          if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON))
-            {
-              if (Storage[0]< -1.5) // first Candidate at all
-                {
-
-                  cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
-                  Opt_Candidate = Candidate;
-                  Storage[0] = sign;
-                  Storage[1] = AlternativeSign;
-                  Storage[2] = BallAngle;
-                  cout << " angle " << Storage[2] << " and Up/Down "
-                    << Storage[0] << endl;
-                }
-              else
-                {
-                  if ( Storage[2] > BallAngle)
-                    {
-                      cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
-                      Opt_Candidate = Candidate;
-                      Storage[0] = sign;
-                      Storage[1] = AlternativeSign;
-                      Storage[2] = BallAngle;
-                      cout << " angle " << Storage[2] << " and Up/Down "
-                        << Storage[0] << endl;
-                    }
-                  else
-                    {
-                      if (DEBUG)
-                        {
-                          cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
-                        }
-                    }
-                }
-              }
-            else
-              {
-              if (DEBUG)
-                {
-                  cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
-                }
-              }
-          }
-        else
-          {
-            if (DEBUG)
-              {
-                cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
-              }
-          }
-      }
-    else
-      {
-        if (DEBUG)
-          {
-            cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
-          }
-      }
-
-
-
-  if (RecursionLevel < 9) // Seven is the recursion level threshold.
-    {
-      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
-        {
-          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
-              Candidate);
-          if (Walker == Parent)
-            { // don't go back the same bond
-              continue;
-            }
-          else
-            {
-              Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel
-                  + 1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS,
-                  mol); //call function again
-            }
-        }
-    }
-  cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+	cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
+	cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
+	cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
+	cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
+	cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
+	cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
+	/* OldNormal is normal vector on the old triangle
+	 * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
+	 * Chord points from b to a!!!
+	 */
+	Vector dif_a; //Vector from a to candidate
+	Vector dif_b; //Vector from b to candidate
+	Vector AngleCheck;
+	Vector TempNormal, Umkreismittelpunkt;
+	Vector Mittelpunkt;
+
+	double CurrentEpsilon = 0.1;
+	double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+	double BallAngle, AlternativeSign;
+	atom *Walker; // variable atom point
+
+	Vector NewUmkreismittelpunkt;
+
+	if (a != Candidate and b != Candidate and c != Candidate) {
+		cout << Verbose(3) << "We have a unique candidate!" << endl;
+		dif_a.CopyVector(&(a->x));
+		dif_a.SubtractVector(&(Candidate->x));
+		dif_b.CopyVector(&(b->x));
+		dif_b.SubtractVector(&(Candidate->x));
+		AngleCheck.CopyVector(&(Candidate->x));
+		AngleCheck.SubtractVector(&(a->x));
+		AngleCheck.ProjectOntoPlane(Chord);
+
+		SideA = dif_b.Norm();
+		SideB = dif_a.Norm();
+		SideC = Chord->Norm();
+		//Chord->Scale(-1);
+
+		alpha = Chord->Angle(&dif_a);
+		beta = M_PI - Chord->Angle(&dif_b);
+		gamma = dif_a.Angle(&dif_b);
+
+		cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+
+		if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
+			cerr << Verbose(0) << "WARNING: sum of angles for base triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+			cout << Verbose(1) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+		}
+
+		if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
+			Umkreisradius = SideA / 2.0 / sin(alpha);
+			//cout << Umkreisradius << endl;
+			//cout << SideB / 2.0 / sin(beta) << endl;
+			//cout << SideC / 2.0 / sin(gamma) << endl;
+
+			if (Umkreisradius < RADIUS) { //Checking whether ball will at least rest on points.
+				cout << Verbose(3) << "Circle of circumference would fit: " << Umkreisradius << " < " << RADIUS << "." << endl;
+				cout << Verbose(2) << "Candidate is "<< *Candidate << endl;
+				sign = AngleCheck.ScalarProduct(direction1);
+				if (fabs(sign)<MYEPSILON) {
+					if (AngleCheck.ScalarProduct(OldNormal)<0) {
+						sign =0;
+						AlternativeSign=1;
+					} else {
+						sign =0;
+						AlternativeSign=-1;
+					}
+				} else {
+					sign /= fabs(sign);
+					cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+				}
+
+				Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+
+				AngleCheck.CopyVector(&ReferencePoint);
+				AngleCheck.Scale(-1);
+				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+				AngleCheck.AddVector(&Mittelpunkt);
+				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+				cout << Verbose(4) << "Reference vector to sphere's center is " << AngleCheck << "." << endl;
+
+				BallAngle = AngleCheck.Angle(OldNormal);
+				cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+
+				//cout << "direction1 is " << direction1->x[0] <<" "<< direction1->x[1] <<" "<< direction1->x[2] <<" " << endl;
+				//cout << "AngleCheck is " << AngleCheck.x[0] << " "<< AngleCheck.x[1] << " "<< AngleCheck.x[2] << " "<< endl;
+
+				cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
+
+				NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
+
+				if ((AngleCheck.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
+					if (Storage[0]< -1.5) { // first Candidate at all
+						if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+							cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
+							Opt_Candidate = Candidate;
+							Storage[0] = sign;
+							Storage[1] = AlternativeSign;
+							Storage[2] = BallAngle;
+							cout << " angle " << Storage[2] << " and Up/Down "
+							<< Storage[0] << endl;
+						} else
+							cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+					} else {
+						if ( Storage[2] > BallAngle) {
+							if (1) { //if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+								cout << Verbose(2) << "Next better candidate is " << *Candidate << " with ";
+								Opt_Candidate = Candidate;
+								Storage[0] = sign;
+								Storage[1] = AlternativeSign;
+								Storage[2] = BallAngle;
+								cout << " angle " << Storage[2] << " and Up/Down "
+								<< Storage[0] << endl;
+							} else
+								cout << "Candidate " << *Candidate << " does not belong to a valid triangle." << endl;
+						} else {
+							if (DEBUG) {
+								cout << Verbose(3) << *Candidate << " looses against better candidate " << *Opt_Candidate << "." << endl;
+							}
+						}
+					}
+				} else {
+					if (DEBUG) {
+						cout << Verbose(3) << *Candidate << " refused due to Up/Down sign which is " << sign << endl;
+					}
+				}
+			} else {
+				if (DEBUG) {
+					cout << Verbose(3) << *Candidate << " would have circumference of " << Umkreisradius << " bigger than ball's radius " << RADIUS << "." << endl;
+				}
+			}
+		} else {
+			if (DEBUG) {
+				cout << Verbose(0) << "Triangle consisting of " << *Candidate << ", " << *a << " and " << *b << " has an angle >150!" << endl;
+			}
+		}
+	} else {
+		if (DEBUG) {
+			cout << Verbose(3) << *Candidate << " is either " << *a << " or " << *b << "." << endl;
+		}
+	}
+
+	if (RecursionLevel < 5) { // Seven is the recursion level threshold.
+		for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+			if (Walker == Parent) { // don't go back the same bond
+				continue;
+			} else {
+				Find_next_suitable_point_via_Angle_of_Sphere(a, b, c, Walker, Candidate, RecursionLevel+1, Chord, direction1, OldNormal, ReferencePoint, Opt_Candidate, Storage, RADIUS, mol); //call function again
+			}
+		}
+	}
+	cout << Verbose(2) << "End of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
 }
 ;
@@ -1709,167 +1765,147 @@
    * @param mol molecule structure with atoms and bonds
    */
-
 void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
     atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
 {
-  /* OldNormal is normal vector on the old triangle
-   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
-   * Chord points from b to a!!!
-   */
-  Vector dif_a; //Vector from a to candidate
-  Vector dif_b; //Vector from b to candidate
-  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
-  Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
-
-  double CurrentEpsilon = 0.1;
-  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-  double BallAngle;
-  atom *Walker; // variable atom point
-
-
-  dif_a.CopyVector(&(a->x));
-  dif_a.SubtractVector(&(Candidate->x));
-  dif_b.CopyVector(&(b->x));
-  dif_b.SubtractVector(&(Candidate->x));
-  DirectionCheckPoint.CopyVector(&dif_a);
-  DirectionCheckPoint.Scale(-1);
-  DirectionCheckPoint.ProjectOntoPlane(Chord);
-
-  SideA = dif_b.Norm();
-  SideB = dif_a.Norm();
-  SideC = Chord->Norm();
-  //Chord->Scale(-1);
-
-  alpha = Chord->Angle(&dif_a);
-  beta = M_PI - Chord->Angle(&dif_b);
-  gamma = dif_a.Angle(&dif_b);
-
-
-  if (DEBUG)
-    {
-      cout << "Atom number" << Candidate->nr << endl;
-      Candidate->x.Output((ofstream *) &cout);
-      cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr]
-          << endl;
-    }
-
-  if (a != Candidate and b != Candidate)
-    {
-      //      alpha = dif_a.Angle(&dif_b) / 2.;
-      //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-      //      SideB = dif_a.Norm();
-      //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-      //          alpha); // note this is squared of center line length
-      //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-      // Those are remains from Freddie. Needed?
-
-
-
-      Umkreisradius = SideA / 2.0 / sin(alpha);
-      //cout << Umkreisradius << endl;
-      //cout << SideB / 2.0 / sin(beta) << endl;
-      //cout << SideC / 2.0 / sin(gamma) << endl;
-
-      if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) //Checking whether ball will at least rest o points.
-        {
-
-          // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
-
-          Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-          Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
-          TempNormal.CopyVector(&dif_a);
-          TempNormal.VectorProduct(&dif_b);
-          if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0)
-            {
-              TempNormal.Scale(-1);
-            }
-          TempNormal.Normalize();
-          Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-          TempNormal.Scale(Restradius);
-
-          Mittelpunkt.CopyVector(&Umkreismittelpunkt);
-          Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-
-          AngleCheck.CopyVector(Chord);
-          AngleCheck.Scale(-0.5);
-          AngleCheck.SubtractVector(&(b->x));
-          AngleCheckReference.CopyVector(&AngleCheck);
-          AngleCheckReference.AddVector(Opt_Mittelpunkt);
-          AngleCheck.AddVector(&Mittelpunkt);
-
-          BallAngle = AngleCheck.Angle(&AngleCheckReference);
-
-          d1->ProjectOntoPlane(&AngleCheckReference);
-          sign = AngleCheck.ScalarProduct(d1);
-          if (fabs(sign) < MYEPSILON)
-            sign = 0;
-          else
-            sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-
-          if (Storage[0]< -1.5) // first Candidate at all
-            {
-
-              cout << "Next better candidate is " << *Candidate << " with ";
-              Opt_Candidate = Candidate;
-              Storage[0] = sign;
-              Storage[1] = BallAngle;
-              Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-              cout << "Angle is " << Storage[1] << ", Halbraum ist "
-                  << Storage[0] << endl;
-
-
-            }
-          else
-            {
-              /*
-               * removed due to change in criterium, now checking angle of ball to old normal.
-              //We will now check for non interference, that is if the new candidate would have the Opt_Candidate
-              //within the ball.
-
-              Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-              //cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
-
-
-              if (Distance >RADIUS) // We have no interference and may now check whether the new point is better.
-                */
-                {
-                  //cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
-
-                  if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon)))  //This will give absolute preference to those in "right-hand" quadrants
-                      //(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
-                    {
-                      cout << "Next better candidate is " << *Candidate << " with ";
-                      Opt_Candidate = Candidate;
-                      Storage[0] = sign;
-                      Storage[1] = BallAngle;
-                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
-                          << Storage[0] << endl;
-
-
-                    }
-                  else
-                    {
-                      if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0
-                          && Storage[1] > BallAngle) ||
-                          (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0
-                          && Storage[1] < BallAngle))
-                      //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
-                        {
-                          cout << "Next better candidate is " << *Candidate << " with ";
-                          Opt_Candidate = Candidate;
-                          Storage[0] = sign;
-                          Storage[1] = BallAngle;
-                          Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
-                          cout << "Angle is " << Storage[1] << ", Halbraum ist "
-                              << Storage[0] << endl;
-                        }
-
-                    }
-                }
-              /*
+	/* OldNormal is normal vector on the old triangle
+	 * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
+	 * Chord points from b to a!!!
+	 */
+	Vector dif_a; //Vector from a to candidate
+	Vector dif_b; //Vector from b to candidate
+	Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
+	Vector TempNormal, Umkreismittelpunkt, Mittelpunkt, helper;
+
+	double CurrentEpsilon = 0.1;
+	double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+	double BallAngle;
+	atom *Walker; // variable atom point
+
+
+	dif_a.CopyVector(&(a->x));
+	dif_a.SubtractVector(&(Candidate->x));
+	dif_b.CopyVector(&(b->x));
+	dif_b.SubtractVector(&(Candidate->x));
+	DirectionCheckPoint.CopyVector(&dif_a);
+	DirectionCheckPoint.Scale(-1);
+	DirectionCheckPoint.ProjectOntoPlane(Chord);
+
+	SideA = dif_b.Norm();
+	SideB = dif_a.Norm();
+	SideC = Chord->Norm();
+	//Chord->Scale(-1);
+
+	alpha = Chord->Angle(&dif_a);
+	beta = M_PI - Chord->Angle(&dif_b);
+	gamma = dif_a.Angle(&dif_b);
+
+
+	if (DEBUG) {
+		cout << "Atom number" << Candidate->nr << endl;
+		Candidate->x.Output((ofstream *) &cout);
+		cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
+	}
+
+	if (a != Candidate and b != Candidate) {
+		//      alpha = dif_a.Angle(&dif_b) / 2.;
+		//      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
+		//      SideB = dif_a.Norm();
+		//      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
+		//          alpha); // note this is squared of center line length
+		//      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
+		// Those are remains from Freddie. Needed?
+
+		Umkreisradius = SideA / 2.0 / sin(alpha);
+		//cout << Umkreisradius << endl;
+		//cout << SideB / 2.0 / sin(beta) << endl;
+		//cout << SideC / 2.0 / sin(gamma) << endl;
+
+		if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
+			// intermediate calculations to aquire centre of sphere, called Mittelpunkt:
+			Umkreismittelpunkt.Zero();
+			helper.CopyVector(&a->x);
+			helper.Scale(sin(2.*alpha));
+			Umkreismittelpunkt.AddVector(&helper);
+			helper.CopyVector(&b->x);
+			helper.Scale(sin(2.*beta));
+			Umkreismittelpunkt.AddVector(&helper);
+			helper.CopyVector(&Candidate->x);
+			helper.Scale(sin(2.*gamma));
+			Umkreismittelpunkt.AddVector(&helper);
+			//Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
+			Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+
+			TempNormal.CopyVector(&dif_a);
+			TempNormal.VectorProduct(&dif_b);
+			if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
+				TempNormal.Scale(-1);
+			}
+			TempNormal.Normalize();
+			Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+			TempNormal.Scale(Restradius);
+
+			Mittelpunkt.CopyVector(&Umkreismittelpunkt);
+			Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
+
+			AngleCheck.CopyVector(Chord);
+			AngleCheck.Scale(-0.5);
+			AngleCheck.SubtractVector(&(b->x));
+			AngleCheckReference.CopyVector(&AngleCheck);
+			AngleCheckReference.AddVector(Opt_Mittelpunkt);
+			AngleCheck.AddVector(&Mittelpunkt);
+
+			BallAngle = AngleCheck.Angle(&AngleCheckReference);
+
+			d1->ProjectOntoPlane(&AngleCheckReference);
+			sign = AngleCheck.ScalarProduct(d1);
+			if (fabs(sign) < MYEPSILON)
+				sign = 0;
+			else
+				sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
+
+
+			if (Storage[0]< -1.5) { // first Candidate at all
+				cout << "Next better candidate is " << *Candidate << " with ";
+				Opt_Candidate = Candidate;
+				Storage[0] = sign;
+				Storage[1] = BallAngle;
+				Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+				cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+			} else {
+				/*
+				 * removed due to change in criterium, now checking angle of ball to old normal.
+				//We will now check for non interference, that is if the new candidate would have the Opt_Candidate
+				//within the ball.
+
+				Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
+				//cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+
+
+				if (Distance >RADIUS) { // We have no interference and may now check whether the new point is better.
+				 */
+					//cout << "Atom " << Candidate << " has distance " << Candidate->x.Distance(Opt_Mittelpunkt) << " to Center of Candidate " << endl;
+
+				if (((Storage[0] < 0 && fabs(sign - Storage[0]) > CurrentEpsilon))) { //This will give absolute preference to those in "right-hand" quadrants
+					//(Candidate->x.Distance(Opt_Mittelpunkt) < RADIUS))    //and those where Candidate would be within old Sphere.
+					cout << "Next better candidate is " << *Candidate << " with ";
+					Opt_Candidate = Candidate;
+					Storage[0] = sign;
+					Storage[1] = BallAngle;
+					Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+					cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+				} else {
+					if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0 && Storage[1] > BallAngle) || (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0 && Storage[1] < BallAngle)) {
+						//Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+						cout << "Next better candidate is " << *Candidate << " with ";
+						Opt_Candidate = Candidate;
+						Storage[0] = sign;
+						Storage[1] = BallAngle;
+						Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+						cout << "Angle is " << Storage[1] << ", Halbraum ist " << Storage[0] << endl;
+					}
+				}
+			}
+/*
                * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
                *
@@ -1913,42 +1949,26 @@
                 }
                 */
-            }
-        }
-      else
-        {
-          if (DEBUG)
-            {
-              cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
-            }
-        }
-    }
-
-  else
-    {
-      if (DEBUG)
-        cout << "identisch mit Ursprungslinie" << endl;
-    }
-
-  if (RecursionLevel < 9) // Five is the recursion level threshold.
-    {
-      for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
-        {
-          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
-              Candidate);
-          if (Walker == Parent)
-            { // don't go back the same bond
-              continue;
-            }
-          else
-            {
-              Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel
-                  + 1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS,
-                  mol); //call function again
-
-            }
-        }
-    }
-}
-;
+		} else {
+			if (DEBUG) {
+				cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
+			}
+		}
+	} else {
+		if (DEBUG) {
+			cout << "identisch mit Ursprungslinie" << endl;
+		}
+	}
+
+	if (RecursionLevel < 9) { // Five is the recursion level threshold.
+		for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) { // go through all bond
+			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+			if (Walker == Parent) { // don't go back the same bond
+				continue;
+			} else {
+				Find_next_suitable_point(a, b, Walker, Candidate, RecursionLevel+1, Chord, d1, OldNormal, Opt_Candidate, Opt_Mittelpunkt, Storage, RADIUS, mol); //call function again
+			}
+		}
+	}
+};
 
 /** This function finds a triangle to a line, adjacent to an existing one.
@@ -1960,176 +1980,219 @@
  * @param RADIUS radius of the rolling ball
  * @param N number of found triangles
+ * @param *filename filename base for intermediate envelopes
  */
-void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
+bool Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
     molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    const double& RADIUS, int N)
-{
-  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
-  Vector direction1;
-  Vector helper;
-  Vector Chord;
-  ofstream *tempstream = NULL;
-  char filename[255];
-  double tmp;
-  //atom* Walker;
-  atom* OldThirdPoint;
-
-  double Storage[3];
-  Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
-  Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
-  Storage[2] = 9999999.;
-  atom* Opt_Candidate = NULL;
-  Vector Opt_Mittelpunkt;
-
-  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
-
-  helper.CopyVector(&(Line.endpoints[0]->node->x));
-  for (int i = 0; i < 3; i++)
-    {
-      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr
-          && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr)
-        {
-          OldThirdPoint = T.endpoints[i]->node;
-          helper.SubtractVector(&T.endpoints[i]->node->x);
-          break;
-        }
-    }
-
-  direction1.CopyVector(&Line.endpoints[0]->node->x);
-  direction1.SubtractVector(&Line.endpoints[1]->node->x);
-  direction1.VectorProduct(&(T.NormalVector));
-
-  if (direction1.ScalarProduct(&helper) < 0)
-    {
-      direction1.Scale(-1);
-    }
-  cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
-
-  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
-  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
-  cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
-
-
-  Vector Umkreismittelpunkt, a, b, c;
-  double alpha, beta, gamma;
-  a.CopyVector(&(T.endpoints[0]->node->x));
-  b.CopyVector(&(T.endpoints[1]->node->x));
-  c.CopyVector(&(T.endpoints[2]->node->x));
-  a.SubtractVector(&(T.endpoints[1]->node->x));
-  b.SubtractVector(&(T.endpoints[2]->node->x));
-  c.SubtractVector(&(T.endpoints[0]->node->x));
-
-//  alpha = a.Angle(&c) - M_PI/2.;
-//  beta = b.Angle(&a);
-//  gamma = c.Angle(&b) - M_PI/2.;
-    alpha = M_PI - a.Angle(&c);
-    beta = M_PI - b.Angle(&a);
-    gamma = M_PI - c.Angle(&b);
-    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
-      cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
-
-    Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-    //cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-    Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-    cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
-    cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
-    cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
-
-  cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-  Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
-  Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-  cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
-  if (DEBUG)
-    cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
-  tmp = 0;
-  for (int i=0;i<NDIM;i++) {
-    helper.CopyVector(&T.endpoints[i]->node->x);
-    helper.SubtractVector(&Umkreismittelpunkt);
-    if (DEBUG)
-      cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
-    if (tmp == 0) // set on first time for comparison against next ones
-      tmp = helper.Norm();
-    if (fabs(helper.Norm() - tmp) > MYEPSILON)
-      cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
-  }
-
-  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
-
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
-      Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
-      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-  Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint,
-      Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
-      &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
-
-
-      cout << "Letzter Winkel bei " << TrianglesOnBoundaryCount << " Winkel ist " << Storage[2] << endl;
-
-
-  if ((TrianglesOnBoundaryCount % 1) == 0)
-    {
-    cout << Verbose(1) << "Writing temporary non convex hull to file testEnvelope-" << TriangleFilesWritten << "d.dat.\n";
-    sprintf(filename, "testEnvelope-%d.dat", TriangleFilesWritten);
-    tempstream = new ofstream(filename, ios::trunc);
-    write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten++);
-    tempstream->close();
-    tempstream->flush();
-    delete(tempstream);
-  }
-
-  if (TrianglesOnBoundaryCount >189 and TrianglesOnBoundaryCount < 200 )
-    {
-      cerr << Verbose(0)
-          << "WARNING: Already more than " << TrianglesOnBoundaryCount-1 << "triangles, construction probably has crashed." << endl;
-      write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
-      cout << Verbose(2) << "This is currently added candidate" << Opt_Candidate << endl;
-    }
-  // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
-
-  cout << Verbose(2) << " Optimal candidate is " << *Opt_Candidate << endl;
-
-  // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
-
-  AddTrianglePoint(Opt_Candidate, 0);
-  LineMap::iterator TryAndError;
-  TryAndError = TPS[0]->lines.find(Line.endpoints[0]->node->nr);
-  bool flag = true;
-  if ((TryAndError != TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
-    flag = false;
-  TryAndError = TPS[0]->lines.find(Line.endpoints[1]->node->nr);
-  if ((TryAndError !=  TPS[0]->lines.end()) && (TryAndError->second->TrianglesCount > 1))
-    flag = false;
-
-  if (flag) { // if so, add
-    AddTrianglePoint(Line.endpoints[0]->node, 1);
-    AddTrianglePoint(Line.endpoints[1]->node, 2);
-
-    AddTriangleLine(TPS[0], TPS[1], 0);
-    AddTriangleLine(TPS[0], TPS[2], 1);
-    AddTriangleLine(TPS[1], TPS[2], 2);
-
-    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-    AddTriangleToLines();
-
-    BTS->GetNormalVector(BTS->NormalVector);
-
-    if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
-        (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) ||
-        (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) )
-
-      {
-        BTS->NormalVector.Scale(-1);
-      };
-    cout << Verbose(2) << "New triangle with " << *BTS << "and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
-    cout << Verbose(2) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
-  } else { // else, yell and do nothing
-    cout << Verbose(2) << "This triangle consisting of ";
-    for (int i=0;i<NDIM;i++)
-      cout << BLS[i] << "(" << BLS[i]->TrianglesCount << "), ";
-    cout << " is invalid!" << endl;
-  }
-  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
-}
-;
+    const double& RADIUS, int N, const char *tempbasename)
+{
+	cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+	Vector direction1;
+	Vector helper;
+	Vector Chord;
+	ofstream *tempstream = NULL;
+	char NumberName[255];
+	double tmp;
+	//atom* Walker;
+	atom* OldThirdPoint;
+
+	double Storage[3];
+	Storage[0] = -2.; // This direction is either +1 or -1 one, so any result will take precedence over initial values
+	Storage[1] = -2.; // This is also lower then any value produced by an eligible atom, which are all positive
+	Storage[2] = 9999999.;
+	atom* Opt_Candidate = NULL;
+	Vector Opt_Mittelpunkt;
+
+	cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+
+	helper.CopyVector(&(Line.endpoints[0]->node->x));
+	for (int i = 0; i < 3; i++) {
+		if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr != Line.endpoints[1]->node->nr) {
+			OldThirdPoint = T.endpoints[i]->node;
+			helper.SubtractVector(&T.endpoints[i]->node->x);
+			break;
+		}
+	}
+
+	direction1.CopyVector(&Line.endpoints[0]->node->x);
+	direction1.SubtractVector(&Line.endpoints[1]->node->x);
+	direction1.VectorProduct(&(T.NormalVector));
+
+	if (direction1.ScalarProduct(&helper) < 0) {
+		direction1.Scale(-1);
+	}
+	cout << Verbose(2) << "Looking in direction " << direction1 << " for candidates.\n";
+
+	Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+	Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+	cout << Verbose(2) << "Baseline vector is " << Chord << ".\n";
+
+
+	Vector Umkreismittelpunkt, a, b, c;
+	double alpha, beta, gamma;
+	a.CopyVector(&(T.endpoints[0]->node->x));
+	b.CopyVector(&(T.endpoints[1]->node->x));
+	c.CopyVector(&(T.endpoints[2]->node->x));
+	a.SubtractVector(&(T.endpoints[1]->node->x));
+	b.SubtractVector(&(T.endpoints[2]->node->x));
+	c.SubtractVector(&(T.endpoints[0]->node->x));
+
+	alpha = M_PI - a.Angle(&c);
+	beta = M_PI - b.Angle(&a);
+	gamma = M_PI - c.Angle(&b);
+	if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON)
+		cerr << Verbose(0) << "WARNING: sum of angles for candidate triangle " << (alpha + beta + gamma)/M_PI*180. << " != 180.\n";
+
+	Umkreismittelpunkt.Zero();
+	helper.CopyVector(&T.endpoints[0]->node->x);
+	helper.Scale(sin(2.*alpha));
+	Umkreismittelpunkt.AddVector(&helper);
+	helper.CopyVector(&T.endpoints[1]->node->x);
+	helper.Scale(sin(2.*beta));
+	Umkreismittelpunkt.AddVector(&helper);
+	helper.CopyVector(&T.endpoints[2]->node->x);
+	helper.Scale(sin(2.*gamma));
+	Umkreismittelpunkt.AddVector(&helper);
+	//Umkreismittelpunkt = (T.endpoints[0]->node->x) * sin(2.*alpha) + T.endpoints[1]->node->x * sin(2.*beta) + (T.endpoints[2]->node->x) * sin(2.*gamma) ;
+	//cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+	Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
+	//cout << "UmkreisMittelpunkt is " << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] << " "<< Umkreismittelpunkt.x[2] << " "<< endl;
+	cout << " We look over line " << Line << " in direction " << direction1.x[0] << " " << direction1.x[1] << " " << direction1.x[2] << " " << endl;
+	cout << " Old Normal is " <<  (T.NormalVector.x)[0] << " " << T.NormalVector.x[1] << " " << (T.NormalVector.x)[2] << " " << endl;
+
+	cout << Verbose(2) << "Base triangle has sides " << a << ", " << b << ", " << c << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
+	cout << Verbose(2) << "Center of circumference is " << Umkreismittelpunkt << "." << endl;
+	if (DEBUG)
+		cout << Verbose(3) << "Check of relative endpoints (same distance, equally spreaded): "<< endl;
+	tmp = 0;
+	for (int i=0;i<NDIM;i++) {
+		helper.CopyVector(&T.endpoints[i]->node->x);
+		helper.SubtractVector(&Umkreismittelpunkt);
+		if (DEBUG)
+			cout << Verbose(3) << "Endpoint[" << i << "]: " << helper << " with length " << helper.Norm() << "." << endl;
+		if (tmp == 0) // set on first time for comparison against next ones
+			tmp = helper.Norm();
+		if (fabs(helper.Norm() - tmp) > MYEPSILON)
+			cerr << Verbose(1) << "WARNING: center of circumference is wrong!" << endl;
+	}
+
+	cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+
+	Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+	Find_next_suitable_point_via_Angle_of_Sphere(Line.endpoints[0]->node, Line.endpoints[1]->node, OldThirdPoint, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Umkreismittelpunkt, Opt_Candidate, Storage, RADIUS, mol);
+	if (Opt_Candidate == NULL) {
+		cerr << "WARNING: Could not find a suitable candidate." << endl;
+		return false;
+	}
+	cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << endl;
+
+	// check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
+	bool flag = CheckPresenceOfTriangle(out, Opt_Candidate, Line.endpoints[0]->node, Line.endpoints[1]->node);
+
+	if (flag) { // if so, add
+		AddTrianglePoint(Opt_Candidate, 0);
+		AddTrianglePoint(Line.endpoints[0]->node, 1);
+		AddTrianglePoint(Line.endpoints[1]->node, 2);
+
+		AddTriangleLine(TPS[0], TPS[1], 0);
+		AddTriangleLine(TPS[0], TPS[2], 1);
+		AddTriangleLine(TPS[1], TPS[2], 2);
+
+		BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+		AddTriangleToLines();
+
+		BTS->GetNormalVector(BTS->NormalVector);
+
+		if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  || (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0) || (fabs(Storage[0]) < MYEPSILON && Storage[1]*BTS->NormalVector.ScalarProduct(&direction1) < 0) ) {
+			BTS->NormalVector.Scale(-1);
+		};
+		cout << Verbose(1) << "New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl;
+		cout << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << Line << "." << endl;
+	} else { // else, yell and do nothing
+		cout << Verbose(1) << "This triangle consisting of ";
+		cout << *Opt_Candidate << ", ";
+		cout << *Line.endpoints[0]->node << " and ";
+		cout << *Line.endpoints[1]->node << " ";
+		cout << "is invalid!" << endl;
+		return false;
+	}
+
+	if ((TrianglesOnBoundaryCount % 10) == 0) {
+		sprintf(NumberName, "-%d", TriangleFilesWritten);
+		if (DoTecplotOutput) {
+			string NameofTempFile(tempbasename);
+			NameofTempFile.append(NumberName);
+			NameofTempFile.append(TecplotSuffix);
+			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+			write_tecplot_file(out, tempstream, this, mol, TriangleFilesWritten);
+			tempstream->close();
+			tempstream->flush();
+			delete(tempstream);
+		}
+
+		if (DoRaster3DOutput) {
+			string NameofTempFile(tempbasename);
+			NameofTempFile.append(NumberName);
+			NameofTempFile.append(Raster3DSuffix);
+			cout << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
+			tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
+			write_raster3d_file(out, tempstream, this, mol);
+			tempstream->close();
+			tempstream->flush();
+			delete(tempstream);
+		}
+		if (DoTecplotOutput || DoRaster3DOutput)
+			TriangleFilesWritten++;
+	}
+
+	cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+	return true;
+};
+
+/** Checks whether the triangle consisting of the three atoms is already present.
+ * Searches for the points in Tesselation::PointsOnBoundary and checks their
+ * lines. If any of the three edges already has two triangles attached, false is
+ * returned.
+ * \param *out output stream for debugging
+ * \param *a first endpoint
+ * \param *b second endpoint
+ * \param *c third endpoint
+ * \return false - triangle invalid due to edge criteria, true - triangle may be added.
+ */
+bool Tesselation::CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c) {
+	LineMap::iterator TryAndError;
+	PointTestPair InsertPoint;
+	bool Present[3];
+
+	*out << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;
+	TPS[0] = new class BoundaryPointSet(a);
+	TPS[1] = new class BoundaryPointSet(b);
+	TPS[2] = new class BoundaryPointSet(c);
+	for (int i=0;i<3;i++) { // check through all endpoints
+		InsertPoint = PointsOnBoundary.insert(PointPair(TPS[i]->node->nr, TPS[i]));
+		Present[i] = !InsertPoint.second;
+		if (Present[i]) { // if new point was not present before, increase counter
+			delete TPS[i];
+			*out << Verbose(2) << "Atom " << *((InsertPoint.first)->second->node) << " gibt's schon in der PointMap." << endl;
+			TPS[i] = (InsertPoint.first)->second;
+		}
+	}
+
+	// check lines
+	for (int i=0;i<3;i++)
+		if (Present[i])
+			for (int j=i;j<3;j++)
+				if (Present[j]) {
+					TryAndError = TPS[i]->lines.find(TPS[j]->node->nr);
+					if ((TryAndError != TPS[i]->lines.end()) && (TryAndError->second->TrianglesCount > 1)) {
+						*out << "WARNING: Line " << *TryAndError->second << " already present with " << TryAndError->second->TrianglesCount << " triangles attached." << endl;
+						*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+						return false;
+					}
+				}
+	*out << Verbose(2) << "End of CheckPresenceOfTriangle" << endl;
+	return true;
+};
+
 
 void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
@@ -2137,217 +2200,192 @@
     molecule* mol, double RADIUS)
 {
-  cout << Verbose(2)
-      << "Begin of Find_second_point_for_Tesselation, recursive level "
-      << RecursionLevel << endl;
-  int i;
-  Vector AngleCheck;
-  atom* Walker;
-  double norm = -1., angle;
-
-  // check if we only have one unique point yet ...
-  if (a != Candidate)
-    {
-    cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
-      AngleCheck.CopyVector(&(Candidate->x));
-      AngleCheck.SubtractVector(&(a->x));
-      norm = AngleCheck.Norm();
-      // second point shall have smallest angle with respect to Oben vector
-      if (norm < RADIUS)
-        {
-        angle = AngleCheck.Angle(&Oben);
-          if (angle < Storage[0])
-            {
-              //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
-              cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
-              Opt_Candidate = Candidate;
-              Storage[0] = AngleCheck.Angle(&Oben);
-              //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
-            }
-          else
-            {
-              cout << "Looses with angle " << angle << " to a better candidate "
-                  << *Opt_Candidate << endl;
-            }
-        }
-      else
-        {
-          cout << "Refused due to Radius " << norm
-              << endl;
-        }
-    }
-
-  // if not recursed to deeply, look at all its bonds
-  if (RecursionLevel < 7)
-    {
-      for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
-        {
-          Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(
-              Candidate);
-          if (Walker == Parent) // don't go back along the bond we came from
-            continue;
-          else
-            Find_second_point_for_Tesselation(a, Walker, Candidate,
-                RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
-        };
-    };
-  cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level "
-      << RecursionLevel << endl;
-}
-;
+	cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+	int i;
+	Vector AngleCheck;
+	atom* Walker;
+	double norm = -1., angle;
+
+	// check if we only have one unique point yet ...
+	if (a != Candidate) {
+		cout << Verbose(3) << "Current candidate is " << *Candidate << ": ";
+		AngleCheck.CopyVector(&(Candidate->x));
+		AngleCheck.SubtractVector(&(a->x));
+		norm = AngleCheck.Norm();
+		// second point shall have smallest angle with respect to Oben vector
+		if (norm < RADIUS) {
+			angle = AngleCheck.Angle(&Oben);
+			if (angle < Storage[0]) {
+				//cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+				cout << "Is a better candidate with distance " << norm << " and " << angle << ".\n";
+				Opt_Candidate = Candidate;
+				Storage[0] = AngleCheck.Angle(&Oben);
+				//cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
+			} else {
+				cout << "Looses with angle " << angle << " to a better candidate " << *Opt_Candidate << endl;
+			}
+		} else {
+			cout << "Refused due to Radius " << norm << endl;
+		}
+	}
+
+	// if not recursed to deeply, look at all its bonds
+	if (RecursionLevel < 7) {
+		for (i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) {
+			Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
+			if (Walker == Parent) // don't go back along the bond we came from
+				continue;
+			else
+				Find_second_point_for_Tesselation(a, Walker, Candidate, RecursionLevel + 1, Oben, Opt_Candidate, Storage, mol, RADIUS);
+		}
+	}
+	cout << Verbose(2) << "End of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
+};
 
 void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
 {
-  cout << Verbose(1) << "Begin of Find_starting_triangle\n";
-  int i = 0;
-  atom* Walker;
-  atom* FirstPoint;
-  atom* SecondPoint;
-  atom* max_index[NDIM];
-  double max_coordinate[NDIM];
-  Vector Oben;
-  Vector helper;
-  Vector Chord;
-  Vector CenterOfFirstLine;
-
-  Oben.Zero();
-
-  for (i = 0; i < 3; i++)
-    {
-      max_index[i] = NULL;
-      max_coordinate[i] = -1;
-    }
-  cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount
-      << " Atoms \n";
-
-  // 1. searching topmost atom with respect to each axis
-  Walker = mol->start;
-  while (Walker->next != mol->end)
-    {
-      Walker = Walker->next;
-      for (i = 0; i < 3; i++)
-        {
-          if (Walker->x.x[i] > max_coordinate[i])
-            {
-              max_coordinate[i] = Walker->x.x[i];
-              max_index[i] = Walker;
-            }
-        }
-    }
-
-  cout << Verbose(2) << "Found maximum coordinates: ";
-  for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *max_index[i] << "\t";
-  cout << endl;
+	cout << Verbose(1) << "Begin of Find_starting_triangle\n";
+	int i = 0;
+	atom* Walker;
+	atom* FirstPoint;
+	atom* SecondPoint;
+	atom* max_index[NDIM];
+	double max_coordinate[NDIM];
+	Vector Oben;
+	Vector helper;
+	Vector Chord;
+	Vector CenterOfFirstLine;
+
+	Oben.Zero();
+
+	for (i = 0; i < 3; i++) {
+		max_index[i] = NULL;
+		max_coordinate[i] = -1;
+	}
+	cout << Verbose(2) << "Molecule mol is there and has " << mol->AtomCount << " Atoms \n";
+
+	// 1. searching topmost atom with respect to each axis
+	Walker = mol->start;
+	while (Walker->next != mol->end) {
+		Walker = Walker->next;
+		for (i = 0; i < 3; i++) {
+			if (Walker->x.x[i] > max_coordinate[i]) {
+				max_coordinate[i] = Walker->x.x[i];
+				max_index[i] = Walker;
+			}
+		}
+	}
+
+	cout << Verbose(2) << "Found maximum coordinates: ";
+	for (int i=0;i<NDIM;i++)
+		cout << i << ": " << *max_index[i] << "\t";
+	cout << endl;
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k = 1;
-  Oben.x[k] = 1.;
-  FirstPoint = max_index[k];
-
-  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
-  double Storage[3];
-  atom* Opt_Candidate = NULL;
-  Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
-  Storage[1] = 999999.; // This will be an angle looking for the third point.
-  Storage[2] = 999999.;
-
-  Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0,
-      Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
-  SecondPoint = Opt_Candidate;
-  cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
-
-  helper.CopyVector(&(FirstPoint->x));
-  helper.SubtractVector(&(SecondPoint->x));
-  helper.Normalize();
-  Oben.ProjectOntoPlane(&helper);
-  Oben.Normalize();
-  helper.VectorProduct(&Oben);
-  Storage[0] = -2.; // This will indicate the quadrant.
-  Storage[1] = 9999999.; // This will be an angle looking for the third point.
-  Storage[2] = 9999999.;
-
-  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
-  Chord.SubtractVector(&(SecondPoint->x));
-  // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
-
-  cout << Verbose(2) << "Looking for third point candidates \n";
-  // look in one direction of baseline for initial candidate
-  Opt_Candidate = NULL;
-  CenterOfFirstLine.CopyVector(&Chord);
-  CenterOfFirstLine.Scale(0.5);
-  CenterOfFirstLine.AddVector(&(SecondPoint->x));
-
-  cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0,
-      &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
-  // look in other direction of baseline for possible better candidate
-  cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
-  Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0,
-      &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
-  cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
-
-  // FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
-
-  // Finally, we only have to add the found points
-  AddTrianglePoint(FirstPoint, 0);
-  AddTrianglePoint(SecondPoint, 1);
-  AddTrianglePoint(Opt_Candidate, 2);
-  // ... and respective lines
-  AddTriangleLine(TPS[0], TPS[1], 0);
-  AddTriangleLine(TPS[1], TPS[2], 1);
-  AddTriangleLine(TPS[0], TPS[2], 2);
-  // ... and triangles to the Maps of the Tesselation class
-  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-  AddTriangleToLines();
-  // ... and calculate its normal vector (with correct orientation)
-  Oben.Scale(-1.);
-  BTS->GetNormalVector(Oben);
-  cout << Verbose(0) << "==> The found starting triangle consists of "
-      << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate
-      << " with normal vector " << BTS->NormalVector << ".\n";
-  cout << Verbose(1) << "End of Find_starting_triangle\n";
-}
-;
-
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
-{
-  int N = 0;
-  struct Tesselation *Tess = new Tesselation;
-  cout << Verbose(1) << "Entering search for non convex hull. " << endl;
-  cout << flush;
-  const double RADIUS = 6.;
-  LineMap::iterator baseline;
-  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
-  bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
-
-  Tess->Find_starting_triangle(mol, RADIUS);
-
-  baseline = Tess->LinesOnBoundary.begin();
-  while ((baseline != Tess->LinesOnBoundary.end()) || (flag))
-    {
-      if (baseline->second->TrianglesCount == 1)
-        {
-          Tess->Find_next_suitable_triangle(out, tecplot, mol,
-              *(baseline->second),
-              *(((baseline->second->triangles.begin()))->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
-          flag = true;
-        }
-      else
-        {
-          cout << Verbose(1) << "Line " << baseline->second << " has "
-              << baseline->second->TrianglesCount << " triangles adjacent"
-              << endl;
-        }
-      N++;
-      baseline++;
-      if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
-        baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-        flag = false;
-      }
-    }
-  cout << Verbose(1) << "Writing final tecplot file\n";
-  write_tecplot_file(out, tecplot, Tess, mol, -1);
-
-  cout << Verbose(0) << "End of Find_non_convex_border\n";
-}
-;
-
+	const int k = 1;
+	Oben.x[k] = 1.;
+	FirstPoint = max_index[k];
+
+	cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << " with " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << " bonds." << endl;
+	double Storage[3];
+	atom* Opt_Candidate = NULL;
+	Storage[0] = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
+	Storage[1] = 999999.; // This will be an angle looking for the third point.
+	Storage[2] = 999999.;
+
+	Find_second_point_for_Tesselation(FirstPoint, FirstPoint, FirstPoint, 0, Oben, Opt_Candidate, Storage, mol, RADIUS); // we give same point as next candidate as its bonds are looked into in find_second_...
+	SecondPoint = Opt_Candidate;
+	cout << Verbose(1) << "Found second point is " << *SecondPoint << " at " << SecondPoint->x << ".\n";
+
+	helper.CopyVector(&(FirstPoint->x));
+	helper.SubtractVector(&(SecondPoint->x));
+	helper.Normalize();
+	Oben.ProjectOntoPlane(&helper);
+	Oben.Normalize();
+	helper.VectorProduct(&Oben);
+	Storage[0] = -2.; // This will indicate the quadrant.
+	Storage[1] = 9999999.; // This will be an angle looking for the third point.
+	Storage[2] = 9999999.;
+
+	Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
+	Chord.SubtractVector(&(SecondPoint->x));
+	// Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
+
+	cout << Verbose(2) << "Looking for third point candidates \n";
+	// look in one direction of baseline for initial candidate
+	Opt_Candidate = NULL;
+	CenterOfFirstLine.CopyVector(&Chord);
+	CenterOfFirstLine.Scale(0.5);
+	CenterOfFirstLine.AddVector(&(SecondPoint->x));
+
+	cout << Verbose(1) << "Looking for third point candidates from " << *FirstPoint << " onward ...\n";
+	Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, SecondPoint, FirstPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine,  Opt_Candidate, Storage, RADIUS, mol);
+	// look in other direction of baseline for possible better candidate
+	cout << Verbose(1) << "Looking for third point candidates from " << *SecondPoint << " onward ...\n";
+	Find_next_suitable_point_via_Angle_of_Sphere(FirstPoint, SecondPoint, SecondPoint, FirstPoint, SecondPoint, 0, &Chord, &helper, &Oben, CenterOfFirstLine, Opt_Candidate, Storage, RADIUS, mol);
+	cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
+
+	// FOUND Starting Triangle: FirstPoint, SecondPoint, Opt_Candidate
+
+	// Finally, we only have to add the found points
+	AddTrianglePoint(FirstPoint, 0);
+	AddTrianglePoint(SecondPoint, 1);
+	AddTrianglePoint(Opt_Candidate, 2);
+	// ... and respective lines
+	AddTriangleLine(TPS[0], TPS[1], 0);
+	AddTriangleLine(TPS[1], TPS[2], 1);
+	AddTriangleLine(TPS[0], TPS[2], 2);
+	// ... and triangles to the Maps of the Tesselation class
+	BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+	AddTriangleToLines();
+	// ... and calculate its normal vector (with correct orientation)
+	Oben.Scale(-1.);
+	BTS->GetNormalVector(Oben);
+	cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+	cout << Verbose(1) << "End of Find_starting_triangle\n";
+};
+
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol, const char *filename, const double RADIUS)
+{
+	int N = 0;
+	struct Tesselation *Tess = new Tesselation;
+	cout << Verbose(1) << "Entering search for non convex hull. " << endl;
+	cout << flush;
+	LineMap::iterator baseline;
+	cout << Verbose(0) << "Begin of Find_non_convex_border\n";
+	bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
+	bool failflag = false;
+	if ((mol->first->next == mol->last) || (mol->last->previous == mol->first))
+		mol->CreateAdjacencyList((ofstream *)&cout, 1.6, true);
+
+	Tess->Find_starting_triangle(mol, RADIUS);
+
+	baseline = Tess->LinesOnBoundary.begin();
+	while ((baseline != Tess->LinesOnBoundary.end()) || (flag)) {
+		if (baseline->second->TrianglesCount == 1) {
+			failflag = Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename); //the line is there, so there is a triangle, but only one.
+			flag = flag || failflag;
+			if (!failflag)
+				cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
+		} else {
+			cout << Verbose(1) << "Line " << *baseline->second << " has " << baseline->second->TrianglesCount << " triangles adjacent" << endl;
+		}
+		N++;
+		baseline++;
+		if ((baseline == Tess->LinesOnBoundary.end()) && (flag)) {
+			baseline = Tess->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
+			flag = false;
+		}
+	}
+	if (failflag) {
+		cout << Verbose(1) << "Writing final tecplot file\n";
+		if (DoTecplotOutput)
+			write_tecplot_file(out, tecplot, Tess, mol, -1);
+		if (DoRaster3DOutput)
+			write_raster3d_file(out, tecplot, Tess, mol);
+	} else {
+		cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
+	}
+
+	cout << Verbose(0) << "End of Find_non_convex_border\n";
+	delete(Tess);
+};
+
