Index: molecuilder/src/border.cpp
===================================================================
--- molecuilder/src/border.cpp	(revision a709c4e760b7e3ec20d189f010224f7e2d9b5096)
+++ molecuilder/src/border.cpp	(revision 6b3826d15e0c823f1e23d4717d70b7955629fe84)
@@ -2,14 +2,10 @@
 #include "boundary.hpp"
 
-inline int round(double x)
-{
-  //brauche ich das hier? Kann ich mit int operieren?
-  return int(x > 0.0 ? x + 0.5 : x - 0.5);
-}
-
-
-
-
-void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, double *Storage, const double RADIUS)
+
+
+
+
+
+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,
@@ -56,18 +52,18 @@
   if (n<5)
     {
-      for(i=0; i<molecule.NumberOfBondsPerAtom[Candidate.nr];i++)
-	{
-	  while (Candidate.nr != molecule->ListOfBonds[Candidate.nr][i])
+      for(i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+	{
+	  while (Candidate.nr != mol->ListOfBonds[Candidate.nr][i])
 	    {
 	      Walker = Walker.next;
 	    }
 
-	  Find_next_suitable_point(a, b, Walker, n+1, d1, d2, RADIUS);
-	}
-    }
-}
-
-
-void Find_next_suitable_triangle(Triangle T, BoundaryLine Line)
+	  Find_next_suitable_point(a, b, Walker, n+1, d1, d2, RADIUS, mol);
+	}
+    }
+}
+
+
+void Find_next_suitable_triangle(molecule mol, BoundaryLineSet Line, BoundaryTriangleSet T)
 {
   Vector CenterOfLine = Line->endpoints.node[0].x;
@@ -78,6 +74,6 @@
   double *Storage[3];
   Storage[0]=-1;   // Id must be positive, we see should nothing be done
-  Storage[1]=-2;   // 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, though due to Storage[1] this is of no concern
+  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
 
   
@@ -106,12 +102,49 @@
   // Konstruiere nun neues Dreieck am Ende der Liste der Dreiecke
   // Next Triangle is Line, atom with number in Storage[0]
-  
-}
-
-
-void Find_starting_triangle()
+
+  Walker=mol->start;
+  while (Walker.nr != (int)Storage[0])
+    {
+      Walker = Walker->next;
+    }
+
+  AddPoint(Walker);
+
+  BPS[0] = BoundaryPointSet(Walker);
+  BPS[1] = Line->enpoints.node[0];
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = Walker;
+  BPS[1] = Line->endpoints.node[1];
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BLS[2] = line;
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
+    {
+      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary->end)
+	{
+	  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+	  LinesOnBoundaryCount++;
+	}
+    }
+  GetNormalVector(BTS.NormalVector);
+
+  if( (BTS.NormalVector.ScalarProduct(T.NormalVecotr)<0 && Storage[1]>0) || \
+      (BTS.NormalVector.ScalarProduct(T.NormalVecotr)>0 && Storage[1]<0))
+    {
+      BTS.NormalVector.Scale(-1);
+    }
+  
+}
+
+
+void Find_starting_triangle(molecule mol)
 {
   atom Walker;
   atom Walker2;
+  atom Walker3;
   int max_index[3];
   double max_coordinate[3];
@@ -128,5 +161,5 @@
     }
 
-  Walker = molecule->start;
+  Walker = mol->start;
   while (Walker->next != NULL)
     {
@@ -145,5 +178,5 @@
 
   Oben.x[k]=1;
-  Walker = molecule->start;
+  Walker = mol->start;
   while (Walker.nr != max_index[k])
     {
@@ -157,8 +190,8 @@
 
 
-  for (i=0; i< molecule->NumberOfBondsPerAtoms[Walker.nr]; i++)
-    {
-      Walker2 = molecule->start;
-      while (Walker2.nr != molecule->ListOfBondsPerAtoms[Walker.nr][i])  // Stimmt die Ueberpruefung $$$
+  for (i=0; i< mol->NumberOfBondsPerAtoms[Walker.nr]; i++)
+    {
+      Walker2 = mol->start;
+      while (Walker2.nr != mol->ListOfBondsPerAtoms[Walker.nr][i])  // Stimmt die Ueberpruefung $$$
 	{
 	  Walker2 =Walker2->next;
@@ -168,5 +201,5 @@
     }
 
-  Walker2=molecule->start;
+  Walker2=mol->start;
 
   while (Walker2.nr != int(Storage[0]))
@@ -181,12 +214,53 @@
 
   Find_next_suitable_point(Walker, Walker2, Candidate, 0, helper, Oben, Storage, Radius);
+  Walker3 = mol->start;
+  while (Walker3.nr != int(Storage[0]))
+    {
+      Walker3 = Walker3->next;
+    }
 
   //Starting Triangle is Walker, Walker2, index Storage[0]
 
-}
-
-
-void Find_non_convex_border()
-{
-
-}
+  AddPoint(Walker);
+  AddPoint(Walker2);
+  AddPoint(Walker3);
+
+  BPS[0] = BoundaryPointSet(Walker);
+  BPS[1] = BoundaryPointSet(Walker2);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = Walker;
+  BPS[1] = Walker3;
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = Walker;
+  BPS[1] = Walker2;
+  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);  
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) 
+    {
+      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+      LinesOnBoundaryCount++;
+    }
+
+}
+
+
+void Find_non_convex_border(Tesselation Tess, molecule mol) // this needs input of type molecule
+{
+
+  Tess.Find_starting_triangle(mol);
+
+  for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++) 
+    if (baseline->second->TrianglesCount == 1) 
+      {
+	Find_next_suitable_triangle(mol, baseline->second, baseline->second->triangles.begin()->second); //the line is there, so there is a triangle, but only one.
+
+      }
+    else
+      {
+	printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+      }
+}
Index: molecuilder/src/boundary.cpp
===================================================================
--- molecuilder/src/boundary.cpp	(revision a709c4e760b7e3ec20d189f010224f7e2d9b5096)
+++ molecuilder/src/boundary.cpp	(revision 6b3826d15e0c823f1e23d4717d70b7955629fe84)
@@ -1,4 +1,10 @@
 #include "molecules.hpp"
 #include "boundary.hpp"
