Changeset cc2ee5 for src

Timestamp:

Dec 29, 2008, 9:46:58 PM (17 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, Candidate_v1.7.0, 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:

a98603

Parents:

02bfde

git-author:

Frederik Heber <heber@…> (12/29/08 21:26:43)

git-committer:

Frederik Heber <heber@…> (12/29/08 21:46:58)

Message:

Just a temporary commit

Location:

src

Files:

• .#border.cpp (deleted)
• .#molecules.hpp (deleted)
• Hbondangle.db (modified) (1 prop)
• Hbonddistance.db (modified) (1 prop)
• Makefile.am (modified) (1 prop)
• analyzer.cpp (modified) (1 prop)
• atom.cpp (modified) (1 prop)
• bond.cpp (modified) (1 prop)
• boundary.cpp (modified) (12 diffs, 1 prop)
• boundary.hpp (modified) (2 diffs, 1 prop)
• builder.cpp (modified) (3 diffs, 1 prop)
• config.cpp (modified) (1 prop)
• datacreator.cpp (modified) (1 prop)
• datacreator.hpp (modified) (1 prop)
• defs.hpp (modified) (1 prop)
• element.cpp (modified) (1 prop)
• elements.db (modified) (1 prop)
• graph.cpp (modified) (1 prop)
• helpers.cpp (modified) (1 prop)
• helpers.hpp (modified) (1 prop)
• joiner.cpp (modified) (1 prop)
• moleculelist.cpp (modified) (1 prop)
• molecules.cpp (modified) (1 prop)
• molecules.hpp (modified) (1 diff, 1 prop)
• orbitals.db (modified) (1 prop)
• parser.cpp (modified) (1 prop)
• parser.hpp (modified) (1 prop)
• periodentafel.cpp (modified) (19 diffs, 1 prop)
• periodentafel.hpp (modified) (1 prop)
• stackclass.hpp (modified) (1 prop)
• valence.db (modified) (1 prop)
• vector.cpp (modified) (1 prop)
• vector.hpp (modified) (1 diff, 1 prop)
• verbose.cpp (modified) (1 prop)

src/boundary.cpp

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

src/boundary.hpp

@@ -86 +86 @@
     void AddTriangleToLines();
     void Find_starting_triangle(molecule* mol, const double RADIUS);
-    void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N);
+    bool Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N, const char *filename);
+    bool CheckPresenceOfTriangle(ofstream *out, atom *a, atom *b, atom *c);
+    void 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);
 
     PointMap PointsOnBoundary;
@@ -110 +112 @@
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem);
 void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
+void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol, const char *tempbasename, const double RADIUS);
 void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule *mol, bool problem);
-void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol);
 
 

src/builder.cpp

@@ -259 +259 @@
 static void CenterAtoms(molecule *mol)
 {
-  Vector x, y;
+  Vector x, y, helper;
   char choice;  // menu choice char
 
@@ -292 +292 @@
       mol->CenterEdge((ofstream *)&cout, &x);  // make every coordinate positive
       mol->Translate(&y); // translate by boundary
-      mol->SetBoxDimension(&(x+y*2));  // update Box of atoms by boundary
+      helper.CopyVector(&y);
+      helper.Scale(2.);
+      helper.AddVector(&x);
+      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
       break;
     case 'd':
@@ -1004 +1007 @@
               }
               break;
-            case 'A':
-                ExitFlag = 1;
-                if ((argptr >= argc) || (argv[argptr][0] == '-')){
-                        ExitFlag =255;
-                        cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
-                }
-                else{
-                        cout << "Parsing bonds from " << argv[argptr] << "." << endl;
-                    ifstream *input = new ifstream(argv[argptr]);
-                        mol->CreateAdjacencyList2((ofstream *)&cout, input);
-                        input->close();
-                }
-                break;
-            case 'N':
-                ExitFlag = 1;
-                if ((argptr >= argc) || (argv[argptr][0] == '-')){
-                        ExitFlag = 255;
-                        cerr << "Not enough or invalid arguments given for non-convex envelope: -o <tecplot output file>" << endl;
-                        }
-                else {
-                        cout << Verbose(0) << "Evaluating npn-convex envelope.";
-                        ofstream *output = new ofstream(argv[argptr], ios::trunc);
-                        cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
-                        Find_non_convex_border((ofstream *)&cout, output, mol);
-                        output->close();
-                        delete(output);
-                        argptr+=1;
-                        }
-                break;
+                                                case 'A':
+                                                        ExitFlag = 1;
+                                                        if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-')) {
+                                                                ExitFlag =255;
+                                                                cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+                                                        } else {
+                                                                cout << "Parsing bonds from " << argv[argptr] << "." << endl;
+                                                                ifstream *input = new ifstream(argv[argptr]);
+                                                                mol->CreateAdjacencyList2((ofstream *)&cout, input);
+                                                                input->close();
+                                                        }
+                                                        break;
+                                                case 'N':
+                                                        ExitFlag = 1;
+                                                        if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
+                                                                ExitFlag = 255;
+                                                                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
+                                                        } else {
+                                                                cout << Verbose(0) << "Evaluating npn-convex envelope.";
+                                                                string TempName(argv[argptr+1]);
+                                                                TempName.append(".r3d");
+                                                                ofstream *output = new ofstream(TempName.c_str(), ios::trunc);
+                                                                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
+                                                                Find_non_convex_border((ofstream *)&cout, output, mol, argv[argptr+1], atof(argv[argptr]));
+                                                                output->close();
+                                                                delete(output);
+                                                                argptr+=2;
+                                                        }
+                                                        break;
             case 'T':
               ExitFlag = 1;

src/molecules.hpp

@@ -60 +60 @@
 #define LineMap map < int, class BoundaryLineSet * >
 #define LinePair pair < int, class BoundaryLineSet * >
-#define LineTestPair pair < LinePair::iterator, bool >
+#define LineTestPair pair < LineMap::iterator, bool >
 
 #define TriangleMap map < int, class BoundaryTriangleSet * >

src/periodentafel.cpp

@@ -1 +1 @@
 /** \file periodentafel.cpp
- * 
+ *
  * Function implementations for the class periodentafel.
- * 
+ *
  */
 
@@ -14 +14 @@
  * Initialises start and end of list and resets periodentafel::checkliste to false.
  */
-periodentafel::periodentafel() 
-{ 
-  start = new element; 
-  end = new element; 
-  start->previous = NULL; 
-  start->next = end; 
-  end->previous = start; 
+periodentafel::periodentafel()
+{
+  start = new element;
+  end = new element;
+  start->previous = NULL;
+  start->next = end;
+  end->previous = start;
   end->next = NULL;
 };
@@ -27 +27 @@
  * Removes every element and afterwards deletes start and end of list.
  */
-periodentafel::~periodentafel() 
-{ 
-  CleanupPeriodtable(); 
-  delete(end); 
-  delete(start); 
-}; 
+periodentafel::~periodentafel()
+{
+  CleanupPeriodtable();
+  delete(end);
+  delete(start);
+};
 
 /** Adds element to period table list
@@ -38 +38 @@
  * \return true - succeeded, false - does not occur
  */
-bool periodentafel::AddElement(element *pointer) 
-{ 
+bool periodentafel::AddElement(element *pointer)
+{
   pointer->sort = &pointer->Z;
   if (pointer->Z < 1 && pointer->Z >= MAX_ELEMENTS)
     cout << Verbose(0) << "Invalid Z number!\n";
-  return add(pointer, end); 
+  return add(pointer, end);
 };
 
@@ -50 +50 @@
  * \return true - succeeded, false - element not found
  */
-bool periodentafel::RemoveElement(element *pointer) 
-{ 
-  return remove(pointer, start, end); 
+bool periodentafel::RemoveElement(element *pointer)
+{
+  return remove(pointer, start, end);
 };
 
@@ -58 +58 @@
  * \return true - succeeded, false - does not occur
  */
-bool periodentafel::CleanupPeriodtable() 
-{ 
-  return cleanup(start,end); 
+bool periodentafel::CleanupPeriodtable()
+{
+  return cleanup(start,end);
 };
 
@@ -76 +76 @@
     cout << Verbose(0) << "Mass: " << endl;
     cin >> walker->mass;
-    walker->Z = Z;  
-    cout << Verbose(0) << "Atomic number: " << walker->Z << endl;  
+    walker->Z = Z;
+    cout << Verbose(0) << "Atomic number: " << walker->Z << endl;
     cout << Verbose(0) << "Name [max 64 chars]: " << endl;
     cin >> walker->name;
@@ -105 +105 @@
 /** Asks for element number and returns pointer to element
  */
-element * periodentafel::AskElement() 
+element * periodentafel::AskElement()
 {
   element *walker = NULL;
@@ -117 +117 @@
 };
 
-
 /** Prints period table to given stream.
  * \param output stream
- */ 
+ */
 bool periodentafel::Output(ofstream *output) const
 {
@@ -131 +130 @@
     }
     return result;
-  } else 
+  } else
     return false;
 };
@@ -138 +137 @@
  * \param *output output stream
  * \param *checkliste elements table for this molecule
- */ 
+ */
 bool periodentafel::Checkout(ofstream *output, const int *checkliste) const
 {
@@ -152 +151 @@
       if ((walker != NULL) && (walker->Z > 0) && (walker->Z < MAX_ELEMENTS) && (checkliste[walker->Z])) {
         walker->No = No;
-        result = result && walker->Checkout(output, No++, checkliste[walker->Z]);      
+        result = result && walker->Checkout(output, No++, checkliste[walker->Z]);
       }
     }
     return result;
-  } else 
+  } else
     return false;
 };
