Changeset e0c5b1 for molecuilder/src/boundary.cpp

Timestamp:

Dec 3, 2008, 2:12:05 PM (17 years ago)

Author:

Christian Neuen <neuen@…>

Children:

3e20fe

Parents:

f5b58e

Message:

Multiple changes to boundary, currently not fully operational.
Molecules has a new routine to create adjacency lists, reading bonds from a dbond file instead of looking for the distances by itself.
Vector function Project onto plane has been updated.

File:

• molecuilder/src/boundary.cpp (modified) (62 diffs)

molecuilder/src/boundary.cpp

@@ -3 +3 @@
 
 
-// =================spaeter gebrauchte
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol);
-void Find_non_convex_border(Tesselation Tess, molecule mol);
 
 
@@ -29 +26 @@
 };
 
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 
+void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
 {
   cout << Verbose(6) << "Adding line " << *line << " to " << *this << "." << endl;
@@ -40 +37 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryPointSet &a) 
+ostream & operator << (ostream &ost, BoundaryPointSet &a)
 {
   ost << "[" << a.Nr << "|" << a.node->Name << "]";
@@ -73 +70 @@
 {
   for (int i=0;i<2;i++) {
-    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 
+    cout << Verbose(5) << "Erasing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;
     endpoints[i]->lines.erase(Nr);
     LineMap::iterator tester = endpoints[i]->lines.begin();
@@ -85 +82 @@
 };
 
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 
+void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
   cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
@@ -92 +89 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryLineSet &a) 
+ostream & operator << (ostream &ost, BoundaryLineSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "]";
@@ -125 +122 @@
     for (int j=0;j<2;j++) { // for both endpoints
       OrderMap.insert ( pair <int, class BoundaryPointSet * >( line[i]->endpoints[j]->Nr, line[i]->endpoints[j]) );
-      // and we don't care whether insertion fails 
+      // and we don't care whether insertion fails
     }
   // set endpoints
@@ -161 +158 @@
   // get normal vector
   NormalVector.MakeNormalVector(&endpoints[0]->node->x, &endpoints[1]->node->x, &endpoints[2]->node->x);
-  
+
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
   if (endpoints[0]->node->x.Projection(&NormalVector) > 0)
@@ -167 +164 @@
 };
 
-ostream & operator << (ostream &ost, BoundaryTriangleSet &a) 
+ostream & operator << (ostream &ost, BoundaryTriangleSet &a)
 {
   ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]";
@@ -213 +210 @@
   LineMap LinesOnBoundary;
   TriangleMap TrianglesOnBoundary;
-  
+
   *out << Verbose(1) << "Finding all boundary points." << endl;
   Boundaries *BoundaryPoints = new Boundaries [NDIM];  // first is alpha, second is (r, nr)
@@ -250 +247 @@
       else
         angle = 0.;  // otherwise it's a vector in Axis Direction and unimportant for boundary issues
-        
+
       //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
       if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
@@ -268 +265 @@
         Walker->x.Output(out);
         *out << endl;
-        double tmp = ProjectedVector.Norm(); 
+        double tmp = ProjectedVector.Norm();
         if (tmp > BoundaryTestPair.first->second.first) {
           BoundaryTestPair.first->second.first = tmp;
-          BoundaryTestPair.first->second.second = Walker; 
+          BoundaryTestPair.first->second.second = Walker;
           *out << Verbose(2) << "Keeping new vector." << endl;
         } else if (tmp == BoundaryTestPair.first->second.first) {
@@ -319 +316 @@
         }
         // check distance
-        
+
         // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
         {
@@ -328 +325 @@
   //          SideA.Output(out);
   //          *out << endl;
-          
+
           SideB.CopyVector(&right->second.second->x);
           SideB.ProjectOntoPlane(&AxisVector);
@@ -334 +331 @@
   //          SideB.Output(out);
   //          *out << endl;
-          
+
           SideC.CopyVector(&left->second.second->x);
           SideC.SubtractVector(&right->second.second->x);
@@ -341 +338 @@
   //          SideC.Output(out);
   //          *out << endl;
-  
+
           SideH.CopyVector(&runner->second.second->x);
           SideH.ProjectOntoPlane(&AxisVector);
@@ -347 +344 @@
   //          SideH.Output(out);
   //          *out << endl;
-          
+
           // calculate each length
           double a = SideA.Norm();
@@ -355 +352 @@
           // calculate the angles
           double alpha = SideA.Angle(&SideH);
-          double beta = SideA.Angle(&SideC); 
+          double beta = SideA.Angle(&SideC);
           double gamma = SideB.Angle(&SideH);
           double delta = SideC.Angle(&SideH);
           double MinDistance = a * sin(beta)/(sin(delta)) * (((alpha < M_PI/2.) || (gamma < M_PI/2.)) ? 1. : -1.);
-  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl; 
+  //          *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
           //*out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
           if ((fabs(h/fabs(h) - MinDistance/fabs(MinDistance)) < MYEPSILON) && (h <  MinDistance)) {
@@ -370 +367 @@
       }
     } while (flag);
-  } 
+  }
   return BoundaryPoints;
 };