+
+
+// =================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);
+
 
 // ======================================== Points on Boundary =================================
@@ -471,8 +477,10 @@
   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);
+  /*
   // 1. calculate center of gravity
   *out << endl;
@@ -518,4 +526,6 @@
 
   *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
+  */
+
 
   // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
@@ -539,5 +549,5 @@
   *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 
-
+  /*
   // 7. translate all points back from CoG
   *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
@@ -548,5 +558,10 @@
     Walker->x.Translate(CenterOfGravity);
   }
-  
+  */
+
+
+
+
+
   // 8. Store triangles in tecplot file
   if (tecplot != NULL) {
@@ -994,2 +1009,313 @@
     PointsOnBoundaryCount++;
 };
+
+
+
+
+
+
+//====================================================================================================================
+//
+// ab hier habe ich das verzapft.
+//
+//====================================================================================================================
+
+
+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.
+   */
+  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);
+  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
+  {
+
+    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);
+      }
+    else
+      {
+	if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]>0 &&  Storage[2]< AngleCheck.Angle(d2)) or \
+	    (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[1] && Storage[1]<0 &&  Storage[2]> AngleCheck.Angle(d2)))
+	  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
+	  {
+	    Storage[0]=(double)Candidate.nr;
+	    Storage[1]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
+	    Storage[2]=AngleCheck.Angle(d2);
+	  }
+      }
+  }
+
+  if (n<5)
+    {
+      for(int i=0; i<mol.NumberOfBondsPerAtom[Candidate.nr];i++)
+	{
+	  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_next_suitable_point(a, b, *Walker, n+1, d1, d2, Storage, RADIUS, mol);
+	}
+    }
+};
+
+
+void Tesselation::Find_next_suitable_triangle(molecule *mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS)
+{
+  Vector CenterOfLine = Line.endpoints[0]->node->x;
+  Vector direction1;
+  Vector direction2;
+  Vector helper;
+  atom* Walker;
+
+  double Storage[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++)
+    {
+      if (T.endpoints[i]->node->nr != Line.endpoints[0]->node->nr && T.endpoints[i]->node->nr!=Line.endpoints[1]->node->nr)
+	{
+	  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);
+    }
+
+  Find_next_suitable_point(*Line.endpoints[0]->node, *Line.endpoints[1]->node, *Line.endpoints[0]->node, 0, &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]
+
+  Walker= mol->start;
+  while (Walker->nr != (int)Storage[0])
+    {
+      Walker = Walker->next;
+    }
+
+  AddPoint(Walker);
+
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[0]->node);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BPS[0] = new class BoundaryPointSet(Walker);
+  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
+  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  BLS[2] = new class BoundaryLineSet(Line);
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
+    {
+      if (LinesOnBoundary.find(BTS->lines[i]) == LinesOnBoundary.end)
+	{
+	  LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+	  LinesOnBoundaryCount++;
+	}
+    }
+  BTS->GetNormalVector(*BTS->NormalVector);
+
+  if( (BTS->NormalVector->ScalarProduct(T.NormalVector)<0 && Storage[1]>0) || \
+      (BTS->NormalVector->ScalarProduct(T.NormalVector)>0 && Storage[1]<0))
+    {
+      BTS->NormalVector->Scale(-1);
+    }
+  
+};
+
+
+void Find_second_point_for_Tesselation(atom a, atom Candidate, int n, Vector Oben, double* Storage, molecule mol)
+{
+  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);
+    };
+
+  if (n<5)
+    {
+      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 = Walker->next;
+	    };
+	  
+	  Find_second_point_for_Tesselation(a, *Walker, n+1, Oben, Storage, mol);
+	    };
+    };
+  
+
+};
+
+
+void Tesselation::Find_starting_triangle(molecule mol, const double RADIUS)
+{
+  int i=0;
+  atom Walker;
+  atom Walker2;
+  atom Walker3;
+  int max_index[3];
+  double max_coordinate[3];
+  Vector Oben;
+  Vector helper;
+
+  Oben.Zero();
+
+
+  for(i =0; i<3; i++)
+    {
+      max_index[i] =-1;
+      max_coordinate[i] =-1;
+    }
+
+  Walker = *mol.start;
+  while (Walker.next != NULL)
+    {
+      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.nr;
+	    }
+	}
+    }
+
+  //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;
+    }
+
+  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) ) 
+	{
+	  Walker2 = *Walker2.next;
+	}
+
+      Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Storage, mol);
+    }
+
+  Walker2 = *mol.start;
+
+  while (Walker2.nr != int(Storage[0]))
+    {
+      Walker = *Walker.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;
+    }
+
+  //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);
+  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
+  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);  
+
+  BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+  TrianglesOnBoundary.insert( TrianglePair(TrianglesOnBoundaryCount, BTS) );
+  TrianglesOnBoundaryCount++;
+
+  for(int i=0;i<NDIM;i++) 
+    {
+      LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]) );
+      LinesOnBoundaryCount++;
+    };
+
+       BTS->GetNormalVector(*BTS->NormalVector);
+
+       if( BTS->NormalVector->ScalarProduct(&Oben)<0)
+	 {
+	   BTS->NormalVector->Scale(-1);
+	 }
+};
+
+
+void Find_non_convex_border(Tesselation* Tess, molecule mol) 
+{
+  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_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
+      {
+	printf("There is a line with %d triangles adjacent", baseline->second->TrianglesCount);
+      }
+};
Index: molecuilder/src/boundary.hpp
===================================================================
--- molecuilder/src/boundary.hpp	(revision a709c4e760b7e3ec20d189f010224f7e2d9b5096)
+++ molecuilder/src/boundary.hpp	(revision 6b3826d15e0c823f1e23d4717d70b7955629fe84)
@@ -69,4 +69,5 @@
     class BoundaryPointSet *endpoints[3];
     class BoundaryLineSet *lines[3];
