Changeset 99c484 for src

Timestamp:

Aug 19, 2009, 12:30:21 PM (15 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

1d9b7aa

Parents:

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

Message:

Merge branch 'VectorUnitTest' into ConcaveHull

Location:

src

Files:

• Makefile.am (modified) (1 diff)
• boundary.cpp (modified) (2 diffs)
• builder.cpp (modified) (2 diffs)
• leastsquaremin.cpp (added)
• leastsquaremin.hpp (added)
• molecules.cpp (modified) (1 diff)
• molecules.hpp (modified) (2 diffs)
• tesselation.cpp (modified) (5 diffs)
• tesselationhelpers.cpp (modified) (2 diffs)
• vector.cpp (modified) (21 diffs)
• vector.hpp (modified) (3 diffs)
• vectorunittest.cpp (added)
• vectorunittest.hpp (added)

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;
       }
@@ -775 +775 @@
       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/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/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 

src/tesselation.cpp

@@ -349 +349 @@
 
   // 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 +737 @@
           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 +897 @@
         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 +957 @@
             TempVector.SubtractVector(Center);
             // make it always point outward
-            if (VirtualNormalVector.Projection(&TempVector) < 0)
+            if (VirtualNormalVector.ScalarProduct(&TempVector) < 0)
               VirtualNormalVector.Scale(-1.);
             // calculate angle
@@ -1530 +1530 @@
         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

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_*/