Changeset 273382

Timestamp:

Apr 13, 2010, 1:22:42 PM (16 years ago)

Author:

Tillmann Crueger <crueger@…>

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, Candidate_v1.7.0, Candidate_v1.7.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:

1bd79e

Parents:

72e7fa

Message:

Prepared interface of Vector Class for transition to VectorComposites

Location:

src

Files:

• Legacy/oldmenu.cpp (modified) (13 diffs)
• Plane.cpp (modified) (7 diffs)
• SingleVector.cpp (modified) (8 diffs)
• analysis_correlation.cpp (modified) (10 diffs)
• atom.cpp (modified) (4 diffs)
• atom_trajectoryparticle.cpp (modified) (1 diff)
• bond.cpp (modified) (2 diffs)
• boundary.cpp (modified) (9 diffs)
• builder.cpp (modified) (3 diffs)
• ellipsoid.cpp (modified) (7 diffs)
• molecule.cpp (modified) (11 diffs)
• molecule.hpp (modified) (1 diff)
• molecule_dynamics.cpp (modified) (4 diffs)
• molecule_fragmentation.cpp (modified) (1 diff)
• molecule_geometry.cpp (modified) (15 diffs)
• molecule_graph.cpp (modified) (1 diff)
• molecule_template.hpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• tesselation.cpp (modified) (63 diffs)
• tesselationhelpers.cpp (modified) (19 diffs)
• test/ActOnAlltest.hpp (modified) (2 diffs)
• unittests/ActOnAllUnitTest.cpp (modified) (1 diff)
• unittests/vectorunittest.cpp (modified) (5 diffs)
• vector.cpp (modified) (33 diffs)
• vector.hpp (modified) (2 diffs)
• vector_ops.cpp (modified) (5 diffs)

src/Legacy/oldmenu.cpp

@@ -103 +103 @@
           dialog->queryVector("Enter reference coordinates.",&x,mol->cell_size,true);
           dialog->queryVector("Enter relative coordinates.",&first->x,mol->cell_size,false);
-          first->x.AddVector(&x);
+          first->x += x;
           dialog->display();
           delete dialog;
@@ -168 +168 @@
           */
           // calc axis vector
-          x.CopyVector(&second->x);
-          x.SubtractVector(&third->x);
+          x= second->x - third->x;
           x.Normalize();
           Log() << Verbose(0) << "x: " << x << endl;
@@ -178 +177 @@
 
           // rotate vector around first angle
-          first->x.CopyVector(&x);
+          first->x = x;
           first->x = RotateVector(first->x,z,b - M_PI);
           Log() << Verbose(0) << "Rotated vector: " << first->x << endl,
           // remove the projection onto the rotation plane of the second angle
-          n.CopyVector(&y);
-          n.Scale(first->x.ScalarProduct(&y));
+          n = y;
+          n.Scale(first->x.ScalarProduct(y));
           Log() << Verbose(0) << "N1: " << n << endl;
-          first->x.SubtractVector(&n);
+          first->x -= n;
           Log() << Verbose(0) << "Subtracted vector: " << first->x << endl;
-          n.CopyVector(&z);
-          n.Scale(first->x.ScalarProduct(&z));
+          n = z;
+          n.Scale(first->x.ScalarProduct(z));
           Log() << Verbose(0) << "N2: " << n << endl;
-          first->x.SubtractVector(&n);
+          first->x -= n;
           Log() << Verbose(0) << "2nd subtracted vector: " << first->x << endl;
 
           // rotate another vector around second angle
-          n.CopyVector(&y);
+          n = y;
           n = RotateVector(n,x,c - M_PI);
           Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 
           // add the two linear independent vectors
-          first->x.AddVector(&n);
+          first->x += n;
           first->x.Normalize();
           first->x.Scale(a);
-          first->x.AddVector(&second->x);
+          first->x += second->x;
 
           Log() << Verbose(0) << "resulting coordinates: " << first->x << endl;
@@ -277 +276 @@
       }
       mol->CenterEdge(&x);  // make every coordinate positive
-      mol->Center.AddVector(&y); // translate by boundary
-      helper.CopyVector(&y);
-      helper.Scale(2.);
-      helper.AddVector(&x);
+      mol->Center += y; // translate by boundary
+      helper = (2*y)+x;
       mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
       break;
@@ -340 +337 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);
-      n.SubtractVector((const Vector *)&second->x);
+      n = first->x - second->x;
       n.Normalize();
       break;
@@ -408 +404 @@
       second = mol->AskAtom("Enter second atom: ");
 
-      n.CopyVector((const Vector *)&first->x);
-      n.SubtractVector((const Vector *)&second->x);
+      n = first->x - second->x;
       n.Normalize();
       break;
@@ -453 +448 @@
         first = second;
         second = first->next;
-        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
+        if (first->x.DistanceSquared(second->x) > tmp1*tmp1) // distance to first above radius ...
           mol->RemoveAtom(first);
       }
@@ -520 +515 @@
         tmp1 = 0.;
         if (first != second) {
-          x.CopyVector((const Vector *)&first->x);
-          x.SubtractVector((const Vector *)&second->x);
+          x = first->x - second->x;
           tmp1 = x.Norm();
         }
@@ -536 +530 @@
       for (int i=NDIM;i--;)
         min[i] = 0.;
-      x.CopyVector((const Vector *)&first->x);
-      x.SubtractVector((const Vector *)&second->x);
+      x = first->x - second->x;
       tmp1 = x.Norm();
       Log() << Verbose(1) << "Distance vector is " << x << "." << "/n"
@@ -549 +542 @@
       third  = mol->AskAtom("Enter last atom: ");
       tmp1 = tmp2 = tmp3 = 0.;
-      x.CopyVector((const Vector *)&first->x);
-      x.SubtractVector((const Vector *)&second->x);
-      y.CopyVector((const Vector *)&third->x);
-      y.SubtractVector((const Vector *)&second->x);
+      x = first->x - second->x;
+      y = third->x - second->x;
       Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