@@ -381 +378 @@
  * \param IsAngstroem whether we have angstroem or atomic units
  * \return NDIM array of the diameters
- */ 
+ */
 double * GetDiametersOfCluster(ofstream *out, Boundaries *BoundaryPtr, molecule *mol, bool IsAngstroem)
 {
@@ -393 +390 @@
     *out << Verbose(1) << "Using given boundary points set." << endl;
   }
-  
+
   // determine biggest "diameter" of cluster for each axis
   Boundaries::iterator Neighbour, OtherNeighbour;
@@ -418 +415 @@
         DistanceVector.SubtractVector(&Neighbour->second.second->x);
         do {  // seek for neighbour pair where it flips
-          OldComponent = DistanceVector.x[Othercomponent]; 
+          OldComponent = DistanceVector.x[Othercomponent];
           Neighbour++;
           if (Neighbour == BoundaryPoints[axis].end())  // make it wrap around
@@ -431 +428 @@
             OtherNeighbour = BoundaryPoints[axis].end();
           OtherNeighbour--;
-          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 
+          //*out << Verbose(2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
           // now we have found the pair: Neighbour and OtherNeighbour
           OtherVector.CopyVector(&runner->second.second->x);
@@ -466 +463 @@
  * \param *BoundaryPoints NDIM set of boundary points on the projected plane per axis, on return if desired
  * \param *mol molecule structure representing the cluster
- * \return determined volume of the cluster in cubed config:GetIsAngstroem() 
+ * \return determined volume of the cluster in cubed config:GetIsAngstroem()
  */
-double VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration, Boundaries *BoundaryPtr, molecule *mol) 
+double VolumeOfConvexEnvelope(ofstream *out, ofstream *tecplot, config *configuration, Boundaries *BoundaryPtr, molecule *mol)
 {
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -477 +474 @@
   double volume = 0.;
   double PyramidVolume = 0.;
-  double G,h; 
+  double G,h;
   Vector x,y;
   double a,b,c;
 
-  Find_non_convex_border(*TesselStruct, *mol);
+  Find_non_convex_border(*TesselStruct, mol);
   /*
   // 1. calculate center of gravity
   *out << endl;
   Vector *CenterOfGravity = mol->DetermineCenterOfGravity(out);
-  
+
   // 2. translate all points into CoG
   *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
@@ -494 +491 @@
     Walker->x.Translate(CenterOfGravity);
   }
-  
+
   // 3. Find all points on the boundary
   if (BoundaryPoints == NULL) {
@@ -502 +499 @@
     *out << Verbose(1) << "Using given boundary points set." << endl;
   }
-  
+
   // 4. fill the boundary point list
   for (int axis=0;axis<NDIM;axis++)
@@ -518 +515 @@
 //  }
 //  *out << endl;
-  
+
   // 5a. guess starting triangle
   TesselStruct->GuessStartingTriangle(out);
-  
+
   // 5b. go through all lines, that are not yet part of two triangles (only of one so far)
   TesselStruct->TesselateOnBoundary(out, configuration, mol);
@@ -545 +542 @@
     PyramidVolume = (1./3.) * G * h;    // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
     *out << Verbose(2) << "Area of triangle is " << G << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-    volume += PyramidVolume; 
+    volume += PyramidVolume;
   }
   *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
@@ -572 +569 @@
     for (int i=0;i<mol->AtomCount;i++)
       LookupList[i] = -1;
-    
+
     // print atom coordinates
     *out << Verbose(2) << "The following triangles were created:";
@@ -586 +583 @@
     for (TriangleMap::iterator runner =  TesselStruct->TrianglesOnBoundary.begin(); runner !=  TesselStruct->TrianglesOnBoundary.end(); runner++) {
       *out << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
-      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 
+      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
     }
     delete[](LookupList);
@@ -595 +592 @@
   if (BoundaryFreeFlag)
     delete[](BoundaryPoints);
-  
+
   return volume;
 };
@@ -612 +609 @@
   // transform to PAS
   mol->PrincipalAxisSystem(out, true);
-  
+
   // some preparations beforehand
   bool IsAngstroem = configuration->GetIsAngstroem();
@@ -619 +616 @@
   if (ClusterVolume == 0)
     clustervolume = VolumeOfConvexEnvelope(out, NULL, configuration, BoundaryPoints, mol);
-  else 
+  else
     clustervolume = ClusterVolume;
   double *GreatestDiameter = GetDiametersOfCluster(out, BoundaryPoints, mol, IsAngstroem);
@@ -636 +633 @@
   }
   *out << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl;
-  
+
   *out << Verbose(0) << "RESULT: The average density is " << setprecision(10) << totalmass/clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  
+
   // solve cubic polynomial
   *out << Verbose(1) << "Solving equidistant suspension in water problem ..." << endl;
@@ -647 +644 @@
     cellvolume = (TotalNoClusters*totalmass/SOLVENTDENSITY_a0 - (totalmass/clustervolume))/(celldensity-1);
   *out << Verbose(1) << "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
-  
+
   double minimumvolume = TotalNoClusters*(GreatestDiameter[0]*GreatestDiameter[1]*GreatestDiameter[2]);
   *out << Verbose(1) << "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
@@ -658 +655 @@
   } else {
     BoxLengths.x[0] = (repetition[0]*GreatestDiameter[0] + repetition[1]*GreatestDiameter[1] + repetition[2]*GreatestDiameter[2]);
-    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1] 
-              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2] 
+    BoxLengths.x[1] = (repetition[0]*repetition[1]*GreatestDiameter[0]*GreatestDiameter[1]
+              + repetition[0]*repetition[2]*GreatestDiameter[0]*GreatestDiameter[2]
               + repetition[1]*repetition[2]*GreatestDiameter[1]*GreatestDiameter[2]);
     BoxLengths.x[2] = minimumvolume - cellvolume;
@@ -669 +666 @@
       x0 = x2;  // sorted in ascending order
     }
-  
+
     cellvolume = 1;
     for(int i=0;i<NDIM;i++) {
@@ -675 +672 @@
       cellvolume *= BoxLengths.x[i];
     }
-  
+
     // set new box dimensions
     *out << Verbose(0) << "Translating to box with these boundaries." << endl;
@@ -692 +689 @@
 Tesselation::Tesselation()
 {
-  PointsOnBoundaryCount = 0; 
-  LinesOnBoundaryCount = 0; 
+  PointsOnBoundaryCount = 0;
+  LinesOnBoundaryCount = 0;
   TrianglesOnBoundaryCount = 0;
 };
@@ -711 +708 @@
  * \param *out output stream for debugging
  * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
- */ 
+ */
 void Tesselation::GuessStartingTriangle(ofstream *out)
 {
@@ -722 +719 @@
   A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
 
-  // with A chosen, take each pair B,C and sort 
+  // with A chosen, take each pair B,C and sort
   if (A != PointsOnBoundary.end()) {
     B = A;
@@ -746 +743 @@
 //    }
 //    *out << endl;
-  
+
   // 4b2. pick three baselines forming a triangle
   // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
@@ -758 +755 @@
     PlaneVector.Output(out);
     *out << endl;
-    // 4. loop over all points 
+    // 4. loop over all points
     double sign = 0.;
     PointMap::iterator checker = PointsOnBoundary.begin();
@@ -785 +782 @@
       tmp = checker->second->node->x.Distance(&A->second->node->x);
       int innerpoint = 0;
-      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+      if ((tmp < A->second->node->x.Distance(&baseline->second.first->second->node->x))
           && (tmp < A->second->node->x.Distance(&baseline->second.second->second->node->x)))
         innerpoint++;
       tmp = checker->second->node->x.Distance(&baseline->second.first->second->node->x);
-      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x)) 
+      if ((tmp < baseline->second.first->second->node->x.Distance(&A->second->node->x))
           && (tmp < baseline->second.first->second->node->x.Distance(&baseline->second.second->second->node->x)))
         innerpoint++;
       tmp = checker->second->node->x.Distance(&baseline->second.second->second->node->x);
-      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x)) 
+      if ((tmp < baseline->second.second->second->node->x.Distance(&baseline->second.first->second->node->x))
           && (tmp < baseline->second.second->second->node->x.Distance(&A->second->node->x)))
         innerpoint++;
