Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision e4ea462e4975fcf1debaf5bbf7e4ced6f9d55b4d)
+++ src/boundary.cpp	(revision 10af0deb228f5e5b25d4fd0c600059654562e5c7)
@@ -1444,5 +1444,5 @@
   /* 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 a to b.
+   * Chord points from b to a!!!
    */
   Vector dif_a; //Vector from a to candidate
@@ -1495,7 +1495,7 @@
 
       Umkreisradius = SideA / 2.0 / sin(alpha);
-      cout << Umkreisradius << endl;
-      cout << SideB / 2.0 / sin(beta) << endl;
-      cout << SideC / 2.0 / sin(gamma) << endl;
+      //cout << Umkreisradius << endl;
+      //cout << SideB / 2.0 / sin(beta) << endl;
+      //cout << SideC / 2.0 / sin(gamma) << endl;
 
       if (Umkreisradius < RADIUS)
@@ -1510,8 +1510,4 @@
           Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
 
-          cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-          cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-          cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-          cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
 
           TempNormal.CopyVector(&dif_a);
@@ -1523,17 +1519,11 @@
           TempNormal.Normalize();
           Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-          cout << "Umkreisradius ist " << Umkreisradius << endl;
-          cout << "Restradius ist " << Restradius << endl;
+
           TempNormal.Scale(Restradius);
-          cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+
           Mittelpunkt.CopyVector(&Umkreismittelpunkt);
           Mittelpunkt.AddVector(&TempNormal);  //this is center of sphere supported by a, b and Candidate
-          cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-          cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-          cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-          cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-          cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+
+
 
           if (Storage[0]< -1.5) // first Candidate at all
@@ -1547,4 +1537,6 @@
               cout << "Angle is " << Storage[1] << ", Halbraum ist "
                   << Storage[0] << endl;
+
+
             }
           else
@@ -1554,13 +1546,13 @@
 
               Distance = Opt_Candidate->x.Distance(&Mittelpunkt);
-              cout << "Opt_Candidate " << Opt_Candidate << " has distance " << Distance << " to Center of Candidate " << endl;
+              //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
-                      //sqrt(Candidate->x.Distance(Opt_Mittelpunkt)) < RADIUS)    //and those where Candidate would be within old Sphere.
+                  //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 ";
@@ -1571,4 +1563,22 @@
                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
                           << Storage[0] << endl;
+/*
+                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
+                                    cout << "Restradius ist " << Restradius << endl;
+                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+                                    */
+
                     }
                   else
@@ -1592,6 +1602,39 @@
               else
                 {
-                  if (DEBUG)
-                    cout << "Looses to better candidate" << endl;
+                  if (sign>0 && AngleCheck.Angle(OldNormal)>0 && Storage[0]<0)
+                    {
+                      cout << "Next better candidate is " << *Candidate << " with ";
+                      Opt_Candidate = Candidate;
+                      Storage[0] = sign;
+                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
+                      cout << "Angle is " << Storage[1] << ", Halbraum ist "
+                      << Storage[0] << endl;
+
+/*
+                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
+                                    cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
+                                    cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
+                                    cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
+                                    cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
+                                    cout << "Umkreisradius ist " << Umkreisradius << endl;
+                                    cout << "Restradius ist " << Restradius << endl;
+                                    cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
+                                    cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
+                                    cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
+                                    cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
+                                    cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
+                                    */
+
+
+                    }
+                  else
+                    {
+                      if (DEBUG)
+                        cout << "Looses to better candidate" << endl;
+                    }
                 }
             }
@@ -1612,5 +1655,5 @@
     }
 
-  if (RecursionLevel < 7) // Five is the recursion level threshold.
+  if (RecursionLevel < 9) // Five is the recursion level threshold.
     {
       for (int i = 0; i < mol->NumberOfBondsPerAtom[Candidate->nr]; i++) // go through all bond
@@ -1685,6 +1728,6 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
-  cout << Verbose(1) << "Looking for third point candidates for triangle ... "
-      << endl;
+  cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+
   Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node,
       Line.endpoints[0]->node, Line.endpoints[1]->node, 0, &Chord, &direction1,
@@ -1693,4 +1736,5 @@
       Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1,
       &(T.NormalVector), Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+
 
   if ((TrianglesOnBoundaryCount % 10) == 0)
@@ -1704,17 +1748,10 @@
   }
 
-  if ((Opt_Candidate == NULL) || (N<12))
+  if (TrianglesOnBoundaryCount >1000 )
     {
       cout << Verbose(1)
           << "No new Atom found, triangle construction will crash" << endl;
       write_tecplot_file(out, tecplot, this, mol, TriangleFilesWritten);
-      Find_next_suitable_point(Line.endpoints[0]->node,
-          Line.endpoints[1]->node, Line.endpoints[0]->node,
-          Line.endpoints[1]->node, 0, &Chord, &direction1, &(T.NormalVector),
-          Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
-      Find_next_suitable_point(Line.endpoints[0]->node,
-          Line.endpoints[1]->node, Line.endpoints[1]->node,
-          Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector),
-          Opt_Candidate, &Opt_Mittelpunkt, Storage, RADIUS, mol);
+      cout << "This is currently added candidate" << Opt_Candidate << endl;
     }
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
@@ -1733,13 +1770,11 @@
   AddTriangleToLines();
   cout << "New triangle with " << *BTS << endl;
-
-  cout << Verbose(1) << "Constructing normal vector for this triangle ... "
-      << endl;
+  cout << "We have "<< TrianglesOnBoundaryCount << endl;
+  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
 
   BTS->GetNormalVector(BTS->NormalVector);
 
-  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)
-      || ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0]
-          < 0)))
+  if ((BTS->NormalVector.ScalarProduct(&(T.NormalVector)) < 0 && Storage[0] > 0)  ||
+      (BTS->NormalVector.ScalarProduct(&(T.NormalVector)) > 0 && Storage[0] < 0))
     {
       BTS->NormalVector.Scale(-1);
