Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/Makefile.am	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -1,4 +1,7 @@
 ATOMSOURCE = atom.cpp atom_atominfo.cpp atom_bondedparticle.cpp atom_bondedparticleinfo.cpp atom_graphnode.cpp atom_graphnodeinfo.cpp atom_particleinfo.cpp atom_trajectoryparticle.cpp atom_trajectoryparticleinfo.cpp
 ATOMHEADER = atom.hpp atom_atominfo.hpp atom_bondedparticle.hpp atom_bondedparticleinfo.hpp atom_graphnode.hpp atom_graphnodeinfo.hpp atom_particleinfo.hpp atom_trajectoryparticle.hpp atom_trajectoryparticleinfo.hpp
+
+LINALGSOURCE = gslmatrix.cpp gslvector.cpp linearsystemofequations.cpp
+LINALGHEADER = gslmatrix.hpp gslvector.hpp linearsystemofequations.hpp
 
 ANALYSISSOURCE = analysis_bonds.cpp analysis_correlation.cpp
@@ -11,12 +14,13 @@
 INCLUDES = -I$(top_srcdir)/src/unittests
 
-noinst_LIBRARIES = libmolecuilder.a
+noinst_LIBRARIES = libmolecuilder.a libgslwrapper.a
 bin_PROGRAMS = molecuilder joiner analyzer
 molecuilderdir = ${bindir}
 libmolecuilder_a_SOURCES = ${SOURCE} ${HEADER}
+libgslwrapper_a_SOURCES = ${LINALGSOURCE} ${LINALGHEADER}
 molecuilder_DATA = elements.db valence.db orbitals.db Hbonddistance.db Hbondangle.db
 molecuilder_LDFLAGS = $(BOOST_LIB)
 molecuilder_SOURCES = builder.cpp
-molecuilder_LDADD = libmolecuilder.a
+molecuilder_LDADD = libmolecuilder.a libgslwrapper.a
 joiner_SOURCES = joiner.cpp datacreator.cpp parser.cpp datacreator.hpp helpers.hpp parser.hpp periodentafel.hpp
 joiner_LDADD = libmolecuilder.a 
@@ -28,5 +32,5 @@
 FORCE:
 $(srcdir)/.git-version: FORCE
-	@if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe --tags HEAD) > .git-version-t 2>/dev/null \
+	@if (test -d $(top_srcdir)/../.git && cd $(srcdir) && git describe HEAD) > .git-version-t 2>/dev/null \
 	  && ! diff .git-version-t $(srcdir)/.git-version >/dev/null 2>&1; then \
 	  mv -f .git-version-t $(srcdir)/.git-version; \
Index: src/analysis_correlation.cpp
===================================================================
--- src/analysis_correlation.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/analysis_correlation.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -10,4 +10,5 @@
 #include "analysis_correlation.hpp"
 #include "element.hpp"