@@ -815 +812 @@
     BPS[0] = baseline->second.first->second;
     BPS[1] = A->second;
-    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
-   
+    BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+
     // 4b3. insert created triangle
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
@@ -825 +822 @@
       LinesOnBoundaryCount++;
     }
-  
+
     *out << Verbose(1) << "Starting triangle is " << *BTS << "." << endl;
   } else {
@@ -857 +854 @@
   do {
     flag = false;
-    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) 
+    for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
       if (baseline->second->TrianglesCount == 1) {
-        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl; 
+        *out << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl;
         // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
         SmallestAngle = M_PI;
@@ -888 +885 @@
         TempVector.CopyVector(&CenterVector);
         TempVector.SubtractVector(&baseline->second->endpoints[0]->node->x);  // TempVector is vector on triangle plane pointing from one baseline egde towards center!
-        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl; 
+        //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
         if (PropagationVector.Projection(&TempVector) > 0)  // make sure normal propagation vector points outward from baseline
           PropagationVector.Scale(-1.);
@@ -895 +892 @@
         *out << endl;
         winner = PointsOnBoundary.end();
-        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) 
+        for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++)
           if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
             *out << Verbose(3) << "Target point is " << *(target->second) << ":";
             bool continueflag = true;
-            
+
             VirtualNormalVector.CopyVector(&baseline->second->endpoints[0]->node->x);
             VirtualNormalVector.AddVector(&baseline->second->endpoints[0]->node->x);