-      Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
+      Log() << Verbose(0) << (acos(x.ScalarProduct(y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
       break;
     case 'd':
@@ -777 +768 @@
       y[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
       for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
-        x.AddVector(&y); // per factor one cell width further
+        x += y; // per factor one cell width further
         for (int k=count;k--;) { // go through every atom of the original cell
           first = World::getInstance().createAtom(); // create a new body
-          first->x.CopyVector(vectors[k]);  // use coordinate of original atom
-          first->x.AddVector(&x);     // translate the coordinates
+          first->x = (*vectors[k]) + x;  // use coordinate of original atom
           first->type = Elements[k];  // insert original element
           mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
@@ -793 +783 @@
       // correct cell size
       if (axis < 0) { // if sign was negative, we have to translate everything
-        x.Zero();
-        x.AddVector(&y);
+        x = y;
         x.Scale(-(faktor-1));
         mol->Translate(&x);
@@ -894 +883 @@
         dialog->display();
         delete dialog;
-        mol->Center.AddVector((const Vector *)&x);
+        mol->Center += x;
      }
      break;

src/Plane.cpp

@@ -20 +20 @@
   normalVector(new Vector())
 {
-  Vector x1, x2;
-
-  x1.CopyVector(&y1);
-  x1.SubtractVector(&y2);
-  x2.CopyVector(&y3);
-  x2.SubtractVector(&y2);
-  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
+  Vector x1 = y1 -y2;
+  Vector x2 = y3 -y2;
+  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(x2)) < MYEPSILON)) {
     throw LinearDependenceException(__FILE__,__LINE__);
   }
@@ -41 +37 @@
   normalVector->Normalize();
 
-  offset=normalVector->ScalarProduct(&y1);
+  offset=normalVector->ScalarProduct(y1);
 }
 /**
@@ -51 +47 @@
     offset(_offset)
 {
-  Vector x1,x2;
-  x1.CopyVector(&y1);
-  x2.CopyVector(&y2);
-  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
+  Vector x1 = y1;
+  Vector x2 = y2;
+  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(x2)) < MYEPSILON)) {
     throw LinearDependenceException(__FILE__,__LINE__);
   }
@@ -84 +79 @@
     normalVector(new Vector(_normalVector))
 {
-  offset = normalVector->ScalarProduct(&_offsetVector);
+  offset = normalVector->ScalarProduct(_offsetVector);
 }
 
@@ -126 +121 @@
   Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
   //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
-  double factor1 = Direction.ScalarProduct(normalVector.get());
+  double factor1 = Direction.ScalarProduct(*normalVector.get());
   if (fabs(factor1) < MYEPSILON) { // Uniqueness: line parallel to plane?
     Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
@@ -132 +127 @@
   }
 
-  double factor2 = Origin.ScalarProduct(normalVector.get());
+  double factor2 = Origin.ScalarProduct(*normalVector.get());
   if (fabs(factor2-offset) < MYEPSILON) { // Origin is in-plane
     Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
@@ -147 +142 @@
 
   // test whether resulting vector really is on plane
-  ASSERT(fabs(res.ScalarProduct(normalVector.get()) - offset) < MYEPSILON,
+  ASSERT(fabs(res.ScalarProduct((*normalVector.get())) - offset) < MYEPSILON,
          "Calculated line-Plane intersection does not lie on plane.");
   return res;

src/SingleVector.cpp

@@ -119 +119 @@
         TranslationVector.MatrixMultiplication(matrix);
         // add onto the original vector to compare with
-        Shiftedy.CopyVector(y);
-        Shiftedy.AddVector(&TranslationVector);
+        Shiftedy = y + TranslationVector;
         // get distance and compare with minimum so far
         tmp = Distance(Shiftedy);
@@ -157 +156 @@
         TranslationVector.MatrixMultiplication(matrix);
         // add onto the original vector to compare with
-        Shiftedy.CopyVector(y);
-        Shiftedy.AddVector(&TranslationVector);
+        Shiftedy = y + TranslationVector;
         // get distance and compare with minimum so far
         tmp = DistanceSquared(Shiftedy);
@@ -180 +178 @@
 //  Output(out);
 //  Log() << Verbose(0) << endl;
-  TestVector.CopyVector(this);
+  TestVector = *this;
   TestVector.InverseMatrixMultiplication(matrix);
   for(int i=NDIM;i--;) { // correct periodically
@@ -190 +188 @@
   }
   TestVector.MatrixMultiplication(matrix);
-  CopyVector(&TestVector);
+  CopyVector(TestVector);
 //  Log() << Verbose(2) << "New corrected vector is: ";
 //  Output(out);
@@ -250 +248 @@
 
   // first create part that is orthonormal to PlaneNormal with withdraw
-  temp = *this - PlaneOffset;
+  temp = (*this )- PlaneOffset;
   temp.MakeNormalTo(PlaneNormal);
   temp.Scale(-1.);
   // then add connecting vector from plane to point
-  temp += *this-PlaneOffset;
-  double sign = temp.ScalarProduct(&PlaneNormal);
+  temp += (*this)-PlaneOffset;
+  double sign = temp.ScalarProduct(PlaneNormal);
   if (fabs(sign) > MYEPSILON)
     sign /= fabs(sign);
@@ -271 +269 @@
   Vector helper = y;
   helper.Scale(-(ScalarProduct(y)));
-  AddVector(&helper);
+  AddVector(helper);
 };
 
@@ -556 +554 @@
 {
   bool result = false;
-  double factor = y1.ScalarProduct(this)/y1.NormSquared();
+  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
   Vector x1;
-  x1.CopyVector(&y1);
-  x1.Scale(factor);
-  SubtractVector(&x1);
+  x1 = factor * y1;
+  SubtractVector(x1);
   for (int i=NDIM;i--;)
     result = result || (fabs(x[i]) > MYEPSILON);
@@ -618 +615 @@
 bool SingleVector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  Vector a;
-  a.CopyVector(this);
-  a.SubtractVector(&offset);
+  Vector a = (*this) - offset;
   a.InverseMatrixMultiplication(parallelepiped);
   bool isInside = true;

src/analysis_correlation.cpp

@@ -55 +55 @@
                 if (Walker->nr < OtherWalker->nr)
                   if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
+                    distance = Walker->node->PeriodicDistance(*(OtherWalker->node), (*MolWalker)->cell_size);
                     //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
@@ -104 +104 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type1 == NULL) || (Walker->type == type1)) {
-          periodicX.CopyVector(Walker->node);
+          periodicX = *(Walker->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -110 +110 @@
             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX.Init(n[0], n[1], n[2]);
-                checkX.AddVector(&periodicX);
+                checkX = Vector(n[0], n[1], n[2]) + periodicX;
                 checkX.MatrixMultiplication(FullMatrix);
                 for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
@@ -122 +121 @@
                       if (Walker->nr < OtherWalker->nr)
                         if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                          periodicOtherX.CopyVector(OtherWalker->node);
+                          periodicOtherX = *(OtherWalker->node);
                           periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
                           // go through every range in xyz and get distance
@@ -128 +127 @@
                             for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                               for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
-                                checkOtherX.Init(Othern[0], Othern[1], Othern[2]);
-                                checkOtherX.AddVector(&periodicOtherX);
+                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
                                 checkOtherX.MatrixMultiplication(FullMatrix);
-                                distance = checkX.Distance(&checkOtherX);
+                                distance = checkX.Distance(checkOtherX);
                                 //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                                 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
@@ -174 +172 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
+          distance = Walker->node->PeriodicDistance(*point, (*MolWalker)->cell_size);
           Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
@@ -216 +214 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          periodicX.CopyVector(Walker->node);
+          periodicX = *(Walker->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -222 +220 @@
             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX.Init(n[0], n[1], n[2]);
-                checkX.AddVector(&periodicX);
+                checkX = Vector(n[0], n[1], n[2]) + periodicX;
                 checkX.MatrixMultiplication(FullMatrix);
-                distance = checkX.Distance(point);
+                distance = checkX.Distance(*point);
                 Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
                 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
@@ -320 +317 @@
         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          periodicX.CopyVector(Walker->node);
+          periodicX = *(Walker->node);
           periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
           // go through every range in xyz and get distance
@@ -327 +324 @@
             for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
               for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX.Init(n[0], n[1], n[2]);
-                checkX.AddVector(&periodicX);
+                checkX = Vector(n[0], n[1], n[2]) + periodicX;
                 checkX.MatrixMultiplication(FullMatrix);
                 triangle = Surface->FindClosestTriangleToVector(&checkX, LC);

src/atom.cpp

@@ -33 +33 @@
 {
   type = pointer->type;  // copy element of atom
-  x.CopyVector(&pointer->x); // copy coordination
-  v.CopyVector(&pointer->v); // copy velocity
+  x = pointer->x; // copy coordination
+  v = pointer->v; // copy velocity
   FixedIon = pointer->FixedIon;
   node = &x;
@@ -46 +46 @@
   res->sort = &nr;
   res->type = type;
-  res->x.CopyVector(&this->x);
-  res->v.CopyVector(&this->v);
+  res->x = this->x;
+  res->v = this->v;
   res->FixedIon = FixedIon;
   res->node = &x;
@@ -252 +252 @@
 double atom::DistanceSquaredToVector(const Vector &origin) const
 {
-  return origin.DistanceSquared(&x);
+  return origin.DistanceSquared(x);
 };
 
@@ -261 +261 @@
 double atom::DistanceToVector(const Vector &origin) const
 {
-  return origin.Distance(&x);
+  return origin.Distance(x);
 };
 

src/atom_trajectoryparticle.cpp

@@ -49 +49 @@
   // set forces
   for (int i=NDIM;i++;)
-    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
+    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(Sprinter->Trajectory.R.at(endstep)));
 };
 

src/bond.cpp

@@ -119 +119 @@
 double bond::GetDistance() const
 {
-  return (leftatom->node->Distance(rightatom->node));
+  return (leftatom->node->Distance(*rightatom->node));
 };
 
@@ -127 +127 @@
 double bond::GetDistanceSquared() const
 {
-  return (leftatom->node->DistanceSquared(rightatom->node));
+  return (leftatom->node->DistanceSquared(*rightatom->node));
 };

src/boundary.cpp

@@ -78 +78 @@
               if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                 Neighbour = BoundaryPoints[axis].begin();
-              DistanceVector.CopyVector(&runner->second.second->x);
-              DistanceVector.SubtractVector(&Neighbour->second.second->x);
+              DistanceVector = runner->second.second->x - Neighbour->second.second->x;
               do { // seek for neighbour pair where it flips
                   OldComponent = DistanceVector[Othercomponent];
@@ -85 +84 @@
                   if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
                     Neighbour = BoundaryPoints[axis].begin();
-                  DistanceVector.CopyVector(&runner->second.second->x);
-                  DistanceVector.SubtractVector(&Neighbour->second.second->x);
+                  DistanceVector = runner->second.second->x - Neighbour->second.second->x;
                   //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
@@ -98 +96 @@
                   //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;
                   // now we have found the pair: Neighbour and OtherNeighbour
-                  OtherVector.CopyVector(&runner->second.second->x);
-                  OtherVector.SubtractVector(&OtherNeighbour->second.second->x);
+                  OtherVector = runner->second.second->x - OtherNeighbour->second.second->x;
                   //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;
                   //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
@@ -170 +167 @@
     while (Walker->next != mol->end) {
       Walker = Walker->next;
-      ProjectedVector.CopyVector(&Walker->x);
-      ProjectedVector.SubtractVector(MolCenter);
-      ProjectedVector.ProjectOntoPlane(&AxisVector);
+      ProjectedVector = Walker->x - (*MolCenter);
+      ProjectedVector.ProjectOntoPlane(AxisVector);
 
       // correct for negative side
       const double radius = ProjectedVector.NormSquared();
       if (fabs(radius) > MYEPSILON)
-        angle = ProjectedVector.Angle(&AngleReferenceVector);
+        angle = ProjectedVector.Angle(AngleReferenceVector);
       else
         angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues
 
       //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
+      if (ProjectedVector.ScalarProduct(AngleReferenceNormalVector) > 0) {
         angle = 2. * M_PI - angle;
       }
@@ -197 +193 @@
           Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl;
         } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
-          helper.CopyVector(&Walker->x);
-          helper.SubtractVector(MolCenter);
+          helper = Walker->x - (*MolCenter);
           const double oldhelperNorm = helper.NormSquared();
-          helper.CopyVector(&BoundaryTestPair.first->second.second->x);
-          helper.SubtractVector(MolCenter);
+          helper = BoundaryTestPair.first->second.second->x - (*MolCenter);
           if (helper.NormSquared() < oldhelperNorm) {
             BoundaryTestPair.first->second.second = Walker;
@@ -251 +245 @@
         {
           Vector SideA, SideB, SideC, SideH;
-          SideA.CopyVector(&left->second.second->x);
-          SideA.SubtractVector(MolCenter);
-          SideA.ProjectOntoPlane(&AxisVector);
+          SideA = left->second.second->x - (*MolCenter);
+          SideA.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideA: " << SideA << endl;
 
-          SideB.CopyVector(&right->second.second->x);
-          SideB.SubtractVector(MolCenter);
-          SideB.ProjectOntoPlane(&AxisVector);
+          SideB = right->second.second->x -(*MolCenter);
+          SideB.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideB: " << SideB << endl;
 
-          SideC.CopyVector(&left->second.second->x);
-          SideC.SubtractVector(&right->second.second->x);
-          SideC.ProjectOntoPlane(&AxisVector);
+          SideC = left->second.second->x - right->second.second->x;
+          SideC.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideC: " << SideC << endl;
 
-          SideH.CopyVector(&runner->second.second->x);
-          SideH.SubtractVector(MolCenter);
-          SideH.ProjectOntoPlane(&AxisVector);
+          SideH = runner->second.second->x -(*MolCenter);
+          SideH.ProjectOntoPlane(AxisVector);
           //          Log() << Verbose(1) << "SideH: " << SideH << endl;
 
@@ -277 +267 @@
           const double h = SideH.Norm();
           // calculate the angles
-          const double alpha = SideA.Angle(&SideH);
-          const double beta = SideA.Angle(&SideC);
-          const double gamma = SideB.Angle(&SideH);
-          const double delta = SideC.Angle(&SideH);
+          const double alpha = SideA.Angle(SideH);
+          const double beta = SideA.Angle(SideC);
+          const double gamma = SideB.Angle(SideH);
+          const double delta = SideC.Angle(SideH);
           const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.);
           //Log() << Verbose(1) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
@@ -629 +619 @@
   for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++)
     { // go through every triangle, calculate volume of its pyramid with CoG as peak
-      x.CopyVector(runner->second->endpoints[0]->node->node);
-      x.SubtractVector(runner->second->endpoints[1]->node->node);
-      y.CopyVector(runner->second->endpoints[0]->node->node);
-      y.SubtractVector(runner->second->endpoints[2]->node->node);
-      const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
-      const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(runner->second->endpoints[2]->node->node));
-      const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(runner->second->endpoints[1]->node->node));
+      x = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[1]->node->node);
+      y = (*runner->second->endpoints[0]->node->node) - (*runner->second->endpoints[2]->node->node);
+      const double a = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
+      const double b = sqrt(runner->second->endpoints[0]->node->node->DistanceSquared(*runner->second->endpoints[2]->node->node));
+      const double c = sqrt(runner->second->endpoints[2]->node->node->DistanceSquared(*runner->second->endpoints[1]->node->node));
       const double G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
       x = Plane(*(runner->second->endpoints[0]->node->node),
                 *(runner->second->endpoints[1]->node->node),
                 *(runner->second->endpoints[2]->node->node)).getNormal();
-      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
+      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(x));
       const double h = x.Norm(); // distance of CoG to triangle
       const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
@@ -917 +905 @@
             if (DoRandomRotation)
               CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
-            CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
-            CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
+            CopyAtoms[Walker->nr]->x += AtomTranslations;
+            CopyAtoms[Walker->nr]->x += FillerTranslations;
+            CopyAtoms[Walker->nr]->x += CurrentPosition;
 
             // insert into Filling

src/builder.cpp

@@ -1817 +1817 @@
                     first = second;
                     second = first->next;
-                    if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
+                    if (first->x.DistanceSquared(third->x) > tmp1*tmp1) // distance to first above radius ...
                       mol->RemoveAtom(first);
                   }
@@ -2096 +2096 @@
                     y[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
                     for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
-                      x.AddVector(&y); // per factor one cell width further
+                      x += y; // per factor one cell width further
                       for (int k=count;k--;) { // go through every atom of the original cell
                         first = World::getInstance().createAtom(); // create a new body
-                        first->x.CopyVector(vectors[k]);  // use coordinate of original atom
-                        first->x.AddVector(&x);      // translate the coordinates
+                        first->x = (*vectors[k]) + x;
                         first->type = Elements[k];  // insert original element
                         mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
@@ -2110 +2109 @@
                     // correct cell size
                     if (axis < 0) { // if sign was negative, we have to translate everything
-                      x.Zero();
-                      x.AddVector(&y);
-                      x.Scale(-(faktor-1));
+                      x =(-(faktor-1)) * y;
                       mol->Translate(&x);
                     }

src/ellipsoid.cpp

@@ -42 +42 @@
 
   // 1. translate coordinate system so that ellipsoid center is in origin
-  helper.CopyVector(&x);
-  helper.SubtractVector(&EllipsoidCenter);
-  RefPoint.CopyVector(&helper);
+  RefPoint = helper = x - EllipsoidCenter;
   //Log() << Verbose(4) << "Translated given point is at " << RefPoint << "." << endl;
 
@@ -86 +84 @@
 
   // 5. determine distance between backtransformed point and x
-  distance = RefPoint.DistanceSquared(&helper);
+  distance = RefPoint.DistanceSquared(helper);
   //Log() << Verbose(4) << "Squared distance between intersection and RefPoint is " << distance << "." << endl;
 
@@ -304 +302 @@
                 Candidate = (*Runner);
                 Log() << Verbose(2) << "Current picked node is " << **Runner << " with index " << index << "." << endl;
-                x[PointsPicked++].CopyVector(Candidate->node);    // we have one more atom picked
+                x[PointsPicked++] = (*Candidate->node);    // we have one more atom picked
                 current++;    // next pre-picked atom
               }
@@ -350 +348 @@
       //Log() << Verbose(3) << "Current node is " << *Runner->second->node << " with " << value << " ... " << threshold << ": ";
       if (value > threshold) {
-        x[PointsPicked].CopyVector(Runner->second->node->node);
+        x[PointsPicked] = (*Runner->second->node->node);
         PointsPicked++;
         //Log() << Verbose(0) << "IN." << endl;
@@ -387 +385 @@
   Center.Zero();
   for (PointMap::iterator Runner = T->PointsOnBoundary.begin(); Runner != T->PointsOnBoundary.end(); Runner++)
-    Center.AddVector(Runner->second->node->node);
+    Center += (*Runner->second->node->node);
   Center.Scale(1./T->PointsOnBoundaryCount);
   Log() << Verbose(1) << "Center is at " << Center << "." << endl;
@@ -405 +403 @@
     // calculate some sensible starting values for parameter fit
     MaxDistance = 0.;
-    MinDistance = x[0].ScalarProduct(&x[0]);
+    MinDistance = x[0].ScalarProduct(x[0]);
     for (int i=0;i<N;i++) {
-      distance = x[i].ScalarProduct(&x[i]);
+      distance = x[i].ScalarProduct(x[i]);
       if (distance > MaxDistance)
         MaxDistance = distance;
@@ -414 +412 @@
     }
     //Log() << Verbose(2) << "MinDistance " << MinDistance << ", MaxDistance " << MaxDistance << "." << endl;
-    EllipsoidCenter.CopyVector(&Center);  // use Center of Gravity as initial center of ellipsoid
+    EllipsoidCenter = Center;  // use Center of Gravity as initial center of ellipsoid
     for (int i=0;i<3;i++)
       EllipsoidAngle[i] = 0.;

src/molecule.cpp

@@ -220 +220 @@
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
-  InBondvector.CopyVector(&TopReplacement->x);
-  InBondvector.SubtractVector(&TopOrigin->x);
+  InBondvector = TopReplacement->x - TopOrigin->x;
   bondlength = InBondvector.Norm();
 
@@ -240 +239 @@
     matrix = ReturnFullMatrixforSymmetric(cell_size);
     Orthovector1.MatrixMultiplication(matrix);
-    InBondvector.SubtractVector(&Orthovector1); // subtract just the additional translation
+    InBondvector.SubtractVector(Orthovector1); // subtract just the additional translation
     Free(&matrix);
     bondlength = InBondvector.Norm();
@@ -265 +264 @@
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
       FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
-      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
       if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
@@ -274 +273 @@
       }
       InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
-      FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ...
-      FirstOtherAtom->x.AddVector(&InBondvector);  // ... and add distance vector to replacement atom
+      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
+      FirstOtherAtom->x = InBondvector;  // ... and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
@@ -316 +315 @@
         }
       } else {
-        Orthovector1.GetOneNormalVector(&InBondvector);
+        Orthovector1.GetOneNormalVector(InBondvector);
       }
       //Log() << Verbose(3)<< "Orthovector1: ";
@@ -331 +330 @@
       FirstOtherAtom->type = elemente->FindElement(1);
       SecondOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      SecondOtherAtom->v = TopReplacement->v; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
@@ -388 +387 @@
       SecondOtherAtom->type = elemente->FindElement(1);
       ThirdOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      SecondOtherAtom->v = TopReplacement->v; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
-      ThirdOtherAtom->v.CopyVector(&TopReplacement->v); // copy velocity
+      ThirdOtherAtom->v = TopReplacement->v; // copy velocity
       ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -399 +398 @@
 
       // we need to vectors orthonormal the InBondvector
-      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector);
+      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(InBondvector);
 //      Log() << Verbose(3) << "Orthovector1: ";
 //      Orthovector1.Output(out);
@@ -426 +425 @@
       factors[1] = f;
       factors[2] = 0.;
-      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
+      FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
-      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
+      SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
-      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
+      ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
 
       // rescale each to correct BondDistance
@@ -439 +438 @@
 
       // and relative to *origin atom
-      FirstOtherAtom->x.AddVector(&TopOrigin->x);
-      SecondOtherAtom->x.AddVector(&TopOrigin->x);
-      ThirdOtherAtom->x.AddVector(&TopOrigin->x);
+      FirstOtherAtom->x += TopOrigin->x;
+      SecondOtherAtom->x += TopOrigin->x;
+      ThirdOtherAtom->x += TopOrigin->x;
 
       // ... and add to molecule
@@ -1022 +1021 @@
     Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl;
     Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl;
-    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
+    if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {
       Log() << Verbose(4) << "Centers of gravity don't match." << endl;
       result = false;

src/molecule.hpp

@@ -145 +145 @@
   template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T), T t ) const;
   template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T) const, T t ) const;
+  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&), T &t ) const;
+  template <typename res, typename T> void ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const;
   template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const;
   template <typename res, typename T, typename U> void ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const;

src/molecule_dynamics.cpp

@@ -39 +39 @@
     // determine normalized trajectories direction vector (n1, n2)
     Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-    trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
-    trajectory1.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
+    trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);
     trajectory1.Normalize();
     Norm1 = trajectory1.Norm();
     Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
-    trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
-    trajectory2.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));
+    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     trajectory2.Normalize();
     Norm2 = trajectory1.Norm();
     // check whether either is zero()
     if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
+      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.startstep));
     } else if (Norm1 < MYEPSILON) {
       Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
-      trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy first offset
-      trajectory1.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));  // subtract second offset
-      trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
-      trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
+      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
+      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
+      trajectory1 -= trajectory2;   // project the part in norm direction away
       tmp = trajectory1.Norm();  // remaining norm is distance
     } else if (Norm2 < MYEPSILON) {
       Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
-      trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy second offset
-      trajectory2.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));  // subtract first offset
-      trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
-      trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
+      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
+      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
+      trajectory2 -= trajectory1;   // project the part in norm direction away
       tmp = trajectory2.Norm();  // remaining norm is distance
-    } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
+    } else if ((fabs(trajectory1.ScalarProduct(trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
   //        Log() << Verbose(0) << trajectory1;
   //        Log() << Verbose(0) << " and ";
   //        Log() << Verbose(0) << trajectory2;
-      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
+      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.startstep));
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
     } else { // determine distance by finding minimum distance
@@ -97 +93 @@
   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
       // test whether we really have the intersection (by checking on c_1 and c_2)
-      TestVector.CopyVector(&Runner->Trajectory.R.at(Params.startstep));
+      trajectory1.Scale(gsl_vector_get(x,0));
       trajectory2.Scale(gsl_vector_get(x,1));
-      TestVector.AddVector(&trajectory2);
       normal.Scale(gsl_vector_get(x,2));
-      TestVector.AddVector(&normal);
-      TestVector.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
-      trajectory1.Scale(gsl_vector_get(x,0));
-      TestVector.SubtractVector(&trajectory1);
+      TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
+                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
       if (TestVector.Norm() < MYEPSILON) {
   //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
@@ -173 +166 @@
     // first term: distance to target
     Runner = Params.PermutationMap[Walker->nr];   // find target point
-    tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)));
+    tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.endstep)));
     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += Params.PenaltyConstants[0] * tmp;
@@ -232 +225 @@
     while(Runner->next != mol->end) {
       Runner = Runner->next;
-      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)), Runner) );
+      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.endstep)), Runner) );
     }
   }

src/molecule_fragmentation.cpp

@@ -1705 +1705 @@
         //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
         ColorList[Walker->nr] = black;    // mark as explored
-        Walker->x.AddVector(&Translationvector); // translate
+        Walker->x += Translationvector; // translate
         for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
           if ((*Runner) != Binder) {

src/molecule_geometry.cpp

@@ -30 +30 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, Center);
+  ActOnAllVectors( &Vector::SubtractVector, *Center);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
@@ -86 +86 @@
 //    Log() << Verbose(0) << endl;
     min->Scale(-1.);
-    max->AddVector(min);
+    (*max) += (*min);
     Translate(min);
     Center.Zero();
@@ -109 +109 @@
       ptr = ptr->next;
       Num++;
-      Center.AddVector(&ptr->x);
+      Center += ptr->x;
     }
     Center.Scale(-1./Num); // divide through total number (and sign for direction)
@@ -133 +133 @@
       ptr = ptr->next;
       Num += 1.;
-      tmp.CopyVector(&ptr->x);
-      a->AddVector(&tmp);
+      tmp = ptr->x;
+      (*a) += tmp;
     }
     a->Scale(1./Num); // divide through total mass (and sign for direction)
@@ -158 +158 @@
       ptr = ptr->next;
       Num += ptr->type->mass;
-      tmp.CopyVector(&ptr->x);
-      tmp.Scale(ptr->type->mass);  // scale by mass
-      a->AddVector(&tmp);
+      tmp = ptr->type->mass * ptr->x;
+      (*a) += tmp;
     }
     a->Scale(-1./Num); // divide through total mass (and sign for direction)
@@ -186 +185 @@
 void molecule::CenterAtVector(Vector *newcenter)
 {
-  Center.CopyVector(newcenter);
+  Center = *newcenter;
 };
 
@@ -215 +214 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Translate(trans);
-    ptr->x.Translate(trans);
+      ptr->Trajectory.R.at(j).Translate(*trans);
+    ptr->x.Translate(*trans);
   }
 };
@@ -230 +229 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, trans);
+  ActOnAllVectors( &Vector::SubtractVector, *trans);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
@@ -243 +242 @@
 void molecule::Mirror(const Vector *n)
 {
-  ActOnAllVectors( &Vector::Mirror, n );
+  ActOnAllVectors( &Vector::Mirror, *n );
 };
 
@@ -266 +265 @@
       if (Walker->type->Z != 1) {
 #endif
-        Testvector.CopyVector(&Walker->x);
+        Testvector = Walker->x;
         Testvector.MatrixMultiplication(inversematrix);
         Translationvector.Zero();
@@ -283 +282 @@
             }
         }
-        Testvector.AddVector(&Translationvector);
+        Testvector += Translationvector;
         Testvector.MatrixMultiplication(matrix);
-        Center.AddVector(&Testvector);
+        Center += Testvector;
         Log() << Verbose(1) << "vector is: " << Testvector << endl;
 #ifdef ADDHYDROGEN
@@ -291 +290 @@
         for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
           if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
-            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
+            Testvector = (*Runner)->GetOtherAtom(Walker)->x;
             Testvector.MatrixMultiplication(inversematrix);
-            Testvector.AddVector(&Translationvector);
+            Testvector += Translationvector;
             Testvector.MatrixMultiplication(matrix);
-            Center.AddVector(&Testvector);
+            Center += Testvector;
             Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
           }
@@ -329 +328 @@
   while (ptr->next != end) {
     ptr = ptr->next;
-    Vector x;
-    x.CopyVector(&ptr->x);
+    Vector x = ptr->x;
     //x.SubtractVector(CenterOfGravity);
     InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
@@ -381 +379 @@
     while (ptr->next != end) {
       ptr = ptr->next;
-      Vector x;
-      x.CopyVector(&ptr->x);
+      Vector x = ptr->x;
       //x.SubtractVector(CenterOfGravity);
       InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
@@ -492 +489 @@
     ptr = ptr->next;
     if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
-      c.CopyVector(&ptr->x);  // copy vector to temporary one
-      c.SubtractVector(&a);   // subtract offset vector
-      t = c.ScalarProduct(&b);           // get direction parameter
-      d.CopyVector(&b);       // and create vector
-      d.Scale(&t);
-      c.SubtractVector(&d);   // ... yielding distance vector
-      res += d.ScalarProduct((const Vector *)&d);        // add squared distance
+      c = ptr->x - a;
+      t = c.ScalarProduct(b);           // get direction parameter
+      d = t*b;       // and create vector
+      c -= d;   // ... yielding distance vector
+      res += d.ScalarProduct(d);        // add squared distance
     }
   }

src/molecule_graph.cpp

@@ -165 +165 @@
                           //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                           (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
-                          const double distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
+                          const double distance = OtherWalker->x.PeriodicDistanceSquared((Walker->x), cell_size);
                           const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
                           if ((OtherWalker->father->nr > Walker->father->nr) && (status)) { // create bond if distance is smaller

src/molecule_template.hpp

@@ -55 +55 @@
   }
 };
+template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&), T &t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)(t);
+  }
+};
+template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    ((Walker->node)->*f)(t);
+  }
+};
 // two arguments
 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const

src/moleculelist.cpp

@@ -161 +161 @@
         Walker = Walker->next;
         counts[Walker->type->getNumber()]++;
-        if (Walker->x.DistanceSquared(&Origin) > size)
-          size = Walker->x.DistanceSquared(&Origin);
+        if (Walker->x.DistanceSquared(Origin) > size)
+          size = Walker->x.DistanceSquared(Origin);
       }
       // output Index, Name, number of atoms, chemical formula
@@ -501 +501 @@
           if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
             // 4. evaluate the morse potential for each matrix component and add up
-            distance = Runner->x.Distance(&Walker->x);
+            distance = Runner->x.Distance(Walker->x);
             //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
             for (int k = 0; k < a; k++) {

src/tesselation.cpp

@@ -234 +234 @@
   // have a normal vector on the base line pointing outwards
   //Log() << Verbose(0) << "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);
-  BaseLineCenter.Scale(1./2.);
-  BaseLine.CopyVector(endpoints[0]->node->node);
-  BaseLine.SubtractVector(endpoints[1]->node->node);
+  BaseLineCenter = (1./2.)*((*endpoints[0]->node->node) + (*endpoints[1]->node->node));
+  BaseLine = (*endpoints[0]->node->node) - (*endpoints[1]->node->node);
   //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;
 
@@ -248 +245 @@
   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);
-    NormalCheck.Scale(sign);
+    NormalCheck += runner->second->NormalVector;
+    NormalCheck *= sign;
     sign = -sign;
     if (runner->second->NormalVector.NormSquared() > MYEPSILON)
-      BaseLineNormal.CopyVector(&runner->second->NormalVector);   // yes, copy second on top of first
+      BaseLineNormal = runner->second->NormalVector;   // yes, copy second on top of first
     else {
       eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl;
@@ -259 +256 @@
     if (node != NULL) {
       //Log() << Verbose(0) << "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] = (*node->node->node) - BaseLineCenter;
       helper[i].MakeNormalTo(BaseLine);  // we want to compare the triangle's heights' angles!
       //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;
@@ -410 +406 @@
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (NormalVector.ScalarProduct(&OtherVector) > 0.)
+  if (NormalVector.ScalarProduct(OtherVector) > 0.)
     NormalVector.Scale(-1.);
   Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl;
@@ -446 +442 @@
   Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl;
 
-  if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) {
+  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) {
+  }   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) {
+  }   else if (Intersection->DistanceSquared(*endpoints[2]->node->node) < MYEPSILON) {
     Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl;
     return true;
@@ -464 +460 @@
                                                     *(endpoints[(i+2)%3]->node->node),
                                                     *Intersection);
-      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();
+      helper = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
+      CrossPoint -= (*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)) {
@@ -510 +505 @@
   }
   catch (LinearDependenceException &excp) {
-    ClosestPoint->CopyVector(x);
+    (*ClosestPoint) = (*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);
+  Vector InPlane = (*x) - (*ClosestPoint); // points from plane intersection to straight-down point
+  InPlane.ProjectOntoPlane(NormalVector);
+  InPlane += *ClosestPoint;
 
   Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl;
@@ -532 +525 @@
   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);
+    Direction = (*endpoints[(i+1)%3]->node->node) - (*endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
 
     CrossPoint[i] = Plane(Direction, InPlane).GetIntersection(*(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();
+    CrossDirection[i] = CrossPoint[i] - InPlane;
+    CrossPoint[i] -= (*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
+      CrossPoint[i] += (*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);
+      const double distance = CrossPoint[i].DistanceSquared(*x);
       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
         ShortestDistance = distance;
-        ClosestPoint->CopyVector(&CrossPoint[i]);
+        (*ClosestPoint) = CrossPoint[i];
       }
     } else
@@ -556 +547 @@
   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]);;
+    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);
+    (*ClosestPoint) = 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);
+      const double distance = x->DistanceSquared(*endpoints[i]->node->node);
       if ((ShortestDistance < 0.) || (ShortestDistance > distance)) {
         ShortestDistance = distance;
-        ClosestPoint->CopyVector(endpoints[i]->node->node);
+        (*ClosestPoint) = (*endpoints[i]->node->node);
       }
     }
@@ -687 +678 @@
   center->Zero();
   for(int i=0;i<3;i++)
-    center->AddVector(endpoints[i]->node->node);
+    (*center) += (*endpoints[i]->node->node);
   center->Scale(1./3.);
   Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl;
@@ -755 +746 @@
     for (int i=0;i<3;i++) // increase each of them
       Runner[i]++;
-    TotalNormal->AddVector(&TemporaryNormal);
+    (*TotalNormal) += TemporaryNormal;
   }
   TotalNormal->Scale(1./(double)counter);
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
-  if (TotalNormal->ScalarProduct(&OtherVector) > 0.)
+  if (TotalNormal->ScalarProduct(OtherVector) > 0.)
     TotalNormal->Scale(-1.);
   Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl;
@@ -777 +768 @@
   int counter = 0;
   for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {
-    center->AddVector((*Runner)->node->node);
+    (*center) += (*(*Runner)->node->node);
     counter++;
   }
@@ -1028 +1019 @@
 {
         Info FunctionInfo(__func__);
-  OptCenter.CopyVector(&OptCandidateCenter);
-  OtherOptCenter.CopyVector(&OtherOptCandidateCenter);
+  OptCenter = OptCandidateCenter;
+  OtherOptCenter = OtherOptCandidateCenter;
 };
 
@@ -1105 +1096 @@
   int num=0;
   for (GoToFirst(); (!IsEnd()); GoToNext()) {
-    Center->AddVector(GetPoint()->node);
+    (*Center) += (*GetPoint()->node);
     num++;
   }
@@ -1217 +1208 @@
           for (; C != PointsOnBoundary.end(); C++)
             {
-              tmp = A->second->node->node->DistanceSquared(B->second->node->node);
+              tmp = A->second->node->node->DistanceSquared(*B->second->node->node);
               distance = tmp * tmp;
-              tmp = A->second->node->node->DistanceSquared(C->second->node->node);
+              tmp = A->second->node->node->DistanceSquared(*C->second->node->node);
               distance += tmp * tmp;
-              tmp = B->second->node->node->DistanceSquared(C->second->node->node);
+              tmp = B->second->node->node->DistanceSquared(*C->second->node->node);
               distance += tmp * tmp;
               DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
@@ -1255 +1246 @@
             continue;
           // 4a. project onto plane vector
-          TrialVector.CopyVector(checker->second->node->node);
-          TrialVector.SubtractVector(A->second->node->node);
-          distance = TrialVector.ScalarProduct(&PlaneVector);
+          TrialVector = (*checker->second->node->node) - (*A->second->node->node);
+          distance = TrialVector.ScalarProduct(PlaneVector);
           if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
             continue;
@@ -1285 +1275 @@
             }
           // 4d. Check whether the point is inside the triangle (check distance to each node
-          tmp = checker->second->node->node->DistanceSquared(A->second->node->node);
+          tmp = checker->second->node->node->DistanceSquared(*A->second->node->node);
           int innerpoint = 0;
           if ((tmp < A->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node)) && (tmp
+              *baseline->second.first->second->node->node)) && (tmp
               < A->second->node->node->DistanceSquared(
-                  baseline->second.second->second->node->node)))
+                  *baseline->second.second->second->node->node)))
             innerpoint++;
           tmp = checker->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node);
+              *baseline->second.first->second->node->node);
           if ((tmp < baseline->second.first->second->node->node->DistanceSquared(
-              A->second->node->node)) && (tmp
+              *A->second->node->node)) && (tmp
               < baseline->second.first->second->node->node->DistanceSquared(
-                  baseline->second.second->second->node->node)))
+                  *baseline->second.second->second->node->node)))
             innerpoint++;
           tmp = checker->second->node->node->DistanceSquared(
-              baseline->second.second->second->node->node);
+              *baseline->second.second->second->node->node);
           if ((tmp < baseline->second.second->second->node->node->DistanceSquared(
-              baseline->second.first->second->node->node)) && (tmp
+              *baseline->second.first->second->node->node)) && (tmp
               < baseline->second.second->second->node->node->DistanceSquared(
-                  A->second->node->node)))
+                  *A->second->node->node)))
             innerpoint++;
           // 4e. If so, break 4. loop and continue with next candidate in 1. loop
@@ -1390 +1380 @@
         // prepare some auxiliary vectors
         Vector BaseLineCenter, BaseLine;
-        BaseLineCenter.CopyVector(baseline->second->endpoints[0]->node->node);
-        BaseLineCenter.AddVector(baseline->second->endpoints[1]->node->node);
-        BaseLineCenter.Scale(1. / 2.); // points now to center of base line
-        BaseLine.CopyVector(baseline->second->endpoints[0]->node->node);
-        BaseLine.SubtractVector(baseline->second->endpoints[1]->node->node);
+        BaseLineCenter = 0.5 * ((*baseline->second->endpoints[0]->node->node) +
+                                (*baseline->second->endpoints[1]->node->node));
+        BaseLine = (*baseline->second->endpoints[0]->node->node) - (*baseline->second->endpoints[1]->node->node);
 
         // offset to center of triangle
         CenterVector.Zero();
         for (int i = 0; i < 3; i++)
-          CenterVector.AddVector(BTS->endpoints[i]->node->node);
+          CenterVector += (*BTS->endpoints[i]->node->node);
         CenterVector.Scale(1. / 3.);
         Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl;
 
         // normal vector of triangle
-        NormalVector.CopyVector(Center);
-        NormalVector.SubtractVector(&CenterVector);
+        NormalVector = (*Center) - CenterVector;
         BTS->GetNormalVector(NormalVector);
-        NormalVector.CopyVector(&BTS->NormalVector);
+        NormalVector = BTS->NormalVector;
         Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl;
 
@@ -1413 +1400 @@
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
         PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
-        TempVector.CopyVector(&CenterVector);
-        TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
+        TempVector = CenterVector - (*baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
         //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
-        if (PropagationVector.ScalarProduct(&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.);
         Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl;
@@ -1427 +1413 @@
 
             // first check direction, so that triangles don't intersect
-            VirtualNormalVector.CopyVector(target->second->node->node);
-            VirtualNormalVector.SubtractVector(&BaseLineCenter); // points from center of base line to target
-            VirtualNormalVector.ProjectOntoPlane(&NormalVector);
-            TempAngle = VirtualNormalVector.Angle(&PropagationVector);
+            VirtualNormalVector = (*target->second->node->node) - BaseLineCenter;
+            TempAngle = VirtualNormalVector.Angle(PropagationVector);
             Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;
             if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees)
@@ -1457 +1441 @@
 
             // check for linear dependence
-            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
-            TempVector.SubtractVector(target->second->node->node);
-            helper.CopyVector(baseline->second->endpoints[1]->node->node);
-            helper.SubtractVector(target->second->node->node);
-            helper.ProjectOntoPlane(&TempVector);
+            TempVector = (*baseline->second->endpoints[0]->node->node) - (*target->second->node->node);
+            helper = (*baseline->second->endpoints[1]->node->node) - (*target->second->node->node);
+            helper.ProjectOntoPlane(TempVector);
             if (fabs(helper.NormSquared()) < MYEPSILON) {
               Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl;
@@ -1472 +1454 @@
                                         *(baseline->second->endpoints[1]->node->node),
                                         *(target->second->node->node)).getNormal();
-            TempVector.CopyVector(baseline->second->endpoints[0]->node->node);
-            TempVector.AddVector(baseline->second->endpoints[1]->node->node);
-            TempVector.AddVector(target->second->node->node);
-            TempVector.Scale(1./3.);
-            TempVector.SubtractVector(Center);
+            TempVector = (1./3.) * ((*baseline->second->endpoints[0]->node->node) +
+                                    (*baseline->second->endpoints[1]->node->node) +
+                                    (*target->second->node->node));
+            TempVector -= (*Center);
             // make it always point outward
-            if (VirtualNormalVector.ScalarProduct(&TempVector) < 0)
+            if (VirtualNormalVector.ScalarProduct(TempVector) < 0)
               VirtualNormalVector.Scale(-1.);
             // calculate angle
-            TempAngle = NormalVector.Angle(&VirtualNormalVector);
+            TempAngle = NormalVector.Angle(VirtualNormalVector);
             Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;
             if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
@@ -1489 +1470 @@
             } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
               // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
-              helper.CopyVector(target->second->node->node);
-              helper.SubtractVector(&BaseLineCenter);
-              helper.ProjectOntoPlane(&BaseLine);
+              helper = (*target->second->node->node) - BaseLineCenter;
+              helper.ProjectOntoPlane(BaseLine);
               // ...the one with the smaller angle is the better candidate
-              TempVector.CopyVector(target->second->node->node);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              TempAngle = TempVector.Angle(&helper);
-              TempVector.CopyVector(winner->second->node->node);
-              TempVector.SubtractVector(&BaseLineCenter);
-              TempVector.ProjectOntoPlane(&VirtualNormalVector);
-              if (TempAngle < TempVector.Angle(&helper)) {
-                TempAngle = NormalVector.Angle(&VirtualNormalVector);
+              TempVector = (*target->second->node->node) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              TempAngle = TempVector.Angle(helper);
+              TempVector = (*winner->second->node->node) - BaseLineCenter;
+              TempVector.ProjectOntoPlane(VirtualNormalVector);
+              if (TempAngle < TempVector.Angle(helper)) {
+                TempAngle = NormalVector.Angle(VirtualNormalVector);
                 SmallestAngle = TempAngle;
                 winner = target;
@@ -1538 +1516 @@
           BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
           BTS->GetCenter(&helper);
-          helper.SubtractVector(Center);
-          helper.Scale(-1);
+          helper -= (*Center);
+          helper *= -1;
           BTS->GetNormalVector(helper);
           TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
@@ -1596 +1574 @@
       Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl;
       // we have the intersection, check whether in- or outside of boundary
-      if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) {
+      if ((Center->DistanceSquared(*Walker->node) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
         // inside, next!
         Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;
@@ -1609 +1587 @@
           OldPoints[i] = BTS->endpoints[i];
         }
-        Normal.CopyVector(&BTS->NormalVector);
+        Normal = BTS->NormalVector;
         // add Walker to boundary points
         Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl;
@@ -2119 +2097 @@
 
     // construct center of circle
-    CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);
-    CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);
-    CircleCenter.Scale(0.5);
+    CircleCenter = 0.5 * ((*BaseLine.endpoints[0]->node->node) + (*BaseLine.endpoints[1]->node->node));
 
     // construct normal vector of circle
-    CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node);
-    CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);
+    CirclePlaneNormal = (*BaseLine.endpoints[0]->node->node) - (*BaseLine.endpoints[1]->node->node);
 
     double radius = CirclePlaneNormal.NormSquared();
     double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.);
 
-    NormalVector.ProjectOntoPlane(&CirclePlaneNormal);
+    NormalVector.ProjectOntoPlane(CirclePlaneNormal);
     NormalVector.Normalize();
     ShortestAngle = 2.*M_PI; // This will indicate the quadrant.
 
-    SphereCenter.CopyVector(&NormalVector);
-    SphereCenter.Scale(CircleRadius);
-    SphereCenter.AddVector(&CircleCenter);
+    SphereCenter = (CircleRadius * NormalVector) +  CircleCenter;
     // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
 
@@ -2337 +2310 @@
 
   // construct center of circle
-  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  CircleCenter.Scale(0.5);
+  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
+                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
 
   // calculate squared radius of circle
-  radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
+  radius = CirclePlaneNormal.ScalarProduct(CirclePlaneNormal);
   if (radius/4. < RADIUS*RADIUS) {
     // construct relative sphere center with now known CircleCenter
-    RelativeSphereCenter.CopyVector(&T.SphereCenter);
-    RelativeSphereCenter.SubtractVector(&CircleCenter);
+    RelativeSphereCenter = T.SphereCenter - CircleCenter;
 
     CircleRadius = RADIUS*RADIUS - radius/4.;
@@ -2360 +2331 @@
     // construct SearchDirection and an "outward pointer"
     SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
-    helper.CopyVector(&CircleCenter);
-    helper.SubtractVector(ThirdNode->node);
-    if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
+    helper = CircleCenter - (*ThirdNode->node);
+    if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
       SearchDirection.Scale(-1.);
     Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-    if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
+    if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
       // rotated the wrong way!
       eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl;
@@ -2512 +2482 @@
     AddTesselationTriangle();
     BTS->GetCenter(&Center);
-    Center.SubtractVector(&CandidateLine.OptCenter);
-    BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);
+    Center -= CandidateLine.OptCenter;
+    BTS->SphereCenter = CandidateLine.OptCenter;
     BTS->GetNormalVector(Center);
 
@@ -2559 +2529 @@
   Vector DistanceToIntersection[2], BaseLine;
   double distance[2];
-  BaseLine.CopyVector(Base->endpoints[1]->node->node);
-  BaseLine.SubtractVector(Base->endpoints[0]->node->node);
+  BaseLine = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
   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]);
+    DistanceToIntersection[i] = (*ClosestPoint) - (*Base->endpoints[i]->node->node);
+    distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
   }
   delete(ClosestPoint);
@@ -2636 +2604 @@
 
   // get the distance vector from Base line to OtherBase line
-  Vector Distance;
-  Distance.CopyVector(ClosestPoint[1]);
-  Distance.SubtractVector(ClosestPoint[0]);
+  Vector Distance = (*ClosestPoint[1]) - (*ClosestPoint[0]);
 
   // calculate volume
@@ -2660 +2626 @@
     for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
       Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
-      BaseLineNormal.AddVector(&(runner->second->NormalVector));
+      BaseLineNormal += (runner->second->NormalVector);
     }
     BaseLineNormal.Scale(1./2.);
 
-    if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
+    if (Distance.ScalarProduct(BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
       Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;
       // calculate volume summand as a general tetraeder
@@ -2699 +2665 @@
   for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
     Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;
-    BaseLineNormal.AddVector(&(runner->second->NormalVector));
+    BaseLineNormal += (runner->second->NormalVector);
   }
   BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
@@ -2840 +2806 @@
               double distance, scaleFactor;
 
-              OrthogonalizedOben.CopyVector(&Oben);
-              aCandidate.CopyVector(a->node);
-              aCandidate.SubtractVector(Candidate->node);
-              OrthogonalizedOben.ProjectOntoPlane(&aCandidate);
+              OrthogonalizedOben = Oben;
+              aCandidate = (*a->node) - (*Candidate->node);
+              OrthogonalizedOben.ProjectOntoPlane(aCandidate);
               OrthogonalizedOben.Normalize();
               distance = 0.5 * aCandidate.Norm();
@@ -2849 +2814 @@
               OrthogonalizedOben.Scale(scaleFactor);
 
-              Center.CopyVector(Candidate->node);
-              Center.AddVector(a->node);
-              Center.Scale(0.5);
-              Center.AddVector(&OrthogonalizedOben);
-
-              AngleCheck.CopyVector(&Center);
-              AngleCheck.SubtractVector(a->node);
+              Center = 0.5 * ((*Candidate->node) + (*a->node));
+              Center += OrthogonalizedOben;
+
+              AngleCheck = Center - (*a->node);
               norm = aCandidate.Norm();
               // second point shall have smallest angle with respect to Oben vector
               if (norm < RADIUS*2.) {
-                angle = AngleCheck.Angle(&Oben);
+                angle = AngleCheck.Angle(Oben);
                 if (angle < Storage[0]) {
                   //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
@@ -2935 +2897 @@
 
   // construct center of circle
-  CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-  CircleCenter.Scale(0.5);
+  CircleCenter = 0.5 * ((*CandidateLine.BaseLine->endpoints[0]->node->node) +
+                        (*CandidateLine.BaseLine->endpoints[1]->node->node));
 
   // construct normal vector of circle
-  CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node);
-  CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node);
-
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
+  CirclePlaneNormal = (*CandidateLine.BaseLine->endpoints[0]->node->node) -
+                      (*CandidateLine.BaseLine->endpoints[1]->node->node);
+
+  RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
 
   // calculate squared radius TesselPoint *ThirdNode,f circle
@@ -2954 +2914 @@
 
     // test whether old center is on the band's plane
-    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);
+    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();
@@ -2964 +2924 @@
       // check SearchDirection
       Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
-      if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
+      if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {  // rotated the wrong way!
         eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl;
       }
@@ -3011 +2971 @@
 
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
-                    radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
+                    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);
+                      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));
+                      NewSphereCenter = NewPlaneCenter;
+                      OtherNewSphereCenter = NewPlaneCenter;
+                      helper = sqrt(RADIUS*RADIUS - radius) * NewNormalVector;
                       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 += helper;
                       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 += helper;
                       Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;
 
@@ -3041 +3000 @@
                       if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
                         if (fabs(alpha - Otheralpha) > MYEPSILON) {
-                          CandidateLine.OptCenter.CopyVector(&NewSphereCenter);
-                          CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);
+                          CandidateLine.OptCenter = NewSphereCenter;
+                          CandidateLine.OtherOptCenter = OtherNewSphereCenter;
                         } else {
-                          CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);
-                          CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);
+                          CandidateLine.OptCenter = OtherNewSphereCenter;
+                          CandidateLine.OtherOptCenter = NewSphereCenter;
                         }
                         // if there is an equal candidate, add it to the list without clearing the list
@@ -3168 +3127 @@
             FindPoint = PointsOnBoundary.find((*Runner)->nr);
             if (FindPoint != PointsOnBoundary.end()) {
-              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) );
+              points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(*x), FindPoint->second) );
               Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl;
             }
@@ -3212 +3171 @@
     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);
+      helper = 0.5 * ((*(LineRunner->second)->endpoints[0]->node->node) +
+                      (*(LineRunner->second)->endpoints[1]->node->node));
+      Center = (*x) - helper;
+      BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
+                 (*(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);
+        helper = (*x) - (*(LineRunner->second)->endpoints[0]->node->node) - Center;
+        const double lengthA = helper.ScalarProduct(BaseLine);
+        helper = (*x) - (*(LineRunner->second)->endpoints[1]->node->node) - Center;
+        const double lengthB = helper.ScalarProduct(BaseLine);
         if (lengthB*lengthA < 0) {  // if have different sign
           ClosestLine = LineRunner->second;
@@ -3280 +3233 @@
     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);
+      BaseLine = (*(LineRunner->second)->endpoints[0]->node->node) -
+                 (*(LineRunner->second)->endpoints[1]->node->node);
       const double lengthBase = BaseLine.NormSquared();
 
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[0]->node->node);
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[0]->node->node);
       const double lengthEndA = BaseLineIntersection.NormSquared();
 
-      BaseLineIntersection.CopyVector(x);
-      BaseLineIntersection.SubtractVector((LineRunner->second)->endpoints[1]->node->node);
+      BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
       const double lengthEndB = BaseLineIntersection.NormSquared();
 
@@ -3307 +3258 @@
       } 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);
+        BaseLineIntersection = (*x) - (*(LineRunner->second)->endpoints[1]->node->node);
+        Center = BaseLineIntersection;
+        Center.ProjectOntoPlane(BaseLine);
+        BaseLineIntersection -= Center;
         const double distance = BaseLineIntersection.NormSquared();
         if (Center.NormSquared() > BaseLine.NormSquared()) {
@@ -3363 +3313 @@
   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);
+    helper = (*x) - Center;
+    const double Alignment = helper.Angle((*Runner)->NormalVector);
     if (Alignment < MinAlignment) {
       result = *Runner;
@@ -3427 +3376 @@
   triangle->GetCenter(&Center);
   Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl;
-  DistanceToCenter.CopyVector(&Center);
-  DistanceToCenter.SubtractVector(&Point);
+  DistanceToCenter = Center - 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(&triangle->NormalVector)) < MYEPSILON) {
+  if (fabs(DistanceToCenter.ScalarProduct(triangle->NormalVector)) < MYEPSILON) {
     // calculate whether inside of triangle
-    DistanceToCenter.CopyVector(&Point);
-    Center.CopyVector(&Point);
-    Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter
-    DistanceToCenter.AddVector(&triangle->NormalVector); // points outside
+    DistanceToCenter = Point + triangle->NormalVector; // points outside
+    Center = Point - triangle->NormalVector; // points towards MolCenter
     Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl;
     if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) {
@@ -3452 +3398 @@
 
     // then check direction to boundary
-    if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) {
+    if (DistanceToCenter.ScalarProduct(triangle->NormalVector) > MYEPSILON) {
       Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl;
       return -distance;
@@ -3568 +3514 @@
   if ((triangles != NULL) && (!triangles->empty())) {
     for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-      PlaneNormal.AddVector(&(*Runner)->NormalVector);
+      PlaneNormal += (*Runner)->NormalVector;
   } else {
     eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl;
@@ -3579 +3525 @@
   // construct one orthogonal vector
   if (Reference != NULL) {
-    AngleZero.CopyVector(Reference);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    AngleZero = (*Reference) - (*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);
+    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*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;
@@ -3602 +3546 @@
   // go through all connected points and calculate angle
   for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
-    helper.CopyVector((*listRunner)->node);
-    helper.SubtractVector(Point->node);
-    helper.ProjectOntoPlane(&PlaneNormal);
+    helper = (*(*listRunner)->node) - (*Point->node);
+    helper.ProjectOntoPlane(PlaneNormal);
     double angle = GetAngle(helper, AngleZero, OrthogonalVector);
     Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;
@@ -3671 +3614 @@
     TesselB++;
     TesselC++;
-    PlaneNormal.AddVector(&helper);
+    PlaneNormal += helper;
   }
   //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
@@ -3686 +3629 @@
   // construct one orthogonal vector
   if (Reference != NULL) {
-    AngleZero.CopyVector(Reference);
-    AngleZero.SubtractVector(Point->node);
-    AngleZero.ProjectOntoPlane(&PlaneNormal);
+    AngleZero = (*Reference) - (*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);
+    AngleZero = (*(*SetOfNeighbours->begin())->node) - (*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;
@@ -3710 +3651 @@
   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);
+    helper = (*(*listRunner)->node) - (*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)
@@ -3959 +3899 @@
 
   // copy old location for the volume
-  OldPoint.CopyVector(point->node->node);
+  OldPoint = (*point->node->node);
 
   // get list of connected points
@@ -3987 +3927 @@
   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
+    NormalVector -= Runner->second->NormalVector; // has to point inward
     RemoveTesselationTriangle(Runner->second);
     count++;
@@ -4023 +3963 @@
           StartNode--;
           //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
-          Point.CopyVector((*StartNode)->node);
-          Point.SubtractVector((*MiddleNode)->node);
+          Point = (*(*StartNode)->node) - (*(*MiddleNode)->node);
           StartNode = MiddleNode;
           StartNode++;
@@ -4030 +3969 @@
             StartNode = connectedPath->begin();
           //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
-          Reference.CopyVector((*StartNode)->node);
-          Reference.SubtractVector((*MiddleNode)->node);
-          OrthogonalVector.CopyVector((*MiddleNode)->node);
-          OrthogonalVector.SubtractVector(&OldPoint);
+          Reference = (*(*StartNode)->node) - (*(*MiddleNode)->node);
+          OrthogonalVector = (*(*MiddleNode)->node) - OldPoint;
           OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
@@ -4502 +4439 @@
   class BoundaryLineSet *BestLine = NULL;
   for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
-    BaseLine.CopyVector(Runner->second->endpoints[0]->node->node);
-    BaseLine.SubtractVector(Runner->second->endpoints[1]->node->node);
-    CenterToPoint.CopyVector(Runner->second->endpoints[0]->node->node);
-    CenterToPoint.AddVector(Runner->second->endpoints[1]->node->node);
-    CenterToPoint.Scale(0.5);
-    CenterToPoint.SubtractVector(point->node);
-    angle = CenterToPoint.Angle(&BaseLine);
+    BaseLine = (*Runner->second->endpoints[0]->node->node) -
+               (*Runner->second->endpoints[1]->node->node);
+    CenterToPoint = 0.5 * ((*Runner->second->endpoints[0]->node->node) +
+                           (*Runner->second->endpoints[1]->node->node));
+    CenterToPoint -= (*point->node);
+    angle = CenterToPoint.Angle(BaseLine);
     if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
       BestAngle = angle;
@@ -4672 +4608 @@
       for (map < int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)
         if (VectorWalker != VectorRunner) { // skip equals
-          const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second);  // ScalarProduct should result in -1. for degenerated triangles
+          const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second);  // ScalarProduct should result in -1. for degenerated triangles
           Log() << Verbose(1) << "Checking " << *VectorWalker->second<< " against " << *VectorRunner->second << ": " << SCP << endl;
           if (fabs(SCP + 1.) < ParallelEpsilon) {
@@ -4758 +4694 @@
     TriangleSet::iterator TriangleWalker = T->begin();  // is the inner iterator
     /// 4a. Get NormalVector for one side (this is "front")
-    NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
+    NormalVector = (*TriangleWalker)->NormalVector;
     Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl;
     TriangleWalker++;
@@ -4769 +4705 @@
       TriangleSprinter++;
       Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl;
-      if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list
+      if (triangle->NormalVector.ScalarProduct(NormalVector) < 0) { // if from other side, then delete and remove from list
         Log() << Verbose(1) << " Removing ... " << endl;
         TriangleNrs.push(triangle->Nr);
@@ -4791 +4727 @@
       AddTesselationTriangle(); // ... and add
       TriangleNrs.pop();
-      BTS->NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);
-      BTS->NormalVector.Scale(-1.);
+      BTS->NormalVector = -1 * (*TriangleWalker)->NormalVector;
       TriangleWalker++;
     }

src/tesselationhelpers.cpp

@@ -88 +88 @@
   center->at(2) =  0.5 * m14/ m11;
 
-  if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
-    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;
+  if (fabs(a.Distance(*center) - RADIUS) > MYEPSILON)
+    eLog() << Verbose(1) << "The given center is further way by " << fabs(a.Distance(*center) - RADIUS) << " from a than RADIUS." << endl;
 
   gsl_matrix_free(A);
@@ -121 +121 @@
   Vector OtherCenter;
   Center->Zero();
-  helper.CopyVector(&a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(&b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(&c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
+  helper = sin(2.*alpha) * a;
+  (*Center) += helper;
+  helper = sin(2.*beta) * b;
+  (*Center) += helper;
+  helper = sin(2.*gamma) * c;
+  (*Center) += helper;
   //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
-  NewUmkreismittelpunkt->CopyVector(Center);
+  (*NewUmkreismittelpunkt) = (*Center);
   Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";
   // Here we calculated center of circumscribing circle, using barycentric coordinates
   Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";
 
-  TempNormal.CopyVector(&a);
-  TempNormal.SubtractVector(&b);
-  helper.CopyVector(&a);
-  helper.SubtractVector(&c);
-  TempNormal.VectorProduct(&helper);
+  TempNormal = a - b;
+  helper = a - c;
+  TempNormal.VectorProduct(helper);
   if (fabs(HalfplaneIndicator) < MYEPSILON)
     {
-      if ((TempNormal.ScalarProduct(AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(AlternativeDirection) >0 && AlternativeIndicator <0))
+      if ((TempNormal.ScalarProduct(*AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(*AlternativeDirection) >0 && AlternativeIndicator <0))
         {
-          TempNormal.Scale(-1);
+          TempNormal *= -1;
         }
     }
   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);
+          TempNormal *= -1;
         }
     }
@@ -163 +158 @@
   Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";
 
-  Center->AddVector(&TempNormal);
+  (*Center) += TempNormal;
   Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n";
   GetSphere(&OtherCenter, a, b, c, RADIUS);
@@ -181 +176 @@
   Vector helper;
   double alpha, beta, gamma;
-  Vector SideA, SideB, SideC;
-  SideA.CopyVector(b);
-  SideA.SubtractVector(&c);
-  SideB.CopyVector(c);
-  SideB.SubtractVector(&a);
-  SideC.CopyVector(a);
-  SideC.SubtractVector(&b);
-  alpha = M_PI - SideB.Angle(&SideC);
-  beta = M_PI - SideC.Angle(&SideA);
-  gamma = M_PI - SideA.Angle(&SideB);
+  Vector SideA = b - c;
+  Vector SideB = c - a;
+  Vector SideC = a - b;
+  alpha = M_PI - SideB.Angle(SideC);
+  beta = M_PI - SideC.Angle(SideA);
+  gamma = M_PI - SideA.Angle(SideB);
   //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
   if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) {
@@ -197 +188 @@
 
   Center->Zero();
-  helper.CopyVector(a);
-  helper.Scale(sin(2.*alpha));
-  Center->AddVector(&helper);
-  helper.CopyVector(b);
-  helper.Scale(sin(2.*beta));
-  Center->AddVector(&helper);
-  helper.CopyVector(c);
-  helper.Scale(sin(2.*gamma));
-  Center->AddVector(&helper);
+  helper = sin(2.*alpha) * a;
+  (*Center) += helper;
+  helper = sin(2.*beta) * b;
+  (*Center) += helper;
+  helper = sin(2.*gamma) * c;
+  (*Center) += helper;
   Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
 };
@@ -227 +215 @@
   Vector helper;
   double radius, alpha;
-  Vector RelativeOldSphereCenter;
-  Vector RelativeNewSphereCenter;
-
-  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
-  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
-  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
-  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
-  helper.CopyVector(&RelativeNewSphereCenter);
+
+  Vector RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
+  Vector RelativeNewSphereCenter = NewSphereCenter - CircleCenter;
+  helper = RelativeNewSphereCenter;
   // test whether new center is on the parameter circle's plane
-  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
-    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal))  << "!" << endl;
-    helper.ProjectOntoPlane(&CirclePlaneNormal);
+  if (fabs(helper.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
+    eLog() << Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(CirclePlaneNormal))  << "!" << endl;
+    helper.ProjectOntoPlane(CirclePlaneNormal);
   }
   radius = helper.NormSquared();
@@ -244 +228 @@
   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(&RelativeOldSphereCenter);
+  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
+  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 " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
-  radius = helper.Distance(&RelativeOldSphereCenter);
-  helper.ProjectOntoPlane(&NormalVector);
+  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)) {
@@ -280 +264 @@
   struct Intersection *I = (struct Intersection *)params;
   Vector intersection;
-  Vector SideA,SideB,HeightA, HeightB;
   for (int i=0;i<NDIM;i++)
     intersection[i] = gsl_vector_get(x, i);
 
-  SideA.CopyVector(&(I->x1));
-  SideA.SubtractVector(&I->x2);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&I->x1);
-  HeightA.ProjectOntoPlane(&SideA);
-
-  SideB.CopyVector(&I->x3);
-  SideB.SubtractVector(&I->x4);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&I->x3);
-  HeightB.ProjectOntoPlane(&SideB);
-
-  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
+  Vector SideA = I->x1 -I->x2 ;
+  Vector HeightA = intersection - I->x1;
+  HeightA.ProjectOntoPlane(SideA);
+
+  Vector SideB = I->x3 - I->x4;
+  Vector HeightB = intersection - I->x3;
+  HeightB.ProjectOntoPlane(SideB);
+
+  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
   //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
 
@@ -321 +300 @@
 
   struct Intersection par;
-    par.x1.CopyVector(&point1);
-    par.x2.CopyVector(&point2);
-    par.x3.CopyVector(&point3);
-    par.x4.CopyVector(&point4);
+    par.x1 = point1;
+    par.x2 = point2;
+    par.x3 = point3;
+    par.x4 = point4;
 
     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
@@ -370 +349 @@
 
     // check whether intersection is in between or not
-  Vector intersection, SideA, SideB, HeightA, HeightB;
+  Vector intersection;
   double t1, t2;
   for (int i = 0; i < NDIM; i++) {
@@ -376 +355 @@
   }
 
-  SideA.CopyVector(&par.x2);
-  SideA.SubtractVector(&par.x1);
-  HeightA.CopyVector(&intersection);
-  HeightA.SubtractVector(&par.x1);
-
-  t1 = HeightA.ScalarProduct(&SideA)/SideA.ScalarProduct(&SideA);
-
-  SideB.CopyVector(&par.x4);
-  SideB.SubtractVector(&par.x3);
-  HeightB.CopyVector(&intersection);
-  HeightB.SubtractVector(&par.x3);
-
-  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
+  Vector SideA = par.x2 - par.x1;
+  Vector HeightA = intersection - par.x1;
+
+  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
+
+  Vector SideB = par.x4 - par.x3;
+  Vector HeightB = intersection - par.x3;
+
+  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
 
   Log() << Verbose(1) << "Intersection " << intersection << " is at "
@@ -428 +403 @@
   if (point.IsZero())
     return M_PI;
-  double phi = point.Angle(&reference);
-  if (OrthogonalVector.ScalarProduct(&point) > 0) {
+  double phi = point.Angle(reference);
+  if (OrthogonalVector.ScalarProduct(point) > 0) {
     phi = 2.*M_PI - phi;
   }
@@ -456 +431 @@
   TetraederVector[2].CopyVector(c);
   for (int j=0;j<3;j++)
-    TetraederVector[j].SubtractVector(&d);
-  Point.CopyVector(&TetraederVector[0]);
-  Point.VectorProduct(&TetraederVector[1]);
-  volume = 1./6. * fabs(Point.ScalarProduct(&TetraederVector[2]));
+    TetraederVector[j].SubtractVector(d);
+  Point.CopyVector(TetraederVector[0]);
+  Point.VectorProduct(TetraederVector[1]);
+  volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
   return volume;
 };
@@ -583 +558 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
-            double currentNorm = helper. Norm();
+            helper = (*Point) - (*(*Runner)->node);
+            double currentNorm = helper.Norm();
             if (currentNorm < distance) {
               // remember second point
@@ -638 +612 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
-            helper.CopyVector(Point);
-            helper.SubtractVector((*Runner)->node);
+            helper = (*Point) - (*(*Runner)->node);
             double currentNorm = helper.NormSquared();
             if (currentNorm < distance) {
@@ -678 +651 @@
         Info FunctionInfo(__func__);
   // construct the plane of the two baselines (i.e. take both their directional vectors)
-  Vector Normal;
-  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);
+  Vector Baseline = (*Base->endpoints[1]->node->node) - (*Base->endpoints[0]->node->node);
+  Vector OtherBaseline = (*OtherBase->endpoints[1]->node->node) - (*OtherBase->endpoints[0]->node->node);
+  Vector Normal = Baseline;
+  Normal.VectorProduct(OtherBaseline);
   Normal.Normalize();
   Log() << Verbose(1) << "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(&NewOffset);
-  NewDirection.AddVector(&OtherBaseline);
+  Vector NewOffset = (*OtherBase->endpoints[0]->node->node) - (*Base->endpoints[0]->node->node);
+  NewOffset.ProjectOntoPlane(Normal);
+  NewOffset += (*Base->endpoints[0]->node->node);
+  Vector NewDirection = NewOffset + OtherBaseline;
 
   // calculate the intersection between this projected baseline and Base
@@ -704 +669 @@
                                                    *(Base->endpoints[1]->node->node),
                                                      NewOffset, NewDirection);
-  Normal.CopyVector(Intersection);
-  Normal.SubtractVector(Base->endpoints[0]->node->node);
-  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
+  Normal = (*Intersection) - (*Base->endpoints[0]->node->node);
+  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl;
 
   return Intersection;
@@ -724 +688 @@
     return -1;
   }
-  distance = x->DistanceToPlane(&triangle->NormalVector, triangle->endpoints[0]->node->node);
+  distance = x->DistanceToPlane(triangle->NormalVector, *triangle->endpoints[0]->node->node);
   return distance;
 };
@@ -787 +751 @@
     Vector *center = cloud->GetCenter();
     // make the circumsphere's center absolute again
-    helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
-    helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
-    helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
-    helper.Scale(1./3.);
-    helper.SubtractVector(center);
+    Vector helper = (1./3.) * ((*Tess->LastTriangle->endpoints[0]->node->node) +
+                               (*Tess->LastTriangle->endpoints[1]->node->node) +
+                               (*Tess->LastTriangle->endpoints[2]->node->node));
+    helper -= (*center);
     // and add to file plus translucency object
     *rasterfile << "# current virtual sphere\n";

src/test/ActOnAlltest.hpp

@@ -27 +27 @@
 
   template <typename klasse, typename res, typename T> void ActOnAll( res (klasse::*f)(T), T t );
+  template <typename klasse, typename res, typename T> void ActOnAll( res (klasse::*f)(T&), T &t );
   template <typename klasse, typename res, typename T, typename U> void ActOnAll( res (klasse::*f)(T, U), T t, U u );
   template <typename klasse, typename res, typename T, typename U, typename V> void ActOnAll( res (klasse::*f)(T, U, V), T t, U u, V v);
@@ -76 +77 @@
 };
 
+template <typename klasse, typename res, typename T> void VectorList::ActOnAll( res (klasse::*f)(T&), T &t )
+{
+  for (ListOfVectors::iterator Runner = Vectors.begin(); Runner != Vectors.end(); Runner++)
+    ((*Runner)->*f)(t);
+};
+
 template <typename klasse, typename res, typename T, typename U> void VectorList::ActOnAll( res (klasse::*f)(T, U), T t, U u )
 {

src/unittests/ActOnAllUnitTest.cpp

@@ -56 +56 @@
 
   // adding, subtracting
-  VL.ActOnAll( &Vector::AddVector, &test );
+  VL.ActOnAll( &Vector::AddVector, test );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , false );
-  VL.ActOnAll( &Vector::SubtractVector, &test );
+  VL.ActOnAll( &Vector::SubtractVector, test );
   CPPUNIT_ASSERT_EQUAL( VL == Ref , true );
 };

src/unittests/vectorunittest.cpp

@@ -72 +72 @@
   CPPUNIT_ASSERT_EQUAL( Vector(2.,3.,4.), fixture );
   // summation and scaling
-  fixture.CopyVector(&zero);
-  fixture.AddVector(&unit);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
-  fixture.CopyVector(&zero);
-  fixture.SubtractVector(&unit);
+  fixture = zero + unit;
+  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
+  fixture = zero - unit;
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
   CPPUNIT_ASSERT_EQUAL( false, fixture.IsZero() );
-  fixture.CopyVector(&zero);
-  fixture.AddVector(&zero);
+  fixture = zero + zero;
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsZero() );
-  fixture.CopyVector(&notunit);
-  fixture.SubtractVector(&otherunit);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
-  fixture.CopyVector(&unit);
-  fixture.AddVector(&otherunit);
+  fixture = notunit - otherunit;
+  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
+  fixture = unit - otherunit;
   CPPUNIT_ASSERT_EQUAL( false, fixture.IsOne() );
-  fixture.CopyVector(&notunit);
-  fixture.SubtractVector(&unit);
-  fixture.SubtractVector(&otherunit);
+  fixture = notunit - unit - otherunit;
   CPPUNIT_ASSERT_EQUAL( false, fixture.IsZero() );
-  fixture.CopyVector(&unit);
-  fixture.Scale(0.98);
+  fixture = 0.98 * unit;
   CPPUNIT_ASSERT_EQUAL( false, fixture.IsOne() );
-  fixture.CopyVector(&unit);
-  fixture.Scale(1.);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
-  fixture.CopyVector(&unit);
+  fixture = 1. * unit;
+  CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
   factor = 0.98;
-  fixture.Scale(factor);
+  fixture = factor * unit;
   CPPUNIT_ASSERT_EQUAL( false, fixture.IsOne() );
-  fixture.CopyVector(&unit);
   factor = 1.;
-  fixture.Scale(factor);
+  fixture = factor * unit;
   CPPUNIT_ASSERT_EQUAL( true, fixture.IsOne() );
 };
@@ -134 +123 @@
 void VectorTest::EuclidianScalarProductTest()
 {
-  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(&zero) );
-  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(&otherunit) );
-  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(&notunit) );
-  CPPUNIT_ASSERT_EQUAL( 1., unit.ScalarProduct(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 0., otherunit.ScalarProduct(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 0., otherunit.ScalarProduct(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 1., otherunit.ScalarProduct(&notunit) );
-  CPPUNIT_ASSERT_EQUAL( 2., two.ScalarProduct(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 1., two.ScalarProduct(&otherunit) );
-  CPPUNIT_ASSERT_EQUAL( 1., two.ScalarProduct(&notunit) );
+  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(zero) );
+  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(unit) );
+  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(otherunit) );
+  CPPUNIT_ASSERT_EQUAL( 0., zero.ScalarProduct(notunit) );
+  CPPUNIT_ASSERT_EQUAL( 1., unit.ScalarProduct(unit) );
+  CPPUNIT_ASSERT_EQUAL( 0., otherunit.ScalarProduct(unit) );
+  CPPUNIT_ASSERT_EQUAL( 0., otherunit.ScalarProduct(unit) );
+  CPPUNIT_ASSERT_EQUAL( 1., otherunit.ScalarProduct(notunit) );
+  CPPUNIT_ASSERT_EQUAL( 2., two.ScalarProduct(unit) );
+  CPPUNIT_ASSERT_EQUAL( 1., two.ScalarProduct(otherunit) );
+  CPPUNIT_ASSERT_EQUAL( 1., two.ScalarProduct(notunit) );
 }
 
@@ -165 +154 @@
 void VectorTest::EuclidianDistancesTest()
 {
-  CPPUNIT_ASSERT_EQUAL( 1., zero.Distance(&unit) );
-  CPPUNIT_ASSERT_EQUAL( sqrt(2.), otherunit.Distance(&unit) );
-  CPPUNIT_ASSERT_EQUAL( sqrt(2.), zero.Distance(&notunit) );
-  CPPUNIT_ASSERT_EQUAL( 1., otherunit.Distance(&notunit) );
-  CPPUNIT_ASSERT_EQUAL( sqrt(5.), two.Distance(&notunit) );
+  CPPUNIT_ASSERT_EQUAL( 1., zero.Distance(unit) );
+  CPPUNIT_ASSERT_EQUAL( sqrt(2.), otherunit.Distance(unit) );
+  CPPUNIT_ASSERT_EQUAL( sqrt(2.), zero.Distance(notunit) );
+  CPPUNIT_ASSERT_EQUAL( 1., otherunit.Distance(notunit) );
+  CPPUNIT_ASSERT_EQUAL( sqrt(5.), two.Distance(notunit) );
 }
 
@@ -176 +165 @@
 void VectorTest::EuclidianAnglesTest()
 {
-  CPPUNIT_ASSERT_EQUAL( M_PI, zero.Angle(&unit) );
-  CPPUNIT_ASSERT_EQUAL( 0., unit.Angle(&unit) );
-  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/2. - otherunit.Angle(&unit)) < MYEPSILON );
-  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/2. - unit.Angle(&notunit)) < MYEPSILON );
-  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/4. - otherunit.Angle(&notunit)) < MYEPSILON );
+  CPPUNIT_ASSERT_EQUAL( M_PI, zero.Angle(unit) );
+  CPPUNIT_ASSERT_EQUAL( 0., unit.Angle(unit) );
+  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/2. - otherunit.Angle(unit)) < MYEPSILON );
+  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/2. - unit.Angle(notunit)) < MYEPSILON );
+  CPPUNIT_ASSERT_EQUAL( true, fabs(M_PI/4. - otherunit.Angle(notunit)) < MYEPSILON );
 };
 
@@ -187 +176 @@
 void VectorTest::ProjectionTest()
 {
-  CPPUNIT_ASSERT_EQUAL( zero, zero.Projection(&unit) );
-  CPPUNIT_ASSERT_EQUAL( zero, otherunit.Projection(&unit) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0.4,0.2,0.),  otherunit.Projection(&two) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0.,1.,0.),  two.Projection(&otherunit) );
+  CPPUNIT_ASSERT_EQUAL( zero, zero.Projection(unit) );
+  CPPUNIT_ASSERT_EQUAL( zero, otherunit.Projection(unit) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0.4,0.2,0.),  otherunit.Projection(two) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0.,1.,0.),  two.Projection(otherunit) );
 };
 

src/vector.cpp

@@ -70 +70 @@
  * \return \f$| x - y |^2\f$
  */
-double Vector::DistanceSquared(const Vector * const y) const
+double Vector::DistanceSquared(const Vector &y) const
 {
   double res = 0.;
   for (int i=NDIM;i--;)
-    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
+    res += (x[i]-y[i])*(x[i]-y[i]);
   return (res);
 };
@@ -82 +82 @@
  * \return \f$| x - y |\f$
  */
-double Vector::Distance(const Vector * const y) const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (sqrt(res));
+double Vector::Distance(const Vector &y) const
+{
+  return (sqrt(DistanceSquared(y)));
 };
 
@@ -95 +92 @@
  * \return \f$| x - y |\f$
  */
-double Vector::PeriodicDistance(const Vector * const y, const double * const cell_size) const
+double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
 {
   double res = Distance(y), tmp, matrix[NDIM*NDIM];
@@ -119 +116 @@
         TranslationVector.MatrixMultiplication(matrix);
         // add onto the original vector to compare with
-        Shiftedy.CopyVector(y);
-        Shiftedy.AddVector(&TranslationVector);
+        Shiftedy = y + TranslationVector;
         // get distance and compare with minimum so far
-        tmp = Distance(&Shiftedy);
+        tmp = Distance(Shiftedy);
         if (tmp < res) res = tmp;
       }
@@ -133 +129 @@
  * \return \f$| x - y |^2\f$
  */
-double Vector::PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const
+double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
 {
   double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
@@ -157 +153 @@
         TranslationVector.MatrixMultiplication(matrix);
         // add onto the original vector to compare with
-        Shiftedy.CopyVector(y);
-        Shiftedy.AddVector(&TranslationVector);
+        Shiftedy = y + TranslationVector;
         // get distance and compare with minimum so far
-        tmp = DistanceSquared(&Shiftedy);
+        tmp = DistanceSquared(Shiftedy);
         if (tmp < res) res = tmp;
       }
@@ -180 +175 @@
 //  Output(out);
 //  Log() << Verbose(0) << endl;
-  TestVector.CopyVector(this);
+  TestVector = (*this);
   TestVector.InverseMatrixMultiplication(matrix);
   for(int i=NDIM;i--;) { // correct periodically
@@ -190 +185 @@
   }
   TestVector.MatrixMultiplication(matrix);
-  CopyVector(&TestVector);
+  (*this) = TestVector;
 //  Log() << Verbose(2) << "New corrected vector is: ";
 //  Output(out);
@@ -201 +196 @@
  * \return \f$\langle x, y \rangle\f$
  */
-double Vector::ScalarProduct(const Vector * const y) const
+double Vector::ScalarProduct(const Vector &y) const
 {
   double res = 0.;
   for (int i=NDIM;i--;)
-    res += x[i]*y->x[i];
+    res += x[i]*y[i];
   return (res);
 };
@@ -216 +211 @@
  *  \return \f$ x \times y \f&
  */
-void Vector::VectorProduct(const Vector * const y)
+void Vector::VectorProduct(const Vector &y)
 {
   Vector tmp;
-  tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
-  tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
-  tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
-  this->CopyVector(&tmp);
+  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
+  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
+  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
+  (*this) = tmp;
 };
 
@@ -230 +225 @@
  * \return \f$\langle x, y \rangle\f$
  */
-void Vector::ProjectOntoPlane(const Vector * const y)
-{
-  Vector tmp;
-  tmp.CopyVector(y);
+void Vector::ProjectOntoPlane(const Vector &y)
+{
+  Vector tmp = y;
   tmp.Normalize();
-  tmp.Scale(ScalarProduct(&tmp));
-  this->SubtractVector(&tmp);
+  tmp *= ScalarProduct(tmp);
+  (*this) -= tmp;
 };
 
@@ -245 +239 @@
  * \return distance to plane
  */
-double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
-{
-  Vector temp;
-
+double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
+{
   // first create part that is orthonormal to PlaneNormal with withdraw
-  temp.CopyVector(this);
-  temp.SubtractVector(PlaneOffset);
-  temp.MakeNormalTo(*PlaneNormal);
-  temp.Scale(-1.);
+  Vector temp = (*this) - PlaneOffset;
+  temp.MakeNormalTo(PlaneNormal);
+  temp *= -1.;
   // then add connecting vector from plane to point
-  temp.AddVector(this);
-  temp.SubtractVector(PlaneOffset);
+  temp += (*this);
+  temp -= PlaneOffset;
   double sign = temp.ScalarProduct(PlaneNormal);
   if (fabs(sign) > MYEPSILON)
@@ -269 +260 @@
  * \param *y array to second vector
  */
-void Vector::ProjectIt(const Vector * const y)
-{
-  Vector helper(*y);
+void Vector::ProjectIt(const Vector &y)
+{
+  Vector helper = y;
   helper.Scale(-(ScalarProduct(y)));
-  AddVector(&helper);
+  AddVector(helper);
 };
 
@@ -280 +271 @@
  * \return Vector
  */
-Vector Vector::Projection(const Vector * const y) const
-{
-  Vector helper(*y);
-  helper.Scale((ScalarProduct(y)/y->NormSquared()));
+Vector Vector::Projection(const Vector &y) const
+{
+  Vector helper = y;
+  helper.Scale((ScalarProduct(y)/y.NormSquared()));
 
   return helper;
@@ -293 +284 @@
 double Vector::Norm() const
 {
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += this->x[i]*this->x[i];
-  return (sqrt(res));
+  return (sqrt(NormSquared()));
 };
 
@@ -304 +292 @@
 double Vector::NormSquared() const
 {
-  return (ScalarProduct(this));
+  return (ScalarProduct(*this));
 };
 
@@ -363 +351 @@
  * @return true - vector is normalized, false - vector is not
  */
-bool Vector::IsNormalTo(const Vector * const normal) const
+bool Vector::IsNormalTo(const Vector &normal) const
 {
   if (ScalarProduct(normal) < MYEPSILON)
@@ -374 +362 @@
  * @return true - vector is normalized, false - vector is not
  */
-bool Vector::IsEqualTo(const Vector * const a) const
+bool Vector::IsEqualTo(const Vector &a) const
 {
   bool status = true;
   for (int i=0;i<NDIM;i++) {
-    if (fabs(x[i] - a->x[i]) > MYEPSILON)
+    if (fabs(x[i] - a[i]) > MYEPSILON)
       status = false;
   }
@@ -388 +376 @@
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
-double Vector::Angle(const Vector * const y) const
-{
-  double norm1 = Norm(), norm2 = y->Norm();
+double Vector::Angle(const Vector &y) const
+{
+  double norm1 = Norm(), norm2 = y.Norm();
   double angle = -1;
   if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
@@ -446 +434 @@
 const Vector& Vector::operator+=(const Vector& b)
 {
-  this->AddVector(&b);
+  this->AddVector(b);
   return *this;
 };
@@ -457 +445 @@
 const Vector& Vector::operator-=(const Vector& b)
 {
-  this->SubtractVector(&b);
+  this->SubtractVector(b);
   return *this;
 };
@@ -480 +468 @@
 {
   Vector x = *this;
-  x.AddVector(&b);
+  x.AddVector(b);
   return x;
 };
@@ -492 +480 @@
 {
   Vector x = *this;
-  x.SubtractVector(&b);
+  x.SubtractVector(b);
   return x;
 };
@@ -556 +544 @@
  * \param trans[] translation vector.
  */
-void Vector::Translate(const Vector * const trans)
-{
-  for (int i=NDIM;i--;)
-    x[i] += trans->x[i];
+void Vector::Translate(const Vector &trans)
+{
+  for (int i=NDIM;i--;)
+    x[i] += trans[i];
 };
 
@@ -639 +627 @@
  * \param *factors three-component vector with the factor for each given vector
  */
-void Vector::LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors)
-{
-  for(int i=NDIM;i--;)
-    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
+void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
+{
+  (*this) = (factors[0]*x1) +
+            (factors[1]*x2) +
+            (factors[2]*x3);
 };
 
@@ -648 +637 @@
  * \param n[] normal vector of mirror plane.
  */
-void Vector::Mirror(const Vector * const n)
+void Vector::Mirror(const Vector &n)
 {
   double projection;
-  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
+  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
   for (int i=NDIM;i--;)
-    x[i] -= 2.*projection*n->x[i];
+    x[i] -= 2.*projection*n[i];
 };
 
@@ -667 +656 @@
 {
   bool result = false;
-  double factor = y1.ScalarProduct(this)/y1.NormSquared();
-  Vector x1;
-  x1.CopyVector(&y1);
-  x1.Scale(factor);
-  SubtractVector(&x1);
+  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
+  Vector x1 = factor * y1 ;
+  SubtractVector(x1);
   for (int i=NDIM;i--;)
     result = result || (fabs(x[i]) > MYEPSILON);
@@ -684 +671 @@
  * \return true - success, false - failure (null vector given)
  */
-bool Vector::GetOneNormalVector(const Vector * const GivenVector)
+bool Vector::GetOneNormalVector(const Vector &GivenVector)
 {
   int Components[NDIM]; // contains indices of non-zero components
@@ -695 +682 @@
   // find two components != 0
   for (j=0;j<NDIM;j++)
-    if (fabs(GivenVector->x[j]) > MYEPSILON)
+    if (fabs(GivenVector[j]) > MYEPSILON)
       Components[Last++] = j;
 
@@ -701 +688 @@
     case 3:  // threecomponent system
     case 2:  // two component system
-      norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
+      norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
       x[Components[2]] = 0.;
       // in skp both remaining parts shall become zero but with opposite sign and third is zero
-      x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
-      x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
+      x[Components[1]] = -1./GivenVector[Components[1]] / norm;
+      x[Components[0]] = 1./GivenVector[Components[0]] / norm;
       return true;
       break;
@@ -720 +707 @@
 };
 
-/** 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
- * \param *C third plane vector
- * \return scaling parameter for this vector
- */
-double Vector::CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const
-{
-//  Log() << Verbose(3) << "For comparison: ";
-//  Log() << Verbose(0) << "A " << A->Projection(this) << "\t";
-//  Log() << Verbose(0) << "B " << B->Projection(this) << "\t";
-//  Log() << Verbose(0) << "C " << C->Projection(this) << "\t";
-//  Log() << Verbose(0) << endl;
-  return A->ScalarProduct(this);
-};
-
-
 /** Adds vector \a *y componentwise.
  * \param *y vector
  */
-void Vector::AddVector(const Vector * const y)
-{
-  for (int i=NDIM;i--;)
-    this->x[i] += y->x[i];
+void Vector::AddVector(const Vector &y)
+{
+  for (int i=NDIM;i--;)
+    this->x[i] += y[i];
 }
 
@@ -749 +719 @@
  * \param *y vector
  */
-void Vector::SubtractVector(const Vector * const y)
-{
-  for (int i=NDIM;i--;)
-    this->x[i] -= y->x[i];
-}
-
-/** Copy vector \a *y componentwise.
- * \param *y vector
- */
-void Vector::CopyVector(const Vector * const y)
-{
-  // check for self assignment
-  if(y!=this){
-    for (int i=NDIM;i--;)
-      this->x[i] = y->x[i];
-  }
+void Vector::SubtractVector(const Vector &y)
+{
+  for (int i=NDIM;i--;)
+    this->x[i] -= y[i];
 }
 
@@ -964 +922 @@
 bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
 {
-  Vector a;
-  a.CopyVector(this);
-  a.SubtractVector(&offset);
+  Vector a = (*this) - offset;
   a.InverseMatrixMultiplication(parallelepiped);
   bool isInside = true;

src/vector.hpp

@@ -33 +33 @@
   virtual ~Vector();
 
-  virtual double Distance(const Vector * const y) const;
-  virtual double DistanceSquared(const Vector * const y) const;
-  virtual double DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const;
-  virtual double PeriodicDistance(const Vector * const y, const double * const cell_size) const;
-  virtual double PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const;
-  virtual double ScalarProduct(const Vector * const y) const;
+  virtual double Distance(const Vector &y) const;
+  virtual double DistanceSquared(const Vector &y) const;
+  virtual double DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const;
+  virtual double PeriodicDistance(const Vector &y, const double * const cell_size) const;
+  virtual double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
+  virtual double ScalarProduct(const Vector &y) const;
   virtual double Norm() const;
   virtual double NormSquared() const;
-  virtual double Angle(const Vector * const y) const;
+  virtual double Angle(const Vector &y) const;
   virtual bool IsZero() const;
   virtual bool IsOne() const;
-  virtual bool IsNormalTo(const Vector * const normal) const;
-  virtual bool IsEqualTo(const Vector * const a) const;
+  virtual bool IsNormalTo(const Vector &normal) const;
+  virtual bool IsEqualTo(const Vector &a) const;
 
-  virtual void AddVector(const Vector * const y);
-  virtual void SubtractVector(const Vector * const y);
-  virtual void CopyVector(const Vector * const y);
+  virtual void AddVector(const Vector &y);
+  virtual void SubtractVector(const Vector &y);
   virtual void CopyVector(const Vector &y);
-  virtual void VectorProduct(const Vector * const y);
-  virtual void ProjectOntoPlane(const Vector * const y);
-  virtual void ProjectIt(const Vector * const y);
-  virtual Vector Projection(const Vector * const y) const;
+  virtual void VectorProduct(const Vector &y);
+  virtual void ProjectOntoPlane(const Vector &y);
+  virtual void ProjectIt(const Vector &y);
+  virtual Vector Projection(const Vector &y) const;
   virtual void Zero();
   virtual void One(const double one);
   virtual void Init(const double x1, const double x2, const double x3);
   virtual void Normalize();
-  virtual void Translate(const Vector * const x);
-  virtual void Mirror(const Vector * const x);
+  virtual void Translate(const Vector &x);
+  virtual void Mirror(const Vector &x);
   virtual void Scale(const double ** const factor);
   virtual void Scale(const double * const factor);
@@ -67 +66 @@
   virtual bool InverseMatrixMultiplication(const double * const M);
   virtual void KeepPeriodic(const double * const matrix);
-  virtual void LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors);
-  virtual double CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const;
-  virtual bool GetOneNormalVector(const Vector * const x1);
+  virtual void LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors);
+  virtual bool GetOneNormalVector(const Vector &x1);
   virtual bool MakeNormalTo(const Vector &y1);
   //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);

src/vector_ops.cpp

@@ -121 +121 @@
   // normalise this vector with respect to axis
   a.CopyVector(vec);
-  a.ProjectOntoPlane(&axis);
+  a.ProjectOntoPlane(axis);
   // construct normal vector
   try {
@@ -135 +135 @@
   y.Scale(sin(alpha));
   a.Scale(cos(alpha));
-  res = vec.Projection(&axis);
+  res = vec.Projection(axis);
   // add scaled normal vector onto this vector
-  res.AddVector(&y);
+  res += y;
   // add part in axis direction
-  res.AddVector(&a);
+  res += a;
   return res;
 };
@@ -197 +197 @@
   Vector parallel;
   double factor = 0.;
-  if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
+  if (fabs(a.ScalarProduct(b)*a.ScalarProduct(b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
     parallel = Line1a - Line2a;
-    factor = parallel.ScalarProduct(&a)/a.Norm();
+    factor = parallel.ScalarProduct(a)/a.Norm();
     if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
       res = Line2a;
@@ -206 +206 @@
     } else {
       parallel = Line1a - Line2b;
-      factor = parallel.ScalarProduct(&a)/a.Norm();
+      factor = parallel.ScalarProduct(a)/a.Norm();
       if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
         res = Line2b;
@@ -221 +221 @@
   double s;
   Vector temp1, temp2;
-  temp1.CopyVector(&c);
-  temp1.VectorProduct(&b);
-  temp2.CopyVector(&a);
-  temp2.VectorProduct(&b);
+  temp1 = c;
+  temp1.VectorProduct(b);
+  temp2 = a;
+  temp2.VectorProduct(b);
   Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
   if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
+    s = temp1.ScalarProduct(temp2)/temp2.NormSquared();
   else
     s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
+  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
 
   // construct intersection