+#include "info.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
@@ -28,4 +29,5 @@
 PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
 {
+  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
@@ -77,4 +79,5 @@
 PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
 {
+  Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
@@ -154,4 +157,5 @@
 CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
 {
+  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
@@ -190,4 +194,5 @@
 CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
 {
+  Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
@@ -243,4 +248,5 @@
 CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
 {
+  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
@@ -261,5 +267,5 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          triangle = Surface->FindClosestTriangleToPoint(Walker->node, LC );
+          triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
           if (triangle != NULL) {
             distance = DistanceToTrianglePlane(Walker->node, triangle);
@@ -288,4 +294,5 @@
 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
 {
+  Info FunctionInfo(__func__);
   CorrelationToSurfaceMap *outmap = NULL;
   double distance = 0;
@@ -320,5 +327,5 @@
                 checkX.AddVector(&periodicX);
                 checkX.MatrixMultiplication(FullMatrix);
-                triangle = Surface->FindClosestTriangleToPoint(&checkX, LC );
+                triangle = Surface->FindClosestTriangleToVector(&checkX, LC );
                 if (triangle != NULL) {
                   distance = DistanceToTrianglePlane(&checkX, triangle);
@@ -342,4 +349,5 @@
 double GetBin ( const double value, const double BinWidth, const double BinStart )
 {
+  Info FunctionInfo(__func__);
   double bin =(double) (floor((value - BinStart)/BinWidth));
   return (bin*BinWidth+BinStart);
@@ -353,4 +361,5 @@
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map )
 {
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tCount" << endl;
   for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
@@ -365,4 +374,5 @@
 void OutputPairCorrelation( ofstream * const file, const PairCorrelationMap * const map )
 {
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom1\tAtom2" << endl;
   for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
@@ -377,4 +387,5 @@
 void OutputCorrelationToPoint( ofstream * const file, const CorrelationToPointMap * const map )
 {
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tAtom::x[i]-point.x[i]" << endl;
   for (CorrelationToPointMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
@@ -392,4 +403,5 @@
 void OutputCorrelationToSurface( ofstream * const file, const CorrelationToSurfaceMap * const map )
 {
+  Info FunctionInfo(__func__);
   *file << "# BinStart\tTriangle" << endl;
   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
Index: src/analyzer.cpp
===================================================================
--- src/analyzer.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/analyzer.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -7,4 +7,6 @@
 
 //============================ INCLUDES ===========================
+
+#include <cstring>
 
 #include "datacreator.hpp"
Index: src/bondgraph.cpp
===================================================================
--- src/bondgraph.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/bondgraph.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -35,5 +35,5 @@
 /** Parses the bond lengths in a given file and puts them int a matrix form.
  * Allocates \a MatrixContainer for BondGraph::BondLengthMatrix, using MatrixContainer::ParseMatrix(),
- * but only if parsing is successfull. Otherwise variable is left as NULL.
+ * but only if parsing is successful. Otherwise variable is left as NULL.
  * \param *out output stream for debugging
  * \param filename file with bond lengths to parse
Index: src/boundary.cpp
===================================================================
--- src/boundary.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/boundary.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -654,9 +654,9 @@
  * \param *out output stream for debugging
  * \param *mol molecule with atoms and bonds
- * \param *&TesselStruct Tesselation with boundary triangles
+ * \param *TesselStruct Tesselation with boundary triangles
  * \param *filename prefix of filename
  * \param *extraSuffix intermediate suffix
  */
-void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix)
+void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix)
 {
 	Info FunctionInfo(__func__);
@@ -789,4 +789,6 @@
  * \param configuration contains box dimensions
  * \param distance[NDIM] distance between filling molecules in each direction
+ * \param boundary length of boundary zone between molecule and filling mollecules
+ * \param epsilon distance to surface which is not filled
  * \param RandAtomDisplacement maximum distance for random displacement per atom
  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
@@ -794,5 +796,5 @@
  * \return *mol pointer to new molecule with filled atoms
  */
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation)
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
 {
 	Info FunctionInfo(__func__);
@@ -817,9 +819,8 @@
   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
-    LCList[i] = new LinkedCell((*ListRunner), 5.); // get linked cell list
-    if (TesselStruct[i] == NULL) {
-      Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-      FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
-    }
+    LCList[i] = new LinkedCell((*ListRunner), 10.); // get linked cell list
+    Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
+    TesselStruct[i] = NULL;
+    FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
     i++;
   }
@@ -835,13 +836,7 @@
   FillerDistance.Init(distance[0], distance[1], distance[2]);
   FillerDistance.InverseMatrixMultiplication(M);
-  Log() << Verbose(1) << "INFO: Grid steps are ";
-  for(int i=0;i<NDIM;i++) {
+  for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance.x[i]);
-    Log() << Verbose(1) << N[i];
-    if (i != NDIM-1)
-      Log() << Verbose(1)<< ", ";
-    else
-      Log() << Verbose(1) << "." << endl;
-  }
+  Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
 
   // go over [0,1]^3 filler grid
@@ -859,8 +854,15 @@
           // get linked cell list
           if (TesselStruct[i] == NULL) {
-            eLog() << Verbose(1) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
+            eLog() << Verbose(0) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
             FillIt = false;
           } else {
-            FillIt = FillIt && (!TesselStruct[i]->IsInnerPoint(CurrentPosition, LCList[i]));
+            const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
+            FillIt = FillIt && (distance > boundary*boundary);
+            if (FillIt) {
+              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
+            } else {
+              Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
+              break;
+            }
             i++;
           }
@@ -931,4 +933,9 @@
       }
   Free(&M);
+
+  // output to file
+  TesselStruct[0]->LastTriangle = NULL;
+  StoreTrianglesinFile(Filling, TesselStruct[0], "Tesselated", ".dat");
+
   for (size_t i=0;i<List->ListOfMolecules.size();i++) {
   	delete(LCList[i]);
@@ -1041,4 +1048,10 @@
 //  TesselStruct->RemoveDegeneratedTriangles();
 
+  // check envelope for consistency
+  status = CheckListOfBaselines(TesselStruct);
+
+  // store before correction
+  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
+
   // correct degenerated polygons
   TesselStruct->CorrectAllDegeneratedPolygons();
Index: src/boundary.hpp
===================================================================
--- src/boundary.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/boundary.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -36,5 +36,5 @@
 #define DEBUG 1
 #define DoSingleStepOutput 0
-#define SingleStepWidth 1
+#define SingleStepWidth 10
 
 #define DistancePair pair < double, atom* >
@@ -49,5 +49,5 @@
 
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandAtomDisplacement, double RandMolDisplacement, bool DoRandomRotation);
+molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
 void FindConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
@@ -58,5 +58,5 @@
 void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity);
 bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
-void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix);
+void StoreTrianglesinFile(const molecule * const mol, const Tesselation * const TesselStruct, const char *filename, const char *extraSuffix);
 double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration);
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/builder.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -49,4 +49,6 @@
 
 using namespace std;
+
+#include <cstring>
 
 #include "analysis_correlation.hpp"
@@ -1410,12 +1412,14 @@
             Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
             Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
-            Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
+            Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
             Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
             Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
             Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
-            Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
+            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
+            Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
             Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
             Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
+            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
             Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
             Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
@@ -1575,16 +1579,4 @@
                 else {
                   Log() << Verbose(2) << "File found and parsed." << endl;
-                  // @TODO rather do the dissection afterwards
-//                  mol->SetNameFromFilename(argv[argptr]);
-//                  molecules->ListOfMolecules.remove(mol);
-//                  molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration);
-//                  delete(mol);
-//                  if (molecules->ListOfMolecules.size() != 0) {
-//                    for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
-//                      if ((*ListRunner)->ActiveFlag) {
-//                        mol = *ListRunner;
-//                        break;
-//                      }
-//                  }
                   configPresent = present;
                 }
@@ -1665,4 +1657,21 @@
               //argptr+=1;
               break;
+            case 'I':
+              Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
+              // @TODO rather do the dissection afterwards
+              molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration);
+              mol = NULL;
+              if (molecules->ListOfMolecules.size() != 0) {
+                for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+                  if ((*ListRunner)->ActiveFlag) {
+                    mol = *ListRunner;
+                    break;
+                  }
+              }
+              if (mol == NULL) {
+                mol = *(molecules->ListOfMolecules.begin());
+                mol->ActiveFlag = true;
+              }
+              break;
             case 'C':
               if (ExitFlag == 0) ExitFlag = 1;
@@ -1672,5 +1681,4 @@
                 performCriticalExit();
               } else {
-                SaveFlag = false;
                 ofstream output(argv[argptr+1]);
                 ofstream binoutput(argv[argptr+2]);
@@ -1692,5 +1700,5 @@
                 counter = 0;
                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
-                  Actives[counter] = (*BigFinder)->ActiveFlag;
+                  Actives[counter++] = (*BigFinder)->ActiveFlag;
                   (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
                 }
@@ -1705,5 +1713,5 @@
                 binoutput.close();
                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
-                  (*BigFinder)->ActiveFlag = Actives[counter];
+                  (*BigFinder)->ActiveFlag = Actives[counter++];
                 Free(&Actives);
                 delete(LCList);
@@ -1728,7 +1736,7 @@
             case 'F':
               if (ExitFlag == 0) ExitFlag = 1;
-              if (argptr+5 >=argc) {
+              if (argptr+6 >=argc) {
                 ExitFlag = 255;
-                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
+                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
                 performCriticalExit();
               } else {
@@ -1738,25 +1746,30 @@
                 molecule *filler = new molecule(periode);;
                 molecule *Filling = NULL;
-                atom *second = NULL, *third = NULL;
+//                atom *second = NULL, *third = NULL;
                 first = new atom();
-                first->type = periode->FindElement(1);
-                first->x.Init(0.441, -0.143, 0.);
+                first->type = periode->FindElement(5);
+                first->x.Zero();
                 filler->AddAtom(first);
-                second = new atom();
-                second->type = periode->FindElement(1);
-                second->x.Init(-0.464, 1.137, 0.0);
-                filler->AddAtom(second);
-                third = new atom();
-                third->type = periode->FindElement(8);
-                third->x.Init(-0.464, 0.177, 0.);
-                filler->AddAtom(third);
-                filler->AddBond(first, third, 1);
-                filler->AddBond(second, third, 1);
+//                first = new atom();
+//                first->type = periode->FindElement(1);
+//                first->x.Init(0.441, -0.143, 0.);
+//                filler->AddAtom(first);
+//                second = new atom();
+//                second->type = periode->FindElement(1);
+//                second->x.Init(-0.464, 1.137, 0.0);
+//                filler->AddAtom(second);
+//                third = new atom();
+//                third->type = periode->FindElement(8);
+//                third->x.Init(-0.464, 0.177, 0.);
+//                filler->AddAtom(third);
+//                filler->AddBond(first, third, 1);
+//                filler->AddBond(second, third, 1);
                 // call routine
                 double distance[NDIM];
                 for (int i=0;i<NDIM;i++)
                   distance[i] = atof(argv[argptr+i]);
-                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
+                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
                 if (Filling != NULL) {
+                  Filling->ActiveFlag = false;
                   molecules->insert(Filling);
                 }
@@ -2107,5 +2120,5 @@
                 if (volume != -1)
                   ExitFlag = 255;
-                  eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
+                  eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
                   performCriticalExit();
               } else {
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/config.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -6,4 +6,5 @@
 
 #include <stdio.h>
+#include <cstring>
 
 #include "atom.hpp"
@@ -27,5 +28,5 @@
     char number1[8];
     char number2[8];
-    char *dummy1, *dummy2;
+    const char *dummy1, *dummy2;
     //Log() << Verbose(0) << s1 << "  " << s2 << endl;
     dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
@@ -1074,6 +1075,6 @@
   // don't do this here ...
   //MolList->DissectMoleculeIntoConnectedSubgraphs(mol,this);
-
-  delete(mol);
+  //delete(mol);
+
   delete(FileBuffer);
 };
@@ -2123,5 +2124,5 @@
               }
               line++;
-            } while (dummy1 != NULL && (dummy1[0] == '#') || (dummy1[0] == '\n'));
+            } while ((dummy1 != NULL) && ((dummy1[0] == '#') || (dummy1[0] == '\n')));
             dummy = dummy1;
           } else { // simple int, strings or doubles start in the same line
Index: src/gslmatrix.cpp
===================================================================
--- src/gslmatrix.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/gslmatrix.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,245 @@
+/*
+ * gslmatrix.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include "gslmatrix.hpp"
+#include "helpers.hpp"
+
+#include <cassert>
+#include <gsl/gsl_linalg.h>
+
+/** Constructor of class GSLMatrix.
+ * Allocates GSL structures
+ * \param m dimension of matrix
+ */
+GSLMatrix::GSLMatrix(size_t m, size_t n) : rows(m), columns(n)
+{
+  matrix = gsl_matrix_calloc(rows, columns);
+};
+
+/** Copy constructor of class GSLMatrix.
+ * Allocates GSL structures and copies components from \a *src.
+ * \param *src source matrix
+ */
+GSLMatrix::GSLMatrix(const GSLMatrix * const src) : rows(src->rows), columns(src->columns)
+{
+  matrix = gsl_matrix_alloc(rows, columns);
+  gsl_matrix_memcpy (matrix, src->matrix);
+};
+
+/** Destructor of class GSLMatrix.
+ * Frees GSL structures
+ */
+GSLMatrix::~GSLMatrix()
+{
+  gsl_matrix_free(matrix);
+  rows = 0;
+  columns = 0;
+};
+
+/** Assignment operator.
+ * \param &rhs right hand side
+ * \return object itself
+ */
+GSLMatrix& GSLMatrix::operator=(const GSLMatrix& rhs)
+{
+  if (this == &rhs)   // not necessary here, but identity assignment check is generally a good idea
+    return *this;
+  assert(rows == rhs.rows && columns == rhs.columns && "Number of rows and columns do not match!");
+
+  gsl_matrix_memcpy (matrix, rhs.matrix);
+
+  return *this;
+};
+
+/* ============================ Accessing =============================== */
+/** This function sets the matrix from a double array.
+ * Creates a matrix view of the array and performs a memcopy.
+ * \param *x array of values (no dimension check is performed)
+ */
+void GSLMatrix::SetFromDoubleArray(double * x)
+{
+  gsl_matrix_view m = gsl_matrix_view_array (x, rows, columns);
+  gsl_matrix_memcpy (matrix, &m.matrix);
+};
+
+/** This function returns the i-th element of a matrix.
+ * If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and 0 is returned.
+ * \param m row index
+ * \param n colum index
+ * \return (m,n)-th element of matrix
+ */
+double GSLMatrix::Get(size_t m, size_t n)
+{
+  return gsl_matrix_get (matrix, m, n);
+};
+
+/** This function sets the value of the \a m -th element of a matrix to \a x.
+ *  If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked.
+ * \param m row index
+ * \param m column index
+ * \param x value to set element (m,n) to
+ */
+void GSLMatrix::Set(size_t m, size_t n, double x)
+{
+  gsl_matrix_set (matrix, m, n, x);
+};
+
+/** These functions return a pointer to the \a m-th element of a matrix.
+ *  If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
+ * \param m index
+ * \return pointer to \a m-th element
+ */
+double *GSLMatrix::Pointer(size_t m, size_t n)
+{
+  return gsl_matrix_ptr (matrix, m, n);
+};
+
+/** These functions return a constant pointer to the \a m-th element of a matrix.
+ *  If \a m or \a n lies outside the allowed range of 0 to GSLMatrix::dimension-1 then the error handler is invoked and a null pointer is returned.
+ * \param m index
+ * \return const pointer to \a m-th element
+ */
+const double *GSLMatrix::const_Pointer(size_t m, size_t n)
+{
+  return gsl_matrix_const_ptr (matrix, m, n);
+};
+
+/* ========================== Initializing =============================== */
+/** This function sets all the elements of the matrix to the value \a x.
+ * \param *x
+ */
+void GSLMatrix::SetAll(double x)
+{
+  gsl_matrix_set_all (matrix, x);
+};
+
+/** This function sets all the elements of the matrix to zero.
+ */
+void GSLMatrix::SetZero()
+{
+  gsl_matrix_set_zero (matrix);
+};
+
+/** This function sets the elements of the matrix to the corresponding elements of the identity matrix, \f$m(i,j) = \delta(i,j)\f$, i.e. a unit diagonal with all off-diagonal elements zero.
+ * This applies to both square and rectangular matrices.
+ */
+void GSLMatrix::SetIdentity()
+{
+  gsl_matrix_set_identity (matrix);
+};
+
+/* ====================== Exchanging elements ============================ */
+/** This function exchanges the \a i-th and \a j-th row of the matrix in-place.
+ * \param i i-th row to swap with ...
+ * \param j ... j-th row to swap against
+ */
+bool GSLMatrix::SwapRows(size_t i, size_t j)
+{
+  return (gsl_matrix_swap_rows (matrix, i, j) == GSL_SUCCESS);
+};
+
+/** This function exchanges the \a i-th and \a j-th column of the matrix in-place.
+ * \param i i-th column to swap with ...
+ * \param j ... j-th column to swap against
+ */
+bool GSLMatrix::SwapColumns(size_t i, size_t j)
+{
+  return (gsl_matrix_swap_columns (matrix, i, j) == GSL_SUCCESS);
+};
+
+/** This function exchanges the \a i-th row and \a j-th column of the matrix in-place.
+ * The matrix must be square for this operation to be possible.
+ * \param i i-th row to swap with ...
+ * \param j ... j-th column to swap against
+ */
+bool GSLMatrix::SwapRowColumn(size_t i, size_t j)
+{
+  assert (rows == columns && "The matrix must be square for swapping row against column to be possible.");
+  return (gsl_matrix_swap_rowcol (matrix, i, j) == GSL_SUCCESS);
+};
+
+/** This function transposes the matrix.
+ * Note that the function is extended to the non-square case, where the matrix is re-allocated and copied.
+ */
+bool GSLMatrix::Transpose()
+{
+  if (rows == columns)// if square, use GSL
+    return (gsl_matrix_transpose (matrix) == GSL_SUCCESS);
+  else { // otherwise we have to copy a bit
+    gsl_matrix *dest = gsl_matrix_alloc(columns, rows);
+    for (size_t i=0;i<rows; i++)
+      for (size_t j=0;j<columns;j++) {
+        gsl_matrix_set(dest, j,i, gsl_matrix_get(matrix, i,j) );
+      }
+    delete(matrix);
+    matrix = dest;
+    flip(rows, columns);
+    return true;
+  }
+};
+
+/* ============================ Properties ============================== */
+/** Checks whether matrix' elements are strictly null.
+ * \return true - is null, false - else
+ */
+bool GSLMatrix::IsNull()
+{
+  return gsl_matrix_isnull (matrix);
+};
+
+/** Checks whether matrix' elements are strictly positive.
+ * \return true - is positive, false - else
+ */
+bool GSLMatrix::IsPositive()
+{
+  return gsl_matrix_ispos (matrix);
+};
+
+/** Checks whether matrix' elements are strictly negative.
+ * \return true - is negative, false - else
+ */
+bool GSLMatrix::IsNegative()
+{
+  return gsl_matrix_isneg (matrix);
+};
+
+/** Checks whether matrix' elements are strictly non-negative.
+ * \return true - is non-negative, false - else
+ */
+bool GSLMatrix::IsNonNegative()
+{
+  return gsl_matrix_isnonneg (matrix);
+};
+
+/** This function performs a Cholesky decomposition to determine whether matrix is positive definite.
+ * We check whether GSL returns GSL_EDOM as error, indicating that decomposition failed due to matrix not being positive-definite.
+ * \return true - matrix is positive-definite, false - else
+ */
+bool GSLMatrix::IsPositiveDefinite()
+{
+  if (rows != columns)  // only possible for square matrices.
+    return false;
+  else
+    return (gsl_linalg_cholesky_decomp (matrix) != GSL_EDOM);
+};
+
+
+/** Calculates the determinant of the matrix.
+ * if matrix is square, uses LU decomposition.
+ */
+double GSLMatrix::Determinant()
+{
+  int signum = 0;
+  assert (rows == columns && "Determinant can only be calculated for square matrices.");
+  gsl_permutation *p = gsl_permutation_alloc(rows);
+  gsl_linalg_LU_decomp(matrix, p, &signum);
+  gsl_permutation_free(p);
+  return gsl_linalg_LU_det(matrix, signum); 
+};
+
Index: src/gslmatrix.hpp
===================================================================
--- src/gslmatrix.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/gslmatrix.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,70 @@
+/*
+ * gslmatrix.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef GSLMATRIX_HPP_
+#define GSLMATRIX_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gsl/gsl_matrix.h>
+
+/****************************************** forward declarations *****************************/
+
+class GSLMatrix;
+
+/********************************************** declarations *******************************/
+
+class GSLMatrix {
+  friend class LinearSystemOfEquations;
+
+public:
+  GSLMatrix(size_t m, size_t n);
+  GSLMatrix(const GSLMatrix * const src);
+  ~GSLMatrix();
+
+  // Accessing
+  void SetFromDoubleArray(double * x);
+  double Get(size_t m, size_t n);
+  void Set(size_t m, size_t n, double x);
+  double *Pointer(size_t m, size_t n);
+  const double *const_Pointer(size_t m, size_t n);
+
+  // Initializing
+  void SetAll(double x);
+  void SetZero();
+  void SetIdentity();
+
+  // Exchanging elements
+  bool SwapRows(size_t i, size_t j);
+  bool SwapColumns(size_t i, size_t j);
+  bool SwapRowColumn(size_t i, size_t j);
+  bool Transpose();
+  bool IsNull();
+  bool IsPositive();
+  bool IsNegative();
+  bool IsNonNegative();
+  bool IsPositiveDefinite();
+  double Determinant();
+
+  GSLMatrix& operator=(const GSLMatrix& rhs);
+
+private:
+  gsl_matrix *matrix;
+
+  size_t rows;
+  size_t columns;
+};
+
+
+#endif /* GSLMATRIX_HPP_ */
Index: src/gslvector.cpp
===================================================================
--- src/gslvector.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/gslvector.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,129 @@
+/*
+ * gslvector.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#include "gslvector.hpp"
+
+/** Constructor of class GSLVector.
+ * Allocates GSL structures
+ * \param m dimension of vector
+ */
+GSLVector::GSLVector(size_t m) : dimension(m)
+{
+  vector = gsl_vector_calloc(dimension);
+};
+
+/** Copy constructor of class GSLVector.
+ * Allocates GSL structures and copies components from \a *src.
+ * \param *src source vector
+ */
+GSLVector::GSLVector(const GSLVector * const src) : dimension(src->dimension)
+{
+  vector = gsl_vector_alloc(dimension);
+  gsl_vector_memcpy (vector, src->vector);
+};
+
+/** Destructor of class GSLVector.
+ * Frees GSL structures
+ */
+GSLVector::~GSLVector()
+{
+  gsl_vector_free(vector);
+  dimension = 0;
+};
+
+/* ============================ Accessing =============================== */
+/** This function sets the vector from a double array.
+ * Creates a vector view of the array and performs a memcopy.
+ * \param *x array of values (no dimension check is performed)
+ */
+void GSLVector::SetFromDoubleArray(double * x)
+{
+  gsl_vector_view m = gsl_vector_view_array (x, dimension);
+  gsl_vector_memcpy (vector, &m.vector);
+};
+
+/** This function returns the i-th element of a vector.
+ * If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked and 0 is returned.
+ * \param m m-th element
+ * \return  m-th element of vector
+ */
+double GSLVector::Get(size_t m)
+{
+  return gsl_vector_get (vector, m);
+};
+
+/** This function sets the value of the \a m -th element of a vector to \a x.
+ *  If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked.
+ * \param m-th element to set
+ * \param x value to set to
+ */
+void GSLVector::Set(size_t m, double x)
+{
+  gsl_vector_set (vector, m, x);
+};
+
+/** These functions return a pointer to the \a m-th element of a vector.
+ *  If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked and a null pointer is returned.
+ * \param m m-th element
+ * \return pointer to \a m-th element
+ */
+double *GSLVector::Pointer(size_t m)
+{
+  return gsl_vector_ptr (vector, m);
+};
+
+/** These functions return a constant pointer to the \a m-th element of a vector.
+ *  If \a m lies outside the allowed range of 0 to GSLVector::dimension-1 then the error handler is invoked and a null pointer is returned.
+ * \param m m-th element
+ * \return const pointer to \a m-th element
+ */
+const double *GSLVector::const_Pointer(size_t m)
+{
+  return gsl_vector_const_ptr (vector, m);
+};
+
+/* ========================== Initializing =============================== */
+/** This function sets all the elements of the vector to the value \a x.
+ * \param *x
+ */
+void GSLVector::SetAll(double x)
+{
+  gsl_vector_set_all (vector, x);
+};
+
+/** This function sets all the elements of the vector to zero.
+ */
+void GSLVector::SetZero()
+{
+  gsl_vector_set_zero (vector);
+};
+
+/** This function makes a basis vector by setting all the elements of the vector to zero except for the i-th element which is set to one.
+ * \param i i-th component to set to unity (all other to zero)
+ * \return vector set
+ */
+int GSLVector::SetBasis(size_t i)
+{
+  return gsl_vector_set_basis (vector, i);
+};
+
+/* ====================== Exchanging elements ============================ */
+/** This function exchanges the \a i-th and \a j-th elements of the vector in-place.
+ * \param i i-th element to swap with ...
+ * \param j ... j-th element to swap against
+ */
+int GSLVector::SwapElements(size_t i, size_t j)
+{
+  return gsl_vector_swap_elements (vector, i, j);
+};
+
+/** This function reverses the order of the elements of the vector.
+ */
+int GSLVector::Reverse()
+{
+  return gsl_vector_reverse (vector);
+};
Index: src/gslvector.hpp
===================================================================
--- src/gslvector.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/gslvector.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,60 @@
+/*
+ * gslvector.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef GSLVECTOR_HPP_
+#define GSLVECTOR_HPP_
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gsl/gsl_vector.h>
+
+/****************************************** forward declarations *****************************/
+
+class GSLVector;
+
+/********************************************** declarations *******************************/
+
+class GSLVector {
+  friend class LinearSystemOfEquations;
+
+public:
+  GSLVector(size_t m);
+  GSLVector(const GSLVector * const src);
+  ~GSLVector();
+
+  // Accessing
+  void SetFromDoubleArray(double *x);
+  double Get(size_t m);
+  void Set(size_t m, double x);
+  double *Pointer(size_t m);
+  const double *const_Pointer(size_t m);
+
+  // Initializing
+  void SetAll(double x);
+  void SetZero();
+  int SetBasis(size_t m);
+
+  // Exchanging elements
+  int SwapElements(size_t i, size_t j);
+  int Reverse();
+
+private:
+  gsl_vector *vector;
+
+  size_t dimension;
+};
+
+
+
+#endif /* GSLVECTOR_HPP_ */
Index: src/helpers.hpp
===================================================================
--- src/helpers.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/helpers.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -74,4 +74,28 @@
   x = y;
   y = tmp;
+};
+
+/** returns greater of the two values.
+ * \param x first value
+ * \param y second value
+ * \return greater of the two (by operator>())
+ */
+template <typename T> T Max(T x, T y)
+{
+  if (x > y)
+    return x;
+  else return y;
+};
+
+/** returns smaller of the two values.
+ * \param x first value
+ * \param y second value
+ * \return smaller of the two (by operator<())
+ */
+template <typename T> T Min(T x, T y)
+{
+  if (x < y)
+    return x;
+  else return y;
 };
 
Index: src/joiner.cpp
===================================================================
--- src/joiner.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/joiner.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -7,4 +7,6 @@
 
 //============================ INCLUDES ===========================
+
+#include <cstring>
 
 #include "datacreator.hpp"
Index: src/linearsystemofequations.cpp
===================================================================
--- src/linearsystemofequations.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/linearsystemofequations.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,129 @@
+/*
+ * linearsystemofequations.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#include "defs.hpp"
+#include "gslmatrix.hpp"
+#include "gslvector.hpp"
+#include "linearsystemofequations.hpp"
+#include "logger.hpp"
+#include "vector.hpp"
+
+#include <cassert>
+#include <gsl/gsl_permutation.h>
+
+/** Constructor for class LinearSystemOfEquations.
+ * Allocates Vector and Matrix classes.
+ * \param m column dimension
+ * \param n row dimension
+ */
+LinearSystemOfEquations::LinearSystemOfEquations(int m, int n) : rows(m), columns(n), IsSymmetric(false)
+{
+  A = new GSLMatrix(m, n);
+  b = new GSLVector(m);
+  x = new GSLVector(n);
+};
+
+/** Destructor for class LinearSystemOfEquations.
+ * Destructs Vector and Matrix classes.
+ */
+LinearSystemOfEquations::~LinearSystemOfEquations()
+{
+  delete(A);
+  delete(b);
+  delete(x);
+};
+
+/** Sets whether matrix is to be regarded as symmetric.
+ * Note that we do not check whether it really is, just take upper diagonal.
+ * \param symmetric true or false
+ */
+bool LinearSystemOfEquations::SetSymmetric(bool symmetric)
+{
+  assert (rows == columns && "Rows and columns don't have equal size! Setting symmetric not possible.");
+  return (IsSymmetric = symmetric);
+};
+
+/** Initializes vector b to the components of the given vector.
+ * \param *x Vector with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::Setb(Vector *x)
+{
+  assert ( columns == NDIM && "Vector class is always three-dimensional, unlike this LEqS!");
+  b->SetFromDoubleArray(x->x);
+};
+
+/** Initializes vector b to the components of the given vector.
+ * \param *x array with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::Setb(double *x)
+{
+  b->SetFromDoubleArray(x);
+};
+
+/** Initializes matrix a to the components of the given array.
+ * note that sort order should be
+ * \param *x array with equal dimension (no check!)
+ */
+void LinearSystemOfEquations::SetA(double *x)
+{
+  A->SetFromDoubleArray(x);
+};
+
+/** Returns the solution vector x \f$A \cdot x = b\f$ as an array of doubles.
+ * \param *array pointer allocated array for solution on return (no bounds check, dimension must match)
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::GetSolutionAsArray(double *&array)
+{
+  bool status = Solve();
+
+  // copy solution
+  for (size_t i=0;i<x->dimension;i++) {
+    array[i] = x->Get(i);
+  }
+  return status;
+};
+
+/** Returns the solution vector x \f$A \cdot x = b\f$ as an array of doubles.
+ * \param &x solution vector on return (must be 3-dimensional)
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::GetSolutionAsVector(Vector &v)
+{
+  assert(rows == NDIM && "Solution can only be returned as vector if number of columns is NDIM.");
+  bool status = Solve();
+
+  // copy solution
+  for (size_t i=0;i<x->dimension;i++)
+    v.x[i] = x->Get(i);
+  return status;
+};
+
+/** Solves the given system of \f$A \cdot x = b\f$.
+ * Use either LU or Householder decomposition.
+ * Solution is stored in LinearSystemOfEquations::x
+ * \return true - solving possible, false - some error occured.
+ */
+bool LinearSystemOfEquations::Solve()
+{
+  // calculate solution
+  int s;
+  if (IsSymmetric) { // use LU
+    gsl_permutation * p = gsl_permutation_alloc (x->dimension);
+    gsl_linalg_LU_decomp (A->matrix, p, &s);
+    gsl_linalg_LU_solve (A->matrix, p, b->vector, x->vector);
+    gsl_permutation_free (p);
+  } else {  // use Householder
+    //GSLMatrix *backup = new GSLMatrix(rows,columns);
+    //*backup = *A;
+    gsl_linalg_HH_solve (A->matrix, b->vector, x->vector);
+    //*A = *backup;
+    //delete(backup);
+  }
+  return true;
+};
+
Index: src/linearsystemofequations.hpp
===================================================================
--- src/linearsystemofequations.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/linearsystemofequations.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,56 @@
+/*
+ * linearsystemofequations.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+/*********************************************** includes ***********************************/
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <gsl/gsl_linalg.h>
+
+/****************************************** forward declarations *****************************/
+
+class Vector;
+class GSLMatrix;
+class GSLVector;
+
+/********************************************** declarations *******************************/
+
+/** Solver class for a linear system of equations in the form \f$A \cdot x = b\F$.
+ *
+ */
+class LinearSystemOfEquations {
+public:
+  LinearSystemOfEquations(int n, int m);
+  ~LinearSystemOfEquations();
+
+  // setting
+  void Setb(Vector *x);
+  void Setb(double *x);
+  void SetA(double *x);
+  bool SetSymmetric(bool symmetric);
+
+  // solving
+  bool GetSolutionAsArray(double *&array);
+  bool GetSolutionAsVector(Vector &v);
+
+  GSLVector *b;
+  GSLVector *x;
+  GSLMatrix *A;
+
+private:
+  bool Solve();
+
+  int rows;
+  int columns;
+
+  bool IsSymmetric;
+};
Index: src/memoryallocator.hpp
===================================================================
--- src/memoryallocator.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/memoryallocator.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -16,4 +16,5 @@
 #endif
 
+#include <cstdlib>
 #include <iostream>
 #include <iomanip>
Index: src/memoryusageobserver.cpp
===================================================================
--- src/memoryusageobserver.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/memoryusageobserver.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -4,4 +4,6 @@
  * This class represents a Singleton for observing memory usage.
  */
+
+#include <cstdlib>
 
 #include "log.hpp"
Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -4,4 +4,6 @@
  *
  */
+
+#include <cstring>
 
 #include "atom.hpp"
@@ -587,5 +589,5 @@
   else
     molname = filename; // contains no slashes
-  char *endname = strchr(molname, '.');
+  const char *endname = strchr(molname, '.');
   if ((endname == NULL) || (endname < molname))
     length = strlen(molname);
Index: src/molecule.hpp
===================================================================
--- src/molecule.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -110,4 +110,5 @@
   TesselPoint *GetPoint() const ;
   TesselPoint *GetTerminalPoint() const ;
+  int GetMaxId() const;
   void GoToNext() const ;
   void GoToPrevious() const ;
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule_dynamics.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -370,9 +370,9 @@
 
   // argument minimise the constrained potential in this injective PermutationMap
-  Log() << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
+  Log() << Verbose(1) << "Argument minimising the PermutationMap." << endl;
   OldPotential = 1e+10;
   round = 0;
   do {
-    Log() << Verbose(2) << "Starting round " << ++round << " ... " << endl;
+    Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl;
     OlderPotential = OldPotential;
     do {
Index: src/molecule_fragmentation.cpp
===================================================================
--- src/molecule_fragmentation.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule_fragmentation.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -5,4 +5,6 @@
  *      Author: heber
  */
+
+#include <cstring>
 
 #include "atom.hpp"
Index: src/molecule_geometry.cpp
===================================================================
--- src/molecule_geometry.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule_geometry.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -101,9 +101,9 @@
 {
   int Num = 0;
-  atom *ptr = start->next;  // start at first in list
+  atom *ptr = start;  // start at first in list
 
   Center.Zero();
 
-  if (ptr != end) {   //list not empty?
+  if (ptr->next != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
       ptr = ptr->next;
Index: src/molecule_graph.cpp
===================================================================
--- src/molecule_graph.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule_graph.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -1117,6 +1117,5 @@
   bool status = true;
   if (ReferenceStack->IsEmpty()) {
-    eLog() << Verbose(0) << "ReferenceStack is empty!" << endl;
-    performCriticalExit();
+    Log() << Verbose(1) << "ReferenceStack is empty!" << endl;
     return false;
   }
Index: src/molecule_pointcloud.cpp
===================================================================
--- src/molecule_pointcloud.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/molecule_pointcloud.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -50,4 +50,12 @@
 };
 
+/** Return the greatest index of all atoms in the list.
+ * \return greatest index
+ */
+int molecule::GetMaxId() const
+{
+  return last_atom;
+};
+
 /** Go to next atom.
  * Stops at last one.
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/moleculelist.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -4,4 +4,6 @@
  *
  */
+
+#include <cstring>
 
 #include "atom.hpp"
@@ -402,6 +404,5 @@
   input.open(line.c_str());
   if (input == NULL) {
-    eLog() << Verbose(0) << endl << "Unable to open " << line << ", is the directory correct?" << endl;
-    performCriticalExit();
+    Log() << Verbose(1) << endl << "Unable to open " << line << ", is the directory correct?" << endl;
     return false;
   }
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/parser.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -6,4 +6,6 @@
 
 // ======================================= INCLUDES ==========================================
+
+#include <cstring>
 
 #include "helpers.hpp"
@@ -158,6 +160,6 @@
   //Log() << Verbose(0) << "Opening " << name << " ... "  << input << endl;
   if (input == NULL) {
-    eLog() << Verbose(0) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
-    performCriticalExit();
+    eLog() << Verbose(1) << endl << "Unable to open " << name << ", is the directory correct?" << endl;
+    //performCriticalExit();
     return false;
   }
Index: src/periodentafel.cpp
===================================================================
--- src/periodentafel.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/periodentafel.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -9,4 +9,5 @@
 #include <iomanip>
 #include <fstream>
+#include <cstring>
 
 #include "element.hpp"
Index: src/tesselation.cpp
===================================================================
--- src/tesselation.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/tesselation.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -35,5 +35,5 @@
  * \param *Walker TesselPoint this boundary point represents
  */
-BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) :
+BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) :
   LinesCount(0),
   node(Walker),
@@ -61,5 +61,5 @@
  * \param *line line to add
  */
-void BoundaryPointSet::AddLine(class BoundaryLineSet *line)
+void BoundaryPointSet::AddLine(BoundaryLineSet * const line)
 {
 	Info FunctionInfo(__func__);
@@ -105,5 +105,5 @@
  * \param number number of the list
  */
-BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number)
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number)
 {
 	Info FunctionInfo(__func__);
@@ -115,4 +115,26 @@
   Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
   Point[1]->AddLine(this); //
+  // set skipped to false
+  skipped = false;
+  // clear triangles list
+  Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl;
+};
+
+/** Constructor of BoundaryLineSet with two endpoints.
+ * Adds line automatically to each endpoints' LineMap
+ * \param *Point1 first boundary point
+ * \param *Point2 second boundary point
+ * \param number number of the list
+ */
+BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number)
+{
+  Info FunctionInfo(__func__);
+  // set number
+  Nr = number;
+  // set endpoints in ascending order
+  SetEndpointsOrdered(endpoints, Point1, Point2);
+  // add this line to the hash maps of both endpoints
+  Point1->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding.
+  Point2->AddLine(this); //
   // set skipped to false
   skipped = false;
@@ -171,5 +193,5 @@
  * \param *triangle to add
  */
-void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
+void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle)
 {
 	Info FunctionInfo(__func__);
@@ -182,5 +204,5 @@
  * \return true - common endpoint present, false - not connected
  */
-bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line)
+bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const
 {
 	Info FunctionInfo(__func__);
@@ -197,5 +219,5 @@
  * \return true - triangles are convex, false - concave or less than two triangles connected
  */
-bool BoundaryLineSet::CheckConvexityCriterion()
+bool BoundaryLineSet::CheckConvexityCriterion() const
 {
 	Info FunctionInfo(__func__);
@@ -221,5 +243,5 @@
   int i=0;
   class BoundaryPointSet *node = NULL;
-  for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
+  for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
     //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
@@ -264,5 +286,5 @@
  * \return true - point is of the line, false - is not
  */
-bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
 	Info FunctionInfo(__func__);
@@ -277,5 +299,5 @@
  * \return NULL - if endpoint not contained in BoundaryLineSet, or pointer to BoundaryPointSet otherwise
  */
-class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point)
+class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const
 {
 	Info FunctionInfo(__func__);
@@ -317,5 +339,5 @@
  * \param number number of triangle
  */
-BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) :
+BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number) :
   Nr(number)
 {
@@ -376,5 +398,5 @@
  * \param &OtherVector direction vector to make normal vector unique.
  */
-void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector)
+void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector)
 {
 	Info FunctionInfo(__func__);
@@ -388,7 +410,7 @@
 };
 
-/** Finds the point on the triangle \a *BTS the line defined by \a *MolCenter and \a *x crosses through.
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
- * This we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
  * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
  * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
@@ -400,7 +422,7 @@
  * \return true - \a *Intersection contains intersection on plane defined by triangle, false - zero vector if outside of triangle.
  */
-bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection)
-{
-	Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const
+{
+  Info FunctionInfo(__func__);
   Vector CrossPoint;
   Vector helper;
@@ -411,4 +433,18 @@
   }
 
+  Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl;
+  Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
+
+  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl;
+    return true;
+  }   else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {
+    Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
+    return true;
+  }
   // Calculate cross point between one baseline and the line from the third endpoint to intersection
   int i=0;