@@ -935 +932 @@
             // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
             continueflag = continueflag && (!(
-                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) 
+                ((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end())
                 && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak))
                ));
@@ -990 +987 @@
           *out << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;
         }
-  
+
         // 5d. If the set of lines is not yet empty, go to 5. and continue
       } else
         *out << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->TrianglesCount << "." << endl;
   } while (flag);
-  
+
 };
 
@@ -1021 +1018 @@
 //====================================================================================================================
 
-
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule mol)
-{
-  /* d2 ist der Normalenvektor auf dem Dreieck,
-   * d1 ist der Vektor, der normal auf der Kante und d2 steht.
+/*!
+ * 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 id, 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.
+ */
+void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, double *Storage, const double RADIUS, molecule* mol)
+{
+  /* d2 is normal vector on the triangle
+   * d1 is normal on the triangle line, from which we come, as well as on d2.
    */
-  Vector dif_a;
-  Vector dif_b;
-  Vector Chord;
-  Vector AngleCheck;
-  atom *Walker;
-
-  dif_a.CopyVector(&a.x);
-  dif_a.SubtractVector(&Candidate.x);
-  dif_b.CopyVector(&b.x);
-  dif_b.SubtractVector(&Candidate.x);
-  Chord.CopyVector(&a.x);
-  Chord.SubtractVector(&b.x);
+  Vector dif_a; //Vector from a to candidate
+  Vector dif_b; //Vector from b to candidate
+  Vector AngleCheck; // Projection of a difference vector on plane orthogonal on triangle side.
+  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));
   AngleCheck.CopyVector(&dif_a);
