Changes in / [7c14ec:34e0592]

Files:

• configure.ac (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/boundary.cpp (modified) (7 diffs)
• src/boundary.hpp (modified) (1 diff)
• src/builder.cpp (modified) (2 diffs)
• src/config.hpp (modified) (1 diff)
• src/leastsquaremin.cpp (added)
• src/leastsquaremin.hpp (added)
• src/molecules.cpp (modified) (1 diff)
• src/molecules.hpp (modified) (3 diffs)
• src/tesselation.cpp (modified) (23 diffs)
• src/tesselation.hpp (modified) (1 diff)
• src/tesselationhelpers.cpp (modified) (2 diffs)
• src/vector.cpp (modified) (21 diffs)
• src/vector.hpp (modified) (3 diffs)
• src/vectorunittest.cpp (added)
• src/vectorunittest.hpp (added)

configure.ac

@@ -11 +11 @@
 
 # Checks for programs.
+AM_PATH_CPPUNIT(1.9.6)
 AC_PROG_CXX
 AC_PROG_CC
+AC_PROG_INSTALL
 AC_CHECK_PROG([LATEX],[latex],[latex],[:])
 AC_CHECK_PROG([BIBTEX],[bibtex],[bibtex],[:])

src/Makefile.am

@@ -1 @@
-SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
-HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp helpers.cpp leastsquaremin.cpp linkedcell.cpp moleculelist.cpp molecules.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
+HEADER = atom.hpp bond.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp helpers.hpp leastsquaremin.hpp linkedcell.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
+
+noinst_PROGRAMS =  VectorUnitTest
 
 bin_PROGRAMS = molecuilder joiner analyzer
 molecuilderdir = ${bindir}
 molecuilder_DATA = elements.db valence.db orbitals.db Hbonddistance.db Hbondangle.db
-molecuilder_SOURCES =  ${SOURCE} ${HEADER}
+molecuilder_SOURCES =  ${SOURCE} builder.cpp  ${HEADER}
 joiner_SOURCES = joiner.cpp datacreator.cpp element.cpp helpers.cpp periodentafel.cpp parser.cpp verbose.cpp datacreator.hpp helpers.hpp parser.hpp periodentafel.hpp 
 analyzer_SOURCES = analyzer.cpp datacreator.cpp element.cpp helpers.cpp periodentafel.cpp parser.cpp verbose.cpp helpers.hpp periodentafel.hpp parser.hpp datacreator.hpp 
 
+TESTS = VectorUnitTest
+check_PROGRAMS = $(TESTS)
+VectorUnitTest_SOURCES = defs.hpp helpers.cpp helpers.hpp leastsquaremin.cpp leastsquaremin.hpp vectorunittest.cpp vectorunittest.hpp vector.cpp vector.hpp verbose.cpp verbose.hpp 
+VectorUnitTest_CXXFLAGS = $(CPPUNIT_CFLAGS)
+VectorUnitTest_LDFLAGS = $(CPPUNIT_LIBS) -ldl
 
 EXTRA_DIST = ${molecuilder_DATA}

src/boundary.cpp

@@ -308 +308 @@
 
       //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+      if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
         angle = 2. * M_PI - angle;
       }
@@ -620 +620 @@
   }
 
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-first");
+
   // First step: RemovePointFromTesselatedSurface
-  PointRunner = TesselStruct->PointsOnBoundary.begin();
-  PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
-  while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
-    PointAdvance++;
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    Concavity = true;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+  do {
+    Concavity = false;
+    PointRunner = TesselStruct->PointsOnBoundary.begin();
+    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
+    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
+      PointAdvance++;
+      point = PointRunner->second;
+      *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+        line = LineRunner->second;
+        *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      }
+      if (!line->CheckConvexityCriterion(out)) {
+        *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+        // remove the point
+        Concavity = true;
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+      }
+      PointRunner = PointAdvance;
+    }
+
+    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+    //StoreTrianglesinFile(out, mol, filename, "-second");
+
+    // second step: PickFarthestofTwoBaselines
+    LineRunner = TesselStruct->LinesOnBoundary.begin();
+    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+      LineAdvance++;
       line = LineRunner->second;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      Concavity = Concavity && (!line->CheckConvexityCriterion(out));