@@ -417,21 +453,107 @@
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
+      CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+      const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared();
+      Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl;
+      if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) {
+        Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl;
+        i=4;
+        break;
+      }
+      i++;
     } else 
-      i++;
-    if (i>2)
       break;
-  } while (CrossPoint.NormSquared() < MYEPSILON);
+  } while (i<3);
   if (i==3) {
-    eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl;
-  }
-  CrossPoint.SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
-
-  // check whether intersection is inside or not by comparing length of intersection and length of cross point
-  if ((CrossPoint.NormSquared() - helper.NormSquared()) < MYEPSILON) { // inside
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl;
     return true;
-  } else { // outside!
-    Intersection->Zero();
+  } else {
+    Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl;
     return false;
   }
+};
+
+/** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
+ * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
+ * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
+ * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
+ * given by the intersection and the third basepoint. Then, we check whether it's on the baseline (i.e. between
+ * the first two basepoints) or not.
+ * \param *x point
+ * \param *ClosestPoint desired closest point inside triangle to \a *x, is absolute vector
+ * \return Distance squared between \a *x and closest point inside triangle
+ */
+double BoundaryTriangleSet::GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const
+{
+  Info FunctionInfo(__func__);
+  Vector Direction;
+
+  // 1. get intersection with plane
+  Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
+  GetCenter(&Direction);
+  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
+    ClosestPoint->CopyVector(x);
+  }
+
+  // 2. Calculate in plane part of line (x, intersection)
+  Vector InPlane;
+  InPlane.CopyVector(x);
+  InPlane.SubtractVector(ClosestPoint);  // points from plane intersection to straight-down point
+  InPlane.ProjectOntoPlane(&NormalVector);
+  InPlane.AddVector(ClosestPoint);
+
+  Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
+  Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl;
+  Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl;
+
+  // Calculate cross point between one baseline and the desired point such that distance is shortest
+  double ShortestDistance = -1.;
+  bool InsideFlag = false;
+  Vector CrossDirection[3];
+  Vector CrossPoint[3];
+  Vector helper;
+  for (int i=0;i<3;i++) {
+    // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point
+    Direction.CopyVector(endpoints[(i+1)%3]->node->node);
+    Direction.SubtractVector(endpoints[i%3]->node->node);
+    // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
+    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
+    CrossDirection[i].CopyVector(&CrossPoint[i]);
+    CrossDirection[i].SubtractVector(&InPlane);
+    CrossPoint[i].SubtractVector(endpoints[i%3]->node->node);  // cross point was returned as absolute vector
+    const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared();
+    Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl;
+    if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) {
+      CrossPoint[i].AddVector(endpoints[i%3]->node->node);  // make cross point absolute again
+      Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl;
+      const double distance = CrossPoint[i].DistanceSquared(x);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(&CrossPoint[i]);
+      }
+    } else
+      CrossPoint[i].Zero();
+  }
+  InsideFlag = true;
+  for (int i=0;i<3;i++) {
+    const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]);
+    const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);;
+    if ((sign > -MYEPSILON) && (othersign > -MYEPSILON))  // have different sign
+      InsideFlag = false;
+  }
+  if (InsideFlag) {
+    ClosestPoint->CopyVector(&InPlane);
+    ShortestDistance = InPlane.DistanceSquared(x);
+  } else {  // also check endnodes
+    for (int i=0;i<3;i++) {
+      const double distance = x->DistanceSquared(endpoints[i]->node->node);
+      if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
+        ShortestDistance = distance;
+        ClosestPoint->CopyVector(endpoints[i]->node->node);
+      }
+    }
+  }
+  Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl;
+  return ShortestDistance;
 };
 
@@ -440,5 +562,5 @@
  * \return true - line is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line)
+bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const
 {
 	Info FunctionInfo(__func__);
@@ -453,5 +575,5 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const
 {
 	Info FunctionInfo(__func__);
@@ -466,5 +588,5 @@
  * \return true - point is of the triangle, false - is not
  */
-bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point)
+bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const
 {
 	Info FunctionInfo(__func__);
@@ -479,7 +601,8 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3])
-{
-	Info FunctionInfo(__func__);
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const
+{
+	Info FunctionInfo(__func__);
+	Log() << Verbose(1) << "INFO: Checking " << Points[0] << ","  << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl;
   return (((endpoints[0] == Points[0])
             || (endpoints[0] == Points[1])
@@ -501,5 +624,5 @@
  * \return true - is the very triangle, false - is not
  */
-bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T)
+bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const
 {
 	Info FunctionInfo(__func__);
@@ -523,5 +646,5 @@
  * \return pointer third endpoint or NULL if line does not belong to triangle.
  */
-class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line)
+class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const
 {
 	Info FunctionInfo(__func__);
@@ -540,5 +663,5 @@
  * \param *center central point on return.
  */
-void BoundaryTriangleSet::GetCenter(Vector *center)
+void BoundaryTriangleSet::GetCenter(Vector * const center) const
 {
 	Info FunctionInfo(__func__);
@@ -547,4 +670,5 @@
     center->AddVector(endpoints[i]->node->node);
   center->Scale(1./3.);
+  Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
 }
 
@@ -1422,5 +1546,5 @@
   TesselPoint *Walker = NULL;
   Vector *Center = cloud->GetCenter();
-  list<BoundaryTriangleSet*> *triangles = NULL;
+  TriangleList *triangles = NULL;
   bool AddFlag = false;
   LinkedCell *BoundaryPoints = NULL;
@@ -1437,5 +1561,5 @@
     Log() << Verbose(0) << "Current point is " << *Walker << "." << endl;
     // get the next triangle
-    triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints);
+    triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints);
     BTS = triangles->front();
     if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
@@ -1587,9 +1711,10 @@
   bool insertNewLine = true;
 
-  if (a->lines.find(b->node->nr) != a->lines.end()) {
-    LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  LineMap::iterator FindLine = a->lines.find(b->node->nr);
+  if (FindLine != a->lines.end()) {
+    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
+
     pair<LineMap::iterator,LineMap::iterator> FindPair;
     FindPair = a->lines.equal_range(b->node->nr);
-    Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;
 
     for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
@@ -1915,10 +2040,13 @@
   double maxCoordinate[NDIM];
   BoundaryLineSet BaseLine;
-  Vector Oben;
   Vector helper;
   Vector Chord;
   Vector SearchDirection;
-
-  Oben.Zero();
+  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 SphereCenter;
+  Vector NormalVector;
+
+  NormalVector.Zero();
 
   for (i = 0; i < 3; i++) {
@@ -1955,6 +2083,6 @@
   BTS = NULL;
   for (int k=0;k<NDIM;k++) {
-    Oben.Zero();
-    Oben.x[k] = 1.;
+    NormalVector.Zero();
+    NormalVector.x[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
@@ -1963,27 +2091,32 @@
     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.
 
-    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, Oben, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
+    FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
     if (Temporary == NULL)  // have we found a second point?
       continue;
     BaseLine.endpoints[1] = new BoundaryPointSet(Temporary);
 
-    helper.CopyVector(BaseLine.endpoints[0]->node->node);
-    helper.SubtractVector(BaseLine.endpoints[1]->node->node);
-    helper.Normalize();
-    Oben.ProjectOntoPlane(&helper);
-    Oben.Normalize();
-    helper.VectorProduct(&Oben);
+    // construct center of circle
+    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
+    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
+    CircleCenter.Scale(0.5);
+
+    // construct normal vector of circle
+    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
+    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+
+    double radius = CirclePlaneNormal.NormSquared();
+    double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
+
+    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.Normalize();
     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
 
-    Chord.CopyVector(BaseLine.endpoints[0]->node->node); // bring into calling function
-    Chord.SubtractVector(BaseLine.endpoints[1]->node->node);
-    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)
+    SphereCenter.CopyVector(&NormalVector);
+    SphereCenter.Scale(CircleRadius);
+    SphereCenter.AddVector(&CircleCenter);
+    // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&Chord, &Oben);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -1993,5 +2126,5 @@
     //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
     CandidateForTesselation OptCandidates(&BaseLine);
-    FindThirdPointForTesselation(Oben, SearchDirection, helper, OptCandidates, NULL, RADIUS, LC);
+    FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
     Log() << Verbose(0) << "List of third Points is:" << endl;
     for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
@@ -2167,5 +2300,5 @@
   Vector CircleCenter;
   Vector CirclePlaneNormal;
-  Vector OldSphereCenter;
+  Vector RelativeSphereCenter;
   Vector SearchDirection;
   Vector helper;
@@ -2174,8 +2307,10 @@
   double radius, CircleRadius;
 
-  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;
   for (int i=0;i<3;i++)
-    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node))
+    if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {
       ThirdNode = T.endpoints[i]->node;
+      break;
+    }
+  Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;
 
   // construct center of circle
@@ -2191,27 +2326,22 @@
   radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
+    // construct relative sphere center with now known CircleCenter
+    RelativeSphereCenter.CopyVector(&T.SphereCenter);
+    RelativeSphereCenter.SubtractVector(&CircleCenter);
+
     CircleRadius = RADIUS*RADIUS - radius/4.;
     CirclePlaneNormal.Normalize();
     Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
 
-    // construct old center
-    GetCenterofCircumcircle(&OldSphereCenter, *T.endpoints[0]->node->node, *T.endpoints[1]->node->node, *T.endpoints[2]->node->node);
-    helper.CopyVector(&T.NormalVector);  // normal vector ensures that this is correct center of the two possible ones
-    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
-    helper.Scale(sqrt(RADIUS*RADIUS - radius));
-    OldSphereCenter.AddVector(&helper);
-    OldSphereCenter.SubtractVector(&CircleCenter);
-    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
-
-    // construct SearchDirection
-    SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal);
-    helper.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
+    Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl;
+
+    // construct SearchDirection and an "outward pointer"
+    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
     if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
-    SearchDirection.ProjectOntoPlane(&OldSphereCenter);
-    SearchDirection.Normalize();
     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -2219,5 +2349,5 @@
 
     // add third point
-    FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
+    FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);
 
   } else {
@@ -2332,13 +2462,14 @@
 
   // fill the set of neighbours
-  Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  Center.SubtractVector(TurningPoint->node);
-  set<TesselPoint*> SetOfNeighbours;
+  TesselPointSet SetOfNeighbours;
   SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
   for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
     SetOfNeighbours.insert(*Runner);
-  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);
+  TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
+  Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl;
+  for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
+    Log() << Verbose(0) << **TesselRunner << endl;
   TesselPointList::iterator Runner = connectedClosestPoints->begin();
   TesselPointList::iterator Sprinter = Runner;
@@ -2350,5 +2481,4 @@
     AddTesselationPoint((*Sprinter), 2);
 
-    Center.CopyVector(&CandidateLine.OptCenter);
     // add the lines
     AddTesselationLine(TPS[0], TPS[1], 0);
@@ -2359,5 +2489,7 @@
     BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
     AddTesselationTriangle();
-    Center.Scale(-1.);
+    BTS->GetCenter(&Center);
+    Center.SubtractVector(&CandidateLine.OptCenter);
+    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
     BTS->GetNormalVector(Center);
 
@@ -2365,4 +2497,7 @@
     Runner = Sprinter;
     Sprinter++;
+    Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl;
+    if (Sprinter != connectedClosestPoints->end())
+      Log() << Verbose(0) << " There are still more triangles to add." << endl;
   }
   delete(connectedClosestPoints);
