Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision 10af0deb228f5e5b25d4fd0c600059654562e5c7)
+++ src/boundary.cpp	(revision e9fa06f36adaf7b5bf0b33b34e7ea747201e8bff)
@@ -1437,6 +1437,5 @@
  * @param mol molecule structure with atoms and bonds
  */
-void
-Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *d1, Vector *OldNormal,
     atom*& Opt_Candidate, Vector *Opt_Mittelpunkt, double *Storage, const double RADIUS, molecule* mol)
@@ -1448,9 +1447,10 @@
   Vector dif_a; //Vector from a to candidate
   Vector dif_b; //Vector from b to candidate
-  Vector AngleCheck; // Projection of a difference vector on plane orthogonal on triangle side.
+  Vector AngleCheck, AngleCheckReference, DirectionCheckPoint;
   Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
 
   double CurrentEpsilon = 0.1;
   double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+  double BallAngle;
   atom *Walker; // variable atom point
 
@@ -1460,7 +1460,7 @@
   dif_b.CopyVector(&(b->x));
   dif_b.SubtractVector(&(Candidate->x));
-  AngleCheck.CopyVector(&dif_a);
-  AngleCheck.Scale(-1);
-  AngleCheck.ProjectOntoPlane(Chord);
+  DirectionCheckPoint.CopyVector(&dif_a);
+  DirectionCheckPoint.Scale(-1);
+  DirectionCheckPoint.ProjectOntoPlane(Chord);
 
   SideA = dif_b.Norm();
@@ -1484,13 +1484,13 @@
   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
-// Those are remains from Freddie. Needed?
-
-//      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-//      Old Chordtest. Replaced by test for Radius of Circumscribing circle
+      //      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);
@@ -1499,15 +1499,11 @@
       //cout << SideC / 2.0 / sin(gamma) << endl;
 
-      if (Umkreisradius < RADIUS)
-        {
-
-          sign = AngleCheck.ScalarProduct(d1);
-          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
-
-          // intermediate calculations:
+      if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) //Checking whether ball will at least rest o points.
+        {
+
+          // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
 
           Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
           Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
 
           TempNormal.CopyVector(&dif_a);
@@ -1519,5 +1515,4 @@
           TempNormal.Normalize();
           Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-
           TempNormal.Scale(Restradius);
 
@@ -1525,4 +1520,16 @@
           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);
+          sign /= fabs(sign); // +1 if in direction of triangle plane, -1 if in other direction...
 
 
@@ -1533,5 +1540,5 @@
               Opt_Candidate = Candidate;
               Storage[0] = sign;
-              Storage[1] = AngleCheck.Angle(OldNormal);
+              Storage[1] = BallAngle;
               Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
               cout << "Angle is " << Storage[1] << ", Halbraum ist "
@@ -1542,4 +1549,6 @@
           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.
@@ -1550,34 +1559,19 @@
 
               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.
+                  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] = AngleCheck.Angle(OldNormal);
+                      Storage[1] = BallAngle;
                       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;
-                                    */
+
 
                     }
@@ -1585,7 +1579,7 @@
                     {
                       if ((fabs(sign - Storage[0]) < CurrentEpsilon && sign > 0
-                          && Storage[1] > AngleCheck.Angle(OldNormal)) ||
+                          && Storage[1] > BallAngle) ||
                           (fabs(sign - Storage[0]) < CurrentEpsilon && sign < 0
-                          && Storage[1] < AngleCheck.Angle(OldNormal)))
+                          && Storage[1] < BallAngle))
                       //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
                         {
@@ -1593,40 +1587,44 @@
                           Opt_Candidate = Candidate;
                           Storage[0] = sign;
-                          Storage[1] = AngleCheck.Angle(OldNormal);
+                          Storage[1] = BallAngle;
                           Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
                           cout << "Angle is " << Storage[1] << ", Halbraum ist "
                               << Storage[0] << endl;
                         }
+
                     }
                 }
-              else
+              /*
+               * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
+               *
+                else
                 {
-                  if (sign>0 && AngleCheck.Angle(OldNormal)>0 && Storage[0]<0)
+                  if (sign>0 && BallAngle>0 && Storage[0]<0)
                     {
                       cout << "Next better candidate is " << *Candidate << " with ";
                       Opt_Candidate = Candidate;
                       Storage[0] = sign;
-                      Storage[1] = AngleCheck.Angle(OldNormal);
+                      Storage[1] = BallAngle;
                       Opt_Mittelpunkt->CopyVector(&Mittelpunkt);
                       cout << "Angle is " << Storage[1] << ", Halbraum ist "
                       << Storage[0] << endl;
 
-/*
+//Debugging purposes only
                       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;
-                                    */
+                      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;
+
 
 
@@ -1638,4 +1636,5 @@
                     }
                 }
+                */
             }
         }
@@ -1728,4 +1727,20 @@
   Chord.SubtractVector(&(Line.endpoints[1]->node->x));
 
+  {
+    Vector TempNormal, Umkreismittelpunkt, Mittelpunkt;
+    double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
+
+    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);
+
+      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)));
+  }
   cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
 