-    }
-    if (Concavity) {
-      *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
-      // remove the point
-      TesselStruct->RemovePointFromTesselatedSurface(out, point);
-    }
-    PointRunner = PointAdvance;
-  }
-
-
-//  // print all lines
-//  LineRunner = TesselStruct->LinesOnBoundary.begin();
-//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary lines for debugging." << endl;
-//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
-//    LineAdvance++;
-//    line = LineRunner->second;
-//    *out << Verbose(2) << "INFO: Current line is " << *line << "." << endl;
-//    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
-//    LineRunner = LineAdvance;
-//  }
-//
-//  // print all triangles
-//  TriangleRunner = TesselStruct->TrianglesOnBoundary.begin();
-//  TriangleAdvance = TriangleRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary triangles for debugging." << endl;
-//  while (TriangleRunner != TesselStruct->TrianglesOnBoundary.end()) {
-//    TriangleAdvance++;
-//    *out << Verbose(2) << "INFO: Current triangle is " << *(TriangleRunner->second) << "." << endl;
-//    if (TriangleAdvance != TesselStruct->TrianglesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next triangle will be " << *(TriangleAdvance->second) << "." << endl;
-//    TriangleRunner = TriangleAdvance;
-//  }
-
-  // second step: PickFarthestofTwoBaselines
+      *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
+      // take highest of both lines
+      if (TesselStruct->IsConvexRectangle(out, line) == NULL) {
+        TesselStruct->PickFarthestofTwoBaselines(out, line);
+        Concavity = true;
+      }
+      LineRunner = LineAdvance;
+    }
+  } while (Concavity);
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-third");
+
+  // third step: IsConvexRectangle
   LineRunner = TesselStruct->LinesOnBoundary.begin();
   LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
@@ -673 +672 @@
     LineAdvance++;
     line = LineRunner->second;
-    *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
-    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-    // take highest of both lines
-    TesselStruct->PickFarthestofTwoBaselines(out, line);
+    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+      //*out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+    if (!line->CheckConvexityCriterion(out)) {
+      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+
+      // take highest of both lines
+      point = TesselStruct->IsConvexRectangle(out, line);
+      if (point != NULL)
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+    }
     LineRunner = LineAdvance;
   }
 
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-fourth");
+
+  // end
+  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
+  return volume;
+};
+
+/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
+ * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
+ * \param *out output stream for debugging
+ * \param *TesselStruct pointer to Tesselation structure
+ */
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+{
+  class BoundaryPointSet *point = NULL;
+  class BoundaryLineSet *line = NULL;
   // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     point = PointRunner->second;
     *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
@@ -692 +715 @@
     }
   }
-
+};
+
+/** Stores triangles to file.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and bonds
+ * \param *filename prefix of filename
+ * \param *extraSuffix intermediate suffix
+ */
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+{
   // 4. Store triangles in tecplot file
   if (filename != NULL) {
     if (DoTecplotOutput) {
       string OutputName(filename);
-      OutputName.append("_intermed");
+      OutputName.append(extraSuffix);
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
@@ -706 +738 @@
     if (DoRaster3DOutput) {
       string OutputName(filename);
-      OutputName.append("_intermed");
+      OutputName.append(extraSuffix);
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
@@ -714 +746 @@
     }
   }
-
-  // third step: IsConvexRectangle
-  LineRunner = TesselStruct->LinesOnBoundary.begin();
-  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
-  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
-    LineAdvance++;
-    line = LineRunner->second;
-    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
-    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-      *out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
-    if (!line->CheckConvexityCriterion(out)) {
-      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
-
-      // take highest of both lines
-      point = TesselStruct->IsConvexRectangle(out, line);
-      if (point != NULL)
-        TesselStruct->RemovePointFromTesselatedSurface(out, point);
-    }
-    LineRunner = LineAdvance;
-  }
-
-  // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    point->value = 0;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
-      line = LineRunner->second;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      if (!line->CheckConvexityCriterion(out))
-        point->value += 1;
-    }
-  }
-
-  // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
-
-  // end
-  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
-  return volume;
 };
 
@@ -804 +779 @@
       G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
       x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
-      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
+      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
       h = x.Norm(); // distance of CoG to triangle
       PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)

src/boundary.hpp

@@ -42 +42 @@
 void FindNextSuitablePoint(class BoundaryTriangleSet *BaseTriangle, class BoundaryLineSet *BaseLine, atom*& OptCandidate, Vector *OptCandidateCenter, double *ShortestAngle, const double RADIUS, LinkedCell *LC);
 Boundaries *GetBoundaryPoints(ofstream *out, molecule *mol);
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct);
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix);
+
 
 #endif /*BOUNDARY_HPP_*/

src/builder.cpp

@@ -188 +188 @@
           // remove the projection onto the rotation plane of the second angle
           n.CopyVector(&y);
-          n.Scale(first->x.Projection(&y));
+          n.Scale(first->x.ScalarProduct(&y));
           cout << "N1: ",
           n.Output((ofstream *)&cout);