@@ -2766,6 +2901,9 @@
   Vector NewNormalVector;   // normal vector of the Candidate's triangle
   Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
+  Vector RelativeOldSphereCenter;
+  Vector NewPlaneCenter;
   double CircleRadius; // radius of this circle
   double radius;
+  double otherradius;
   double alpha, Otheralpha; // angles (i.e. parameter for the circle).
   int N[NDIM], Nlower[NDIM], Nupper[NDIM];
@@ -2783,23 +2921,26 @@
   CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
 
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+
   // calculate squared radius TesselPoint *ThirdNode,f circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
-  if (radius/4. < RADIUS*RADIUS) {
-    CircleRadius = RADIUS*RADIUS - radius/4.;
+  radius = CirclePlaneNormal.NormSquared()/4.;
+  if (radius < RADIUS*RADIUS) {
+    CircleRadius = RADIUS*RADIUS - radius;
     CirclePlaneNormal.Normalize();
-    //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+    Log() << Verbose(1) << "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) {
-      eLog() << Verbose(1) << "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(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
+      eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;
+      RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);
+    }
+    radius = RelativeOldSphereCenter.NormSquared();
     if (fabs(radius - CircleRadius) < HULLEPSILON) {
-      //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
+      Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;
 
       // check SearchDirection
-      //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -2832,28 +2973,35 @@
 
                 // check for three unique points
-                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl;
+                Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl;
                 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){
 
-                  // construct both new centers
-                  GetCenterofCircumcircle(&NewSphereCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
-                  OtherNewSphereCenter.CopyVector(&NewSphereCenter);
-
-                  if ((NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node))
-                  && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON)
+                  // find center on the plane
+                  GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);
+                  Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
+
+                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
+                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
                   ) {
-                    helper.CopyVector(&NewNormalVector);
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);
+                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
+                    Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
+                    Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl;
                     if (radius < RADIUS*RADIUS) {
+                      otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);
+                      if (fabs(radius - otherradius) > HULLEPSILON) {
+                        eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;
+                      }
+                      // construct both new centers
+                      NewSphereCenter.CopyVector(&NewPlaneCenter);
+                      OtherNewSphereCenter.CopyVector(&NewPlaneCenter);
+                      helper.CopyVector(&NewNormalVector);
                       helper.Scale(sqrt(RADIUS*RADIUS - radius));
-                      Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;
+                      Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl;
                       NewSphereCenter.AddVector(&helper);
-                      NewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;
-
                       // OtherNewSphereCenter is created by the same vector just in the other direction
                       helper.Scale(-1.);
                       OtherNewSphereCenter.AddVector(&helper);
-                      OtherNewSphereCenter.SubtractVector(&CircleCenter);
                       Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 
@@ -2861,4 +3009,5 @@
                       Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection);
                       alpha = min(alpha, Otheralpha);
+
                       // if there is a better candidate, drop the current list and add the new candidate
                       // otherwise ignore the new candidate and keep the list
@@ -2892,5 +3041,4 @@
                         }
                       }
-
                     } else {
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
@@ -2936,6 +3084,6 @@
   const BoundaryLineSet * lines[2] = { line1, line2 };
   class BoundaryPointSet *node = NULL;
-  map<int, class BoundaryPointSet *> OrderMap;
-  pair<map<int, class BoundaryPointSet *>::iterator, bool> OrderTest;
+  PointMap OrderMap;
+  PointTestPair OrderTest;
   for (int i = 0; i < 2; i++)
     // for both lines
@@ -2957,85 +3105,209 @@
 };
 
+/** Finds the boundary points that are closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
+ */
+DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+  PointMap::const_iterator FindPoint;
+  int N[NDIM], Nlower[NDIM], Nupper[NDIM];
+
+  if (LinesOnBoundary.empty()) {
+    eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl;
+    return NULL;
+  }
+
+  // gather all points close to the desired one
+  LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
+  for(int i=0;i<NDIM;i++) // store indices of this cell
+    N[i] = LC->n[i];
+  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
+
+  DistanceToPointMap * points = new DistanceToPointMap;
+  LC->GetNeighbourBounds(Nlower, Nupper);
+  //Log() << Verbose(1) << 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]++) {
+        const LinkedNodes *List = LC->GetCurrentCell();
+        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
+        if (List != NULL) {
+          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
+            FindPoint = PointsOnBoundary.find((*Runner)->nr);
+            if (FindPoint != PointsOnBoundary.end()) {
+              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
+            }
+          }
+        } else {
+          eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl;
+        }
+      }
+
+  // check whether we found some points
+  if (points->empty()) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    delete(points);
+    return NULL;
+  }
+  return points;
+};
+
+/** Finds the boundary line that is closest to a given Vector \a *x.
+ * \param *out output stream for debugging
+ * \param *x Vector to look from
+ * \return closest BoundaryLineSet or NULL in degenerate case.
+ */
+BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const
+{
+  Info FunctionInfo(__func__);
+
+  // get closest points
+  DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
+    return NULL;
+  }
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl;
+  BoundaryLineSet *ClosestLine = NULL;
+  double MinDistance = -1.;
+  Vector helper;
+  Vector Center;
+  Vector BaseLine;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+      // calculate closest point on line to desired point
+      helper.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      helper.AddVector((LineRunner->second)->endpoints[1]->node->node);
+      helper.Scale(0.5);
+      Center.CopyVector(x);
+      Center.SubtractVector(&helper);
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      Center.ProjectOntoPlane(&BaseLine);
+      const double distance = Center.NormSquared();
+      if ((ClosestLine == NULL) || (distance < MinDistance)) {
+        // additionally calculate intersection on line (whether it's on the line section or not)
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthA = helper.ScalarProduct(&BaseLine);
+        helper.CopyVector(x);
+        helper.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        helper.SubtractVector(&Center);
+        const double lengthB = helper.ScalarProduct(&BaseLine);
+        if (lengthB*lengthA < 0) {  // if have different sign
+          ClosestLine = LineRunner->second;
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl;
+        }
+      } else {
+        Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl;
+      }
+    }
+  }
+  delete(points);
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLine == NULL) {
+    Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
+    return NULL;
+  }
+  return ClosestLine;
+};
+
+
 /** Finds the triangle that is closest to a given Vector \a *x.
  * \param *out output stream for debugging
  * \param *x Vector to look from
- * \return list of BoundaryTriangleSet of nearest triangles or NULL in degenerate case.
- */
-list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const
-{
-	Info FunctionInfo(__func__);
-  TesselPoint *trianglePoints[3];
-  TesselPoint *SecondPoint = NULL;
-  list<BoundaryTriangleSet*> *triangles = NULL;
-
-  if (LinesOnBoundary.empty()) {
-    eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";
+ * \return BoundaryTriangleSet of nearest triangle or NULL.
+ */
+TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const
+{
+	Info FunctionInfo(__func__);
+
+	// get closest points
+	DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC);
+  if (points == NULL) {
+    eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl;
     return NULL;
   }
-  Log() << Verbose(1) << "Finding closest Tesselpoint to " << *x << " ... " << endl;
-  trianglePoints[0] = FindClosestPoint(x, SecondPoint, LC);
-  
-  // check whether closest point is "too close" :), then it's inside
-  if (trianglePoints[0] == NULL) {
+
+  // for each point, check its lines, remember closest
+  Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl;
+  LineSet ClosestLines;
+  double MinDistance = 1e+16;
+  Vector BaseLineIntersection;
+  Vector Center;
+  Vector BaseLine;
+  Vector BaseLineCenter;
+  for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
+    for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
+
+      BaseLine.CopyVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLine.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthBase = BaseLine.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      const double lengthEndA = BaseLineIntersection.NormSquared();
+
+      BaseLineIntersection.CopyVector(x);
+      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      const double lengthEndB = BaseLineIntersection.NormSquared();
+
+      if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) {  // intersection would be outside, take closer endpoint
+        const double lengthEnd = Min(lengthEndA, lengthEndB);
+        if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
+          ClosestLines.clear();
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = lengthEnd;
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl;
+        } else if  (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
+          ClosestLines.insert(LineRunner->second);
+          Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl;
+        } else { // line is worse
+          Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl;
+        }
+      } else { // intersection is closer, calculate
+        // calculate closest point on line to desired point
+        BaseLineIntersection.CopyVector(x);
+        BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+        Center.CopyVector(&BaseLineIntersection);
+        Center.ProjectOntoPlane(&BaseLine);
+        BaseLineIntersection.SubtractVector(&Center);
+        const double distance = BaseLineIntersection.NormSquared();
+        if (Center.NormSquared() > BaseLine.NormSquared()) {
+          eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl;
+        }
+        if ((ClosestLines.empty()) || (distance < MinDistance)) {
+          ClosestLines.insert(LineRunner->second);
+          MinDistance = distance;
+          Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl;
+        } else {
+          Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl;
+        }
+      }
+    }
+  }
+  delete(points);
+
+  // check whether closest line is "too close" :), then it's inside
+  if (ClosestLines.empty()) {
     Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl;
     return NULL;
   }
-  if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) {
-    Log() << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl;
-    PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr);
-    triangles = new list<BoundaryTriangleSet*>;
-    if (PointRunner != PointsOnBoundary.end()) {
-      for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-        for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-          triangles->push_back(TriangleRunner->second);
-      triangles->sort();
-      triangles->unique();
-    } else {
-      PointRunner = PointsOnBoundary.find(SecondPoint->nr);
-      trianglePoints[0] = SecondPoint;
-      if (PointRunner != PointsOnBoundary.end()) {
-        for(LineMap::iterator LineRunner = PointRunner->second->lines.begin(); LineRunner != PointRunner->second->lines.end(); LineRunner++)
-          for(TriangleMap::iterator TriangleRunner = LineRunner->second->triangles.begin(); TriangleRunner != LineRunner->second->triangles.end(); TriangleRunner++)
-            triangles->push_back(TriangleRunner->second);
-        triangles->sort();
-        triangles->unique();
-      } else {
-        eLog() << Verbose(1) << "I cannot find a boundary point to the tessel point " << *trianglePoints[0] << "." << endl;
-        return NULL;
-      }
-    }
-  } else {
-    set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(trianglePoints[0]);
-    TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(connectedPoints, trianglePoints[0], x);
-    delete(connectedPoints);
-    if (connectedClosestPoints != NULL) {
-      trianglePoints[1] = connectedClosestPoints->front();
-      trianglePoints[2] = connectedClosestPoints->back();
-      for (int i=0;i<3;i++) {
-        if (trianglePoints[i] == NULL) {
-          eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl;
-        }
-        //Log() << Verbose(1) << "List of triangle points:" << endl;
-        //Log() << Verbose(2) << *trianglePoints[i] << endl;
-      }
-
-      triangles = FindTriangles(trianglePoints);
-      Log() << Verbose(1) << "List of possible triangles:" << endl;
-      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-        Log() << Verbose(2) << **Runner << endl;
-
-      delete(connectedClosestPoints);
-    } else {
-      triangles = NULL;
-      eLog() << Verbose(2) << "There is no circle of connected points!" << endl;
-    }
-  }
-  
-  if ((triangles == NULL) || (triangles->empty())) {
-    eLog() << Verbose(1) << "There is no nearest triangle. Please check the tesselation structure.";
-    delete(triangles);
-    return NULL;
-  } else
-    return triangles;
+  TriangleList * candidates = new TriangleList;
+  for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
+    for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
+    candidates->push_back(Runner->second);
+  }
+  return candidates;
 };
 
@@ -3046,77 +3318,120 @@
  * \return list of BoundaryTriangleSet of nearest triangles or NULL.
  */
-class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const
+class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const
 {
 	Info FunctionInfo(__func__);
   class BoundaryTriangleSet *result = NULL;
-  list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC);
+  TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
+  TriangleList candidates;
   Vector Center;
-
-  if (triangles == NULL)
+  Vector helper;
+
+  if ((triangles == NULL) || (triangles->empty()))
     return NULL;
 
-  if (triangles->size() == 1) { // there is no degenerate case
-    result = triangles->front();
-    Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;
-  } else {
-    result = triangles->front();
-    result->GetCenter(&Center);
-    Center.SubtractVector(x);
-    Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;
-    if (Center.ScalarProduct(&result->NormalVector) < 0) {
-      result = triangles->back();
-      Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;
-      if (Center.ScalarProduct(&result->NormalVector) < 0) {
-        eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl;
-      }
+  // go through all and pick the one with the best alignment to x
+  double MinAlignment = 2.*M_PI;
+  for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
+    (*Runner)->GetCenter(&Center);
+    helper.CopyVector(x);
+    helper.SubtractVector(&Center);
+    const double Alignment = helper.Angle(&(*Runner)->NormalVector);
+    if (Alignment < MinAlignment) {
+      result = *Runner;
+      MinAlignment = Alignment;
+      Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl;
+    } else {
+      Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl;
     }
   }
   delete(triangles);
+
   return result;
 };
 
-/** Checks whether the provided Vector is within the tesselation structure.
+/** Checks whether the provided Vector is within the Tesselation structure.
+ * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
+ * @param point of which to check the position
+ * @param *LC LinkedCell structure
+ *
+ * @return true if the point is inside the Tesselation structure, false otherwise
+ */
+bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
+{
+  return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
+}
+
+/** Returns the distance to the surface given by the tesselation.
+ * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
+ * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
+ * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
+ * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
+ * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
+ *  -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
+ *  -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
+ *  -# If inside, take it to calculate closest distance
+ *  -# If not, take intersection with BoundaryLine as distance
+ *
+ * @note distance is squared despite it still contains a sign to determine in-/outside!
  *
  * @param point of which to check the position
  * @param *LC LinkedCell structure
  *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
-{
-	Info FunctionInfo(__func__);
-  class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC);
+ * @return >0 if outside, ==0 if on surface, <0 if inside (Note that distance can be at most LinkedCell::RADIUS.)
+ */
+double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
+{
+  Info FunctionInfo(__func__);
+  class BoundaryTriangleSet *result = FindClosestTriangleToVector(&Point, LC);
   Vector Center;
+  Vector helper;
+  Vector DistanceToCenter;
+  Vector Intersection;
+  double distance = 0.;
 
   if (result == NULL) {// is boundary point or only point in point cloud?
     Log() << Verbose(1) << Point << " is the only point in vicinity." << endl;
-    return false;
+    return LC->RADIUS;
+  } else {
+    Log() << Verbose(1) << "INFO: Closest triangle found is " << *result << " with normal vector " << result->NormalVector << "." << endl;
   }
 
   result->GetCenter(&Center);
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  Center.SubtractVector(&Point);
-  Log() << Verbose(2) << "INFO: Vector from center to point to test is " << Center << "." << endl;
-  if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) {
-    Log() << Verbose(1) << Point << " is an inner point." << endl;
-    return true;
+  DistanceToCenter.CopyVector(&Center);
+  DistanceToCenter.SubtractVector(&Point);
+  Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl;
+
+  // check whether we are on boundary
+  if (fabs(DistanceToCenter.ScalarProduct(&result->NormalVector)) < MYEPSILON) {
+    // calculate whether inside of triangle
+    DistanceToCenter.CopyVector(&Point);
+    Center.CopyVector(&Point);
+    Center.SubtractVector(&result->NormalVector); // points towards MolCenter
+    DistanceToCenter.AddVector(&result->NormalVector); // points outside
+    Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
+    if (result->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
+      Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl;
+      return 0.;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl;
+      return false;
+    }
   } else {
-    Log() << Verbose(1) << Point << " is NOT an inner point." << endl;
-    return false;
-  }
-}
-
-/** Checks whether the provided TesselPoint is within the tesselation structure.
- *
- * @param *Point of which to check the position
- * @param *LC Linked Cell structure
- *
- * @return true if the point is inside the tesselation structure, false otherwise
- */
-bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const
-{
-	Info FunctionInfo(__func__);
-  return IsInnerPoint(*(Point->node), LC);
-}
+    // calculate smallest distance
+    distance = result->GetClosestPointInsideTriangle(&Point, &Intersection);
+    Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl;
+    distance = Min(distance, (LC->RADIUS*LC->RADIUS));
+
+    // then check direction to boundary
+    if (DistanceToCenter.ScalarProduct(&result->NormalVector) > MYEPSILON) {
+      Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
+      return -distance;
+    } else {
+      Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl;
+      return +distance;
+    }
+  }
+};
 
 /** Gets all points connected to the provided point by triangulation lines.
@@ -3126,8 +3441,8 @@
  * @return set of the all points linked to the provided one
  */