-  AngleCheck.ProjectOntoPlane(&Chord);
-
-  //Storage eintrag fuer aktuelles Atom
-  if (Chord.Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
+  AngleCheck.ProjectOntoPlane(Chord);
+
+  if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
   {
 
     if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[1]) //This will give absolute preference to those in "right-hand" quadrants
       {
-        Storage[0]=(double)Candidate.nr;
-        Storage[1]=1;
-        Storage[2]=AngleCheck.Angle(d2);
+        Storage[0]=(double)Candidate->nr;
+        Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+        Storage[2]=AngleCheck.Angle(d2);
       }
     else
@@ -1058 +1060 @@
           //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
           {
-            Storage[0]=(double)Candidate.nr;
+            Storage[0]=(double)Candidate->nr;
             Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
             Storage[2]=AngleCheck.Angle(d2);
@@ -1065 +1067 @@
   }
 
-  if (n<5)
-    {
-      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+  if (n<5) // Five is the recursion level threshold.
+    {
+      for(int i=0; i<mol->NumberOfBondsPerAtom[Candidate->nr];i++) // go through all bond
         {
-          while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
-            {
+                Walker= mol->start; // go through all atoms
+
+          while (Walker->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr))
+            { // until atom found which belongs to bond
               Walker = Walker->next;
             }
 
-          Find_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
+
+          Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Storage, RADIUS, mol);  //call function again
         }
     }
 };
 
-
-void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
-{
-  Vector CenterOfLine = Line.endpoints[0]->node->x;
+/*!
+ * this function fins a triangle to a line, adjacent to an existing one.
+ */
+void Tesselation::Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
+{
+        printf("Looking for next suitable triangle \n");
   Vector direction1;
-  Vector direction2;
   Vector helper;
+  Vector Chord;
   atom* Walker;
 
-  double Storage[3];
+  double *Storage;
+  Storage = new double[3];
   Storage[0]=-1.;   // Id must be positive, we see should nothing be done
   Storage[1]=-1.;   // This direction is either +1 or -1 one, so any result will take precedence over initial values
   Storage[2]=-10.;  // This is also lower then any value produced by an eligible atom, which are all positive
 
-  
+
   helper.CopyVector(&(Line.endpoints[0]->node->x));
   for (int i =0; i<3; i++)
@@ -1114 +1122 @@
     }
 
-  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &direction1, T.NormalVector, Storage, RADIUS, *mol);
-  
+  Chord.CopyVector(&(Line.endpoints[0]->node->x)); // bring into calling function
+  Chord.SubtractVector(&(Line.endpoints[1]->node->x));
+
+printf("Looking for third point candidates for triangle \n");
+  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, T.NormalVector, Storage, RADIUS, mol);
+
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
   // Next Triangle is Line, atom with number in Storage[0]
@@ -1141 +1153 @@
   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
     {
-      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
+
+
+  if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
         {
-          LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
-          LinesOnBoundaryCount++;
+                LinesOnBoundaryCount++;
         }
     }
@@ -1154 +1167 @@
       BTS->NormalVector->Scale(-1);
     }