@@ -197 +197 @@
           cout << endl;
           n.CopyVector(&z);
-          n.Scale(first->x.Projection(&z));
+          n.Scale(first->x.ScalarProduct(&z));
           cout << "N2: ",
           n.Output((ofstream *)&cout);

src/config.hpp

@@ -18 +18 @@
 #include "molecules.hpp"
 #include "periodentafel.hpp"
+
+class ConfigFileBuffer {
+  public:
+    char **buffer;
+    int *LineMapping;
+    int CurrentLine;
+    int NoLines;
+
+    ConfigFileBuffer();
+    ConfigFileBuffer(char *filename);
+    ~ConfigFileBuffer();
+
+    void InitMapping();
+    void MapIonTypesInBuffer(int NoAtoms);
+};
 
 /** The config file.

src/molecules.cpp

@@ -7 +7 @@
 #include "config.hpp"
 #include "molecules.hpp"
-
-/************************************* Other Functions *************************************/
-
-/** Determines sum of squared distances of \a X to all \a **vectors.
- * \param *x reference vector
- * \param *params
- * \return sum of square distances
- */
-double LSQ (const gsl_vector * x, void * params)
-{
-  double sum = 0.;
-  struct LSQ_params *par = (struct LSQ_params *)params;
-  Vector **vectors = par->vectors;
-  int num = par->num;
-
-  for (int i=num;i--;) {
-    for(int j=NDIM;j--;)
-      sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
-  }
-
-  return sum;
-};
 
 /************************************* Functions for class molecule *********************************/

src/molecules.hpp

@@ -28 +28 @@
 #include "bond.hpp"
 #include "element.hpp"
+#include "leastsquaremin.hpp"
 #include "linkedcell.hpp"
 #include "parser.hpp"
@@ -83 +84 @@
 
 
-// some algebraic matrix stuff
-#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
-#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix
-
-
-/** Parameter structure for least square minimsation.
- */
-struct LSQ_params {
-  Vector **vectors;
-  int num;
-};
-
-double LSQ(const gsl_vector * x, void * params);
-
-/** Parameter structure for least square minimsation.
- */
-struct lsq_params {
-  gsl_vector *x;
-  const molecule *mol;
-  element *type;
-};
 
 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented 
@@ -323 +303 @@
 };
 
-class ConfigFileBuffer {
-  public:
-    char **buffer;
-    int *LineMapping;
-    int CurrentLine;
-    int NoLines;
-
-    ConfigFileBuffer();
-    ConfigFileBuffer(char *filename);
-    ~ConfigFileBuffer();
-
-    void InitMapping();
-    void MapIonTypesInBuffer(int NoAtoms);
-};
-
 
 #endif /*MOLECULES_HPP_*/

src/tesselation.cpp

@@ -16 +16 @@
   LinesCount = 0;
   Nr = -1;
+  value = 0.;
 };
 
@@ -26 +27 @@
   LinesCount = 0;
   Nr = Walker->nr;
+  value = 0.;
 };
 
@@ -148 +150 @@
 void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle)
 {
-  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "."
-      << endl;
+  cout << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl;
   triangles.insert(TrianglePair(triangle->Nr, triangle));
 };
@@ -177 +178 @@
   if (triangles.size() != 2) {
     *out << Verbose(1) << "ERROR: Baseline " << *this << " is connect to less than two triangles, Tesselation incomplete!" << endl;
-    return false;
+    return true;
   }
   // check normal vectors
   // have a normal vector on the base line pointing outwards
-  *out << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
+  //*out << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;
   BaseLineCenter.CopyVector(endpoints[0]->node->node);
   BaseLineCenter.AddVector(endpoints[1]->node->node);
@@ -187 +188 @@
   BaseLine.CopyVector(endpoints[0]->node->node);
   BaseLine.SubtractVector(endpoints[1]->node->node);
-  *out << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
+  //*out << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
 
   BaseLineNormal.Zero();
@@ -195 +196 @@
   class BoundaryPointSet *node = NULL;
   for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) {
-    *out << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
+    //*out << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;
     NormalCheck.AddVector(&runner->second->NormalVector);
     NormalCheck.Scale(sign);
@@ -202 +203 @@
     node = runner->second->GetThirdEndpoint(this);
     if (node != NULL) {
-      *out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
+      //*out << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
       helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
-      *out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
+      //*out << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
       i++;
     } else {
-      *out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
+      //*out << Verbose(2) << "WARNING: I cannot find third node in triangle, something's wrong." << endl;
       return true;
     }
   }