-set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
-{
-	Info FunctionInfo(__func__);
-  set<TesselPoint*> *connectedPoints = new set<TesselPoint*>;
+TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
+{
+	Info FunctionInfo(__func__);
+	TesselPointSet *connectedPoints = new TesselPointSet;
   class BoundaryPointSet *ReferencePoint = NULL;
   TesselPoint* current;
@@ -3170,5 +3485,5 @@
   }
 
-  if (connectedPoints->size() == 0) { // if have not found any points
+  if (connectedPoints->empty()) { // if have not found any points
     eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl;
     return NULL;
@@ -3191,14 +3506,15 @@
  * @return list of the all points linked to the provided one
  */
-list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
 {
 	Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<TesselPoint*> *connectedCircle = new list<TesselPoint*>;
-  Vector center;
+  TesselPointList *connectedCircle = new TesselPointList;
   Vector PlaneNormal;
   Vector AngleZero;
   Vector OrthogonalVector;
   Vector helper;
+  const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL};
+  TriangleList *triangles = NULL;
 
   if (SetOfNeighbours == NULL) {
@@ -3209,16 +3525,15 @@
 
   // calculate central point
-  for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
-    center.AddVector((*TesselRunner)->node);
-  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
-  //  << "; scale factor " << 1.0/connectedPoints.size();
-  center.Scale(1.0/SetOfNeighbours->size());
-  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
-
-  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
-  PlaneNormal.CopyVector(Point->node);
-  PlaneNormal.SubtractVector(&center);
+  triangles = FindTriangles(TrianglePoints);
+  if ((triangles != NULL) && (!triangles->empty())) {
+    for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+      PlaneNormal.AddVector(&(*Runner)->NormalVector);
+  } else {
+    eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
+    performCriticalExit();
+  }
+  PlaneNormal.Scale(1.0/triangles->size());
+  Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl;
   PlaneNormal.Normalize();
-  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
 
   // construct one orthogonal vector
@@ -3246,5 +3561,5 @@
 
   // go through all connected points and calculate angle
-  for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
     helper.CopyVector((*listRunner)->node);
     helper.SubtractVector(Point->node);
@@ -3262,4 +3577,116 @@
 }
 
+/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
+ * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
+ * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
+ * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
+ * triangle we are looking for.
+ *
+ * @param *SetOfNeighbours all points for which the angle should be calculated
+ * @param *Point of which get all connected points
+ * @param *Reference Reference vector for zero angle or NULL for no preference
+ * @return list of the all points linked to the provided one
+ */
+TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const
+{
+  Info FunctionInfo(__func__);
+  map<double, TesselPoint*> anglesOfPoints;
+  TesselPointList *connectedCircle = new TesselPointList;
+  Vector center;
+  Vector PlaneNormal;
+  Vector AngleZero;
+  Vector OrthogonalVector;
+  Vector helper;
+
+  if (SetOfNeighbours == NULL) {
+    eLog() << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
+    return NULL;
+  }
+
+  // check whether there's something to do
+  if (SetOfNeighbours->size() < 3) {
+    for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
+      connectedCircle->push_back(*TesselRunner);
+    return connectedCircle;
+  }
+
+  Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl;
+  // calculate central point
+
+  TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
+  TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
+  TesselB++;
+  TesselC++;
+  TesselC++;
+  int counter = 0;
+  while (TesselC != SetOfNeighbours->end()) {
+    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
+    counter++;
+    TesselA++;
+    TesselB++;
+    TesselC++;
+    PlaneNormal.AddVector(&helper);
+  }
+  //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
+  //  << "; scale factor " << counter;
+  PlaneNormal.Scale(1.0/(double)counter);
+//  Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
+//
+//  // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
+//  PlaneNormal.CopyVector(Point->node);
+//  PlaneNormal.SubtractVector(&center);
+//  PlaneNormal.Normalize();
+  Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;
+
+  // construct one orthogonal vector
+  if (Reference != NULL) {
+    AngleZero.CopyVector(Reference);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+  }
+  if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {
+    Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;
+    AngleZero.CopyVector((*SetOfNeighbours->begin())->node);
+    AngleZero.SubtractVector(Point->node);
+    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    if (AngleZero.NormSquared() < MYEPSILON) {
+      eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl;
+      performCriticalExit();
+    }
+  }
+  Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
+  if (AngleZero.NormSquared() > MYEPSILON)
+    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+  else
+    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+  Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
+
+  // go through all connected points and calculate angle
+  pair <map<double, TesselPoint*>::iterator, bool > InserterTest;
+  for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
+    helper.CopyVector((*listRunner)->node);
+    helper.SubtractVector(Point->node);
+    helper.ProjectOntoPlane(&PlaneNormal);
+    double angle = GetAngle(helper, AngleZero, OrthogonalVector);
+    if (angle > M_PI) // the correction is of no use here (and not desired)
+      angle = 2.*M_PI - angle;
+    Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl;
+    InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner)));
+    if (!InserterTest.second) {
+      eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl;
+      performCriticalExit();
+    }
+  }
+
+  for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
+    connectedCircle->push_back(AngleRunner->second);
+  }
+
+  return connectedCircle;
+}
+
 /** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
  *
@@ -3268,10 +3695,10 @@
  * @return list of the all points linked to the provided one
  */