+    Vector *NormalVector;
     int Nr;
 };
@@ -81,4 +82,6 @@
     void GuessStartingTriangle(ofstream *out);
     void AddPoint(atom * Walker);
+    void Find_starting_triangle(molecule mol, const double RADIUS);
+    void Find_next_suitable_triangle(molecule* mol, BoundaryLineSet Line, BoundaryTriangleSet T, const double& RADIUS);
     
     PointMap PointsOnBoundary;
Index: molecuilder/src/vector.cpp
===================================================================
--- molecuilder/src/vector.cpp	(revision a709c4e760b7e3ec20d189f010224f7e2d9b5096)
+++ molecuilder/src/vector.cpp	(revision 6b3826d15e0c823f1e23d4717d70b7955629fe84)
@@ -115,4 +115,22 @@
 };
 
+
+/** Calculates VectorProduct between this and another vector.
+ *  -# returns the Product in place of vector from which it was initiated
+ *  -# ATTENTION: Only three dim.
+ *  \param *y array to vector with which to calculate crossproduct
+ *  \return \f$ x \times y \f&
+ */ 
+void Vector::VectorProduct(const Vector *y)
+{
+  Vector tmp;
+  tmp[0] = x[1]*y->x[2] - x[2]*y->x[1];
+  tmp[1] = x[2]*y->x[0] - x[0]*y->x[2];
+  tmp[2] = x[0]*y->x[1] - x[1]*Y->x[0];
+  this->CopyVector(&tmp);
+
+};
+
+
 /** projects this vector onto plane defined by \a *y.
  * \param *y array to normal vector of plane
@@ -184,5 +202,5 @@
 };
 
-/** Calculates the angle between this and another vector.
+/** Calculates the angle between this and another vector. 
  * \param *y array to second vector
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
Index: molecuilder/src/vector.hpp
===================================================================
--- molecuilder/src/vector.hpp	(revision a709c4e760b7e3ec20d189f010224f7e2d9b5096)
+++ molecuilder/src/vector.hpp	(revision 6b3826d15e0c823f1e23d4717d70b7955629fe84)
@@ -26,4 +26,5 @@
   void CopyVector(const Vector *y);
   void RotateVector(const Vector *y, const double alpha);
+  void VectorProduct(const Vector *y);
   void ProjectOntoPlane(const Vector *y);
   void Zero(); 