-  
-};
-
-
-void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
-{
+
+};
+
+
+void Find_second_point_for_Tesselation(atom* a, atom* Candidate, int n, Vector Oben, double Storage[3], molecule* mol)
+{
+        printf("Looking for second point of starting triangle \n");
   int i;
   Vector *AngleCheck;
   atom* Walker;
 
-  AngleCheck->CopyVector(&Candidate.x);
-  AngleCheck->SubtractVector(&a.x);
-  if (AngleCheck->ScalarProduct(&Oben) < Storage[1])
-    {
-      Storage[0]=(double)(Candidate.nr);
-      Storage[1]=AngleCheck->ScalarProduct(&Oben);
+  AngleCheck->CopyVector(&(Candidate->x));
+  AngleCheck->SubtractVector(&(a->x));
+  if (AngleCheck->Angle(&Oben) < Storage[1])
+    {
+          //printf("Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
+      Storage[0]=(double)(Candidate->nr);
+      Storage[1]=AngleCheck->Angle(&Oben);
+      //printf("Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     };
-
+  printf("%d \n", n);
   if (n<5)
     {
-      for (i = 0; i< mol.NumberOfBondsPerAtom[Candidate.nr]; i++)
+      for (i = 0; i< mol->NumberOfBondsPerAtom[Candidate->nr]; i++)
         {
-          Walker = mol.start;
-          while (Candidate.nr != (mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr ==Candidate.nr ? mol.ListOfBondsPerAtom[Candidate.nr][i]->leftatom->nr : mol.ListOfBondsPerAtom[Candidate.nr][i]->rightatom->nr))
+          Walker = mol->start;
+          while (Candidate->nr != (mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr ==Candidate->nr ? mol->ListOfBondsPerAtom[Candidate->nr][i]->leftatom->nr : mol->ListOfBondsPerAtom[Candidate->nr][i]->rightatom->nr))
             {
               Walker = Walker->next;
             };
-          
-          Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
+
+          Find_second_point_for_Tesselation(a, Walker, n+1, Oben, Storage, mol);
             };
     };
-  
-
-};
-
-
-void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
-{
+
+
+};
+
+
+void Tesselation::Find_starting_triangle(molecule* mol, const double RADIUS)
+{
+        printf("Looking for starting triangle \n");
   int i=0;
-  atom Walker;
-  atom Walker2;
-  atom Walker3;
+  atom* Walker;
+  atom* Walker2;
+  atom* Walker3;
   int max_index[3];
   double max_coordinate[3];
   Vector Oben;
   Vector helper;
+  Vector Chord;
 
   Oben.Zero();
@@ -1209 +1227 @@
       max_coordinate[i] =-1;
     }
-
-  Walker = *mol.start;
-  while (Walker.next != NULL)
-    {
+printf("Molecule mol is there and has %d Atoms \n", mol->AtomCount);
+  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])
+      {
+          if (Walker->x.x[i] > max_coordinate[i])
             {
-              max_coordinate[i]=Walker.x.x[i];
-              max_index[i]=Walker.nr;
+              max_coordinate[i]=Walker->x.x[i];
+              max_index[i]=Walker->nr;
             }
-        }
-    }
-
+      }
+    }
+
+printf("Found starting atom \n");
   //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
-  const int k=0;
-
-  Oben.x[k]=1;
-  Walker = *mol.start;
-  while (Walker.nr != max_index[k])
-    {
-      Walker = *Walker.next;
-    }
-
+  const int k=2;
+
+  Oben.x[k]=1.;
+  Walker = mol->start;
+  Walker = Walker->next;
+  while (Walker->nr != max_index[k])
+    {
+      Walker = Walker->next;
+    }
+printf("%d \n", Walker->nr);
   double Storage[3];
   Storage[0]=-1.;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2.;   // 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[2]=-10.;  // This will be an angle looking for the third point.
-
-
-  for (i=0; i< mol.NumberOfBondsPerAtom[Walker.nr]; i++)
-    {
-      Walker2 = *mol.start;
-      while (Walker2.nr != (mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom->nr : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr) ) 
+  Storage[1]=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[2]=999999.;  // This will be an angle looking for the third point.
+printf("%d \n", mol->NumberOfBondsPerAtom[Walker->nr]);
+
+  for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
+    {
+      Walker2 = mol->start;
+      Walker2 = Walker2->next;
+      while (Walker2->nr != (mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][i]->rightatom->nr : mol->ListOfBondsPerAtom[Walker->nr][i]->leftatom->nr) )
         {
-          Walker2 = *Walker2.next;
+          Walker2 = Walker2->next;
         }
 
@@ -1250 +1272 @@
     }
 