-list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
+list< TesselPointList *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
 {
 	Info FunctionInfo(__func__);
   map<double, TesselPoint*> anglesOfPoints;
-  list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
+  list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
   Vector center;
   Vector PlaneNormal;
@@ -3310,5 +3737,5 @@
       } else if (!LineRunner->second) {
         LineRunner->second = true;
-        connectedPath = new list<TesselPoint*>;
+        connectedPath = new TesselPointList;
         triangle = NULL;
         CurrentLine = runner->second;
@@ -3384,18 +3811,18 @@
  * @return list of the closed paths
  */
-list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
-{
-	Info FunctionInfo(__func__);
-  list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point);
-  list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>;
-  list<TesselPoint*> *connectedPath = NULL;
-  list<TesselPoint*> *newPath = NULL;
+list<TesselPointList *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
+{
+	Info FunctionInfo(__func__);
+  list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
+  list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>;
+  TesselPointList *connectedPath = NULL;
+  TesselPointList *newPath = NULL;
   int count = 0;
 
 
-  list<TesselPoint*>::iterator CircleRunner;
-  list<TesselPoint*>::iterator CircleStart;
-
-  for(list<list<TesselPoint*> *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
+  TesselPointList::iterator CircleRunner;
+  TesselPointList::iterator CircleStart;
+
+  for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
     connectedPath = *ListRunner;
 
@@ -3406,11 +3833,11 @@
 
     // go through list, look for reappearance of starting Point and create list
-    list<TesselPoint*>::iterator Marker = CircleStart;
+    TesselPointList::iterator Marker = CircleStart;
     for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
       if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
         // we have a closed circle from Marker to new Marker
         Log() << Verbose(1) << count+1 << ". closed path consists of: ";
-        newPath = new list<TesselPoint*>;
-        list<TesselPoint*>::iterator CircleSprinter = Marker;
+        newPath = new TesselPointList;
+        TesselPointList::iterator CircleSprinter = Marker;
         for (; CircleSprinter != CircleRunner; CircleSprinter++) {
           newPath->push_back(*CircleSprinter);
@@ -3446,8 +3873,8 @@
  * \return pointer to allocated list of triangles
  */
-set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
-{
-	Info FunctionInfo(__func__);
-  set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>;
+TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
+{
+	Info FunctionInfo(__func__);
+	TriangleSet *connectedTriangles = new TriangleSet;
 
   if (Point == NULL) {
@@ -3498,16 +3925,16 @@
   }
 
-  list<list<TesselPoint*> *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  list<TesselPoint*> *connectedPath = NULL;
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
+  TesselPointList *connectedPath = NULL;
 
   // gather all triangles
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
     count+=LineRunner->second->triangles.size();
-  map<class BoundaryTriangleSet *, int> Candidates;
+  TriangleMap Candidates;
   for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
     line = LineRunner->second;
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      Candidates.insert( pair<class BoundaryTriangleSet *, int> (triangle, triangle->Nr) );
+      Candidates.insert( TrianglePair (triangle->Nr, triangle) );
     }
   }
@@ -3516,16 +3943,16 @@
   count=0;
   NormalVector.Zero();
-  for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;
-    NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward
-    RemoveTesselationTriangle(Runner->first);
+  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
+    Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl;
+    NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward
+    RemoveTesselationTriangle(Runner->second);
     count++;
   }
   Log() << Verbose(1) << count << " triangles were removed." << endl;
 
-  list<list<TesselPoint*> *>::iterator ListAdvance = ListOfClosedPaths->begin();
-  list<list<TesselPoint*> *>::iterator ListRunner = ListAdvance;
-  map<class BoundaryTriangleSet *, int>::iterator NumberRunner = Candidates.begin();
-  list<TesselPoint*>::iterator StartNode, MiddleNode, EndNode;
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+  TriangleMap::iterator NumberRunner = Candidates.begin();
+  TesselPointList::iterator StartNode, MiddleNode, EndNode;
   double angle;
   double smallestangle;
@@ -3541,5 +3968,5 @@
 
       // re-create all triangles by going through connected points list
-      list<class BoundaryLineSet *> NewLines;
+      LineList NewLines;
       for (;!connectedPath->empty();) {
         // search middle node with widest angle to next neighbours
@@ -3647,10 +4074,10 @@
       // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
       if (NewLines.size() > 1) {
-        list<class BoundaryLineSet *>::iterator Candidate;
+        LineList::iterator Candidate;
         class BoundaryLineSet *OtherBase = NULL;
         double tmp, maxgain;
         do {
           maxgain = 0;
-          for(list<class BoundaryLineSet *>::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
+          for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
             tmp = PickFarthestofTwoBaselines(*Runner);
             if (maxgain < tmp) {
@@ -3694,47 +4121,98 @@
  * Finds triangles belonging to the three provided points.
  *
- * @param *Points[3] list, is expected to contain three points
+ * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
  *
  * @return triangles which belong to the provided points, will be empty if there are none,
  *         will usually be one, in case of degeneration, there will be two
  */
-list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
-{
-	Info FunctionInfo(__func__);
-  list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>;
+TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
+{
+	Info FunctionInfo(__func__);
+	TriangleList *result = new TriangleList;
   LineMap::const_iterator FindLine;
   TriangleMap::const_iterator FindTriangle;
   class BoundaryPointSet *TrianglePoints[3];
+  size_t NoOfWildcards = 0;
 
   for (int i = 0; i < 3; i++) {
-    PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
-    if (FindPoint != PointsOnBoundary.end()) {
-      TrianglePoints[i] = FindPoint->second;
+    if (Points[i] == NULL) {
+      NoOfWildcards++;
+      TrianglePoints[i] = NULL;
     } else {
-      TrianglePoints[i] = NULL;
-    }
-  }
-
-  // checks lines between the points in the Points for their adjacent triangles
-  for (int i = 0; i < 3; i++) {
-    if (TrianglePoints[i] != NULL) {
-      for (int j = i+1; j < 3; j++) {
-        if (TrianglePoints[j] != NULL) {
-          for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
-              (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
-              FindLine++) {
-            for (FindTriangle = FindLine->second->triangles.begin();
-                FindTriangle != FindLine->second->triangles.end();
-                FindTriangle++) {
-              if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
-                result->push_back(FindTriangle->second);
+      PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
+      if (FindPoint != PointsOnBoundary.end()) {
+        TrianglePoints[i] = FindPoint->second;
+      } else {
+        TrianglePoints[i] = NULL;
+      }
+    }
+  }
+
+  switch (NoOfWildcards) {
+    case 0: // checks lines between the points in the Points for their adjacent triangles
+      for (int i = 0; i < 3; i++) {
+        if (TrianglePoints[i] != NULL) {
+          for (int j = i+1; j < 3; j++) {
+            if (TrianglePoints[j] != NULL) {
+              for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
+                  (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr);
+                  FindLine++) {
+                for (FindTriangle = FindLine->second->triangles.begin();
+                    FindTriangle != FindLine->second->triangles.end();
+                    FindTriangle++) {
+                  if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+                    result->push_back(FindTriangle->second);
+                  }
+                }
               }
+              // Is it sufficient to consider one of the triangle lines for this.
+              return result;
             }
           }
-          // Is it sufficient to consider one of the triangle lines for this.
-          return result;
         }
       }
-    }
+      break;
+    case 1: // copy all triangles of the respective line
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] == NULL)
+          break;
+      for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap!
+          (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr);
+          FindLine++) {
+        for (FindTriangle = FindLine->second->triangles.begin();
+            FindTriangle != FindLine->second->triangles.end();
+            FindTriangle++) {
+          if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
+            result->push_back(FindTriangle->second);
+          }
+        }
+      }
+      break;
+    }
+    case 2: // copy all triangles of the respective point
+    {
+      int i=0;
+      for (; i < 3; i++)
+        if (TrianglePoints[i] != NULL)
+          break;
+      for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
+        for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
+          result->push_back(triangle->second);
+      result->sort();
+      result->unique();
+      break;
+    }
+    case 3: // copy all triangles
+    {
+      for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
+        result->push_back(triangle->second);
+      break;
+    }
+    default:
+      eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl;
+      performCriticalExit();
+      break;
   }
 
@@ -3779,9 +4257,9 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedLines()
+IndexToIndex * Tesselation::FindAllDegeneratedLines()
 {
 	Info FunctionInfo(__func__);
 	UniqueLines AllLines;
-  map<int, int> * DegeneratedLines = new map<int, int>;
+  IndexToIndex * DegeneratedLines = new IndexToIndex;
 
   // sanity check
@@ -3804,5 +4282,5 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
     const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
@@ -3823,9 +4301,9 @@
  *         in the list, once as key and once as value
  */
-map<int, int> * Tesselation::FindAllDegeneratedTriangles()
-{
-	Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedLines = FindAllDegeneratedLines();
-  map<int, int> * DegeneratedTriangles = new map<int, int>;
+IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
+{
+	Info FunctionInfo(__func__);
+  IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
+  IndexToIndex * DegeneratedTriangles = new IndexToIndex;
 
   TriangleMap::iterator TriangleRunner1, TriangleRunner2;
@@ -3833,5 +4311,5 @@
   class BoundaryLineSet *line1 = NULL, *line2 = NULL;
 
-  for (map<int, int>::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
+  for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
     // run over both lines' triangles
     Liner = LinesOnBoundary.find(LineRunner->first);
@@ -3854,5 +4332,5 @@
 
   Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
@@ -3868,10 +4346,10 @@
 {
 	Info FunctionInfo(__func__);
-  map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
   TriangleMap::iterator finder;
   BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
   int count  = 0;
 
-  for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
+  for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin();
     TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner
   ) {
@@ -3961,5 +4439,5 @@
   // find nearest boundary point
   class TesselPoint *BackupPoint = NULL;
-  class TesselPoint *NearestPoint = FindClosestPoint(point->node, BackupPoint, LC);
+  class TesselPoint *NearestPoint = FindClosestTesselPoint(point->node, BackupPoint, LC);
   class BoundaryPointSet *NearestBoundaryPoint = NULL;
   PointMap::iterator PointRunner;
@@ -4128,4 +4606,5 @@
 
   /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
+  IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
   set < BoundaryPointSet *> EndpointCandidateList;
   pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;
@@ -4138,7 +4617,11 @@
     for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
       for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
-        TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
-        if (TriangleInsertionTester.second)
-          Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
+          TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );
+          if (TriangleInsertionTester.second)
+            Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;
+        } else {
+          Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;
+        }
       }
     // check whether there are two that are parallel
@@ -4213,4 +4696,20 @@
     /// 4a. Gather all triangles of this polygon
     TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
+
+    // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.
+    if (T->size() == 2) {
+      Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;
+      delete(T);
+      continue;
+    }
+
+    // check whether number is even
+    // If this case occurs, we have to think about it!
+    // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
+    // connections to either polygon ...
+    if (T->size() % 2 != 0) {
+      eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;
+      performCriticalExit();
+    }
 
     TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
@@ -4259,7 +4758,7 @@
   }
 
-  map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
+  IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
   Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;
-  map<int,int>::iterator it;
+  IndexToIndex::iterator it;
   for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
       Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;
Index: src/tesselation.hpp
===================================================================
--- src/tesselation.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/tesselation.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -52,4 +52,6 @@
 // ======================================================= some template functions =========================================
 
+#define IndexToIndex map <int, int>
+
 #define PointMap map < int, class BoundaryPointSet * >
 #define PointSet set < class BoundaryPointSet * >
@@ -76,4 +78,7 @@
 #define PolygonSet set < class BoundaryPolygonSet * >
 #define PolygonList list < class BoundaryPolygonSet * >
+
+#define DistanceToPointMap multimap <double, class BoundaryPointSet * >
+#define DistanceToPointPair pair <double, class BoundaryPointSet * >
 
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
@@ -101,8 +106,8 @@
   public:
     BoundaryPointSet();
-    BoundaryPointSet(TesselPoint * Walker);
+    BoundaryPointSet(TesselPoint * const Walker);
     ~BoundaryPointSet();
 
-    void AddLine(class BoundaryLineSet *line);
+    void AddLine(BoundaryLineSet * const line);
 
     LineMap lines;
@@ -120,12 +125,13 @@
   public:
     BoundaryLineSet();
-    BoundaryLineSet(class BoundaryPointSet *Point[2], const int number);
+    BoundaryLineSet(BoundaryPointSet * const Point[2], const int number);
+    BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number);
     ~BoundaryLineSet();
 
-    void AddTriangle(class BoundaryTriangleSet *triangle);
-    bool IsConnectedTo(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-    bool CheckConvexityCriterion();
-    class BoundaryPointSet *GetOtherEndpoint(class BoundaryPointSet *);
+    void AddTriangle(BoundaryTriangleSet * const triangle);
+    bool IsConnectedTo(const BoundaryLineSet * const line) const;
+    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
+    bool CheckConvexityCriterion() const;
+    class BoundaryPointSet *GetOtherEndpoint(const BoundaryPointSet * const point) const;
 
     class BoundaryPointSet *endpoints[2];
@@ -142,20 +148,22 @@
   public:
     BoundaryTriangleSet();
-    BoundaryTriangleSet(class BoundaryLineSet *line[3], int number);
+    BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number);
     ~BoundaryTriangleSet();
 
-    void GetNormalVector(Vector &NormalVector);
-    void GetCenter(Vector *center);
-    bool GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection);
-    bool ContainsBoundaryLine(class BoundaryLineSet *line);
-    bool ContainsBoundaryPoint(class BoundaryPointSet *point);
-    bool ContainsBoundaryPoint(class TesselPoint *point);
-    class BoundaryPointSet *GetThirdEndpoint(class BoundaryLineSet *line);
-    bool IsPresentTupel(class BoundaryPointSet *Points[3]);
-    bool IsPresentTupel(class BoundaryTriangleSet *T);
+    void GetNormalVector(const Vector &NormalVector);
+    void GetCenter(Vector * const center) const;
+    bool GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const;
+    double GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const;
+    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
+    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
+    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
+    class BoundaryPointSet *GetThirdEndpoint(const BoundaryLineSet * const line) const;
+    bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
+    bool IsPresentTupel(const BoundaryTriangleSet * const T) const;
 
     class BoundaryPointSet *endpoints[3];
     class BoundaryLineSet *lines[3];
     Vector NormalVector;
+    Vector SphereCenter;
     int Nr;
 };
@@ -225,4 +233,5 @@
   virtual TesselPoint *GetPoint() const { return NULL; };
   virtual TesselPoint *GetTerminalPoint() const { return NULL; };
+  virtual int GetMaxId() const { return 0; };
   virtual void GoToNext() const {};
   virtual void GoToPrevious() const {};
@@ -300,4 +309,5 @@
     list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
     list<TesselPoint*> * GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
+    list<TesselPoint*> * GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
     class BoundaryPointSet *GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
     list<BoundaryTriangleSet*> *FindTriangles(const TesselPoint* const Points[3]) const;
@@ -305,6 +315,10 @@
     class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
     bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
-    bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const;
+    double GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const;
     bool AddBoundaryPoint(TesselPoint * Walker, const int n);
+    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;
+    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
+    BoundaryTriangleSet* FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
 
     // print for debugging
Index: src/tesselationhelpers.cpp
===================================================================
--- src/tesselationhelpers.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/tesselationhelpers.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -143,5 +143,5 @@
   if (fabs(HalfplaneIndicator) < MYEPSILON)
     {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 and AlternativeIndicator >0) or (TempNormal.ScalarProduct(AlternativeDirection) >0 and AlternativeIndicator <0))
+      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
         {
           TempNormal.Scale(-1);
@@ -150,5 +150,5 @@
   else
     {
-      if (TempNormal.ScalarProduct(Direction)<0 && HalfplaneIndicator >0 || TempNormal.ScalarProduct(Direction)>0 && HalfplaneIndicator<0)
+      if (((TempNormal.ScalarProduct(Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(Direction)>0) && (HalfplaneIndicator<0)))
         {
           TempNormal.Scale(-1);
@@ -226,6 +226,12 @@
   Vector helper;
   double radius, alpha;
-
-  helper.CopyVector(&NewSphereCenter);
+  Vector RelativeOldSphereCenter;
+  Vector RelativeNewSphereCenter;
+
+  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
+  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
+  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
+  helper.CopyVector(&RelativeNewSphereCenter);
   // test whether new center is on the parameter circle's plane
   if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
@@ -233,21 +239,21 @@
     helper.ProjectOntoPlane(&CirclePlaneNormal);
   }
-  radius = helper.ScalarProduct(&helper);
+  radius = helper.NormSquared();
   // test whether the new center vector has length of CircleRadius
   if (fabs(radius - CircleRadius) > HULLEPSILON)
     eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
-  alpha = helper.Angle(&OldSphereCenter);
+  alpha = helper.Angle(&RelativeOldSphereCenter);
   // 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;
-  //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&OldSphereCenter);
+  Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
+  radius = helper.Distance(&RelativeOldSphereCenter);
   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)) {
-    //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
+    Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     return alpha;
   } else {
-    //Log() << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
+    Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl;
     return 2.*M_PI;
   }
@@ -552,5 +558,5 @@
  * @return point which is second closest to the provided one
  */
-TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
+TesselPoint* FindSecondClosestTesselPoint(const Vector* Point, const LinkedCell* const LC)
 {
 	Info FunctionInfo(__func__);
@@ -607,5 +613,5 @@
  * @return point which is closest to the provided one, NULL if none found
  */
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
+TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
 {
 	Info FunctionInfo(__func__);
@@ -633,5 +639,5 @@
             helper.CopyVector(Point);
             helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
+            double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
               secondDistance = distance;
@@ -860,6 +866,5 @@
     }
     *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
-    int i=0;
-    for (cloud->GoToFirst(); !cloud->IsEnd(); cloud->GoToNext(), i++);
+    int i=cloud->GetMaxId();
     int *LookupList = new int[i];
     for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
Index: src/tesselationhelpers.hpp
===================================================================
--- src/tesselationhelpers.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/tesselationhelpers.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -59,6 +59,6 @@
 bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3]);
 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation *candidate2);
-TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
-TesselPoint* FindSecondClosestPoint(const Vector*, const LinkedCell* const LC);
+TesselPoint* FindClosestTesselPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC);
+TesselPoint* FindSecondClosestTesselPoint(const Vector*, const LinkedCell* const LC);
 Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase);
 
Index: src/unittests/AnalysisCorrelationToPointUnitTest.cpp
===================================================================
--- src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/AnalysisCorrelationToPointUnitTest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -11,4 +11,6 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
 
 #include "analysis_correlation.hpp"
Index: src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
===================================================================
--- src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -11,4 +11,6 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
 
 #include "analysis_correlation.hpp"
Index: src/unittests/AnalysisPairCorrelationUnitTest.cpp
===================================================================
--- src/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/AnalysisPairCorrelationUnitTest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -11,4 +11,6 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
 
 #include "analysis_correlation.hpp"
Index: src/unittests/Makefile.am
===================================================================
--- src/unittests/Makefile.am	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/Makefile.am	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -4,46 +4,88 @@
 AM_CXXFLAGS = $(CPPUNIT_CFLAGS)
 
-TESTS = ActOnAllUnitTest AnalysisBondsUnitTests AnalysisCorrelationToPointUnitTest AnalysisCorrelationToSurfaceUnitTest AnalysisPairCorrelationUnitTest BondGraphUnitTest InfoUnitTest ListOfBondsUnitTest LogUnitTest MemoryUsageObserverUnitTest MemoryAllocatorUnitTest StackClassUnitTest VectorUnitTest
+TESTS = \
+  ActOnAllUnitTest \
+  AnalysisBondsUnitTests \
+  AnalysisCorrelationToPointUnitTest \
+  AnalysisCorrelationToSurfaceUnitTest \
+  AnalysisPairCorrelationUnitTest \
+  BondGraphUnitTest \
+  GSLMatrixSymmetricUnitTest \
+  GSLMatrixUnitTest \
+  GSLVectorUnitTest \
+  InfoUnitTest \
+  LinearSystemOfEquationsUnitTest \
+  ListOfBondsUnitTest \
+  LogUnitTest \
+  MemoryUsageObserverUnitTest \
+  MemoryAllocatorUnitTest \
+  StackClassUnitTest \
+  TesselationUnitTest \
+  Tesselation_BoundaryTriangleUnitTest \
+  Tesselation_InOutsideUnitTest \
+  VectorUnitTest
+  
 check_PROGRAMS = $(TESTS)
 noinst_PROGRAMS = $(TESTS)
 
 ActOnAllUnitTest_SOURCES = ../test/ActOnAllTest.hpp ActOnAllUnitTest.cpp ActOnAllUnitTest.hpp
-ActOnAllUnitTest_LDADD = ../libmolecuilder.a
+ActOnAllUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 AnalysisBondsUnitTests_SOURCES = analysisbondsunittest.cpp analysisbondsunittest.hpp
-AnalysisBondsUnitTests_LDADD = ../libmolecuilder.a
+AnalysisBondsUnitTests_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 AnalysisCorrelationToPointUnitTest_SOURCES = analysis_correlation.hpp AnalysisCorrelationToPointUnitTest.cpp AnalysisCorrelationToPointUnitTest.hpp
-AnalysisCorrelationToPointUnitTest_LDADD = ../libmolecuilder.a
+AnalysisCorrelationToPointUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 AnalysisCorrelationToSurfaceUnitTest_SOURCES = analysis_correlation.hpp AnalysisCorrelationToSurfaceUnitTest.cpp AnalysisCorrelationToSurfaceUnitTest.hpp
-AnalysisCorrelationToSurfaceUnitTest_LDADD = ../libmolecuilder.a
+AnalysisCorrelationToSurfaceUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 AnalysisPairCorrelationUnitTest_SOURCES = analysis_correlation.hpp AnalysisPairCorrelationUnitTest.cpp AnalysisPairCorrelationUnitTest.hpp
-AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a
+AnalysisPairCorrelationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 BondGraphUnitTest_SOURCES = bondgraphunittest.cpp bondgraphunittest.hpp
-BondGraphUnitTest_LDADD = ../libmolecuilder.a
+BondGraphUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+
+GSLMatrixSymmetricUnitTest_SOURCES = gslmatrixsymmetricunittest.cpp gslmatrixsymmetricunittest.hpp
+GSLMatrixSymmetricUnitTest_LDADD = ../libgslwrapper.a
+
+GSLMatrixUnitTest_SOURCES = gslmatrixunittest.cpp gslmatrixunittest.hpp
+GSLMatrixUnitTest_LDADD = ../libgslwrapper.a
+
+GSLVectorUnitTest_SOURCES = gslvectorunittest.cpp gslvectorunittest.hpp
+GSLVectorUnitTest_LDADD = ../libgslwrapper.a
 
 InfoUnitTest_SOURCES = infounittest.cpp infounittest.hpp
-InfoUnitTest_LDADD = ../libmolecuilder.a
+InfoUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+
+LinearSystemOfEquationsUnitTest_SOURCES = linearsystemofequationsunittest.cpp linearsystemofequationsunittest.hpp
+LinearSystemOfEquationsUnitTest_LDADD = ../libgslwrapper.a ../libmolecuilder.a
 
 ListOfBondsUnitTest_SOURCES = listofbondsunittest.cpp listofbondsunittest.hpp
-ListOfBondsUnitTest_LDADD = ../libmolecuilder.a
+ListOfBondsUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 LogUnitTest_SOURCES = logunittest.cpp logunittest.hpp 
-LogUnitTest_LDADD = ../libmolecuilder.a
+LogUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 MemoryAllocatorUnitTest_SOURCES = memoryallocatorunittest.cpp memoryallocatorunittest.hpp 
-MemoryAllocatorUnitTest_LDADD = ../libmolecuilder.a
+MemoryAllocatorUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 MemoryUsageObserverUnitTest_SOURCES = memoryusageobserverunittest.cpp memoryusageobserverunittest.hpp 
-MemoryUsageObserverUnitTest_LDADD = ../libmolecuilder.a
+MemoryUsageObserverUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 StackClassUnitTest_SOURCES = stackclassunittest.cpp stackclassunittest.hpp
-StackClassUnitTest_LDADD = ../libmolecuilder.a
+StackClassUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+
+TesselationUnitTest_SOURCES = tesselationunittest.cpp tesselationunittest.hpp
+TesselationUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+
+Tesselation_BoundaryTriangleUnitTest_SOURCES = tesselation_boundarytriangleunittest.cpp tesselation_boundarytriangleunittest.hpp
+Tesselation_BoundaryTriangleUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
+
+Tesselation_InOutsideUnitTest_SOURCES = tesselation_insideoutsideunittest.cpp tesselation_insideoutsideunittest.hpp
+Tesselation_InOutsideUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 VectorUnitTest_SOURCES = vectorunittest.cpp vectorunittest.hpp 
-VectorUnitTest_LDADD = ../libmolecuilder.a
+VectorUnitTest_LDADD = ../libmolecuilder.a ../libgslwrapper.a
 
 
Index: src/unittests/analysisbondsunittest.cpp
===================================================================
--- src/unittests/analysisbondsunittest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/analysisbondsunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -14,4 +14,5 @@
 #include <iostream>
 #include <stdio.h>
+#include <cstring>
 
 #include "analysis_bonds.hpp"
Index: src/unittests/bondgraphunittest.cpp
===================================================================
--- src/unittests/bondgraphunittest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/bondgraphunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -14,4 +14,5 @@
 #include <iostream>
 #include <stdio.h>
+#include <cstring>
 
 #include "atom.hpp"
Index: src/unittests/gslmatrixsymmetricunittest.cpp
===================================================================
--- src/unittests/gslmatrixsymmetricunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslmatrixsymmetricunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,288 @@
+/*
+ * gslmatrixunittest.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "gslmatrixsymmetricunittest.hpp"
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( GSLMatrixSymmetricTest );
+
+
+void GSLMatrixSymmetricTest::setUp()
+{
+  m = new GSLMatrix(3,3);
+};
+
+void GSLMatrixSymmetricTest::tearDown()
+{
+  delete(m);
+};
+
+/** Unit Test for accessing matrix elements.
+ *
+ */
+void GSLMatrixSymmetricTest::AccessTest()
+{
+  // check whether all elements are initially zero
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 0., m->Get(i,j) );
+
+  // set
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i*3+j );
+
+  // and check
+  double *ptr = NULL;
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++) {
+      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), m->Get(i,j) );
+      ptr = m->Pointer(i,j);
+      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), *ptr );
+    }
+
+  // assignment
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i*3+j );
+  GSLMatrix *dest = new GSLMatrix(3,3);
+  *dest = *m;
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( dest->Get(i,j), m->Get(i,j) );
+  delete(dest);
+
+  // out of bounds
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(4,2) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(2,17) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(1024,140040) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(-1,0) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(0,-1) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(-1,-1) );
+};
+
+/** Unit Test for initializating matrices.
+ *
+ */
+void GSLMatrixSymmetricTest::InitializationTest()
+{
+  // set zero
+  m->SetZero();
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 0., m->Get(i,j) );
+
+  // set all
+  m->SetAll(1.5);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 1.5, m->Get(i,j) );
+
+  // set basis
+  m->SetIdentity();
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->Get(i,j) );
+
+  // set from array
+  double array[] = { 1., 0., 0.,
+                     0., 1., 0.,
+                     0., 0., 1. };
+  m->SetFromDoubleArray(array);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->Get(i,j) );
+
+};
+
+/** Unit Test for copying matrices.
+ *
+ */
+void GSLMatrixSymmetricTest::CopyTest()
+{
+  // set basis
+  GSLMatrix *dest = NULL;
+  for (int i=0;i<3;i++) {
+    m->SetAll(i);
+    dest = new GSLMatrix(m);
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( m->Get(i,j) , dest->Get(i,j) );
+
+    delete(dest);
+  }
+};
+
+/** Unit Test for exchanging rows and columns.
+ *
+ */
+void GSLMatrixSymmetricTest::ExchangeTest()
+{
+  // set to 1,1,1,2, ...
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i+1 );
+
+  // swap such that nothing happens
+  m->SwapColumns(1,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->Get(i,j) );
+
+  // swap two rows
+  m->SwapRows(1,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      switch (j) {
+        case 0:
+          CPPUNIT_ASSERT_EQUAL( 1., m->Get(j,i) );
+          break;
+        case 1:
+          CPPUNIT_ASSERT_EQUAL( 3., m->Get(j,i) );
+          break;
+        case 2:
+          CPPUNIT_ASSERT_EQUAL( 2., m->Get(j,i) );
+          break;
+        default:
+          CPPUNIT_ASSERT_EQUAL( -1., m->Get(i,j) );
+      }
+  // check that op is reversable
+  m->SwapRows(1,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->Get(i,j) );
+
+  // set to 1,2,3,1, ...
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, j+1. );
+
+  // swap such that nothing happens
+  m->SwapRows(0,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->Get(i,j) );
+
+  // swap two columns
+  m->SwapColumns(0,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      switch (j) {
+        case 0:
+          CPPUNIT_ASSERT_EQUAL( 3., m->Get(i,j) );
+          break;
+        case 1:
+          CPPUNIT_ASSERT_EQUAL( 2., m->Get(i,j) );
+          break;
+        case 2:
+          CPPUNIT_ASSERT_EQUAL( 1., m->Get(i,j) );
+          break;
+        default:
+          CPPUNIT_ASSERT_EQUAL( -1., m->Get(i,j) );
+      }
+  // check that op is reversable
+  m->SwapColumns(0,2);
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->Get(i,j) );
+
+
+  // set to 1,2,3,4, ...
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, 3*i+j+1 );
+  // transpose
+  m->Transpose();
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(j,i) );
+  // second transpose
+  m->Transpose();
+  for (int i=0;i<3;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(i,j) );
+};
+
+/** Unit Test for matrix properties.
+ *
+ */
+void GSLMatrixSymmetricTest::PropertiesTest()
+{
+  // is zero
+  m->SetZero();
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
+
+  // is positive
+  m->SetAll(0.5);
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
+
+  // is negative
+  m->SetAll(-0.1);
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNonNegative() );
+
+  // is positive definite
+  double array[] = { 1., 0., 0.,
+                     0., 1., 1.,
+                     0., 0., 1. };
+  m->SetFromDoubleArray(array);
+  CPPUNIT_ASSERT_EQUAL( true, m->IsPositiveDefinite() );
+
+  //determinant
+  m->SetIdentity();
+  CPPUNIT_ASSERT_EQUAL( 1., m->Determinant() );
+
+  m->SetZero();
+  CPPUNIT_ASSERT_EQUAL( 0., m->Determinant() );
+
+  m->Set( 0, 0, 1.);
+  m->Set( 1, 1, 1.);
+  m->Set( 2, 1, 1.);
+  CPPUNIT_ASSERT_EQUAL( 0., m->Determinant() );
+  
+  double array2[] = { 2., 0., 1.,
+                     -3., 1., 1.,
+                     1., 5.5, 1. };
+  m->SetFromDoubleArray(array2);
+  CPPUNIT_ASSERT_EQUAL( -26.5, m->Determinant() );
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/gslmatrixsymmetricunittest.hpp
===================================================================
--- src/unittests/gslmatrixsymmetricunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslmatrixsymmetricunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,42 @@
+/*
+ * gslmatrixunittest.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef GSLMATRIXUNITTEST_HPP_
+#define GSLMATRIXUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "gslmatrix.hpp"
+
+/********************************************** Test classes **************************************/
+
+class GSLMatrixSymmetricTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( GSLMatrixSymmetricTest) ;
+    CPPUNIT_TEST (AccessTest );
+    CPPUNIT_TEST (InitializationTest );
+    CPPUNIT_TEST (CopyTest );
+    CPPUNIT_TEST (ExchangeTest );
+    CPPUNIT_TEST (PropertiesTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+    void setUp();
+    void tearDown();
+
+    void AccessTest();
+    void InitializationTest();
+    void CopyTest();
+    void ExchangeTest();
+    void PropertiesTest();
+
+private:
+    GSLMatrix *m;
+};
+
+
+#endif /* GSLMATRIXUNITTEST_HPP_ */
Index: src/unittests/gslmatrixunittest.cpp
===================================================================
--- src/unittests/gslmatrixunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslmatrixunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,270 @@
+/*
+ * gslmatrixunittest.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "gslmatrixunittest.hpp"
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( GSLMatrixTest );
+
+
+void GSLMatrixTest::setUp()
+{
+  m = new GSLMatrix(4,3);
+};
+
+void GSLMatrixTest::tearDown()
+{
+  delete(m);
+};
+
+/** Unit Test for accessing matrix elements.
+ *
+ */
+void GSLMatrixTest::AccessTest()
+{
+  // check whether all elements are initially zero
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 0., m->Get(i,j) );
+
+  // set
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i*3+j );
+
+  // and check
+  double *ptr = NULL;
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++) {
+      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), m->Get(i,j) );
+      ptr = m->Pointer(i,j);
+      CPPUNIT_ASSERT_EQUAL( (double)(i*3+j), *ptr );
+    }
+
+  // assignment
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i*3+j );
+  GSLMatrix *dest = new GSLMatrix(4,3);
+  *dest = *m;
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( dest->Get(i,j), m->Get(i,j) );
+  delete(dest);
+
+  // out of bounds
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(5,2) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(2,17) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(1024,140040) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(-1,0) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(0,-1) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(-1,-1) );
+};
+
+/** Unit Test for initializating matrices.
+ *
+ */
+void GSLMatrixTest::InitializationTest()
+{
+  // set zero
+  m->SetZero();
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 0., m->Get(i,j) );
+
+  // set all
+  m->SetAll(1.5);
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( 1.5, m->Get(i,j) );
+
+  // set basis
+  m->SetIdentity();
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->Get(i,j) );
+
+  // set from array
+  double array[] = { 1., 0., 0.,
+                     0., 1., 0.,
+                     0., 0., 1.,
+                     0., 0., 0. };
+  m->SetFromDoubleArray(array);
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , m->Get(i,j) );
+
+};
+
+/** Unit Test for copying matrices.
+ *
+ */
+void GSLMatrixTest::CopyTest()
+{
+  // set basis
+  GSLMatrix *dest = NULL;
+  for (int i=0;i<4;i++) {
+    m->SetAll(i);
+    dest = new GSLMatrix(m);
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( m->Get(i,j) , dest->Get(i,j) );
+
+    delete(dest);
+  }
+};
+
+/** Unit Test for exchanging rows and columns.
+ *
+ */
+void GSLMatrixTest::ExchangeTest()
+{
+  // set to 1,1,1,2, ...
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, i+1 );
+
+  // swap such that nothing happens
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(1,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->Get(i,j) );
+
+  // swap two rows
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      switch (i) {
+        case 0:
+          CPPUNIT_ASSERT_EQUAL( 1., m->Get(i,j) );
+          break;
+        case 1:
+          CPPUNIT_ASSERT_EQUAL( 3., m->Get(i,j) );
+          break;
+        case 2:
+          CPPUNIT_ASSERT_EQUAL( 2., m->Get(i,j) );
+          break;
+        case 3:
+          CPPUNIT_ASSERT_EQUAL( 4., m->Get(i,j) );
+          break;
+        default:
+          CPPUNIT_ASSERT_EQUAL( -1., m->Get(i,j) );
+      }
+  // check that op is reversable
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(1,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)i+1., m->Get(i,j) );
+
+  // set to 1,2,3,1, ...
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, j+1. );
+
+  // swap such that nothing happens
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapRows(0,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->Get(i,j) );
+
+  // swap two columns
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      switch (j) {
+        case 0:
+          CPPUNIT_ASSERT_EQUAL( 3., m->Get(i,j) );
+          break;
+        case 1:
+          CPPUNIT_ASSERT_EQUAL( 2., m->Get(i,j) );
+          break;
+        case 2:
+          CPPUNIT_ASSERT_EQUAL( 1., m->Get(i,j) );
+          break;
+        default:
+          CPPUNIT_ASSERT_EQUAL( -1., m->Get(i,j) );
+      }
+  // check that op is reversable
+  CPPUNIT_ASSERT_EQUAL( true, m->SwapColumns(0,2) );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)j+1., m->Get(i,j) );
+
+
+  // set to 1,2,3,4, ...
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      m->Set(i,j, 3*i+j+1 );
+  // transpose
+  CPPUNIT_ASSERT_EQUAL( true, m->Transpose() );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(j,i) );
+  // second transpose
+  CPPUNIT_ASSERT_EQUAL( true, m->Transpose() );
+  for (int i=0;i<4;i++)
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (double)(3*i+j+1), m->Get(i,j) );
+};
+
+/** Unit Test for matrix properties.
+ *
+ */
+void GSLMatrixTest::PropertiesTest()
+{
+  // is zero
+  m->SetZero();
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
+
+  // is positive
+  m->SetAll(0.5);
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNonNegative() );
+
+  // is negative
+  m->SetAll(-0.1);
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNull() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsPositive() );
+  CPPUNIT_ASSERT_EQUAL( true, m->IsNegative() );
+  CPPUNIT_ASSERT_EQUAL( false, m->IsNonNegative() );
+
+  // is positive definite
+  CPPUNIT_ASSERT_EQUAL( false, m->IsPositiveDefinite() );
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/gslmatrixunittest.hpp
===================================================================
--- src/unittests/gslmatrixunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslmatrixunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,42 @@
+/*
+ * gslmatrixunittest.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef GSLMATRIXUNITTEST_HPP_
+#define GSLMATRIXUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "gslmatrix.hpp"
+
+/********************************************** Test classes **************************************/
+
+class GSLMatrixTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( GSLMatrixTest) ;
+    CPPUNIT_TEST (AccessTest );
+    CPPUNIT_TEST (InitializationTest );
+    CPPUNIT_TEST (CopyTest );
+    CPPUNIT_TEST (ExchangeTest );
+    CPPUNIT_TEST (PropertiesTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+    void setUp();
+    void tearDown();
+
+    void AccessTest();
+    void InitializationTest();
+    void CopyTest();
+    void ExchangeTest();
+    void PropertiesTest();
+
+private:
+    GSLMatrix *m;
+};
+
+
+#endif /* GSLMATRIXUNITTEST_HPP_ */
Index: src/unittests/gslvectorunittest.cpp
===================================================================
--- src/unittests/gslvectorunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslvectorunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,136 @@
+/*
+ * gslvectorunittest.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "gslvectorunittest.hpp"
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( GSLVectorTest );
+
+
+void GSLVectorTest::setUp()
+{
+  v = new GSLVector(3);
+};
+
+void GSLVectorTest::tearDown()
+{
+  delete(v);
+};
+
+/** Unit test for accessing elements in the vector.
+ *
+ */
+void GSLVectorTest::AccessTest()
+{
+  // set
+  for (int i=0;i<3;i++)
+    v->Set(i,i);
+
+  // and check
+  double *ptr = NULL;
+  for (int i=0;i<3;i++) {
+    CPPUNIT_ASSERT_EQUAL( (double)i, v->Get(i) );
+    ptr = v->Pointer(i);
+    CPPUNIT_ASSERT_EQUAL( (double)i, *ptr );
+  }
+
+  // out of bounds
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(3) );
+  //CPPUNIT_ASSERT_EQUAL(0., v->Get(-1) );
+};
+
+
+/** Unit test for initializing the vector.
+ *
+ */
+void GSLVectorTest::InitializaionTest()
+{
+  // set zero
+  v->SetZero();
+  for (int i=0;i<3;i++)
+    CPPUNIT_ASSERT_EQUAL( 0., v->Get(i) );
+
+  // set basis
+  for (int i=0;i<3;i++) {
+    v->SetBasis(i);
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( i == j ? 1. : 0. , v->Get(j) );
+  }
+
+  // set all
+  v->SetAll(1.5);
+  for (int i=0;i<3;i++)
+    CPPUNIT_ASSERT_EQUAL( 1.5, v->Get(i) );
+};
+
+/** Unit test for copying vectors.
+ *
+ */
+void GSLVectorTest::CopyTest()
+{
+  // set basis
+  GSLVector *dest = NULL;
+  for (int i=0;i<3;i++) {
+    v->SetBasis(i);
+    dest = new GSLVector(v);
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( v->Get(j) , dest->Get(j) );
+
+    delete(dest);
+  }
+};
+
+/** Unit test for exchanging elements of a vector.
+ *
+ */
+void GSLVectorTest::ExchangeTest()
+{
+  // set basis
+  for (int i=0;i<3;i++) {
+    v->SetBasis(i);
+    v->SwapElements((i)%3,(i+1)%3);
+    for (int j=0;j<3;j++)
+      CPPUNIT_ASSERT_EQUAL( (i+1)%3 == j ? 1. : 0. , v->Get(j) );
+  }
+
+  // set to 1,2,3 and reverse
+  for (int i=0;i<3;i++)
+    v->Set(i,i+1);
+  v->Reverse();
+  for (int j=0;j<3;j++)
+    CPPUNIT_ASSERT_EQUAL( (double)(3-j), v->Get(j) );
+};
+
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/gslvectorunittest.hpp
===================================================================
--- src/unittests/gslvectorunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/gslvectorunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,41 @@
+/*
+ * gslvectorunittest.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef GSLVECTORUNITTEST_HPP_
+#define GSLVECTORUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "gslvector.hpp"
+
+/********************************************** Test classes **************************************/
+
+class GSLVectorTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( GSLVectorTest) ;
+    CPPUNIT_TEST (AccessTest );
+    CPPUNIT_TEST (InitializaionTest );
+    CPPUNIT_TEST (CopyTest );
+    CPPUNIT_TEST (ExchangeTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+    void setUp();
+    void tearDown();
+
+    void AccessTest();
+    void InitializaionTest();
+    void CopyTest();
+    void ExchangeTest();
+
+private:
+
+    GSLVector *v;
+};
+
+
+#endif /* GSLVECTORUNITTEST_HPP_ */
Index: src/unittests/linearsystemofequationsunittest.cpp
===================================================================
--- src/unittests/linearsystemofequationsunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/linearsystemofequationsunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,125 @@
+/*
+ * linearsystemofequationsunittest.cpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <iomanip>
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include "linearsystemofequationsunittest.hpp"
+#include "vector.hpp"
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( LinearSystemOfEquationsTest );
+
+
+void LinearSystemOfEquationsTest::setUp()
+{
+  s = new LinearSystemOfEquations(4,4);
+};
+
+void LinearSystemOfEquationsTest::tearDown()
+{
+  delete(s);
+};
+
+/** Unit test for initialization.
+ *
+ */
+void LinearSystemOfEquationsTest::InitializationTest()
+{
+  // Setting A
+  double array[] = { 1., 0., 0., 0.,
+                     0., 2., 0., 0.,
+                     0., 0., 3., 0.,
+                     0., 0., 0., 4. };
+  s->SetA(array);
+
+  // Setting b
+};
+
+/** Unit test for symmetry property.
+ *
+ */
+void LinearSystemOfEquationsTest::SymmetricTest()
+{
+  CPPUNIT_ASSERT_EQUAL( true, s->SetSymmetric(true) );
+  //LinearSystemOfEquations *t = new LinearSystemOfEquations(4,3);
+  //CPPUNIT_ASSERT_EQUAL( true, t->SetSymmetric(true) );
+  //delete(t);
+};
+
+/** Unit test for simple Solve.
+ *
+ */
+void LinearSystemOfEquationsTest::SolveSimpleTest()
+{
+  // set up A and b
+  double m_array[] = { 1., 0., 0., 0.,
+                       0., 2., 0., 0.,
+                       0., 0., 3., 0.,
+                       0., 0., 0., 4. };
+  double v_array[] = { 1., 2., 3., 4. };
+  s->SetA(m_array);
+  s->Setb(v_array);
+
+  // solve
+  double *array = new double[4];
+  s->GetSolutionAsArray(array);
+  for (int i=0;i<4;i++)
+    CPPUNIT_ASSERT_EQUAL( 1., array[i]);
+  delete[](array);
+};
+
+/** Unit test for advanced Solve.
+ *
+ */
+void LinearSystemOfEquationsTest::SolveAdvancedTest()
+{
+  // set up A and b
+  double m_array[] = { 0.18, 0.60, 0.57, 0.96,
+                       0.41, 0.24, 0.99, 0.58,
+                       0.14, 0.30, 0.97, 0.66,
+                       0.51, 0.13, 0.19, 0.85 };
+  double v_array[] = { 1.0, 2.0, 3.0, 4.0 };
+  double x_array[] = { -4.05205022957398, -12.60561139590692, 1.660911626708844, 8.693766928795233 };
+  s->SetA(m_array);
+  s->Setb(v_array);
+
+  // solve
+  double *array = new double[4];
+  s->GetSolutionAsArray(array);
+  for (int i=0;i<4;i++) {
+    CPPUNIT_ASSERT( fabs(x_array[i] - array[i]) < MYEPSILON );
+  }
+  delete[](array);
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/linearsystemofequationsunittest.hpp
===================================================================
--- src/unittests/linearsystemofequationsunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/linearsystemofequationsunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,40 @@
+/*
+ * linearsystemofequationsunittest.hpp
+ *
+ *  Created on: Jan 8, 2010
+ *      Author: heber
+ */
+
+#ifndef LINEARSYSTEMOFEQUATIONSUNITTEST_HPP_
+#define LINEARSYSTEMOFEQUATIONSUNITTEST_HPP_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "linearsystemofequations.hpp"
+
+/********************************************** Test classes **************************************/
+
+class LinearSystemOfEquationsTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( LinearSystemOfEquationsTest) ;
+    CPPUNIT_TEST ( InitializationTest );
+    CPPUNIT_TEST ( SymmetricTest );
+    CPPUNIT_TEST ( SolveSimpleTest );
+    CPPUNIT_TEST ( SolveAdvancedTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+    void setUp();
+    void tearDown();
+
+    void InitializationTest();
+    void SymmetricTest();
+    void SolveSimpleTest();
+    void SolveAdvancedTest();
+
+private:
+    LinearSystemOfEquations *s;
+};
+
+
+#endif /* LINEARSYSTEMOFEQUATIONSUNITTEST_HPP_ */
Index: src/unittests/listofbondsunittest.cpp
===================================================================
--- src/unittests/listofbondsunittest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/listofbondsunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -11,4 +11,6 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
 
 #include "listofbondsunittest.hpp"
Index: src/unittests/tesselation_boundarytriangleunittest.cpp
===================================================================
--- src/unittests/tesselation_boundarytriangleunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/tesselation_boundarytriangleunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,214 @@
+/*
+ * tesselation_boundarytriangleunittest.cpp
+ *
+ *  Created on: Jan 13, 2010
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
+
+#include "defs.hpp"
+#include "tesselation.hpp"
+#include "tesselation_boundarytriangleunittest.hpp"
+
+#define SPHERERADIUS 2.
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( TesselationBoundaryTriangleTest );
+
+
+void TesselationBoundaryTriangleTest::setUp()
+{
+  setVerbosity(5);
+
+  // create nodes
+  TesselPoint *Walker = NULL;
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationBoundaryTriangleTest::setUp - *Name");
+  strcpy(Walker->Name, "1");
+  Walker->nr = 1;
+  points[0] = new BoundaryPointSet(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 1., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationBoundaryTriangleTest::setUp - *Name");
+  strcpy(Walker->Name, "2");
+  Walker->nr = 2;
+  points[1] = new BoundaryPointSet(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationBoundaryTriangleTest::setUp - *Name");
+  strcpy(Walker->Name, "3");
+  Walker->nr = 3;
+  points[2] = new BoundaryPointSet(Walker);
+
+  // create line
+  lines[0] = new BoundaryLineSet(points[0], points[1], 0);
+  lines[1] = new BoundaryLineSet(points[1], points[2], 1);
+  lines[2] = new BoundaryLineSet(points[0], points[2], 2);
+
+  // create triangle
+  triangle = new BoundaryTriangleSet(lines, 0);
+  triangle->GetNormalVector(Vector(0.,0.,1.));
+};
+
+
+void TesselationBoundaryTriangleTest::tearDown()
+{
+  delete(triangle);
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
+};
+
+/** UnitTest for Tesselation::IsInnerPoint()
+ */
+void TesselationBoundaryTriangleTest::GetClosestPointOnPlaneTest()
+{
+  Vector TestIntersection;
+  Vector Point;
+
+  // simple test on y line
+  Point.Init(-1.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(-4.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( 16., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on x line
+  Point.Init(0.5,-1.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(0.5,-6.,0.);
+  CPPUNIT_ASSERT_EQUAL( 36., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on slanted line
+  Point.Init(1.,1.,0.);
+  CPPUNIT_ASSERT_EQUAL( 0.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(5.,5.,0.);
+  CPPUNIT_ASSERT_EQUAL( 40.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on first node
+  Point.Init(-1.,-1.,0.);
+  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on second node
+  Point.Init(0.,2.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,1.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on third node
+  Point.Init(2.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( 1., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(1.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+};
+
+/** UnitTest for Tesselation::IsInnerPoint()
+ */
+void TesselationBoundaryTriangleTest::GetClosestPointOffPlaneTest()
+{
+  Vector TestIntersection;
+  Vector Point;
+
+  // straight down/up
+  Point.Init(1./3.,1./3.,+5.);
+  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(1./3.,1./3.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(1./3.,1./3.,-5.);
+  CPPUNIT_ASSERT_EQUAL( 25. , triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(1./3.,1./3.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on y line
+  Point.Init(-1.,0.5,+2.);
+  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(-1.,0.5,-3.);
+  CPPUNIT_ASSERT_EQUAL( 10., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on x line
+  Point.Init(0.5,-1.,+1.);
+  CPPUNIT_ASSERT_EQUAL( 2., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(0.5,-1.,-2.);
+  CPPUNIT_ASSERT_EQUAL( 5., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on slanted line
+  Point.Init(1.,1.,+3.);
+  CPPUNIT_ASSERT_EQUAL( 9.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+  Point.Init(1.,1.,-4.);
+  CPPUNIT_ASSERT_EQUAL( 16.5, triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.5,0.5,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on first node
+  Point.Init(-1.,-1.,5.);
+  CPPUNIT_ASSERT_EQUAL( 27., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on second node
+  Point.Init(0.,2.,5.);
+  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(0.,1.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+
+  // simple test on third node
+  Point.Init(2.,0.,5.);
+  CPPUNIT_ASSERT_EQUAL( 26., triangle->GetClosestPointInsideTriangle(&Point, &TestIntersection) );
+  Point.Init(1.,0.,0.);
+  CPPUNIT_ASSERT_EQUAL( true , Point == TestIntersection );
+};
+
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/tesselation_boundarytriangleunittest.hpp
===================================================================
--- src/unittests/tesselation_boundarytriangleunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/tesselation_boundarytriangleunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,41 @@
+/*
+ * tesselation_boundarytriangleunittest.hpp
+ *
+ *  Created on: Jan 13, 2010
+ *      Author: heber
+ */
+
+#ifndef TESSELATION_BOUNDARYTRIANGLEUNITTEST_HPP_
+#define TESSELATION_BOUNDARYTRIANGLEUNITTEST_HPP_
+
+/*********************************************** includes ***********************************/
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "linkedcell.hpp"
+#include "tesselation.hpp"
+
+/********************************************** Test classes **************************************/
+
+class TesselationBoundaryTriangleTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( TesselationBoundaryTriangleTest) ;
+    CPPUNIT_TEST ( GetClosestPointOnPlaneTest );
+    CPPUNIT_TEST ( GetClosestPointOffPlaneTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void GetClosestPointOnPlaneTest();
+      void GetClosestPointOffPlaneTest();
+
+private:
+      class BoundaryTriangleSet *triangle;
+      class BoundaryLineSet *lines[3];
+      class BoundaryPointSet *points[3];
+      LinkedNodes Corners;
+};
+
+
+#endif /* TESSELATION_BOUNDARYTRIANGLEUNITTEST_HPP_ */
Index: src/unittests/tesselation_insideoutsideunittest.cpp
===================================================================
--- src/unittests/tesselation_insideoutsideunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/tesselation_insideoutsideunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,183 @@
+/*
+ * tesselation_insideoutsideunittest.cpp
+ *
+ *  Created on: Dec 28, 2009
+ *      Author: heber
+ */
+
+using namespace std;
+
+#include <cppunit/CompilerOutputter.h>
+#include <cppunit/extensions/TestFactoryRegistry.h>
+#include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
+
+#include "defs.hpp"
+#include "tesselation.hpp"
+#include "tesselation_insideoutsideunittest.hpp"
+
+#define SPHERERADIUS 2.
+
+/********************************************** Test classes **************************************/
+
+// Registers the fixture into the 'registry'
+CPPUNIT_TEST_SUITE_REGISTRATION( TesselationInOutsideTest );
+
+
+void TesselationInOutsideTest::setUp()
+{
+  setVerbosity(2);
+
+  // create corners
+  class TesselPoint *Walker;
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "1");
+  Walker->nr = 1;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 1., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "2");
+  Walker->nr = 2;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 0., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "3");
+  Walker->nr = 3;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 1., 0.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "4");
+  Walker->nr = 4;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 0., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "5");
+  Walker->nr = 5;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(0., 1., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "6");
+  Walker->nr = 6;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 0., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "7");
+  Walker->nr = 7;
+  Corners.push_back(Walker);
+  Walker = new TesselPoint;
+  Walker->node = new Vector(1., 1., 1.);
+  Walker->Name = Malloc<char>(3, "TesselationInOutsideTest::setUp - *Name");
+  strcpy(Walker->Name, "8");
+  Walker->nr = 8;
+  Corners.push_back(Walker);
+
+  // create linkedcell
+  LinkedList = new LinkedCell(&Corners, 2.*SPHERERADIUS);
+
+  // create tesselation
+  TesselStruct = new Tesselation;
+  TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
+
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
+      Log() << Verbose(2) << *(Runner->second) << endl;
+
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl;
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
+    }
+
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++)
+      Log() << Verbose(2) << *(Runner->second) << endl;
+
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+        //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
+    }
+  }
+};
+
+
+void TesselationInOutsideTest::tearDown()
+{
+  delete(LinkedList);
+  delete(TesselStruct);
+  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
+    delete((*Runner)->node);
+    delete(*Runner);
+  }
+  Corners.clear();
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
+};
+
+/** UnitTest for Tesselation::IsInnerPoint()
+ */
+void TesselationInOutsideTest::IsInnerPointTest()
+{
+  double n[3];
+  const double boundary = 2.;
+  const double step = 1.;
+
+  // go through the mesh and check each point
+  for (n[0] = -boundary; n[0] <= boundary; n[0]+=step)
+    for (n[1] = -boundary; n[1] <= boundary; n[1]+=step)
+      for (n[2] = -boundary; n[2] <= boundary; n[2]+=step) {
+        if ( ((n[0] >= 0.) && (n[1] >= 0.) && (n[2] >= 0.)) && ((n[0] <= 1.) && (n[1] <= 1.) && (n[2] <= 1.)))
+          CPPUNIT_ASSERT_EQUAL( true , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
+        else
+          CPPUNIT_ASSERT_EQUAL( false , TesselStruct->IsInnerPoint(Vector(n[0], n[1], n[2]), LinkedList) );
+      }
+};
+
+/********************************************** Main routine **************************************/
+
+int main(int argc, char **argv)
+{
+  // Get the top level suite from the registry
+  CppUnit::Test *suite = CppUnit::TestFactoryRegistry::getRegistry().makeTest();
+
+  // Adds the test to the list of test to run
+  CppUnit::TextUi::TestRunner runner;
+  runner.addTest( suite );
+
+  // Change the default outputter to a compiler error format outputter
+  runner.setOutputter( new CppUnit::CompilerOutputter( &runner.result(),
+                                                       std::cerr ) );
+  // Run the tests.
+  bool wasSucessful = runner.run();
+
+  // Return error code 1 if the one of test failed.
+  return wasSucessful ? 0 : 1;
+};
Index: src/unittests/tesselation_insideoutsideunittest.hpp
===================================================================
--- src/unittests/tesselation_insideoutsideunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
+++ src/unittests/tesselation_insideoutsideunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -0,0 +1,38 @@
+/*
+ * tesselation_insideoutsideunittest.hpp
+ *
+ *  Created on: Dec 28, 2009
+ *      Author: heber
+ */
+
+#ifndef TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_
+#define TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_
+
+/*********************************************** includes ***********************************/
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "linkedcell.hpp"
+#include "tesselation.hpp"
+
+/********************************************** Test classes **************************************/
+
+class TesselationInOutsideTest : public CppUnit::TestFixture
+{
+    CPPUNIT_TEST_SUITE( TesselationInOutsideTest) ;
+    CPPUNIT_TEST ( IsInnerPointTest );
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+      void setUp();
+      void tearDown();
+      void IsInnerPointTest();
+
+private:
+      class Tesselation *TesselStruct;
+      LinkedNodes Corners;
+      class LinkedCell *LinkedList;
+};
+
+
+#endif /* TESSELATION_INSIDEOUTSIDEUNITTEST_HPP_ */
Index: src/unittests/tesselationunittest.cpp
===================================================================
--- src/unittests/tesselationunittest.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/tesselationunittest.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -12,4 +12,6 @@
 #include <cppunit/extensions/TestFactoryRegistry.h>
 #include <cppunit/ui/text/TestRunner.h>
+
+#include <cstring>
 
 #include "defs.hpp"
@@ -30,18 +32,18 @@
   class TesselPoint *Walker;
   Walker = new TesselPoint;
-  Walker->node = new Vector(1., 0., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(1., 0., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "1");
   Walker->nr = 1;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., 1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., 1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "2");
   Walker->nr = 2;
   Corners.push_back(Walker);
   Walker = new TesselPoint;
-  Walker->node = new Vector(-1., -1., 0.);
-  Walker->Name = new char[3];
+  Walker->node = new Vector(-1., -1., -1.);
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "3");
   Walker->nr = 3;
@@ -49,5 +51,5 @@
   Walker = new TesselPoint;
   Walker->node = new Vector(-1., 0., 1.);
-  Walker->Name = new char[3];
+  Walker->Name = Malloc<char>(3, "TesselationTest::setUp");
   strcpy(Walker->Name, "4");
   Walker->nr = 4;
@@ -59,19 +61,35 @@
   // create tesselation
   TesselStruct = new Tesselation;
-  TesselStruct->PointsOnBoundary.clear();
-  TesselStruct->LinesOnBoundary.clear();
-  TesselStruct->TrianglesOnBoundary.clear();
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->PointsOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->LinesOnBoundary.empty() );
+  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->TrianglesOnBoundary.empty() );
   TesselStruct->FindStartingTriangle(SPHERERADIUS, LinkedList);
-  bool flag = false;
 
-  LineMap::iterator baseline = TesselStruct->LinesOnBoundary.begin();
-  while (baseline != TesselStruct->LinesOnBoundary.end()) {
-    if (baseline->second->triangles.size() == 1) {
-      flag = TesselStruct->FindNextSuitableTriangle(*(baseline->second), *(((baseline->second->triangles.begin()))->second), SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+  CandidateForTesselation *baseline = NULL;
+  BoundaryTriangleSet *T = NULL;
+  bool OneLoopWithoutSuccessFlag = true;
+  bool TesselationFailFlag = false;
+  while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) {
+    // 2a. fill all new OpenLines
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      baseline = Runner->second;
+      if (baseline->pointlist.empty()) {
+        T = (((baseline->BaseLine->triangles.begin()))->second);
+        TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); //the line is there, so there is a triangle, but only one.
+      }
     }
-    baseline++;
-    if ((baseline == TesselStruct->LinesOnBoundary.end()) && (flag)) {
-      baseline = TesselStruct->LinesOnBoundary.begin();   // restart if we reach end due to newly inserted lines
-      flag = false;
+
+    // 2b. search for smallest ShortestAngle among all candidates
+    double ShortestAngle = 4.*M_PI;
+    for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) {
+      if (Runner->second->ShortestAngle < ShortestAngle) {
+        baseline = Runner->second;
+        ShortestAngle = baseline->ShortestAngle;
+      }
+    }
+    if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty()))
+      OneLoopWithoutSuccessFlag = false;
+    else {
+      TesselStruct->AddCandidateTriangle(*baseline);
     }
   }
@@ -84,33 +102,11 @@
   delete(TesselStruct);
   for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++) {
-    delete[]((*Runner)->Name);
     delete((*Runner)->node);
     delete(*Runner);
   }
   Corners.clear();
-};
-
-/** UnitTest for Tesselation::IsInnerPoint()
- */
-void TesselationTest::IsInnerPointTest()
-{
-  // true inside points
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.5,0.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.5,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint(Vector(0.,0.,0.5), LinkedList) );
-
-  // corners
-  for (LinkedNodes::iterator Runner = Corners.begin(); Runner != Corners.end(); Runner++)
-    CPPUNIT_ASSERT_EQUAL( true, TesselStruct->IsInnerPoint((*Runner), LinkedList) );
-
-  // true outside points
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,5.,0.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(0.,0.,5.), LinkedList) );
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(1.,1.,1.), LinkedList) );
-
-  // tricky point, there are three equally close triangles
-  CPPUNIT_ASSERT_EQUAL( false, TesselStruct->IsInnerPoint(Vector(5.,0.,0.), LinkedList) );
-
+  MemoryUsageObserver::purgeInstance();
+  logger::purgeInstance();
+  errorLogger::purgeInstance();
 };
 