-  *out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
+  //*out << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;
   if (NormalCheck.NormSquared() < MYEPSILON) {
-    *out << Verbose(3) << "Normalvectors of both triangles are the same: convex." << endl;
+    *out << Verbose(2) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;
     return true;
   }
   double angle = GetAngle(helper[0], helper[1], BaseLineNormal);
-  if ((angle - M_PI) > -MYEPSILON)
+  if ((angle - M_PI) > -MYEPSILON) {
+    *out << Verbose(2) << "ACCEPT: Angle is greater than pi: convex." << endl;
     return true;
-  else
+  } else {
+    *out << Verbose(2) << "REJECT: Angle is less than pi: concave." << endl;
     return false;
+  }
 }
 
@@ -349 +353 @@
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (NormalVector.Projection(&OtherVector) > 0)
+  if (NormalVector.ScalarProduct(&OtherVector) > 0.)
     NormalVector.Scale(-1.);
 };
@@ -737 +741 @@
           TrialVector.CopyVector(checker->second->node->node);
           TrialVector.SubtractVector(A->second->node->node);
-          distance = TrialVector.Projection(&PlaneVector);
+          distance = TrialVector.ScalarProduct(&PlaneVector);
           if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
             continue;
@@ -897 +901 @@
         TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //*out << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
-        if (PropagationVector.Projection(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
+        if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline
           PropagationVector.Scale(-1.);
         *out << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl;
@@ -957 +961 @@
             TempVector.SubtractVector(Center);
             // make it always point outward
-            if (VirtualNormalVector.Projection(&TempVector) < 0)
+            if (VirtualNormalVector.ScalarProduct(&TempVector) < 0)
               VirtualNormalVector.Scale(-1.);
             // calculate angle
@@ -1530 +1534 @@
         AddTesselationLine(TPS[0], TPS[1], 0);
       }
-      cout << Verbose(2) << "Projection is " << BTS->NormalVector.Projection(&Oben) << "." << endl;
+      cout << Verbose(2) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl;
     }
     if (BTS != NULL) // we have created one starting triangle
@@ -1751 +1755 @@
   class BoundaryPointSet *Spot = NULL;
   class BoundaryLineSet *OtherBase;
-  Vector *ClosestPoint[2];
+  Vector *ClosestPoint;
 
   int m=0;
@@ -1764 +1768 @@
 
   // get the closest point on each line to the other line
-  ClosestPoint[0] = GetClosestPointBetweenLine(out, Base, OtherBase);
-  ClosestPoint[1] = GetClosestPointBetweenLine(out, OtherBase, Base);
+  ClosestPoint = GetClosestPointBetweenLine(out, Base, OtherBase);
 
   // delete the temporary other base line
@@ -1771 +1774 @@
 
   // get the distance vector from Base line to OtherBase line
-  Vector DistanceToIntersection, BaseLine;
+  Vector DistanceToIntersection[2], BaseLine;
+  double distance[2];
   BaseLine.CopyVector(Base->endpoints[1]->node->node);
   BaseLine.SubtractVector(Base->endpoints[0]->node->node);
-  DistanceToIntersection.CopyVector(ClosestPoint[0]);
-  DistanceToIntersection.SubtractVector(Base->endpoints[0]->node->node);
-  Spot = Base->endpoints[1];
-  if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
-    DistanceToIntersection.CopyVector(ClosestPoint[0]);
-    DistanceToIntersection.SubtractVector(Base->endpoints[1]->node->node);
-    Spot = Base->endpoints[0];
-    if (BaseLine.ScalarProduct(&DistanceToIntersection) < -MYEPSILON) {
-      *out << Verbose(3) << "ERROR: Something's very wrong here, both SKPs return negative?!" << endl;
-      return NULL;
-    }
-  }
-  if ((BaseLine.NormSquared() - DistanceToIntersection.NormSquared()) < -MYEPSILON) {  // distance is smaller: concave
-    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave: Base line squared length " << BaseLine.NormSquared() << " against Distance to intersection squared " << DistanceToIntersection.NormSquared() << "." << endl;
+  for (int i=0;i<2;i++) {
+    DistanceToIntersection[i].CopyVector(ClosestPoint);
+    DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node);
+    distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]);
+  }
+  delete(ClosestPoint);
+  if ((distance[0] * distance[1]) > 0)  { // have same sign?
+    *out << Verbose(3) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1]  << ". " << *Base << "' rectangle is concave." << endl;
+    if (distance[0] < distance[1]) {
+      Spot = Base->endpoints[0];
+    } else {
+      Spot = Base->endpoints[1];
+    }
     return Spot;