-  Walker2 = *mol.start;
-
-  while (Walker2.nr != int(Storage[0]))
-    {
-      Walker = *Walker.next;
-    }
-  
-  helper.CopyVector(&Walker.x);
-  helper.SubtractVector(&Walker2.x);
+  Walker2 = mol->start;
+
+  while (Walker2->nr != int(Storage[0]))
+    {
+      Walker2 = Walker2->next;
+    }
+
+  helper.CopyVector(&(Walker->x));
+  helper.SubtractVector(&(Walker2->x));
   Oben.ProjectOntoPlane(&helper);
   helper.VectorProduct(&Oben);
-
-       Find_next_suitable_point(Walker, Walker2, *(mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom->nr == Walker.nr ? mol.ListOfBondsPerAtom[Walker.nr][i]->rightatom : mol.ListOfBondsPerAtom[Walker.nr][i]->leftatom), 0, &helper, &Oben, Storage, RADIUS, mol);
-  Walker3 = *mol.start;
-  while (Walker3.nr != int(Storage[0]))
-    {
-      Walker3 = *Walker3.next;
+  Storage[0]=-1.;   // Id must be positive, we see should nothing be done
+  Storage[1]=-2.;   // 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[2]= -10.;  // This will be an angle looking for the third point.
+
+  Chord.CopyVector(&(Walker->x)); // bring into calling function
+  Chord.SubtractVector(&(Walker2->x));
+
+printf("Looking for third point candidates \n");
+       Find_next_suitable_point(Walker, Walker2, (mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom->nr == Walker->nr ? mol->ListOfBondsPerAtom[Walker->nr][0]->rightatom : mol->ListOfBondsPerAtom[Walker->nr][0]->leftatom), 0, &Chord, &helper, &Oben, Storage, RADIUS, mol);
+  Walker3 = mol->start;
+  while (Walker3->nr != int(Storage[0]))
+    {
+      Walker3 = Walker3->next;
     }
 
   //Starting Triangle is Walker, Walker2, index Storage[0]
 
-  AddPoint(&Walker);
-  AddPoint(&Walker2);
-  AddPoint(&Walker3);
-
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
+  AddPoint(Walker);
+  AddPoint(Walker2);
+  AddPoint(Walker3);
+
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker2);
   BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker3);
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker3);
   BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
-  BPS[0] = new class BoundaryPointSet(&Walker);
-  BPS[1] = new class BoundaryPointSet(&Walker2);
-  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Walker2);
+  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
 
   BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
@@ -1289 +1318 @@
   TrianglesOnBoundaryCount++;
 
-  for(int i=0;i<NDIM;i++) 
+  for(int i=0;i<NDIM;i++)
     {
       LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
@@ -1304 +1333 @@
 
 
-void Find_non_convex_border(Tesselation* Tess, molecule mol) 
-{
+void Find_non_convex_border(Tesselation Tess, molecule* mol)
+{
+        printf("Entering finding of non convex hull. \n");
   const double RADIUS =6;
-  Tess->Find_starting_triangle(mol, RADIUS);
-
-  for (LineMap::iterator baseline = Tess->LinesOnBoundary.begin(); baseline != Tess->LinesOnBoundary.end(); baseline++) 
-    if (baseline->second->TrianglesCount == 1) 
+  Tess.Find_starting_triangle(mol, RADIUS);
+
+  for (LineMap::iterator baseline = Tess.LinesOnBoundary.begin(); baseline != Tess.LinesOnBoundary.end(); baseline++)
+    if (baseline->second->TrianglesCount == 1)
       {
-        Tess->Find_next_suitable_triangle(&mol, baseline->second, baseline->second->triangles.begin()->second, RADIUS); //the line is there, so there is a triangle, but only one.
-
+                Tess.Find_next_suitable_triangle(mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
       }
     else