-
 
 /** Loads element list from file.
@@ -171 +169 @@
   bool status = true;
   bool otherstatus = true;
-  char *filename = new char[MAXSTRINGSIZE];
-  
+  char filename[255];
+
   // fill elements DB
   strncpy(filename, path, MAXSTRINGSIZE);
@@ -225 +223 @@
   if (infile != NULL) {
     while (!infile.eof()) {
-        infile >> tmp;
-        infile >> ws;
-        infile >> FindElement((int)tmp)->Valence;
-        infile >> ws;
-        //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
+      infile >> tmp;
+      infile >> ws;
+      infile >> FindElement((int)tmp)->Valence;
+      infile >> ws;
+      //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->Valence << " valence electrons." << endl;
     }
     infile.close();
@@ -253 +251 @@
   } else
     otherstatus = false;
-  
+
   // fill H-BondDistance DB per element
   strncpy(filename, path, MAXSTRINGSIZE);
@@ -261 +259 @@
   if (infile != NULL) {
     while (!infile.eof()) {
-        infile >> tmp;
+      infile >> tmp;
       ptr = FindElement((int)tmp);
-        infile >> ws;
+      infile >> ws;
       infile >> ptr->HBondDistance[0];
       infile >> ptr->HBondDistance[1];
       infile >> ptr->HBondDistance[2];
-        infile >> ws;
+      infile >> ws;
       //cout << Verbose(3) << "Element " << (int)tmp << " has " << FindElement((int)tmp)->HBondDistance[0] << " Angstrom typical distance to hydrogen." << endl;
     }
@@ -274 +272 @@
   } else
     otherstatus = false;
-  
+
   // fill H-BondAngle DB per element
   strncpy(filename, path, MAXSTRINGSIZE);
@@ -294 +292 @@
   } else
     otherstatus = false;
-  
+
   if (!otherstatus)
     cerr << "WARNING: Something went wrong while parsing the other databases!" << endl;
-  
+
   return status;
 };
@@ -308 +306 @@
   ofstream f;
   char filename[MAXSTRINGSIZE];
-  
+
   strncpy(filename, path, MAXSTRINGSIZE);
   strncat(filename, "/", MAXSTRINGSIZE-strlen(filename));

src/vector.hpp

@@ -56 +56 @@
 
 ostream & operator << (ostream& ost, Vector &m);
-Vector& operator+=(Vector& a, const Vector& b);
-Vector& operator*=(Vector& a, const double m);
-Vector& operator*(const Vector& a, const double m);
-Vector& operator+(const Vector& a, const Vector& b);
+//Vector& operator+=(Vector& a, const Vector& b);
+//Vector& operator*=(Vector& a, const double m);
+//Vector& operator*(const Vector& a, const double m);
+//Vector& operator+(const Vector& a, const Vector& b);
 
 #endif /*VECTOR_HPP_*/