Changeset 542ab3

Timestamp:

Jul 9, 2009, 10:22:05 AM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

4777e9

Parents:

86234b

Message:

Removed unnecessary code

• Choose_preferable_third_point()
• already commented out Find_next_suitable_point()
• Tesselation::Find_next_suitable_point_via_Angle_of_Sphere()

File:

• src/boundary.cpp (modified) (3 diffs)

src/boundary.cpp

@@ -1710 +1710 @@
 
 
-
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
- *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
- *  upon which we operate.
- *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
- *  direction and angle into Storage.
- *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
- *  with all neighbours of the candidate.
- * @param a first point
- * @param b second point
- * *param c atom old third point of old triangle
- * @param Candidate base point along whose bonds to start looking from
- * @param Parent point to avoid during search as its wrong direction
- * @param RecursionLevel contains current recursion depth
- * @param Chord baseline vector of first and second point
- * @param direction1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
- * @param OldNormal normal of the triangle which the baseline belongs to
- * @param ReferencePoint Vector of center of circumscribing circle from which we look towards center of sphere
- * @param Opt_Candidate candidate reference to return
- * @param Storage array containing two angles of current Opt_Candidate
- * @param RADIUS radius of ball
- * @param mol molecule structure with atoms and bonds
- */
-void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
-    atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol)
-{
-  cout << Verbose(2) << "Begin of Find_next_suitable_point_via_Angle_of_Sphere, recursion level " << RecursionLevel << ".\n";
-  cout << Verbose(3) << "Candidate is "<< *Candidate << endl;
-  cout << Verbose(4) << "Baseline vector is " << *Chord << "." << endl;
-  cout << Verbose(4) << "ReferencePoint is " << ReferencePoint << "." << endl;
-  cout << Verbose(4) << "Normal of base triangle is " << *OldNormal << "." << endl;
-  cout << Verbose(4) << "Search direction is " << *direction1 << "." << endl;
-  /* OldNormal is normal vector on the old triangle
-   * direction1 is normal on the triangle line, from which we come, as well as on OldNormal.
-   * Chord points from b to a!!!
-   */
-  Vector dif_a; //Vector from a to candidate
-  Vector dif_b; //Vector from b to candidate
-  Vector AngleCheck;
-  Vector TempNormal, Umkreismittelpunkt;
-  Vector Mittelpunkt;
-
-  double CurrentEpsilon = 0.1;
-  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
-  double BallAngle, AlternativeSign;
-  atom *Walker; // variable atom point
-
-  Vector NewUmkreismittelpunkt;
-
-  if (a != Candidate and b != Candidate and c != Candidate) {
-    cout << Verbose(3) << "We have a unique candidate!" << endl;
-    dif_a.CopyVector(&(a->x));
-    dif_a.SubtractVector(&(Candidate->x));
-    dif_b.CopyVector(&(b->x));
-    dif_b.SubtractVector(&(Candidate->x));
-    AngleCheck.CopyVector(&(Candidate->x));
-    AngleCheck.SubtractVector(&(a->x));
-    AngleCheck.ProjectOntoPlane(Chord);
-
-    SideA = dif_b.Norm();
-    SideB = dif_a.Norm();
-    SideC = Chord->Norm();
-    //Chord->Scale(-1);
-
-    alpha = Chord->Angle(&dif_a);
-    beta = M_PI - Chord->Angle(&dif_b);
-    gamma = dif_a.Angle(&dif_b);
-
-    cout << Verbose(2) << "Base triangle has sides " << dif_a << ", " << dif_b << ", " << *Chord << " with angles " << alpha/M_PI*180. << ", " << beta/M_PI*180. << ", " << gamma/M_PI*180. << "." << endl;
-
-    if (fabs(M_PI - alpha - beta - gamma) > MYEPSILON) {
-      cerr << Verbose(0) << "WARNING: sum of angles for 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*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
-      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);
-        }
-        if (sign >= 0) {
-          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);
-
-          Mittelpunkt.SubtractVector(&ReferencePoint);
-          cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
-
-          BallAngle = Mittelpunkt.Angle(OldNormal);
-          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
-
-          //cout << "direction1 is " << *direction1 << "." << endl;
-          //cout << "Mittelpunkt is " << Mittelpunkt << "."<< endl;
-
-          //cout << Verbose(3) << "BallAngle is " << BallAngle << " Sign is " << sign << endl;
-
-          NewUmkreismittelpunkt.SubtractVector(&ReferencePoint);
-
-          if ((Mittelpunkt.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 << " is not in search direction." << 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";
-}
-;
-
-
 /** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
  * \param *Center new center on return
@@ -1978 +1791 @@
   }
 };
-
-
-/** This recursive function finds a third point, to form a triangle with two given ones.
- * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
- * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
- * the center of the sphere is still fixed up to a single parameter. The band of possible values
- * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
- * us the "null" on this circle, the new center of the candidate point will be some way along this
- * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
- * by the normal vector of the base triangle that always points outwards by construction.
- * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
- * We construct the normal vector that defines the plane this circle lies in, it is just in the
- * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
- * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
- * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
- * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
- * both.
- * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
- * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
- * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
- * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
- * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
- * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
- * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points
- * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line
- * @param OptCandidate candidate reference on return
- * @param OptCandidateCenter candidate's sphere center on return
- * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
- * @param RADIUS radius of sphere
- * @param *LC LinkedCell structure with neighbouring atoms
- */
-// void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
-// {
-//   atom *Walker = NULL;
-//   Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
-//   Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
-//   Vector OldSphereCenter;   // center of the sphere defined by the three points of BaseTriangle
-//   Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
-//   Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
-//   Vector NewNormalVector;   // normal vector of the Candidate's triangle
-//   Vector SearchDirection;   // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
-//   Vector helper;
-//   LinkedAtoms *List = NULL;
-//   double CircleRadius; // radius of this circle
-//   double radius;
-//   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
-//   double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
-//   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
-//   atom *Candidate = NULL;
-//
-//   cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
-//
-//   cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
-//
-//   // construct center of circle
-//   CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
-//   CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
-//   CircleCenter.Scale(0.5);
-//
-//   // construct normal vector of circle
-//   CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
-//   CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
-//
-//   // calculate squared radius of circle
-//   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-//   if (radius/4. < RADIUS*RADIUS) {
-//     CircleRadius = RADIUS*RADIUS - radius/4.;
-//     CirclePlaneNormal.Normalize();
-//     cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
-//
-//     // construct old center
-//     GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
-//     helper.CopyVector(&BaseTriangle->NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-//     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
-//     helper.Scale(sqrt(RADIUS*RADIUS - radius));
-//     OldSphereCenter.AddVector(&helper);
-//     OldSphereCenter.SubtractVector(&CircleCenter);
-//     cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-//
-//     // test whether old center is on the band's plane
-//     if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-//       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
-//       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
-//     }
-//     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
-//     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-//
-//       // construct SearchDirection
-//       SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
-//       helper.CopyVector(&BaseLine->endpoints[0]->node->x);
-//       for(int i=0;i<3;i++)  // just take next different endpoint
-//         if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
-//           helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
-//         }
-//       if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // ohoh, SearchDirection points inwards!
-//         SearchDirection.Scale(-1.);
-//       SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-//       SearchDirection.Normalize();
-//       cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-//       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
-//         cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
-//       }
-//
-//       if (LC->SetIndexToVector(&CircleCenter)) {  // get cell for the starting atom
-//         for(int i=0;i<NDIM;i++) // store indices of this cell
-//           N[i] = LC->n[i];
-//         cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
-//       } else {
-//         cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
-//         return;
-//       }
-//       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
-//       cout << Verbose(2) << "LC Intervals:";
-//       for (int i=0;i<NDIM;i++) {
-//         Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
-//         Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
-//         cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
-//       }
-//       cout << endl;
-//       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
-//         for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
-//           for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
-//             List = LC->GetCurrentCell();
-//             cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
-//             if (List != NULL) {
-//               for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-//                 Candidate = (*Runner);
-//
-//                 // check for three unique points
-//                 if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
-//                   cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
-//
-//                   // construct both new centers
-//                   GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
-//                   OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-//
-//                   if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
-//                     helper.CopyVector(&NewNormalVector);
-//                     cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-//                     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
-//                     if (radius < RADIUS*RADIUS) {
-//                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
-//                       cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
-//                       NewSphereCenter.AddVector(&helper);
-//                       NewSphereCenter.SubtractVector(&CircleCenter);
-//                       cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-//
-//                       helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
-//                       OtherNewSphereCenter.AddVector(&helper);
-//                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
-//                       cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
-//
-//                       // check both possible centers
-//                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
-//                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
-//                       alpha = min(alpha, Otheralpha);
-//                       if (*ShortestAngle > alpha) {
-//                           OptCandidate = Candidate;
-//                           *ShortestAngle = alpha;
-//                           if (alpha != Otheralpha)
-//                             OptCandidateCenter->CopyVector(&NewSphereCenter);
-//                           else
-//                             OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
-//                           cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
-//                       } else {
-//                         if (OptCandidate != NULL)
-//                           cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
-//                         else
-//                           cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
-//                       }
-//
-//                     } else {
-//                       cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
-//                     }
-//                   } else {
-//                     cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
-//                   }
-//                 } else {
-//                   cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
-//                 }
-//               }
-//             }
-//           }
-//     } else {
-//       cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-//     }
-//   } else {
-//     cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
-//   }
-//
-//   cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
-// };
 
 
@@ -2450 +2071 @@
 };
 
-/**
- * Finds the preferable out of two third-point candidates with equal angles.
- *
- * @param Candidate - this and the second parameter are evaluated
- * @param OptCandidate - this and the second parameter are evaluated
- * @param current base line
- * @param third node of the base triangle
- * @param tesselation object
- *
- * @return true if Candidate should be taken, false if OptCandidate should be kept
- */
-bool Choose_preferable_third_point(
-        atom *Candidate, atom *OptCandidate, class BoundaryLineSet *BaseLine,
-        atom *ThirdNode, Tesselation *Tess
-) {
-        bool takeNewCandidate;
-
-        ofstream *out = new ofstream();
-        atom *Atoms[3];
-        bool optCandidateAndBaseLineFormTriangle = (ThirdNode != NULL) && (OptCandidate == ThirdNode);
-        bool candidateAndBaseLineFormTriangle = (ThirdNode != NULL) && (Candidate == ThirdNode);
-        Atoms[0] = Candidate;
-        Atoms[1] = OptCandidate;
-        Atoms[2] = BaseLine->endpoints[0]->node;
-        bool candidatesAndBaseLineNode0FormTriangle = (Tess->CheckPresenceOfTriangle(out, Atoms) > 0);
-        Atoms[0] = Candidate;
-        Atoms[1] = OptCandidate;
-        Atoms[2] = BaseLine->endpoints[1]->node;
-        bool candidatesAndBaseLineNode1FormTriangle = (Tess->CheckPresenceOfTriangle(out, Atoms) > 0);
-        Vector halfBaseLine;
-        halfBaseLine.CopyVector(&BaseLine->endpoints[0]->node->x);
-        halfBaseLine.AddVector(&BaseLine->endpoints[1]->node->x);
-        halfBaseLine.Scale(0.5);
-
-        if (optCandidateAndBaseLineFormTriangle) {
-                takeNewCandidate = (!existsIntersection(Candidate->x, halfBaseLine, OptCandidate->x, BaseLine->endpoints[0]->node->x)
-                        && !existsIntersection(Candidate->x, halfBaseLine, OptCandidate->x, BaseLine->endpoints[1]->node->x));
-        } else if (candidateAndBaseLineFormTriangle) {
-                takeNewCandidate = (existsIntersection(OptCandidate->x, halfBaseLine, Candidate->x, BaseLine->endpoints[0]->node->x)
-                        || existsIntersection(OptCandidate->x, halfBaseLine, Candidate->x, BaseLine->endpoints[1]->node->x));
-        } else if (candidatesAndBaseLineNode0FormTriangle) {
-                takeNewCandidate = !existsIntersection(OptCandidate->x, BaseLine->endpoints[0]->node->x, Candidate->x, halfBaseLine);
-        } else if (candidatesAndBaseLineNode1FormTriangle) {
-                takeNewCandidate = !existsIntersection(OptCandidate->x, BaseLine->endpoints[1]->node->x, Candidate->x, halfBaseLine);
-        } else {
-                takeNewCandidate = (ThirdNode == NULL)
-                        || ((!existsIntersection(Candidate->x, halfBaseLine, ThirdNode->x, BaseLine->endpoints[0]->node->x)
-                                && !existsIntersection(Candidate->x, halfBaseLine, ThirdNode->x, BaseLine->endpoints[1]->node->x)));
-        }
-
-        return takeNewCandidate;
-};
-
 
 struct Intersection {