Context Navigation

-              r986c80
+              rce5ac3
 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;
                         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;
+        }
+  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;
+  }
+}
+;
 …
 BoundaryTriangleSet::~BoundaryTriangleSet()
+{
         for (int i = 0; i < 3; i++) {
                 cout << Verbose(5) << "Erasing triangle Nr." << Nr << endl;
                 lines[i]->triangles.erase(Nr);
                 if (lines[i]->triangles.empty()) {
                         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;
+        }
+  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;
+  }
+}
+;
 …
 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);
+  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);
 };
 …
 Tesselation::~Tesselation()
+{
         cout << Verbose(1) << "Free'ing TesselStruct ... " << endl;
         for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
                 if (runner->second != NULL) {
                         delete (runner->second);
                         runner->second = NULL;
                 } else
                         cerr << "ERROR: The triangle " << runner->first << " has already been free'd." << endl;
+        }
         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;
+        }
+  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;
+  }
+}
+;
 …
     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 alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
         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";
+  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 alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
+  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";
+}
+;
 …
     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, helper;
         double CurrentEpsilon = 0.1;
         double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
         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;
+                                        }
+                                }
+                        }
+  /* 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;
+  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.
 …
+                }
                 */
                 } 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
+      }
+    }
+  }
 };
 …
     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;
+  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;
 };
 …
  */
 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;
+  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;
 };
 …
     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";
+  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);
+  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);
 };

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset ce5ac3 for src/boundary.cpp

Legend:

src/boundary.cpp

Download in other formats: