Changeset d067d45 for src/boundary.cpp

Timestamp:

Jul 23, 2009, 1:45:24 PM (16 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:

51c910

Parents:

ce5ac3 (diff), 437922 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Frederik Heber <heber@…> (07/23/09 12:34:47)

git-committer:

Frederik Heber <heber@…> (07/23/09 13:45:24)

Message:

Merge branch 'MultipleMolecules'

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/atom.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/helpers.cpp
molecuilder/src/joiner.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp
molecuilder/src/vector.cpp
molecuilder/src/verbose.cpp

merges:

• analyzer.cpp: all from HEAD, Hessian stuff
• atom.cpp: ostream & operator << (ostream &ost, const atom &a), 2nd param is const
• boundary.cpp: DoSingleStepOutput, VRMLSUffix, HULLEPSILON, Get_center_of_sphere(), Find_next_suitable_point_via_Angle_of_Sphere() and all of tesselation stuff from MultipleMolecules; DoRaster3DOutput, ~BoundaryLineSet(), ~Tesselation() from HEAD; write_raster3d_file with some changes from MultipleMolecules
• boundary.hpp: all from MultipleMolecules
• builder.cpp: all in ParseCommandLineOptions() from MultipleMolecules
• config.cpp: all from HEAD, is stuff from ThermoStat and ConstrainedMD
• datacreator.hpp: all from HEAD, Hessian stuff
• helpers.cpp: all from HEAD, ReAlloc() differed just by indentation
• moleculelist.cpp: all from MultipleMolecules, has to do with introduced molecule centers
• molecules.cpp: all from HEAD, was ConstrainedMD stuff
• molecules.hpp: ConstrainedMD and VerletForceIntegration from HEAD, moleculelist stuff from MultipleMolecules
• parser.cpp: all from HEAD, ColumnCounter is *ColumnCounter there
• parser.hpp: all from HEAD
• vector.cpp: completely from MultipleMolecules (more lines)
• verbose.cpp: all from HEAD

compilation fixes:

• CheckPresenceOfTriangle() from MultipleMolecules
• MoleculeListClass::ListOfMolecules has changed to STL list, hence use insert(), ReturnIndex() and MoleculeListClass::OutputConfigForListOfFragments
• MoleculeListClass::OutputConfigForListOfFragments() had fragmentprefix instead of FRAGMENTPREFIX
• parser.hpp: EnergyMatrix::ParseIndices() declaration was missing

File:

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

src/boundary.cpp

@@ -1 @@
-#include "molecules.hpp"
 #include "boundary.hpp"
 
 #define DEBUG 1
-#define DoTecplotOutput 0
+#define DoSingleStepOutput 0
+#define DoTecplotOutput 1
 #define DoRaster3DOutput 1
+#define DoVRMLOutput 1
 #define TecplotSuffix ".dat"
 #define Raster3DSuffix ".r3d"
+#define VRMLSUffix ".wrl"
+#define HULLEPSILON MYEPSILON
 
 // ======================================== Points on Boundary =================================
@@ -28 +31 @@
 {
   cout << Verbose(5) << "Erasing point nr. " << Nr << "." << endl;
+  if (!lines.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some lines." << endl;
   node = NULL;
   lines.clear();
@@ -33 +38 @@
 ;
 
-void
-BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
   cout << Verbose(6) << "Adding " << *this << " to line " << *line << "."
@@ -101 +105 @@
       cout << Verbose(5) << *endpoints[i] << " has still lines it's attached to." << endl;
   }
+  if (!triangles.empty())
+    cerr << "WARNING: Memory Leak! I " << *this << " am still connected to some triangles." << endl;
 }
 ;
@@ -109 +115 @@
   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
       << endl;
-  triangles.insert(TrianglePair(TrianglesCount, triangle));
+  triangles.insert(TrianglePair(triangle->Nr, triangle));
   TrianglesCount++;
 }
@@ -205 +211 @@
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (endpoints[0]->node->x.Projection(&OtherVector) > 0)
+  if (NormalVector.Projection(&OtherVector) > 0)
     NormalVector.Scale(-1.);
 }
@@ -565 +571 @@
 ;
 
-/** Creates the objects in a raster3d file (renderable with a header.r3d)
+/** Creates the objects in a VRML file.
  * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
+ * \param *vrmlfile output stream for tecplot data
+ * \param *Tess Tesselation structure with constructed triangles
+ * \param *mol molecule structure with atom positions
+ */
+void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
+{
+  atom *Walker = mol->start;
+  bond *Binder = mol->first;
+  int i;
+  Vector *center = mol->DetermineCenterOfAll(out);
+  if (vrmlfile != NULL) {
+    //cout << Verbose(1) << "Writing Raster3D file ... ";
+    *vrmlfile << "#VRML V2.0 utf8" << endl;
+    *vrmlfile << "#Created by molecuilder" << endl;
+    *vrmlfile << "#All atoms as spheres" << endl;
+    while (Walker->next != mol->end) {
+      Walker = Walker->next;
+      *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Walker->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
+    }
+
+    *vrmlfile << "# All bonds as vertices" << endl;
+    while (Binder->next != mol->last) {
+      Binder = Binder->next;
+      *vrmlfile << "3" << endl << "  "; // 2 is round-ended cylinder type
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Binder->leftatom->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.03\t";
+      for (i=0;i<NDIM;i++)
+        *vrmlfile << Binder->rightatom->x.x[i]+center->x[i] << " ";
+      *vrmlfile << "\t0.03\t0. 0. 1." << endl; // radius 0.05 and blue as colour
+    }
+
+    *vrmlfile << "# All tesselation triangles" << endl;
+    for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
+      *vrmlfile << "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
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->x.x[j]+center->x[j] << " ";
+        *vrmlfile << "\t";
+      }
+      *vrmlfile << "1. 0. 0." << endl;  // red as colour
+      *vrmlfile << "18" << endl << "  0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
+    }
+  } else {
+    cerr << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl;
+  }
+  delete(center);
+};
+
+/** Creates the objects in a raster3d file (renderable with a header.r3d).
+ * \param *out output stream for debugging
+ * \param *rasterfile output stream for tecplot data
  * \param *Tess Tesselation structure with constructed triangles
  * \param *mol molecule structure with atom positions
@@ -603 +663 @@
 
     *rasterfile << "# All tesselation triangles" << endl;
+    *rasterfile << "8\n  25. -1.   1. 1. 1.   0.0    0 0 0 2\n  SOLID     1.0 0.0 0.0\n  BACKFACE  0.3 0.3 1.0   0 0\n";
     for (TriangleMap::iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
       *rasterfile << "1" << endl << "  ";  // 1 is triangle type
@@ -613 +674 @@
       *rasterfile << "18" << endl << "  0.5 0.5 0.5" << endl;  // 18 is transparency type for previous object
     }
+    *rasterfile << "9\n  terminating special property\n";
   } else {
     cerr << "ERROR: Given rasterfile is " << rasterfile << "." << endl;
@@ -673 +735 @@
  * Determines first the convex envelope, then tesselates it and calculates its volume.
  * \param *out output stream for debugging
- * \param *tecplot output stream for tecplot data
+ * \param *filename filename prefix for output of vertex data
  * \param *configuration needed for path to store convex envelope file
  * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
@@ -680 +742 @@
  */
 double
-VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration,
+VolumeOfConvexEnvelope(ofstream *out, const char *filename, config *configuration,
     Boundaries *BoundaryPtr, molecule *mol)
 {
@@ -761 +823 @@
       y.CopyVector(&runner->second->endpoints[0]->node->x);
       y.SubtractVector(&runner->second->endpoints[2]->node->x);
-      a = sqrt(runner->second->endpoints[0]->node->x.Distance(
+      a = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
           &runner->second->endpoints[1]->node->x));
-      b = sqrt(runner->second->endpoints[0]->node->x.Distance(
+      b = sqrt(runner->second->endpoints[0]->node->x.DistanceSquared(
           &runner->second->endpoints[2]->node->x));
-      c = sqrt(runner->second->endpoints[2]->node->x.Distance(
+      c = sqrt(runner->second->endpoints[2]->node->x.DistanceSquared(
           &runner->second->endpoints[1]->node->x));
-      G = sqrt(((a * a + b * b + c * c) * (a * a + b * b + c * c) - 2 * (a * a
-          * a * a + b * b * b * b + c * c * c * c)) / 16.); // area of tesselated triangle
+      G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
       x.MakeNormalVector(&runner->second->endpoints[0]->node->x,
           &runner->second->endpoints[1]->node->x,
@@ -797 +858 @@
 
   // 8. Store triangles in tecplot file
+  string OutputName(filename);
+  OutputName.append(TecplotSuffix);
+  ofstream *tecplot = new ofstream(OutputName.c_str());
   write_tecplot_file(out, tecplot, TesselStruct, mol, 0);
+  tecplot->close();
+  delete(tecplot);
 
   // free reference lists
@@ -998 +1064 @@
           for (; C != PointsOnBoundary.end(); C++)
             {
-              tmp = A->second->node->x.Distance(&B->second->node->x);
+              tmp = A->second->node->x.DistanceSquared(&B->second->node->x);
               distance = tmp * tmp;
-              tmp = A->second->node->x.Distance(&C->second->node->x);
+              tmp = A->second->node->x.DistanceSquared(&C->second->node->x);
               distance += tmp * tmp;
-              tmp = B->second->node->x.Distance(&C->second->node->x);
+              tmp = B->second->node->x.DistanceSquared(&C->second->node->x);
               distance += tmp * tmp;
               DistanceMMap.insert(DistanceMultiMapPair(distance, pair<
@@ -1070 +1136 @@
             }
           // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->x.Distance(&A->second->node->x);
+          tmp = checker->second->node->x.DistanceSquared(&A->second->node->x);
           int innerpoint = 0;
-          if ((tmp < A->second->node->x.Distance(
+          if ((tmp < A->second->node->x.DistanceSquared(
               &baseline->second.first->second->node->x)) && (tmp
-              < A->second->node->x.Distance(
+              < A->second->node->x.DistanceSquared(
                   &baseline->second.second->second->node->x)))
             innerpoint++;
-          tmp = checker->second->node->x.Distance(
+          tmp = checker->second->node->x.DistanceSquared(
               &baseline->second.first->second->node->x);
-          if ((tmp < baseline->second.first->second->node->x.Distance(
+          if ((tmp < baseline->second.first->second->node->x.DistanceSquared(
               &A->second->node->x)) && (tmp
-              < baseline->second.first->second->node->x.Distance(
+              < baseline->second.first->second->node->x.DistanceSquared(
                   &baseline->second.second->second->node->x)))
             innerpoint++;
-          tmp = checker->second->node->x.Distance(
+          tmp = checker->second->node->x.DistanceSquared(
               &baseline->second.second->second->node->x);
-          if ((tmp < baseline->second.second->second->node->x.Distance(
+          if ((tmp < baseline->second.second->second->node->x.DistanceSquared(
               &baseline->second.first->second->node->x)) && (tmp
-              < baseline->second.second->second->node->x.Distance(
+              < baseline->second.second->second->node->x.DistanceSquared(
                   &A->second->node->x)))
             innerpoint++;
@@ -1483 +1549 @@
 ;
 
+
+double det_get(gsl_matrix *A, int inPlace) {
+  /*
+  inPlace = 1 => A is replaced with the LU decomposed copy.
+  inPlace = 0 => A is retained, and a copy is used for LU.
+  */
+
+  double det;
+  int signum;
+  gsl_permutation *p = gsl_permutation_alloc(A->size1);
+  gsl_matrix *tmpA;
+
+  if (inPlace)
+  tmpA = A;
+  else {
+  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+  gsl_matrix_memcpy(tmpA , A);
+  }
+
+
+  gsl_linalg_LU_decomp(tmpA , p , &signum);
+  det = gsl_linalg_LU_det(tmpA , signum);
+  gsl_permutation_free(p);
+  if (! inPlace)
+  gsl_matrix_free(tmpA);
+
+  return det;
+};
+
+void get_sphere(Vector *center, Vector &a, Vector &b, Vector &c, double RADIUS)
+{
+  gsl_matrix *A = gsl_matrix_calloc(3,3);
+  double m11, m12, m13, m14;
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m11 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, b.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m12 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, c.x[i]);
+  }
+  m13 = det_get(A, 1);
+
+  for(int i=0;i<3;i++) {
+    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
+    gsl_matrix_set(A, i, 1, a.x[i]);
+    gsl_matrix_set(A, i, 2, b.x[i]);
+  }
+  m14 = det_get(A, 1);
+
+  if (fabs(m11) < MYEPSILON)
+    cerr << "ERROR: three points are colinear." << endl;
+
+  center->x[0] =  0.5 * m12/ m11;
+  center->x[1] = -0.5 * m13/ m11;
+  center->x[2] =  0.5 * m14/ m11;
+
+  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
+    cerr << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+
+  gsl_matrix_free(A);
+};
+
+
+
 /**
  * Function returns center of sphere with RADIUS, which rests on points a, b, c
@@ -1489 +1632 @@
  * @param b vector second point of triangle
  * @param c vector third point of triangle
- * @param *Umkreismittelpunkt new cneter point of circumference
  * @param Direction vector indicates up/down
  * @param AlternativeDirection vecotr, needed in case the triangles have 90 deg angle
@@ -1500 +1642 @@
  * @param Umkreisradius double radius of circumscribing circle
  */
-
-  void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
-      double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
-  {
-    Vector TempNormal, helper;
-    double Restradius;
-    cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
-    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);
-    cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
-    // Here we calculated center of circumscribing circle, using barycentric coordinates
-    cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
-
-    TempNormal.CopyVector(&a);
-    TempNormal.SubtractVector(&b);
-    helper.CopyVector(&a);
-    helper.SubtractVector(&c);
-    TempNormal.VectorProduct(&helper);
-    if (fabs(HalfplaneIndicator) < MYEPSILON)
-      {
-        if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
-          {
-            TempNormal.Scale(-1);
-          }
-      }
-    else
-      {
-        if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
-          {
-            TempNormal.Scale(-1);
-          }
-      }
-
-    TempNormal.Normalize();
-    Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
-    cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
-    TempNormal.Scale(Restradius);
-    cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
-
-    Center->AddVector(&TempNormal);
-    cout << Verbose(4) << "Center of sphere of circumference is " << *Center << ".\n";
-    cout << Verbose(3) << "End of Get_center_of_sphere.\n";
-  }
-  ;
-
+void Get_center_of_sphere(Vector* Center, Vector a, Vector b, Vector c, Vector *NewUmkreismittelpunkt, Vector* Direction, Vector* AlternativeDirection,
+    double HalfplaneIndicator, double AlternativeIndicator, double alpha, double beta, double gamma, double RADIUS, double Umkreisradius)
+{
+  Vector TempNormal, helper;
+  double Restradius;
+  Vector OtherCenter;
+  cout << Verbose(3) << "Begin of Get_center_of_sphere.\n";
+  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);
+  cout << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
+  // Here we calculated center of circumscribing circle, using barycentric coordinates
+  cout << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
+
+  TempNormal.CopyVector(&a);
+  TempNormal.SubtractVector(&b);
+  helper.CopyVector(&a);
+  helper.SubtractVector(&c);
+  TempNormal.VectorProduct(&helper);
+  if (fabs(HalfplaneIndicator) < MYEPSILON)
+    {
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  else
+    {
+      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+        {
+          TempNormal.Scale(-1);
+        }
+    }
+  TempNormal.Normalize();
+  Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
+  cout << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";
+  TempNormal.Scale(Restradius);
+  cout << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
+
+  Center->AddVector(&TempNormal);
+  cout << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n";
+  get_sphere(&OtherCenter, a, b, c, RADIUS);
+  cout << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n";
+  cout << Verbose(3) << "End of Get_center_of_sphere.\n";
+};
 
 /** This recursive function finds a third point, to form a triangle with two given ones.
@@ -1580 +1721 @@
  * @param mol molecule structure with atoms and bonds
  */
-
 void Tesselation::Find_next_suitable_point_via_Angle_of_Sphere(atom* a, atom* b, atom* c, atom* Candidate, atom* Parent,
     int RecursionLevel, Vector *Chord, Vector *direction1, Vector *OldNormal, Vector ReferencePoint,
@@ -1633 +1773 @@
     }
 
-    if ((M_PI*4. > alpha*5.) && (M_PI*4. > beta*5.) && (M_PI*4 > gamma*5.)) {
+    if ((M_PI*179./180. > alpha) && (M_PI*179./180. > beta) && (M_PI*179./180. > gamma)) {
       Umkreisradius = SideA / 2.0 / sin(alpha);
       //cout << Umkreisradius << endl;
@@ -1655 +1795 @@
           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;
+        if (sign >= 0) {
+          cout << Verbose(3) << "Candidate is in search direction: " << sign << "." << endl;
+          Get_center_of_sphere(&Mittelpunkt, (a->x), (b->x), (Candidate->x), &NewUmkreismittelpunkt, OldNormal, direction1, sign, AlternativeSign, alpha, beta, gamma, RADIUS, Umkreisradius);
+          Mittelpunkt.SubtractVector(&ReferencePoint);
+          cout << Verbose(3) << "Reference vector to sphere's center is " << Mittelpunkt << "." << endl;
+          BallAngle = Mittelpunkt.Angle(OldNormal);
+          cout << Verbose(3) << "Angle between normal of base triangle and center of ball sphere is :" << BallAngle << "." << endl;
+
 
         cout << 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 ";
+          if ((Mittelpunkt.ScalarProduct(direction1) >=0) || (fabs(NewUmkreismittelpunkt.Norm()) < MYEPSILON)) {
+            if (Storage[0]< -1.5) { // first Candidate at all
+              if (1) {//if (CheckPresenceOfTriangle((ofstream *)&cout,a,b,Candidate)) {
+                cout << Verbose(2) << "First good candidate is " << *Candidate << " with ";
                 Opt_Candidate = Candidate;
                 Storage[0] = sign;
                 Storage[1] = AlternativeSign;
                 Storage[2] = BallAngle;
-                cout << " angle " << Storage[2] << " and Up/Down "
-                << Storage[0] << endl;
+                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;
+              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 << " refused due to Up/Down sign which is " << sign << endl;
+            cout << Verbose(3) << *Candidate << " is not in search direction." << endl;
           }
         }
@@ -1741 +1875 @@
 
 
-  /** This recursive function finds a third point, to form a triangle with two given ones.
-   * Two atoms are fixed, a candidate is supplied, additionally two vectors for direction distinction, a Storage area to \
-   *  supply results to the calling function, the radius of the sphere which the triangle shall support and the molecule \
-   *  upon which we operate.
-   *  If the candidate is more fitting to support the sphere than the already stored atom is, then we write its general \
-   *  direction and angle into Storage.
-   *  We the determine the recursive level we have reached and if this is not on the threshold yet, call this function again, \
-   *  with all neighbours of the candidate.
-   * @param a first point
-   * @param b second point
-   * @param Candidate base point along whose bonds to start looking from
-   * @param Parent point to avoid during search as its wrong direction
-   * @param RecursionLevel contains current recursion depth
-   * @param Chord baseline vector of first and second point
-   * @param d1 second in plane vector (along with \a Chord) of the triangle the baseline belongs to
-   * @param OldNormal normal of the triangle which the baseline belongs to
-   * @param Opt_Candidate candidate reference to return
-   * @param Opt_Mittelpunkt Centerpoint of ball, when resting on Opt_Candidate
-   * @param Storage array containing two angles of current Opt_Candidate
-   * @param RADIUS radius of ball
-   * @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, helper;
-
-  double CurrentEpsilon = 0.1;
-  double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
-  double BallAngle;
-  atom *Walker; // variable atom point
-
-
-  dif_a.CopyVector(&(a->x));
-  dif_a.SubtractVector(&(Candidate->x));
-  dif_b.CopyVector(&(b->x));
-  dif_b.SubtractVector(&(Candidate->x));
-  DirectionCheckPoint.CopyVector(&dif_a);
-  DirectionCheckPoint.Scale(-1);
-  DirectionCheckPoint.ProjectOntoPlane(Chord);
-
-  SideA = dif_b.Norm();
-  SideB = dif_a.Norm();
-  SideC = Chord->Norm();
-  //Chord->Scale(-1);
-
-  alpha = Chord->Angle(&dif_a);
-  beta = M_PI - Chord->Angle(&dif_b);
-  gamma = dif_a.Angle(&dif_b);
-
-
-  if (DEBUG) {
-    cout << "Atom number" << Candidate->nr << endl;
-    Candidate->x.Output((ofstream *) &cout);
-    cout << "number of bonds " << mol->NumberOfBondsPerAtom[Candidate->nr] << endl;
-  }
-
-  if (a != Candidate and b != Candidate) {
-    //      alpha = dif_a.Angle(&dif_b) / 2.;
-    //      SideA = Chord->Norm() / 2.;// (Chord->Norm()/2.) / sin(0.5*alpha);
-    //      SideB = dif_a.Norm();
-    //      centerline = SideA * SideA + SideB * SideB - 2. * SideA * SideB * cos(
-    //          alpha); // note this is squared of center line length
-    //      centerline = (Chord->Norm()/2.) / sin(0.5*alpha);
-    // Those are remains from Freddie. Needed?
-
-    Umkreisradius = SideA / 2.0 / sin(alpha);
-    //cout << Umkreisradius << endl;
-    //cout << SideB / 2.0 / sin(beta) << endl;
-    //cout << SideC / 2.0 / sin(gamma) << endl;
-
-    if (Umkreisradius < RADIUS && DirectionCheckPoint.ScalarProduct(&(Candidate->x))>0) { //Checking whether ball will at least rest o points.
-      // intermediate calculations to aquire centre of sphere, called Mittelpunkt:
-      Umkreismittelpunkt.Zero();
-      helper.CopyVector(&a->x);
-      helper.Scale(sin(2.*alpha));
-      Umkreismittelpunkt.AddVector(&helper);
-      helper.CopyVector(&b->x);
-      helper.Scale(sin(2.*beta));
-      Umkreismittelpunkt.AddVector(&helper);
-      helper.CopyVector(&Candidate->x);
-      helper.Scale(sin(2.*gamma));
-      Umkreismittelpunkt.AddVector(&helper);
-      //Umkreismittelpunkt = (a->x) * sin(2.*alpha) + b->x * sin(2.*beta) + (Candidate->x) * sin(2.*gamma) ;
-      Umkreismittelpunkt.Scale(1/(sin(2*alpha) + sin(2*beta) + sin(2*gamma)));
-
-      TempNormal.CopyVector(&dif_a);
-      TempNormal.VectorProduct(&dif_b);
-      if (TempNormal.ScalarProduct(OldNormal)<0 && sign>0 || TempNormal.ScalarProduct(OldNormal)>0 && sign<0) {
-        TempNormal.Scale(-1);
+/** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
+ * \param *Center new center on return
+ * \param *a first point
+ * \param *b second point
+ * \param *c third point
+ */
+void GetCenterofCircumcircle(Vector *Center, Vector *a, Vector *b, Vector *c)
+{
+  Vector helper;
+  double alpha, beta, gamma;
+  Vector SideA, SideB, SideC;
+  SideA.CopyVector(b);
+  SideA.SubtractVector(c);
+  SideB.CopyVector(c);
+  SideB.SubtractVector(a);
+  SideC.CopyVector(a);
+  SideC.SubtractVector(b);
+  alpha = M_PI - SideB.Angle(&SideC);
+  beta = M_PI - SideC.Angle(&SideA);
+  gamma = M_PI - SideA.Angle(&SideB);
+  cout << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
+  if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON)
+    cerr << "Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl;
+
+  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->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
+};
+
+/** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
+ * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
+ * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
+ * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
+ * \param CircleCenter Center of the parameter circle
+ * \param CirclePlaneNormal normal vector to plane of the parameter circle
+ * \param CircleRadius radius of the parameter circle
+ * \param NewSphereCenter new center of a circumcircle
+ * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
+ * \param NormalVector normal vector
+ * \param SearchDirection search direction to make angle unique on return.
+ * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
+ */
+double GetPathLengthonCircumCircle(Vector &CircleCenter, Vector &CirclePlaneNormal, double CircleRadius, Vector &NewSphereCenter, Vector &OldSphereCenter, Vector &NormalVector, Vector &SearchDirection)
+{
+  Vector helper;
+  double radius, alpha;
+
+  helper.CopyVector(&NewSphereCenter);
+  // test whether new center is on the parameter circle's plane
+  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+    cerr << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  }
+  radius = helper.ScalarProduct(&helper);
+  // test whether the new center vector has length of CircleRadius
+  if (fabs(radius - CircleRadius) > HULLEPSILON)
+    cerr << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+  alpha = helper.Angle(&OldSphereCenter);
+  // make the angle unique by checking the halfplanes/search direction
+  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
+    alpha = 2.*M_PI - alpha;
+  cout << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&OldSphereCenter);
+  helper.ProjectOntoPlane(&NormalVector);
+  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
+  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
+    cout << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    return alpha;
+  } else {
+    cout << Verbose(1) << "ERROR: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    return 2.*M_PI;
+  }
+};
+
+
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param BaseTriangle BoundaryTriangleSet of the current base triangle with all three points
+ * @param BaseLine BoundaryLineSet of BaseTriangle with the current base line
+ * @param OptCandidate candidate reference on return
+ * @param OptCandidateCenter candidate's sphere center on return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+// void Find_next_suitable_point(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+// {
+//   atom *Walker = NULL;
+//   Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+//   Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+//   Vector OldSphereCenter;   // center of the sphere defined by the three points of BaseTriangle
+//   Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+//   Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+//   Vector NewNormalVector;   // normal vector of the Candidate's triangle
+//   Vector SearchDirection;   // vector that points out of BaseTriangle and is orthonormal to its NormalVector (i.e. the desired direction for the best Candidate)
+//   Vector helper;
+//   LinkedAtoms *List = NULL;
+//   double CircleRadius; // radius of this circle
+//   double radius;
+//   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+//   double Nullalpha; // angle between OldSphereCenter and NormalVector of base triangle
+//   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+//   atom *Candidate = NULL;
+//
+//   cout << Verbose(1) << "Begin of Find_next_suitable_point" << endl;
+//
+//   cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << BaseTriangle->NormalVector << "." << endl;
+//
+//   // construct center of circle
+//   CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+//   CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+//   CircleCenter.Scale(0.5);
+//
+//   // construct normal vector of circle
+//   CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+//   CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+//
+//   // calculate squared radius of circle
+//   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+//   if (radius/4. < RADIUS*RADIUS) {
+//     CircleRadius = RADIUS*RADIUS - radius/4.;
+//     CirclePlaneNormal.Normalize();
+//     cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+//
+//     // construct old center
+//     GetCenterofCircumcircle(&OldSphereCenter, &(BaseTriangle->endpoints[0]->node->x), &(BaseTriangle->endpoints[1]->node->x), &(BaseTriangle->endpoints[2]->node->x));
+//     helper.CopyVector(&BaseTriangle->NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+//     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+//     helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//     OldSphereCenter.AddVector(&helper);
+//     OldSphereCenter.SubtractVector(&CircleCenter);
+//     cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+//
+//     // test whether old center is on the band's plane
+//     if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+//       cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+//       OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+//     }
+//     radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+//     if (fabs(radius - CircleRadius) < HULLEPSILON) {
+//
+//       // construct SearchDirection
+//       SearchDirection.MakeNormalVector(&BaseTriangle->NormalVector, &CirclePlaneNormal);
+//       helper.CopyVector(&BaseLine->endpoints[0]->node->x);
+//       for(int i=0;i<3;i++)  // just take next different endpoint
+//         if ((BaseTriangle->endpoints[i]->node != BaseLine->endpoints[0]->node) && (BaseTriangle->endpoints[i]->node != BaseLine->endpoints[1]->node)) {
+//           helper.SubtractVector(&BaseTriangle->endpoints[i]->node->x);
+//         }
+//       if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // ohoh, SearchDirection points inwards!
+//         SearchDirection.Scale(-1.);
+//       SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+//       SearchDirection.Normalize();
+//       cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+//       if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+//         cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+//       }
+//
+//       if (LC->SetIndexToVector(&CircleCenter)) {  // get cell for the starting atom
+//         for(int i=0;i<NDIM;i++) // store indices of this cell
+//           N[i] = LC->n[i];
+//         cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+//       } else {
+//         cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+//         return;
+//       }
+//       // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+//       cout << Verbose(2) << "LC Intervals:";
+//       for (int i=0;i<NDIM;i++) {
+//         Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+//         Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+//         cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+//       }
+//       cout << endl;
+//       for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+//         for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+//           for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+//             List = LC->GetCurrentCell();
+//             cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+//             if (List != NULL) {
+//               for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+//                 Candidate = (*Runner);
+//
+//                 // check for three unique points
+//                 if ((Candidate != BaseTriangle->endpoints[0]->node) && (Candidate != BaseTriangle->endpoints[1]->node) && (Candidate != BaseTriangle->endpoints[2]->node)) {
+//                   cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+//
+//                   // construct both new centers
+//                   GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+//                   OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+//
+//                   if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+//                     helper.CopyVector(&NewNormalVector);
+//                     cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+//                     radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+//                     if (radius < RADIUS*RADIUS) {
+//                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
+//                       cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+//                       NewSphereCenter.AddVector(&helper);
+//                       NewSphereCenter.SubtractVector(&CircleCenter);
+//                       cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+//
+//                       helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+//                       OtherNewSphereCenter.AddVector(&helper);
+//                       OtherNewSphereCenter.SubtractVector(&CircleCenter);
+//                       cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+//
+//                       // check both possible centers
+//                       alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, BaseTriangle->NormalVector, SearchDirection);
+//                       alpha = min(alpha, Otheralpha);
+//                       if (*ShortestAngle > alpha) {
+//                           OptCandidate = Candidate;
+//                           *ShortestAngle = alpha;
+//                           if (alpha != Otheralpha)
+//                             OptCandidateCenter->CopyVector(&NewSphereCenter);
+//                           else
+//                             OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+//                           cout << Verbose(1) << "We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+//                       } else {
+//                         if (OptCandidate != NULL)
+//                           cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+//                         else
+//                           cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+//                       }
+//
+//                     } else {
+//                       cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+//                     }
+//                   } else {
+//                     cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+//                   }
+//                 } else {
+//                   cout << Verbose(1) << "REJECT: Base triangle " << *BaseTriangle << " contains Candidate " << *Candidate << "." << endl;
+//                 }
+//               }
+//             }
+//           }
+//     } else {
+//       cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
+//     }
+//   } else {
+//     cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and base triangle " << *BaseTriangle << " is too big!" << endl;
+//   }
+//
+//   cout << Verbose(1) << "End of Find_next_suitable_point" << endl;
+// };
+
+
+/** This recursive function finds a third point, to form a triangle with two given ones.
+ * Note that this function is for the starting triangle.
+ * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
+ * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
+ * the center of the sphere is still fixed up to a single parameter. The band of possible values
+ * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
+ * us the "null" on this circle, the new center of the candidate point will be some way along this
+ * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
+ * by the normal vector of the base triangle that always points outwards by construction.
+ * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
+ * We construct the normal vector that defines the plane this circle lies in, it is just in the
+ * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
+ * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
+ * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
+ * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
+ * both.
+ * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
+ * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
+ * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
+ * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
+ * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
+ * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
+ * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa Find_starting_triangle())
+ * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
+ * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
+ * @param BaseLine BoundaryLineSet with the current base line
+ * @param ThirdNode third atom to avoid in search
+ * @param OptCandidate candidate reference on return
+ * @param OptCandidateCenter candidate's sphere center on return
+ * @param ShortestAngle the current path length on this circle band for the current Opt_Candidate
+ * @param RADIUS radius of sphere
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_third_point_for_Tesselation(Vector NormalVector, Vector SearchDirection, Vector OldSphereCenter, class BoundaryLineSet *BaseLine, atom *ThirdNode, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC)
+{
+  Vector CircleCenter;  // center of the circle, i.e. of the band of sphere's centers
+  Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
+  Vector NewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
+  Vector OtherNewSphereCenter;   // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
+  Vector NewNormalVector;   // normal vector of the Candidate's triangle
+  Vector helper;
+  LinkedAtoms *List = NULL;
+  double CircleRadius; // radius of this circle
+  double radius;
+  double alpha, Otheralpha; // angles (i.e. parameter for the circle).
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+  atom *Candidate = NULL;
+
+  cout << Verbose(1) << "Begin of Find_third_point_for_Tesselation" << endl;
+
+  cout << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl;
+
+  // construct center of circle
+  CircleCenter.CopyVector(&(BaseLine->endpoints[0]->node->x));
+  CircleCenter.AddVector(&BaseLine->endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&BaseLine->endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&BaseLine->endpoints[1]->node->x);
+
+  // calculate squared radius of circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    // test whether old center is on the band's plane
+    if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      cerr << "ERROR: Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = OldSphereCenter.ScalarProduct(&OldSphereCenter);
+    if (fabs(radius - CircleRadius) < HULLEPSILON) {
+
+      // check SearchDirection
+      cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+        cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
-      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;
+      // get cell for the starting atom
+      if (LC->SetIndexToVector(&CircleCenter)) {
+          for(int i=0;i<NDIM;i++) // store indices of this cell
+          N[i] = LC->n[i];
+        cout << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
       } else {
-        /*
-         * 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;
+        cerr << "ERROR: Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl;
+        return;
+      }
+      // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+      cout << Verbose(2) << "LC Intervals:";
+      for (int i=0;i<NDIM;i++) {
+        Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+        Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+        cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+      }
+      cout << endl;
+      for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+        for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+          for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+            List = LC->GetCurrentCell();
+            //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+            if (List != NULL) {
+              for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+                Candidate = (*Runner);
+
+                // check for three unique points
+                if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) && (Candidate != ThirdNode)) {
+                  cout << Verbose(1) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->x << "." << endl;
+
+                  // construct both new centers
+                  GetCenterofCircumcircle(&NewSphereCenter, &(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x));
+                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
+
+                  if ((NewNormalVector.MakeNormalVector(&(BaseLine->endpoints[0]->node->x), &(BaseLine->endpoints[1]->node->x), &(Candidate->x))) && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)) {
+                    helper.CopyVector(&NewNormalVector);
+                    cout << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
+                    radius = BaseLine->endpoints[0]->node->x.DistanceSquared(&NewSphereCenter);
+                    if (radius < RADIUS*RADIUS) {
+                      helper.Scale(sqrt(RADIUS*RADIUS - radius));
+                      cout << Verbose(3) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << "." << endl;
+                      NewSphereCenter.AddVector(&helper);
+                      NewSphereCenter.SubtractVector(&CircleCenter);
+                      cout << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
+
+                      helper.Scale(-1.); // OtherNewSphereCenter is created by the same vector just in the other direction
+                      OtherNewSphereCenter.AddVector(&helper);
+                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
+                      cout << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
+
+                      // check both possible centers
+                      alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
+                      alpha = min(alpha, Otheralpha);
+                      if (*ShortestAngle > alpha) {
+                          OptCandidate = Candidate;
+                          *ShortestAngle = alpha;
+                          if (alpha != Otheralpha)
+                            OptCandidateCenter->CopyVector(&NewSphereCenter);
+                          else
+                            OptCandidateCenter->CopyVector(&OtherNewSphereCenter);
+                          cout << Verbose(1) << "ACCEPT: We have found a better candidate: " << *OptCandidate << " with " << alpha << " and circumsphere's center at " << *OptCandidateCenter << "." << endl;
+                      } else {
+                        if (OptCandidate != NULL)
+                          cout << Verbose(1) << "REJECT: Old candidate: " << *OptCandidate << " is better than " << alpha << " with " << *ShortestAngle << "." << endl;
+                        else
+                          cout << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;
+                      }
+
+                    } else {
+                      cout << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " is too far away: " << radius << "." << endl;
+                    }
+                  } else {
+                    cout << Verbose(1) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
+                  }
+                } else {
+                  if (ThirdNode != NULL)
+                    cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;
+                  else
+                    cout << Verbose(1) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;
+                }
+              }
+            }
           }
-        }
-      }
-/*
-               * This is for checking point-angle and presence of Candidates in Ball, currently we are checking the ball Angle.
-               *
-                else
-                {
-                  if (sign>0 && BallAngle>0 && Storage[0]<0)
-                    {
-                      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;
-
-//Debugging purposes only
-                      cout << "Umkreismittelpunkt has coordinates" << Umkreismittelpunkt.x[0] << " "<< Umkreismittelpunkt.x[1] <<" "<<Umkreismittelpunkt.x[2] << endl;
-                      cout << "Candidate has coordinates" << Candidate->x.x[0]<< " " << Candidate->x.x[1] << " " << Candidate->x.x[2] << endl;
-                      cout << "a has coordinates" << a->x.x[0]<< " " << a->x.x[1] << " " << a->x.x[2] << endl;
-                      cout << "b has coordinates" << b->x.x[0]<< " " << b->x.x[1] << " " << b->x.x[2] << endl;
-                      cout << "Mittelpunkt has coordinates" << Mittelpunkt.x[0] << " " << Mittelpunkt.x[1]<< " "  <<Mittelpunkt.x[2] << endl;
-                      cout << "Umkreisradius ist " << Umkreisradius << endl;
-                      cout << "Restradius ist " << Restradius << endl;
-                      cout << "TempNormal has coordinates " << TempNormal.x[0] << " " << TempNormal.x[1] << " " << TempNormal.x[2] << " " << endl;
-                      cout << "OldNormal has coordinates " << OldNormal->x[0] << " " << OldNormal->x[1] << " " << OldNormal->x[2] << " " << endl;
-                      cout << "Dist a to UmkreisMittelpunkt " << a->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist b to UmkreisMittelpunkt " << b->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist Candidate to UmkreisMittelpunkt " << Candidate->x.Distance(&Umkreismittelpunkt) << endl;
-                      cout << "Dist a to Mittelpunkt " << a->x.Distance(&Mittelpunkt) << endl;
-                      cout << "Dist b to Mittelpunkt " << b->x.Distance(&Mittelpunkt) << endl;
-                      cout << "Dist Candidate to Mittelpunkt " << Candidate->x.Distance(&Mittelpunkt) << endl;
-
-
-
-                    }
-                  else
-                    {
-                      if (DEBUG)
-                        cout << "Looses to better candidate" << endl;
-                    }
-                }
-                */
     } else {
-      if (DEBUG) {
-        cout << "Doesn't satisfy requirements for circumscribing circle" << endl;
-      }
+      cerr << Verbose(1) << "ERROR: The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << 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.
- * @param out   output stream for debugging
- * @param tecplot output stream for writing found triangles in TecPlot format
- * @param mol molecule structure with all atoms and bonds
- * @param Line current baseline to search from
- * @param T current triangle which \a Line is edge of
- * @param RADIUS radius of the rolling ball
- * @param N number of found triangles
- * @param *filename filename base for intermediate envelopes
- */
-bool Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot,
-    molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
-    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;
+    if (ThirdNode != NULL)
+      cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;
+    else
+      cout << Verbose(1) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl;
+  }
+
+
+  cout << Verbose(1) << "End of Find_third_point_for_Tesselation" << endl;
 };
 
@@ -2194 +2375 @@
 
 
-void Find_second_point_for_Tesselation(atom* a, atom* Candidate, atom* Parent,
-    int RecursionLevel, Vector Oben, atom*& Opt_Candidate, double Storage[3],
-    molecule* mol, double RADIUS)
-{
-  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation, recursive level " << RecursionLevel << endl;
-  int i;
+/** Finds the second point of starting triangle.
+ * \param *a first atom
+ * \param *Candidate pointer to candidate atom on return
+ * \param Oben vector indicating the outside
+ * \param Opt_Candidate reference to recommended candidate on return
+ * \param Storage[3] array storing angles and other candidate information
+ * \param RADIUS radius of virtual sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, Vector Oben, atom*& Opt_Candidate, double Storage[3], double RADIUS, LinkedCell *LC)
+{
+  cout << Verbose(2) << "Begin of Find_second_point_for_Tesselation" << endl;
   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;
+  LinkedAtoms *List = NULL;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LC->SetIndexToAtom(a)) {  // get cell for the starting atom
+    for(int i=0;i<NDIM;i++) // store indices of this cell
+      N[i] = LC->n[i];
+  } else {
+    cerr << "ERROR: Atom " << *a << " is not found in cell " << LC->index << "." << endl;
+    return;
+  }
+  // then go through the current and all neighbouring cells and check the contained atoms for possible candidates
+  cout << Verbose(2) << "LC Intervals:";
+  for (int i=0;i<NDIM;i++) {
+    Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0;
+    Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1;
+    cout << " [" << Nlower[i] << "," << Nupper[i] << "] ";
+  }
+  cout << endl;
+  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
+    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
+      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            Candidate = (*Runner);
+                // check 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;
+              }
+            }
+          }
+        }
       }
-    } 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) << "End of Find_second_point_for_Tesselation" << endl;
 };
 
-void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
+/** Finds the starting triangle for find_non_convex_border().
+ * Looks at the outermost atom per axis, then Find_second_point_for_Tesselation()
+ * for the second and Find_next_suitable_point_via_Angle_of_Sphere() for the third
+ * point are called.
+ * \param RADIUS radius of virtual rolling sphere
+ * \param *LC LinkedCell structure with neighbouring atoms
+ */
+void Tesselation::Find_starting_triangle(ofstream *out, molecule *mol, const double RADIUS, LinkedCell *LC)
 {
   cout << Verbose(1) << "Begin of Find_starting_triangle\n";
   int i = 0;
-  atom* Walker;
+  LinkedAtoms *List = NULL;
   atom* FirstPoint;
   atom* SecondPoint;
-  atom* max_index[NDIM];
+  atom* MaxAtom[NDIM];
   double max_coordinate[NDIM];
   Vector Oben;
   Vector helper;
   Vector Chord;
-  Vector CenterOfFirstLine;
+  Vector SearchDirection;
+  Vector OptCandidateCenter;
 
   Oben.Zero();
 
   for (i = 0; i < 3; i++) {
-    max_index[i] = NULL;
+    MaxAtom[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;
+  for (int i=0;i<NDIM;i++) { // each axis
+    LC->n[i] = LC->N[i]-1; // current axis is topmost cell
+    for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++)
+      for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) {
+        List = LC->GetCurrentCell();
+        //cout << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
+        if (List != NULL) {
+          for (LinkedAtoms::iterator Runner = List->begin();Runner != List->end();Runner++) {
+            cout << Verbose(2) << "Current atom is " << *(*Runner) << "." << endl;
+            if ((*Runner)->x.x[i] > max_coordinate[i]) {
+              max_coordinate[i] = (*Runner)->x.x[i];
+              MaxAtom[i] = (*Runner);
+            }
+          }
+        } else {
+          cerr << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
       }
-    }
   }
 
   cout << Verbose(2) << "Found maximum coordinates: ";
   for (int i=0;i<NDIM;i++)
-    cout << i << ": " << *max_index[i] << "\t";
+    cout << i << ": " << *MaxAtom[i] << "\t";
   cout << endl;
-  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k = 1;
+  const int k = 1;    // arbitrary choice
   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];
+  FirstPoint = MaxAtom[k];
+  cout << Verbose(1) << "Coordinates of start atom " << *FirstPoint << " at " << FirstPoint->x << "." << endl;
+
+  // add first point
+  AddTrianglePoint(FirstPoint, 0);
+
+  double ShortestAngle;
   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_...
+  ShortestAngle = 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.
+
+  Find_second_point_for_Tesselation(FirstPoint, NULL, Oben, Opt_Candidate, &ShortestAngle, RADIUS, LC); // 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";
+
+  // add second point and first baseline
+  AddTrianglePoint(SecondPoint, 1);
+  AddTriangleLine(TPS[0], TPS[1], 0);
 
   helper.CopyVector(&(FirstPoint->x));
@@ -2300 +2523 @@
   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.;
+  ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
 
   Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
   Chord.SubtractVector(&(SecondPoint->x));
+  double radius = Chord.ScalarProduct(&Chord);
+  double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+  helper.CopyVector(&Oben);
+  helper.Scale(CircleRadius);
   // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized)
 
@@ -2311 +2536 @@
   // 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);
+  SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+
+  cout << Verbose(1) << "Looking for third point candidates ...\n";
+  Find_third_point_for_Tesselation(Oben, SearchDirection, helper, BLS[0], NULL, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
   cout << Verbose(1) << "Third Point is " << *Opt_Candidate << endl;
 
+  // add third point
+  AddTrianglePoint(Opt_Candidate, 2);
+
   // 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);
+  // Finally, we only have to add the found further lines
   AddTriangleLine(TPS[1], TPS[2], 1);
   AddTriangleLine(TPS[0], TPS[2], 2);
@@ -2336 +2554 @@
   AddTriangleToLines();
   // ... and calculate its normal vector (with correct orientation)
-  Oben.Scale(-1.);
-  BTS->GetNormalVector(Oben);
+  OptCandidateCenter.Scale(-1.);
+  cout << Verbose(2) << "Oben is currently " << OptCandidateCenter << "." << endl;
+  BTS->GetNormalVector(OptCandidateCenter);
   cout << Verbose(0) << "==> The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << " with normal vector " << BTS->NormalVector << ".\n";
+  cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
   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)
+/** This function finds a triangle to a line, adjacent to an existing one.
+ * @param out output stream for debugging
+ * @param *mol molecule with Atom's and Bond's
+ * @param Line current baseline to search from
+ * @param T current triangle which \a Line is edge of
+ * @param RADIUS radius of the rolling ball
+ * @param N number of found triangles
+ * @param *filename filename base for intermediate envelopes
+ * @param *LC LinkedCell structure with neighbouring atoms
+ */
+bool Tesselation::Find_next_suitable_triangle(ofstream *out,
+    molecule *mol, BoundaryLineSet &Line, BoundaryTriangleSet &T,
+    const double& RADIUS, int N, const char *tempbasename, LinkedCell *LC)
+{
+  cout << Verbose(1) << "Begin of Find_next_suitable_triangle\n";
+  ofstream *tempstream = NULL;
+  char NumberName[255];
+
+  atom* Opt_Candidate = NULL;
+  Vector OptCandidateCenter;
+
+  Vector CircleCenter;
+  Vector CirclePlaneNormal;
+  Vector OldSphereCenter;
+  Vector SearchDirection;
+  Vector helper;
+  atom *ThirdNode = NULL;
+  double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
+  double radius, CircleRadius;
+
+  cout << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl;
+  for (int i=0;i<3;i++)
+    if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node))
+      ThirdNode = T.endpoints[i]->node;
+
+  // construct center of circle
+  CircleCenter.CopyVector(&Line.endpoints[0]->node->x);
+  CircleCenter.AddVector(&Line.endpoints[1]->node->x);
+  CircleCenter.Scale(0.5);
+
+  // construct normal vector of circle
+  CirclePlaneNormal.CopyVector(&Line.endpoints[0]->node->x);
+  CirclePlaneNormal.SubtractVector(&Line.endpoints[1]->node->x);
+
+  // calculate squared radius of circle
+  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  if (radius/4. < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius/4.;
+    CirclePlaneNormal.Normalize();
+    cout << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+
+    // construct old center
+    GetCenterofCircumcircle(&OldSphereCenter, &(T.endpoints[0]->node->x), &(T.endpoints[1]->node->x), &(T.endpoints[2]->node->x));
+    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
+    radius = Line.endpoints[0]->node->x.DistanceSquared(&OldSphereCenter);
+    helper.Scale(sqrt(RADIUS*RADIUS - radius));
+    OldSphereCenter.AddVector(&helper);
+    OldSphereCenter.SubtractVector(&CircleCenter);
+    cout << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+
+    // construct SearchDirection
+    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
+    helper.CopyVector(&Line.endpoints[0]->node->x);
+    helper.SubtractVector(&ThirdNode->x);
+    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // ohoh, SearchDirection points inwards!
+      SearchDirection.Scale(-1.);
+    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
+    SearchDirection.Normalize();
+    cout << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      cerr << "ERROR: SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
+    }
+
+    // add third point
+    cout << Verbose(1) << "Looking for third point candidates for triangle ... " << endl;
+    Find_third_point_for_Tesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, Opt_Candidate, &OptCandidateCenter, &ShortestAngle, RADIUS, LC);
+
+  } else {
+    cout << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl;
+  }
+
+  if (Opt_Candidate == NULL) {
+    cerr << "WARNING: Could not find a suitable candidate." << endl;
+    return false;
+  }
+  cout << Verbose(1) << " Optimal candidate is " << *Opt_Candidate << " with circumsphere's center at " << OptCandidateCenter << "." << 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();
+
+    OptCandidateCenter.Scale(-1.);
+    BTS->GetNormalVector(OptCandidateCenter);
+    OptCandidateCenter.Scale(-1.);
+
+    cout << "--> 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 (flag && (DoSingleStepOutput && (TrianglesOnBoundaryCount % 10 == 0))) { // if we have a new triangle and want to output each new triangle configuration
+    sprintf(NumberName, "-%04d-%s_%s_%s", TriangleFilesWritten, BTS->endpoints[0]->node->Name, BTS->endpoints[1]->node->Name, BTS->endpoints[2]->node->Name);
+    if (DoTecplotOutput) {
+      string NameofTempFile(tempbasename);
+      NameofTempFile.append(NumberName);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+      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);
+      for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
+        NameofTempFile.erase(npos, 1);
+      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);
+      // include the current position of the virtual sphere in the temporary raster3d file
+      // make the circumsphere's center absolute again
+      helper.CopyVector(&Line.endpoints[0]->node->x);
+      helper.AddVector(&Line.endpoints[1]->node->x);
+      helper.Scale(0.5);
+      OptCandidateCenter.AddVector(&helper);
+      Vector *center = mol->DetermineCenterOfAll(out);
+      OptCandidateCenter.AddVector(center);
+      delete(center);
+      // and add to file plus translucency object
+      *tempstream << "# current virtual sphere\n";
+      *tempstream << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
+      *tempstream << "2\n  " << OptCandidateCenter.x[0] << " " << OptCandidateCenter.x[1] << " " << OptCandidateCenter.x[2] << "\t" << RADIUS << "\t1 0 0\n";
+      *tempstream << "9\n  terminating special property\n";
+      tempstream->close();
+      tempstream->flush();
+      delete(tempstream);
+    }
+    if (DoTecplotOutput || DoRaster3DOutput)
+      TriangleFilesWritten++;
+  }
+
+  cout << Verbose(1) << "End of Find_next_suitable_triangle\n";
+  return true;
+};
+
+/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
+ * \param *out output stream for debugging
+ * \param *mol molecule structure with Atom's and Bond's
+ * \param *Tess Tesselation filled with points, lines and triangles on boundary on return
+ * \param *LCList linked cell list of all atoms
+ * \param *filename filename prefix for output of vertex data
+ * \para RADIUS radius of the virtual sphere
+ */
+void Find_non_convex_border(ofstream *out, molecule* mol, class Tesselation *Tess, class LinkedCell *LCList, 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;
+  bool freeTess = false;
+  *out << Verbose(1) << "Entering search for non convex hull. " << endl;
+  if (Tess == NULL) {
+    *out << Verbose(1) << "Allocating Tesselation struct ..." << endl;
+    Tess = new Tesselation;
+    freeTess = true;
+  }
+  bool freeLC = false;
   LineMap::iterator baseline;
-  cout << Verbose(0) << "Begin of Find_non_convex_border\n";
+  *out << 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);
+
+  if (LCList == NULL) {
+    LCList = new LinkedCell(mol, 2.*RADIUS);
+    freeLC = true;
+  }
+
+  Tess->Find_starting_triangle(out, mol, RADIUS, LCList);
 
   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.
+      failflag = Tess->Find_next_suitable_triangle(out, mol, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, filename, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
@@ -2374 +2778 @@
     }
   }
-  if (failflag) {
-    cout << Verbose(1) << "Writing final tecplot file\n";
-    if (DoTecplotOutput)
+  if (1) { //failflag) {
+    *out << Verbose(1) << "Writing final tecplot file\n";
+    if (DoTecplotOutput) {
+      string OutputName(filename);
+      OutputName.append(TecplotSuffix);
+      ofstream *tecplot = new ofstream(OutputName.c_str());
       write_tecplot_file(out, tecplot, Tess, mol, -1);
-    if (DoRaster3DOutput)
-      write_raster3d_file(out, tecplot, Tess, mol);
+      tecplot->close();
+      delete(tecplot);
+    }
+    if (DoRaster3DOutput) {
+      string OutputName(filename);
+      OutputName.append(Raster3DSuffix);
+      ofstream *raster = new ofstream(OutputName.c_str());
+      write_raster3d_file(out, raster, Tess, mol);
+      raster->close();
+      delete(raster);
+    }
   } else {
     cerr << "ERROR: Could definately not find all necessary triangles!" << endl;
   }
-
-  cout << Verbose(0) << "End of Find_non_convex_border\n";
-  delete(Tess);
+  if (freeTess)
+    delete(Tess);
+  if (freeLC)
+    delete(LCList);
+  *out << Verbose(0) << "End of Find_non_convex_border\n";
 };
 
+/** Finds a hole of sufficient size in \a this molecule to embed \a *srcmol into it.
+ * \param *out output stream for debugging
+ * \param *srcmol molecule to embed into
+ * \return *Vector new center of \a *srcmol for embedding relative to \a this
+ */
+Vector* molecule::FindEmbeddingHole(ofstream *out, molecule *srcmol)
+{
+  Vector *Center = new Vector;
+  Center->Zero();
+  // calculate volume/shape of \a *srcmol
+
+  // find embedding holes
+
+  // if more than one, let user choose
+
+  // return embedding center
+  return Center;
+};
+