-  } else {  // base line is longer: convex
-    *out << Verbose(3) << "INFO: Rectangle of triangles of base line " << *Base << " is concave." << endl;
+  } else {  // different sign, i.e. we are in between
+    *out << Verbose(3) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;
     return NULL;
   }
@@ -1796 +1799 @@
 };
 
+void Tesselation::PrintAllBoundaryPoints(ofstream *out)
+{
+  // print all lines
+  *out << Verbose(1) << "Printing all boundary points for debugging:" << endl;
+  for (PointMap::iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++)
+    *out << Verbose(2) << *(PointRunner->second) << endl;
+};
+
+void Tesselation::PrintAllBoundaryLines(ofstream *out)
+{
+  // print all lines
+  *out << Verbose(1) << "Printing all boundary lines for debugging:" << endl;
+  for (LineMap::iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
+    *out << Verbose(2) << *(LineRunner->second) << endl;
+};
+
+void Tesselation::PrintAllBoundaryTriangles(ofstream *out)
+{
+  // print all triangles
+  *out << Verbose(1) << "Printing all boundary triangles for debugging:" << endl;
+  for (TriangleMap::iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
+    *out << Verbose(2) << *(TriangleRunner->second) << endl;
+};
 
 /** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
@@ -1826 +1852 @@
   Distance.SubtractVector(ClosestPoint[0]);
 
-  // delete the temporary other base line
+  // delete the temporary other base line and the closest points
+  delete(ClosestPoint[0]);
+  delete(ClosestPoint[1]);
   delete(OtherBase);
 
@@ -1843 +1871 @@
       BaseLineNormal.AddVector(&(runner->second->NormalVector));
     }
-    BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
+    BaseLineNormal.Scale(1./2.);
 
     if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
@@ -1866 +1894 @@
   // construct the plane of the two baselines (i.e. take both their directional vectors)
   Vector Normal;
-  Vector OtherBaseline;
-  Normal.CopyVector(Base->endpoints[1]->node->node);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  Vector Baseline, OtherBaseline;
+  Baseline.CopyVector(Base->endpoints[1]->node->node);
+  Baseline.SubtractVector(Base->endpoints[0]->node->node);
   OtherBaseline.CopyVector(OtherBase->endpoints[1]->node->node);
   OtherBaseline.SubtractVector(OtherBase->endpoints[0]->node->node);
+  Normal.CopyVector(&Baseline);
   Normal.VectorProduct(&OtherBaseline);
   Normal.Normalize();
-
-  // project one offset point of OtherBase onto this plane
+  *out << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
+
+  // project one offset point of OtherBase onto this plane (and add plane offset vector)
   Vector NewOffset;
   NewOffset.CopyVector(OtherBase->endpoints[0]->node->node);
+  NewOffset.SubtractVector(Base->endpoints[0]->node->node);
   NewOffset.ProjectOntoPlane(&Normal);
+  NewOffset.AddVector(Base->endpoints[0]->node->node);
   Vector NewDirection;
-  NewDirection.CopyVector(OtherBase->endpoints[1]->node->node);
-  NewDirection.ProjectOntoPlane(&Normal);
+  NewDirection.CopyVector(&NewOffset);
+  NewDirection.AddVector(&OtherBaseline);
 
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
   Intersection->GetIntersectionOfTwoLinesOnPlane(out, Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  Normal.CopyVector(Intersection);
+  Normal.SubtractVector(Base->endpoints[0]->node->node);
+  *out << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
 
   return Intersection;
@@ -2595 +2630 @@
   int i;
 
+  if (point == NULL) {
+    *out << Verbose(1) << "ERROR: Cannot remove the point " << point << ", it's NULL!" << endl;
+    return 0.;
+  } else
+    *out << Verbose(2) << "Removing point " << *point << " from tesselated boundary ..." << endl;
+
   // copy old location for the volume
   OldPoint.CopyVector(point->node->node);
@@ -2609 +2650 @@
     count+=LineRunner->second->triangles.size();
   numbers = new int[count];
+  class BoundaryTriangleSet **Candidates = new BoundaryTriangleSet*[count];
   i=0;
   for (LineMap::iterator LineRunner = point->lines.begin(); (point != NULL) && (LineRunner != point->lines.end()); LineRunner++) {
@@ -2614 +2656 @@
     for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
       triangle = TriangleRunner->second;
-      *out << Verbose(2) << "Erasing triangle " << *triangle << "." << endl;
+      Candidates[i] = triangle;
       numbers[i++] = triangle->Nr;
-      RemoveTesselationTriangle(triangle);
-      triangle = NULL;
-    }
-  }
+    }
+  }
+  for (int j=0;j<i;j++) {
+    RemoveTesselationTriangle(Candidates[j]);
+  }
+  delete[](Candidates);
   *out << Verbose(1) << i << " triangles were removed." << endl;
 
@@ -3023 +3067 @@
 double GetAngle(const Vector &point, const Vector &reference, const Vector OrthogonalVector)
 {
-  if (reference.IsNull())
+  if (reference.IsZero())
     return M_PI;
 
   // calculate both angles and correct with in-plane vector
-  if (point.IsNull())
+  if (point.IsZero())
     return M_PI;
   double phi = point.Angle(&reference);

src/tesselation.hpp

@@ -230 +230 @@
     bool AddBoundaryPoint(TesselPoint *Walker, int n);
 
+    // print for debugging
+    void PrintAllBoundaryPoints(ofstream *out);
+    void PrintAllBoundaryLines(ofstream *out);
+    void PrintAllBoundaryTriangles(ofstream *out);
+
+
     PointMap PointsOnBoundary;
     LineMap LinesOnBoundary;

src/tesselationhelpers.cpp

@@ -359 +359 @@
   HeightA.SubtractVector(&par.x1);
 
-  t1 = HeightA.Projection(&SideA)/SideA.ScalarProduct(&SideA);
+  t1 = HeightA.ScalarProduct(&SideA)/SideA.ScalarProduct(&SideA);
 
   SideB.CopyVector(&par.x4);
@@ -366 +366 @@
   HeightB.SubtractVector(&par.x3);
 
-  t2 = HeightB.Projection(&SideB)/SideB.ScalarProduct(&SideB);
+  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
 
   cout << Verbose(2) << "Intersection " << intersection << " is at "

src/vector.cpp

@@ -6 +6 @@
 
 
-#include "molecules.hpp"
-
+#include "defs.hpp"
+#include "helpers.hpp"
+#include "leastsquaremin.hpp"
+#include "vector.hpp"
+#include "verbose.hpp"
 
 /************************************ Functions for class vector ************************************/
@@ -209 +212 @@
  * \param *PlaneNormal Plane's normal vector
  * \param *PlaneOffset Plane's offset vector
- * \param *LineVector first vector of line
- * \param *LineVector2 second vector of line
+ * \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
  */
@@ -221 +224 @@
   Direction.CopyVector(LineVector);
   Direction.SubtractVector(Origin);
+  //*out << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
   if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
@@ -227 +231 @@
   }
   helper.CopyVector(PlaneOffset);
-  helper.SubtractVector(LineVector);
+  helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
   //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
   Direction.Scale(factor);
-  CopyVector(LineVector);
+  CopyVector(Origin);
+  //*out << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 
@@ -238 +243 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    *out << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
+    //*out << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
   } else {
@@ -247 +252 @@
 
 /** Calculates the intersection of the two lines that are both on the same plane.
- * Note that we do not check whether they are 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.
  * \param *out output stream for debugging
  * \param *Line1a first vector of first line
@@ -258 +265 @@
 bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *PlaneNormal)
 {
-  double factor1, factor2;
-  Vector helper, Line, LineNormal, *OtherNormal = NULL;
-  const Vector *Normal;
-  bool result = false;
-
-  // create Plane normal vector
-  if (PlaneNormal == NULL) {
-    OtherNormal = new Vector(0.,0.,0.);
-    if (!OtherNormal->MakeNormalVector(Line1a, Line1b, Line2a))
-      if (!OtherNormal->MakeNormalVector(Line1a, Line1b, Line2b)) {
-        *out << Verbose(1) << "ERROR: GetIntersectionOfTwoLinesOnPlane() cannot create a normal of the plane, everything is linear dependent." << endl;
-        return false;
-      }
-    Normal = OtherNormal;
-  } else
-    Normal = PlaneNormal;
-  *out << Verbose(3) << "INFO: Normal of plane is " << *Normal << "." << endl;
-
-  // create normal vector to one line
-  Line.CopyVector(Line1b);
-  Line.SubtractVector(Line1a);
-  LineNormal.MakeNormalVector(&Line, Normal);
-  *out << Verbose(3) << "INFO: Normal of first line is " << LineNormal << "." << endl;
-
-  // check if lines are parallel
-  helper.CopyVector(Line2b);
-  helper.SubtractVector(Line2a);
-  if (fabs(helper.ScalarProduct(&LineNormal)) < MYEPSILON) {
-    *out << Verbose(1) << "Lines " << helper << " and " << Line << " are parallel, no cross point!" << endl;
-    result = false;
-  } else {  
-    helper.CopyVector(Line2a);
-    helper.SubtractVector(Line1a);
-    factor1 = helper.ScalarProduct(&LineNormal);
-    helper.CopyVector(Line2b);
-    helper.SubtractVector(Line1a);
-    factor2 = helper.ScalarProduct(&LineNormal);
-    if (fabs(factor2) > MYEPSILON) {
-      CopyVector(Line2a);
-      helper.Scale(factor1/factor2);
-      AddVector(&helper);
+  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())
+    return false;
+  OtherDirection.CopyVector(Line2b);
+  OtherDirection.SubtractVector(Line2a);
+  if (OtherDirection.IsZero())
+    return false;
+
+  Direction.Normalize();
+  OtherDirection.Normalize();
+
+  //*out << 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;
+    *out << 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);
+    //*out << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
+
+    Distance.CopyVector(Line2a);
+    Distance.SubtractVector(Line1a);
+    //*out << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
+    if (Distance.IsZero()) {
+      // offsets are equal, match found
+      CopyVector(Line1a);
       result = true;
     } else {
-      Zero();
-      result = false;
+      CopyVector(Distance.Projection(&AuxiliaryNormal));
+      //*out << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
+      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
+      //*out << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
+      Scale(1./(factor*factor));
+      //*out << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
+      CopyVector(Projection(&Direction));
+      //*out << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
+      if (this->IsZero())
+        result = false;
+      else
+        result = true;
+      AddVector(Line1a);
     }
-  }
-
-  if (OtherNormal != NULL)
-    delete(OtherNormal);
+
+    if (FreeNormal)
+      delete(ConstructedNormal);
+  }
+  if (result)
+    *out << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
 
   return result;
@@ -314 +345 @@
 /** Calculates the projection of a vector onto another \a *y.
  * \param *y array to second vector
- * \return \f$\langle x, y \rangle\f$
- */
-double Vector::Projection(const Vector *y) const
-{
-  return (ScalarProduct(y));
+ */
+void Vector::ProjectIt(const Vector *y)
+{
+  Vector helper(*y);
+  helper.Scale(-(ScalarProduct(y)));
+  AddVector(&helper);
+};
+
+/** Calculates the projection of a vector onto another \a *y.
+ * \param *y array to second vector
+ * \return Vector
+ */
+Vector Vector::Projection(const Vector *y) const
+{
+  Vector helper(*y);
+  helper.Scale((ScalarProduct(y)/y->NormSquared()));
+
+  return helper;
 };
 
@@ -380 +424 @@
  * @return true - vector is zero, false - vector is not
  */
-bool Vector::IsNull() const
-{
-  return (fabs(x[0]+x[1]+x[2]) < MYEPSILON);
+bool Vector::IsZero() const
+{
+  return (fabs(x[0])+fabs(x[1])+fabs(x[2]) < MYEPSILON);
+};
+
+/** Checks whether vector has length of 1.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsOne() const
+{
+  return (fabs(Norm() - 1.) < MYEPSILON);
+};
+
+/** Checks whether vector is normal to \a *normal.
+ * @return true - vector is normalized, false - vector is not
+ */
+bool Vector::IsNormalTo(const Vector *normal) const
+{
+  if (ScalarProduct(normal) < MYEPSILON)
+    return true;
+  else
+    return false;
 };
 
@@ -392 +455 @@
 {
   double norm1 = Norm(), norm2 = y->Norm();
-  double angle = 1;
+  double angle = -1;
   if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
     angle = this->ScalarProduct(y)/norm1/norm2;
@@ -413 +476 @@
   // normalise this vector with respect to axis
   a.CopyVector(this);
-  a.Scale(Projection(axis));
-  SubtractVector(&a);
+  a.ProjectOntoPlane(axis);
   // construct normal vector
   y.MakeNormalVector(axis,this);
@@ -427 +489 @@
 };
 
+/** Compares vector \a to vector \a b component-wise.
+ * \param a base vector
+ * \param b vector components to add
+ * \return a == b
+ */
+bool operator==(const Vector& a, const Vector& b)
+{
+  bool status = true;
+  for (int i=0;i<NDIM;i++)
+    status = status && (fabs(a.x[i] - b.x[i]) < MYEPSILON);
+  return status;
+};
+
 /** Sums vector \a to this lhs component-wise.
  * \param a base vector
@@ -437 +512 @@
   return a;
 };
+
+/** Subtracts vector \a from this lhs component-wise.
+ * \param a base vector
+ * \param b vector components to add
+ * \return lhs - a
+ */
+Vector& operator-=(Vector& a, const Vector& b)
+{
+  a.SubtractVector(&b);
+  return a;
+};
+
 /** factor each component of \a a times a double \a m.
  * \param a base vector
@@ -461 +548 @@
 };
 
+/** Subtracts vector \a from \b component-wise.
+ * \param a first vector
+ * \param b second vector
+ * \return a - b
+ */
+Vector& operator-(const Vector& a, const Vector& b)
+{
+  Vector *x = new Vector;
+  x->CopyVector(&a);
+  x->SubtractVector(&b);
+  return *x;
+};
+
 /** Factors given vector \a a times \a m.
  * \param a vector
  * \param m factor
- * \return a + b
+ * \return m * a
  */
 Vector& operator*(const Vector& a, const double m)
+{
+  Vector *x = new Vector;
+  x->CopyVector(&a);
+  x->Scale(m);
+  return *x;
+};
+
+/** Factors given vector \a a times \a m.
+ * \param m factor
+ * \param a vector
+ * \return m * a
+ */
+Vector& operator*(const double m, const Vector& a )
 {
   Vector *x = new Vector;
@@ -660 +773 @@
   x2.SubtractVector(y2);
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    cout << Verbose(4) << "Given vectors are linear dependent." << endl;
+    cout << Verbose(4) << "WARNING: Given vectors are linear dependent." << endl;
     return false;
   }
@@ -694 +807 @@
   Zero();
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    cout << Verbose(4) << "Given vectors are linear dependent." << endl;
+    cout << Verbose(4) << "WARNING: Given vectors are linear dependent." << endl;
     return false;
   }
@@ -714 +827 @@
 /** Calculates orthonormal vector to one given vectors.
  * Just subtracts the projection onto the given vector from this vector.
+ * The removed part of the vector is Vector::Projection()
  * \param *x1 vector
  * \return true - success, false - vector is zero
@@ -720 +834 @@
 {
   bool result = false;
-  double factor = y1->Projection(this)/y1->Norm()/y1->Norm();
+  double factor = y1->ScalarProduct(this)/y1->NormSquared();
   Vector x1;
   x1.CopyVector(y1);
@@ -777 +891 @@
 };
 
-/** Determines paramter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
+/** Determines parameter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
  * \param *A first plane vector
  * \param *B second plane vector
@@ -790 +904 @@
 //  cout << "C " << C->Projection(this) << "\t";
 //  cout << endl;
-  return A->Projection(this);
+  return A->ScalarProduct(this);
 };
 
@@ -902 +1016 @@
   for (int i=NDIM;i--;)
     this->x[i] = y->x[i];
+}
+
+/** Copy vector \a y componentwise.
+ * \param y vector
+ */
+void Vector::CopyVector(const Vector y)
+{
+  for (int i=NDIM;i--;)
+    this->x[i] = y.x[i];
 }
 

src/vector.hpp

@@ -5 +5 @@
 
 #include "helpers.hpp"
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_multimin.h>
 
 class Vector;
@@ -24 +27 @@
   double PeriodicDistanceSquared(const Vector *y, const double *cell_size) const;
   double ScalarProduct(const Vector *y) const;
-  double Projection(const Vector *y) const;
   double Norm() const;
   double NormSquared() const;
   double Angle(const Vector *y) const;
-  bool IsNull() const;
+  bool IsZero() const;
+  bool IsOne() const;
+  bool IsNormalTo(const Vector *normal) const;
 
   void AddVector(const Vector *y);
   void SubtractVector(const Vector *y);
   void CopyVector(const Vector *y);
+  void CopyVector(const Vector y);
   void RotateVector(const Vector *y, const double alpha);
   void VectorProduct(const Vector *y);
   void ProjectOntoPlane(const Vector *y);
+  void ProjectIt(const Vector *y);
+  Vector Projection(const Vector *y) const;
   void Zero();
   void One(double one);
@@ -64 +71 @@
 
 ostream & operator << (ostream& ost, const Vector &m);
-//Vector& operator+=(Vector& a, const Vector& b);
-//Vector& operator*=(Vector& a, const double m);
-//Vector& operator*(const Vector& a, const double m);
-//Vector& operator+(const Vector& a, const Vector& b);
+bool operator==(const Vector& a, const Vector& b);
+Vector& operator+=(Vector& a, const Vector& b);
+Vector& operator-=(Vector& a, const Vector& b);
+Vector& operator*=(Vector& a, const double m);
+Vector& operator*(const Vector& a, const double m);
+Vector& operator*(const double m, const Vector& a);
+Vector& operator+(const Vector& a, const Vector& b);
+Vector& operator-(const Vector& a, const Vector& b);
+
+// some algebraic matrix stuff
+#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
+#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix
+
+
 
 #endif /*VECTOR_HPP_*/