Index: src/unittests/tesselationunittest.hpp
===================================================================
--- src/unittests/tesselationunittest.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/unittests/tesselationunittest.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -21,5 +21,4 @@
 {
     CPPUNIT_TEST_SUITE( TesselationTest) ;
-    CPPUNIT_TEST ( IsInnerPointTest );
     CPPUNIT_TEST ( GetAllTrianglesTest );
     CPPUNIT_TEST ( ContainmentTest );
@@ -29,5 +28,4 @@
       void setUp();
       void tearDown();
-      void IsInnerPointTest();
       void GetAllTrianglesTest();
       void ContainmentTest();
Index: src/vector.cpp
===================================================================
--- src/vector.cpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/vector.cpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -8,9 +8,16 @@
 #include "defs.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
+#include "info.hpp"
+#include "gslmatrix.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
+#include "memoryallocator.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
+
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
 
 /************************************ Functions for class vector ************************************/
@@ -215,8 +222,9 @@
  * \param *Origin first vector of line
  * \param *LineVector second vector of line
- * \return true -  \a this contains intersection point on return, false - line is parallel to plane
+ * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
  */
 bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
 {
+  Info FunctionInfo(__func__);
   double factor;
   Vector Direction, helper;
@@ -226,8 +234,9 @@
   Direction.SubtractVector(Origin);
   Direction.Normalize();
-  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
+  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
-  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
-    eLog() << Verbose(2) << "Line is parallel to plane, no intersection." << endl;
+  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
+    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
     return false;
   }
@@ -235,6 +244,6 @@
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (factor < MYEPSILON) { // Origin is in-plane
-    //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
+  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
+    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
     CopyVector(Origin);
     return true;
@@ -243,5 +252,5 @@
   Direction.Scale(factor);
   CopyVector(Origin);
-  //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
+  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 
@@ -250,5 +259,5 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
     return true;
   } else {
@@ -286,7 +295,5 @@
 
 /** Calculates the intersection of the two lines that are both on the same plane.
- * We construct auxiliary plane with its vector normal to one line direction and the PlaneNormal, then a vector
- * from the first line's offset onto the plane. Finally, scale by factor is 1/cos(angle(line1,line2..)) = 1/SP(...), and
- * project onto the first line's direction and add its offset.
+ * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html 
  * \param *out output stream for debugging
  * \param *Line1a first vector of first line
@@ -299,80 +306,96 @@
 bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
 {
-  bool result = true;
-  Vector Direction, OtherDirection;
-  Vector AuxiliaryNormal;
-  Vector Distance;
-  const Vector *Normal = NULL;
-  Vector *ConstructedNormal = NULL;
-  bool FreeNormal = false;
-
-  // construct both direction vectors
-  Zero();
-  Direction.CopyVector(Line1b);
-  Direction.SubtractVector(Line1a);
-  if (Direction.IsZero())
+  Info FunctionInfo(__func__);
+
+  GSLMatrix *M = new GSLMatrix(4,4);
+
+  M->SetAll(1.);
+  for (int i=0;i<3;i++) {
+    M->Set(0, i, Line1a->x[i]);
+    M->Set(1, i, Line1b->x[i]);
+    M->Set(2, i, Line2a->x[i]);
+    M->Set(3, i, Line2b->x[i]);
+  }
+  
+  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
+  //for (int i=0;i<4;i++) {
+  //  for (int j=0;j<4;j++)
+  //    cout << "\t" << M->Get(i,j);
+  //  cout << endl;
+  //}
+  if (fabs(M->Determinant()) > MYEPSILON) {
+    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
     return false;
-  OtherDirection.CopyVector(Line2b);
-  OtherDirection.SubtractVector(Line2a);
-  if (OtherDirection.IsZero())
+  }
+  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
+
+
+  // constuct a,b,c
+  Vector a;
+  Vector b;
+  Vector c;
+  Vector d;
+  a.CopyVector(Line1b);
+  a.SubtractVector(Line1a);
+  b.CopyVector(Line2b);
+  b.SubtractVector(Line2a);
+  c.CopyVector(Line2a);
+  c.SubtractVector(Line1a);
+  d.CopyVector(Line2b);
+  d.SubtractVector(Line1b);
+  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
+  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
+   Zero();
+   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
+   return false;
+  }
+
+  // check for parallelity
+  Vector parallel;
+  double factor = 0.;
+  if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+    parallel.CopyVector(Line1a);
+    parallel.SubtractVector(Line2a);
+    factor = parallel.ScalarProduct(&a)/a.Norm();
+    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+      CopyVector(Line2a);
+      Log() << Verbose(1) << "Lines conincide." << endl;
+      return true;
+    } else {
+      parallel.CopyVector(Line1a);
+      parallel.SubtractVector(Line2b);
+      factor = parallel.ScalarProduct(&a)/a.Norm();
+      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
+        CopyVector(Line2b);
+        Log() << Verbose(1) << "Lines conincide." << endl;
+        return true;
+      }
+    }
+    Log() << Verbose(1) << "Lines are parallel." << endl;
+    Zero();
     return false;
-
-  Direction.Normalize();
-  OtherDirection.Normalize();
-
-  //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
-
-  if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
-    if ((Line1a == Line2a) || (Line1a == Line2b))
-      CopyVector(Line1a);
-    else if ((Line1b == Line2b) || (Line1b == Line2b))
-        CopyVector(Line1b);
-    else
-      return false;
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
-    return true;
-  } else {
-    // check whether we have a plane normal vector
-    if (PlaneNormal == NULL) {
-      ConstructedNormal = new Vector;
-      ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
-      Normal = ConstructedNormal;
-      FreeNormal = true;
-    } else
-      Normal = PlaneNormal;
-
-    AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
-    //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
-
-    Distance.CopyVector(Line2a);
-    Distance.SubtractVector(Line1a);
-    //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
-    if (Distance.IsZero()) {
-      // offsets are equal, match found
-      CopyVector(Line1a);
-      result = true;
-    } else {
-      CopyVector(Distance.Projection(&AuxiliaryNormal));
-      //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
-      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
-      //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
-      Scale(1./(factor*factor));
-      //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
-      CopyVector(Projection(&Direction));
-      //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
-      if (this->IsZero())
-        result = false;
-      else
-        result = true;
-      AddVector(Line1a);
-    }
-
-    if (FreeNormal)
-      delete(ConstructedNormal);
-  }
-  if (result)
-    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
-
-  return result;
+  }
+
+  // obtain s
+  double s;
+  Vector temp1, temp2;
+  temp1.CopyVector(&c);
+  temp1.VectorProduct(&b);
+  temp2.CopyVector(&a);
+  temp2.VectorProduct(&b);
+  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
+  if (fabs(temp2.NormSquared()) > MYEPSILON)
+    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+  else
+    s = 0.;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+
+  // construct intersection
+  CopyVector(&a);
+  Scale(s);
+  AddVector(Line1a);
+  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
+
+  return true;
 };
 
@@ -480,4 +503,17 @@
   else
     return false;
+};
+
+/** Checks whether vector is normal to \a *normal.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsEqualTo(const Vector * const a) const
+{
+  bool status = true;
+  for (int i=0;i<NDIM;i++) {
+    if (fabs(x[i] - a->x[i]) > MYEPSILON)
+      status = false;
+  }
+  return status;
 };
 
Index: src/vector.hpp
===================================================================
--- src/vector.hpp	(revision f34def7cdde1404608ef61377d53e30b1c5d0bdb)
+++ src/vector.hpp	(revision 315164b6c8e5a1302b98d6dd1aefc8ede4ec56b5)
@@ -42,4 +42,5 @@
   bool IsOne() const;
   bool IsNormalTo(const Vector * const normal) const;
+  bool IsEqualTo(const Vector * const a) const;
 
   void AddVector(const Vector * const y);
