Changeset 0a4f7f for src

Timestamp:

Apr 7, 2010, 3:45:38 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:

72e7fa

Parents:

f9352d

Message:

Made data internal data-structure of vector class private

• Replaced occurences of access to internals with operator
• moved Vector-class into LinAlg-Module
• Reworked Vector to allow clean modularization
• Added Plane class to describe arbitrary planes in 3d space

Location:

src

Files:

• Exceptions/CustomException.cpp (added)
• Exceptions/CustomException.hpp (added)
• Exceptions/LinearDependenceException.cpp (added)
• Exceptions/LinearDependenceException.hpp (added)
• Helpers/fast_functions.hpp (added)
• Legacy/oldmenu.cpp (modified) (20 diffs)
• Makefile.am (modified) (5 diffs)
• Plane.cpp (added)
• Plane.hpp (added)
• UIElements/TextDialog.cpp (modified) (1 diff)
• analysis_correlation.cpp (modified) (1 diff)
• atom.cpp (modified) (6 diffs)
• atom_trajectoryparticle.cpp (modified) (12 diffs)
• boundary.cpp (modified) (12 diffs)
• builder.cpp (modified) (6 diffs)
• config.cpp (modified) (9 diffs)
• ellipsoid.cpp (modified) (4 diffs)
• gslmatrix.cpp (modified) (1 diff)
• helpers.cpp (modified) (3 diffs)
• helpers.hpp (modified) (3 diffs)
• leastsquaremin.cpp (modified) (1 diff)
• linearsystemofequations.cpp (modified) (2 diffs)
• linkedcell.cpp (modified) (9 diffs)
• molecule.cpp (modified) (11 diffs)
• molecule_dynamics.cpp (modified) (5 diffs)
• molecule_fragmentation.cpp (modified) (3 diffs)
• molecule_geometry.cpp (modified) (9 diffs)
• molecule_graph.cpp (modified) (1 diff)
• moleculelist.cpp (modified) (2 diffs)
• tesselation.cpp (modified) (22 diffs)
• tesselationhelpers.cpp (modified) (13 diffs)
• unittests/ActOnAllUnitTest.cpp (modified) (2 diffs)
• unittests/ActOnAllUnitTest.hpp (modified) (1 diff)
• unittests/vectorunittest.cpp (modified) (2 diffs)
• vector.cpp (modified) (17 diffs)
• vector.hpp (modified) (3 diffs)
• vector_ops.cpp (added)
• vector_ops.hpp (added)

src/Legacy/oldmenu.cpp

@@ -24 +24 @@
 #include "molecule.hpp"
 #include "periodentafel.hpp"
+#include "vector_ops.hpp"
+#include "Plane.hpp"
 
 #include "UIElements/UIFactory.hpp"
@@ -77 +79 @@
       break;
       case 'a': // absolute coordinates of atom
-        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
+      {
         first = World::getInstance().createAtom();
-        first->x.AskPosition(mol->cell_size, false);
-        first->type = periode->AskElement();  // give type
-        mol->AddAtom(first);  // add to molecule
-        break;
+        Dialog *dialog = UIFactory::getInstance().makeDialog();
+        dialog->queryVector("Enter absolute coordinates.",&first->x,mol->cell_size, false);
+        dialog->queryElement("Choose element for this atom",&first->type);
+        if(dialog->display()){
+          mol->AddAtom(first);  // add to molecule
+        }
+        else{
+          // dialog was canceled... destroy the atom that was used
+          World::getInstance().destroyAtom(first);
+        }
+        delete dialog;
+      }
+      break;
 
       case 'b': // relative coordinates of atom wrt to reference point
@@ -88 +99 @@
         valid = true;
         do {
+          Dialog *dialog = UIFactory::getInstance().makeDialog();
           if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
-          Log() << Verbose(0) << "Enter reference coordinates." << endl;
-          x.AskPosition(mol->cell_size, true);
-          Log() << Verbose(0) << "Enter relative coordinates." << endl;
-          first->x.AskPosition(mol->cell_size, false);
-          first->x.AddVector((const Vector *)&x);
-          Log() << Verbose(0) << "\n";
-        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
+          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);
+          dialog->display();
+          delete dialog;
+        } while (!(valid = mol->CheckBounds(&first->x)));
         first->type = periode->AskElement();  // give type
         mol->AddAtom(first);  // add to molecule
@@ -101 +112 @@
 
       case 'c': // relative coordinates of atom wrt to already placed atom
+      {
         first = World::getInstance().createAtom();
         valid = true;
@@ -106 +118 @@
           if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
           second = mol->AskAtom("Enter atom number: ");
-          Log() << Verbose(0) << "Enter relative coordinates." << endl;
-          first->x.AskPosition(mol->cell_size, false);
+          Dialog *dialog = UIFactory::getInstance().makeDialog();
+          dialog->queryVector("Enter relative coordinates.",&first->x,mol->cell_size,false);
+          dialog->display();
           for (int i=NDIM;i--;) {
-            first->x.x[i] += second->x.x[i];
+            first->x[i] += second->x[i];
           }
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
         first->type = periode->AskElement();  // give type
         mol->AddAtom(first);  // add to molecule
-        break;
+      }
+      break;
 
     case 'd': // two atoms, two angles and a distance
@@ -157 +171 @@
           x.SubtractVector(&third->x);
           x.Normalize();
-          Log() << Verbose(0) << "x: ",
-          x.Output();
-          Log() << Verbose(0) << endl;
-          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
-          Log() << Verbose(0) << "z: ",
-          z.Output();
-          Log() << Verbose(0) << endl;
-          y.MakeNormalVector(&x,&z);
-          Log() << Verbose(0) << "y: ",
-          y.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "x: " << x << endl;
+          z = Plane(second->x,third->x,fourth->x).getNormal();
+          Log() << Verbose(0) << "z: " << z << endl;
+          y = Plane(x,z,0).getNormal();
+          Log() << Verbose(0) << "y: " << y << endl;
 
           // rotate vector around first angle
           first->x.CopyVector(&x);
-          first->x.RotateVector(&z,b - M_PI);
-          Log() << Verbose(0) << "Rotated vector: ",
-          first->x.Output();
-          Log() << Verbose(0) << endl;
+          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));
-          Log() << Verbose(0) << "N1: ",
-          n.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "N1: " << n << endl;
           first->x.SubtractVector(&n);
-          Log() << Verbose(0) << "Subtracted vector: ",
-          first->x.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "Subtracted vector: " << first->x << endl;
           n.CopyVector(&z);
           n.Scale(first->x.ScalarProduct(&z));
-          Log() << Verbose(0) << "N2: ",
-          n.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "N2: " << n << endl;
           first->x.SubtractVector(&n);
-          Log() << Verbose(0) << "2nd subtracted vector: ",
-          first->x.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "2nd subtracted vector: " << first->x << endl;
 
           // rotate another vector around second angle
           n.CopyVector(&y);
-          n.RotateVector(&x,c - M_PI);
-          Log() << Verbose(0) << "2nd Rotated vector: ",
-          n.Output();
-          Log() << Verbose(0) << endl;
+          n = RotateVector(n,x,c - M_PI);
+          Log() << Verbose(0) << "2nd Rotated vector: " << n << endl;
 
           // add the two linear independent vectors
@@ -208 +204 @@
           first->x.AddVector(&second->x);
 
-          Log() << Verbose(0) << "resulting coordinates: ";
-          first->x.Output();
-          Log() << Verbose(0) << endl;
+          Log() << Verbose(0) << "resulting coordinates: " << first->x << endl;
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
         first->type = periode->AskElement();  // give type
@@ -233 +227 @@
         } while ((j != -1) && (i<128));
         if (i >= 2) {
-          first->x.LSQdistance((const Vector **)atoms, i);
-
-          first->x.Output();
+          LSQdistance(first->x,(const Vector **)atoms, i);
+
+          Log() << Verbose(0) << first->x;
           first->type = periode->AskElement();  // give type
           mol->AddAtom(first);  // add to molecule
@@ -280 +274 @@
       for (int i=0;i<NDIM;i++) {
         Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
-        cin >> y.x[i];
+        cin >> y[i];
       }
       mol->CenterEdge(&x);  // make every coordinate positive
@@ -293 +287 @@
       for (int i=0;i<NDIM;i++) {
         Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
-        cin >> x.x[i];
+        cin >> x[i];
       }
       // update Box of atoms by boundary
@@ -330 +324 @@
       third = mol->AskAtom("Enter third atom: ");
 
-      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
+      n = Plane(first->x,second->x,third->x).getNormal();
       break;
     case 'b': // normal vector of mirror plane
-      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-      n.AskPosition(mol->cell_size,0);
+    {
+      Dialog *dialog = UIFactory::getInstance().makeDialog();
+      dialog->queryVector("Enter normal vector of mirror plane.",&n,mol->cell_size,false);
+      dialog->display();
+      delete dialog;
       n.Normalize();
-      break;
+    }
+    break;
+
     case 'c': // three atoms defining mirror plane
       first = mol->AskAtom("Enter first atom: ");
@@ -356 +355 @@
       mol->GetAlignvector(&param);
       for (int i=NDIM;i--;) {
-        x.x[i] = gsl_vector_get(param.x,i);
-        n.x[i] = gsl_vector_get(param.x,i+NDIM);
+        x[i] = gsl_vector_get(param.x,i);
+        n[i] = gsl_vector_get(param.x,i+NDIM);
       }
       gsl_vector_free(param.x);
-      Log() << Verbose(0) << "Offset vector: ";
-      x.Output();
-      Log() << Verbose(0) << endl;
+      Log() << Verbose(0) << "Offset vector: " << x << endl;
       n.Normalize();
       break;
   };
-  Log() << Verbose(0) << "Alignment vector: ";
-  n.Output();
-  Log() << Verbose(0) << endl;
+  Log() << Verbose(0) << "Alignment vector: " << n << endl;
   mol->Align(&n);
 };
@@ -397 +392 @@
       third = mol->AskAtom("Enter third atom: ");
 
-      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
+      n = Plane(first->x,second->x,third->x).getNormal();
       break;
     case 'b': // normal vector of mirror plane
-      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-      n.AskPosition(mol->cell_size,0);
+    {
+      Dialog *dialog = UIFactory::getInstance().makeDialog();
+      dialog->queryVector("Enter normal vector of mirror plane.",&n,mol->cell_size,false);
+      dialog->display();
+      delete dialog;
       n.Normalize();
-      break;
+    }
+    break;
+
     case 'c': // three atoms defining mirror plane
       first = mol->AskAtom("Enter first atom: ");
@@ -413 +413 @@
       break;
   };
-  Log() << Verbose(0) << "Normal vector: ";
-  n.Output();
-  Log() << Verbose(0) << endl;
+  Log() << Verbose(0) << "Normal vector: " << n << endl;
   mol->Mirror((const Vector *)&n);
 };
@@ -471 +469 @@
         first = second;
         second = first->next;
-        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
+        if ((first->x[axis] < tmp1) || (first->x[axis] > tmp2)) {// out of boundary ...
           //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
           mol->RemoveAtom(first);
@@ -541 +539 @@
       x.SubtractVector((const Vector *)&second->x);
       tmp1 = x.Norm();
-      Log() << Verbose(1) << "Distance vector is ";
-      x.Output();
-      Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
+      Log() << Verbose(1) << "Distance vector is " << x << "." << "/n"
+            << "Norm of distance is " << tmp1 << "." << endl;
       break;
 
@@ -676 +673 @@
         minBond = 0.;
         for (int i=NDIM;i--;)
-          minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
+          minBond += (first->x[i]-second->x[i])*(first->x[i] - second->x[i]);
         minBond = sqrt(minBond);
         Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
@@ -682 +679 @@
         cin >> bond;
         for (int i=NDIM;i--;) {
-          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
+          second->x[i] -= (second->x[i]-first->x[i])/minBond*(minBond-bond);
         }
         //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
@@ -778 +775 @@
       x.Zero();
       y.Zero();
-      y.x[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
+      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
@@ -893 +890 @@
         mol = *ListRunner;
         Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
-        Log() << Verbose(0) << "Enter translation vector." << endl;
-        x.AskPosition(mol->cell_size,0);
+        Dialog *dialog = UIFactory::getInstance().makeDialog();
+        dialog->queryVector("Enter translation vector.",&x,mol->cell_size,false);
+        dialog->display();
+        delete dialog;
         mol->Center.AddVector((const Vector *)&x);
      }

src/Makefile.am

@@ +1 @@
+# this includes source files that need to be present at multiple points
+HELPERSOURCE =  Helpers/Assert.cpp
+                        
 ATOMSOURCE = atom.cpp atom_atominfo.cpp atom_bondedparticle.cpp atom_bondedparticleinfo.cpp atom_graphnode.cpp atom_graphnodeinfo.cpp atom_particleinfo.cpp atom_trajectoryparticle.cpp atom_trajectoryparticleinfo.cpp
 ATOMHEADER = atom.hpp atom_atominfo.hpp atom_bondedparticle.hpp atom_bondedparticleinfo.hpp atom_graphnode.hpp atom_graphnodeinfo.hpp atom_particleinfo.hpp atom_trajectoryparticle.hpp atom_trajectoryparticleinfo.hpp
 
-LINALGSOURCE = gslmatrix.cpp gslvector.cpp linearsystemofequations.cpp
-LINALGHEADER = gslmatrix.hpp gslvector.hpp linearsystemofequations.hpp
+LINALGSOURCE = ${HELPERSOURCE} \
+               gslmatrix.cpp \
+                           gslvector.cpp \
+                           linearsystemofequations.cpp \
+                           vector.cpp
+                           
+LINALGHEADER = gslmatrix.hpp \
+                           gslvector.hpp \
+                           linearsystemofequations.hpp \
+                           vector.hpp
+                           
 
 ANALYSISSOURCE = analysis_bonds.cpp analysis_correlation.cpp
@@ -70 +82 @@
 
 
+EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
+                                  Exceptions/LinearDependenceException.cpp
+                                  
+EXCEPTIONHEADER = Exceptions/CustomException.hpp \
+                                  Exceptions/LinearDependenceException.hpp
+
 SOURCE = ${ANALYSISSOURCE} \
                  ${ATOMSOURCE} \
@@ -75 +93 @@
                  ${UISOURCE} \
                  ${DESCRIPTORSOURCE} \
+                 ${HELPERSOURCE} \
                  ${LEGACYSOURCE} \
+                 ${EXCEPTIONSOURCE} \
                  bond.cpp \
                  bondgraph.cpp \
@@ -85 +105 @@
                  graph.cpp \
                  helpers.cpp \
-                 Helpers/Assert.cpp \
                  info.cpp \
                  leastsquaremin.cpp \
@@ -102 +121 @@
                  parser.cpp \
                  periodentafel.cpp \
+                 Plane.cpp \
                  tesselation.cpp \
                  tesselationhelpers.cpp \
-                 vector.cpp \
                  verbose.cpp \
+                 vector_ops.cpp \
                  World.cpp
-HEADER = ${ANALYSISHEADER} ${ATOMHEADER} ${PATTERNHEADER} ${UIHEADER} ${DESCRIPTORHEADER} ${LEGACYHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp World.hpp 
+                 
+HEADER = ${ANALYSISHEADER} \
+         ${ATOMHEADER} \
+         ${PATTERNHEADER} \
+         ${UIHEADER} \
+         ${DESCRIPTORHEADER} \
+         ${LEGACYHEADER} \
+         ${EXCEPTIONHEADER} \
+         bond.hpp \
+         bondgraph.hpp \
+         boundary.hpp \
+         config.hpp \
+         defs.hpp \
+         element.hpp \
+         ellipsoid.hpp \
+         errorlogger.hpp \
+         graph.hpp \
+         info.hpp \
+         leastsquaremin.hpp \
+         linkedcell.hpp \
+         lists.hpp \
+         log.hpp \
+         logger.hpp \
+         memoryallocator.hpp \
+         memoryusageobserver.hpp \
+         molecule.hpp \
+         molecule_template.hpp \
+         parser.hpp \
+         periodentafel.hpp \
+         Plane.hpp \
+         stackclass.hpp \
+         tesselation.hpp \
+         tesselationhelpers.hpp \
+         verbose.hpp \
+         vector_ops.hpp \
+         World.hpp 
 
 BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)

src/UIElements/TextDialog.cpp

@@ -119 +119 @@
 
 bool TextDialog::VectorTextQuery::handle() {
- Log() << Verbose(0) << getTitle();
- tmp->AskPosition(cellSize,check);
- return true;
+  Log() << Verbose(0) << getTitle();
+
+  char coords[3] = {'x','y','z'};
+  int j = -1;
+  for (int i=0;i<3;i++) {
+    j += i+1;
+    do {
+      Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
+      cin >> (*tmp)[i];
+    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
+  }
+  return true;
 }
 

src/analysis_correlation.cpp

@@ -400 +400 @@
     *file << runner->first;
     for (int i=0;i<NDIM;i++)
-      *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
+      *file << "\t" << (runner->second.first->node->at(i) - runner->second.second->at(i));
     *file << endl;
   }

src/atom.cpp

@@ -130 +130 @@
   if (out != NULL) {
     *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint;
-    *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
+    *out << x[0] << "\t" << x[1] << "\t" << x[2];
     *out << "\t" << FixedIon;
     if (v.Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
+      *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[2] << "\t";
     if (comment != NULL)
       *out << " # " << comment << endl;
@@ -155 +155 @@
   if (out != NULL) {
     *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
-    *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
+    *out << x[0] << "\t" << x[1] << "\t" << x[2];
     *out << "\t" << FixedIon;
     if (v.Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
+      *out << "\t" << scientific << setprecision(6) << v[0] << "\t" << v[1] << "\t" << v[2] << "\t";
     if (comment != NULL)
       *out << " # " << comment << endl;
@@ -175 +175 @@
 {
   if (out != NULL) {
-    *out << type->symbol << "\t" << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2] << "\t" << endl;
+    *out << type->symbol << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << "\t" << endl;
     return true;
   } else
@@ -193 +193 @@
   if (out != NULL) {
     *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint;
-    *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
+    *out << Trajectory.R.at(step)[0] << "\t" << Trajectory.R.at(step)[1] << "\t" << Trajectory.R.at(step)[2];
     *out << "\t" << FixedIon;
     if (Trajectory.U.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t";
+      *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step)[0] << "\t" << Trajectory.U.at(step)[1] << "\t" << Trajectory.U.at(step)[2] << "\t";
     if (Trajectory.F.at(step).Norm() > MYEPSILON)
-      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t";
+      *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step)[0] << "\t" << Trajectory.F.at(step)[1] << "\t" << Trajectory.F.at(step)[2] << "\t";
     *out << "\t# Number in molecule " << nr << endl;
     return true;
@@ -214 +214 @@
   if (out != NULL) {
     *out << type->symbol << "\t";
-    *out << Trajectory.R.at(step).x[0] << "\t";
-    *out << Trajectory.R.at(step).x[1] << "\t";
-    *out << Trajectory.R.at(step).x[2] << endl;
+    *out << Trajectory.R.at(step)[0] << "\t";
+    *out << Trajectory.R.at(step)[1] << "\t";
+    *out << Trajectory.R.at(step)[2] << endl;
     return true;
   } else
@@ -229 +229 @@
 void atom::OutputMPQCLine(ofstream * const out, const Vector *center, int *AtomNo = NULL) const
 {
-  *out << "\t\t" << type->symbol << " [ " << x.x[0]-center->x[0] << "\t" << x.x[1]-center->x[1] << "\t" << x.x[2]-center->x[2] << " ]" << endl;
+  *out << "\t\t" << type->symbol << " [ " << x[0]-center->at(0) << "\t" << x[1]-center->at(1) << "\t" << x[2]-center->at(2) << " ]" << endl;
   if (AtomNo != NULL)
     *AtomNo++;

src/atom_trajectoryparticle.cpp

@@ -34 +34 @@
 {
   for (int i=NDIM;i--;)
-    *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
+    *temperature += type->mass * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i];
 };
 
@@ -60 +60 @@
 {
   for(int d=0;d<NDIM;d++) {
-    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
-    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
+    Trajectory.U.at(Step)[d] -= CoGVelocity->at(d);
+    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step)[d] * Trajectory.U.at(Step)[d];
   }
 };
@@ -89 +89 @@
 
   for (int n=NDIM;n--;) {
-    Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
-    Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
-    Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[n];
+    Trajectory.R.at(dest)[n] = Trajectory.R.at(src)[n];
+    Trajectory.U.at(dest)[n] = Trajectory.U.at(src)[n];
+    Trajectory.F.at(dest)[n] = Trajectory.F.at(src)[n];
   }
 };
@@ -105 +105 @@
   //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
   for (int d=0; d<NDIM; d++) {
-    Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-    Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
-    Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-    Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass);     // F = m * a and s =
+    Trajectory.F.at(NextStep)[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    Trajectory.R.at(NextStep)[d] = Trajectory.R.at(NextStep-1)[d];
+    Trajectory.R.at(NextStep)[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/type->mass);     // F = m * a and s =
   }
   // Update U
   for (int d=0; d<NDIM; d++) {
-    Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
-    Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
+    Trajectory.U.at(NextStep)[d] = Trajectory.U.at(NextStep-1)[d];
+    Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/type->mass); // v = F/m * t
   }
   // Update R (and F)
@@ -134 +134 @@
   *TotalMass += type->mass;  // sum up total mass
   for(int d=0;d<NDIM;d++) {
-    TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
+    TotalVelocity->at(d) += Trajectory.U.at(Step)[d]*type->mass;
   }
 };
@@ -145 +145 @@
 void TrajectoryParticle::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
 {
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
@@ -160 +160 @@
 void TrajectoryParticle::Thermostat_Gaussian_init(int Step, double *G, double *E)
 {
-  double *U = Trajectory.U.at(Step).x;
-  double *F = Trajectory.F.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
+  Vector &F = Trajectory.F.at(Step);
   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
@@ -177 +177 @@
 void TrajectoryParticle::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
 {
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
@@ -194 +194 @@
 {
   double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     // throw a dice to determine whether it gets hit by a heat bath particle
@@ -218 +218 @@
 void TrajectoryParticle::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
 {
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
@@ -233 +233 @@
 void TrajectoryParticle::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
 {
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {
@@ -248 +248 @@
 void TrajectoryParticle::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
 {
-  double *U = Trajectory.U.at(Step).x;
+  Vector &U = Trajectory.U.at(Step);
   if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
     for (int d=0; d<NDIM; d++) {

src/boundary.cpp

@@ -18 +18 @@
 #include "tesselation.hpp"
 #include "tesselationhelpers.hpp"
+#include "Plane.hpp"
 
 #include<gsl/gsl_poly.h>
@@ -80 +81 @@
               DistanceVector.SubtractVector(&Neighbour->second.second->x);
               do { // seek for neighbour pair where it flips
-                  OldComponent = DistanceVector.x[Othercomponent];
+                  OldComponent = DistanceVector[Othercomponent];
                   Neighbour++;
                   if (Neighbour == BoundaryPoints[axis].end()) // make it wrap around
@@ -88 +89 @@
                   //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;
                 } while ((runner != Neighbour) && (fabs(OldComponent / fabs(
-                  OldComponent) - DistanceVector.x[Othercomponent] / fabs(
-                  DistanceVector.x[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
+                  OldComponent) - DistanceVector[Othercomponent] / fabs(
+                  DistanceVector[Othercomponent])) < MYEPSILON)); // as long as sign does not flip
               if (runner != Neighbour) {
                   OtherNeighbour = Neighbour;
@@ -102 +103 @@
                   //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;
                   // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour
-                  w1 = fabs(OtherVector.x[Othercomponent]);
-                  w2 = fabs(DistanceVector.x[Othercomponent]);
-                  tmp = fabs((w1 * DistanceVector.x[component] + w2
-                      * OtherVector.x[component]) / (w1 + w2));
+                  w1 = fabs(OtherVector[Othercomponent]);
+                  w2 = fabs(DistanceVector[Othercomponent]);
+                  tmp = fabs((w1 * DistanceVector[component] + w2
+                      * OtherVector[component]) / (w1 + w2));
                   // mark if it has greater diameter
                   //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;
@@ -159 +160 @@
     AngleReferenceVector.Zero();
     AngleReferenceNormalVector.Zero();
-    AxisVector.x[axis] = 1.;
-    AngleReferenceVector.x[(axis + 1) % NDIM] = 1.;
-    AngleReferenceNormalVector.x[(axis + 2) % NDIM] = 1.;
+    AxisVector[axis] = 1.;
+    AngleReferenceVector[(axis + 1) % NDIM] = 1.;
+    AngleReferenceNormalVector[(axis + 2) % NDIM] = 1.;
 
     Log() << Verbose(1) << "Axisvector is " << AxisVector << " and AngleReferenceVector is " << AngleReferenceVector << ", and AngleReferenceNormalVector is " << AngleReferenceNormalVector << "." << endl;
@@ -636 +637 @@
       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.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
+      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));
       const double h = x.Norm(); // distance of CoG to triangle
@@ -751 +754 @@
     Log() << Verbose(0) << "Setting Box dimensions to minimum possible, the greatest diameters." << endl;
     for (int i = 0; i < NDIM; i++)
-      BoxLengths.x[i] = GreatestDiameter[i];
+      BoxLengths[i] = GreatestDiameter[i];
     mol->CenterEdge(&BoxLengths);
   } else {
-    BoxLengths.x[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
-    BoxLengths.x[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
-    BoxLengths.x[2] = minimumvolume - cellvolume;
+    BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
+    BoxLengths[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
+    BoxLengths[2] = minimumvolume - cellvolume;
     double x0 = 0.;
     double x1 = 0.;
     double x2 = 0.;
-    if (gsl_poly_solve_cubic(BoxLengths.x[0], BoxLengths.x[1], BoxLengths.x[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
+    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
       Log() << Verbose(0) << "RESULT: The resulting spacing is: " << x0 << " ." << endl;
     else {
@@ -769 +772 @@
     cellvolume = 1.;
     for (int i = 0; i < NDIM; i++) {
-      BoxLengths.x[i] = repetition[i] * (x0 + GreatestDiameter[i]);
-      cellvolume *= BoxLengths.x[i];
+      BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
+      cellvolume *= BoxLengths[i];
     }
 
@@ -780 +783 @@
   // update Box of atoms by boundary
   mol->SetBoxDimension(&BoxLengths);
-  Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths.x[0] << " and " << BoxLengths.x[1] << " and " << BoxLengths.x[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
+  Log() << Verbose(0) << "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
 };
 
@@ -839 +842 @@
   FillerDistance.InverseMatrixMultiplication(M);
   for(int i=0;i<NDIM;i++)
-    N[i] = (int) ceil(1./FillerDistance.x[i]);
+    N[i] = static_cast<int>(ceil(1./FillerDistance[i]));
   Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl;
 
@@ -880 +883 @@
           // create molecule random translation vector ...
           for (int i=0;i<NDIM;i++)
-            FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+            FillerTranslations[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
           Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
 
@@ -892 +895 @@
             // create atomic random translation vector ...
             for (int i=0;i<NDIM;i++)
-              AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+              AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
 
             // ... and rotation matrix

src/builder.cpp

@@ -1487 +1487 @@
                   Log() << Verbose(2) << "found element " << first->type->name << endl;
                 for (int i=NDIM;i--;)
-                  first->x.x[i] = atof(argv[argptr+1+i]);
+                  first->x[i] = atof(argv[argptr+1+i]);
                 if (first->type != NULL) {
                   mol->AddAtom(first);  // add to molecule
@@ -1837 +1837 @@
                 Log() << Verbose(1) << "Translating all ions by given vector." << endl;
                 for (int i=NDIM;i--;)
-                  x.x[i] = atof(argv[argptr+i]);
+                  x[i] = atof(argv[argptr+i]);
                 mol->Translate((const Vector *)&x);
                 argptr+=3;
@@ -1853 +1853 @@
                 Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl;
                 for (int i=NDIM;i--;)
-                  x.x[i] = atof(argv[argptr+i]);
+                  x[i] = atof(argv[argptr+i]);
                 mol->TranslatePeriodically((const Vector *)&x);
                 argptr+=3;
@@ -1875 +1875 @@
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
-                  x.x[i] = atof(argv[NDIM+i]);
+                  x[i] = atof(argv[NDIM+i]);
                   mol->cell_size[j]*=factor[i];
                 }
@@ -1936 +1936 @@
                 for (int i=0;i<NDIM;i++) {
                   j += i+1;
-                  x.x[i] = atof(argv[argptr+i]);
-                  mol->cell_size[j] += x.x[i]*2.;
+                  x[i] = atof(argv[argptr+i]);
+                  mol->cell_size[j] += x[i]*2.;
                 }
                 mol->Translate((const Vector *)&x);
@@ -2094 +2094 @@
                     x.Zero();
                     y.Zero();
-                    y.x[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
+                    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

src/config.cpp

@@ -740 +740 @@
               neues = AtomList[i][j];
             status = (status &&
-                    ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
-                    ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
-                    ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x[0], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x[1], 1, (repetition == 0) ? critical : optional) &&
+                    ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x[2], 1, (repetition == 0) ? critical : optional) &&
                     ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
             if (!status) break;
@@ -756 +756 @@
             // put into trajectories list
             for (int d=0;d<NDIM;d++)
-              neues->Trajectory.R.at(repetition).x[d] = neues->x.x[d];
+              neues->Trajectory.R.at(repetition)[d] = neues->x[d];
 
             // parse velocities if present
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
-              neues->v.x[0] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
-              neues->v.x[1] = 0.;
-            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
-              neues->v.x[2] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v[0], 1,optional))
+              neues->v[0] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v[1], 1,optional))
+              neues->v[1] = 0.;
+            if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v[2], 1,optional))
+              neues->v[2] = 0.;
             for (int d=0;d<NDIM;d++)
-              neues->Trajectory.U.at(repetition).x[d] = neues->v.x[d];
+              neues->Trajectory.U.at(repetition)[d] = neues->v[d];
 
             // parse forces if present
@@ -776 +776 @@
               value[2] = 0.;
             for (int d=0;d<NDIM;d++)
-              neues->Trajectory.F.at(repetition).x[d] = value[d];
+              neues->Trajectory.F.at(repetition)[d] = value[d];
 
   //            Log() << Verbose(0) << "Parsed position of step " << (repetition) << ": (";
@@ -819 +819 @@
             neues = AtomList[i][j];
           // then parse for each atom the coordinates as often as present
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
-          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &neues->x[0], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &neues->x[1], repetition,critical);
+          ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &neues->x[2], repetition,critical);
           ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
-            neues->v.x[0] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
-            neues->v.x[1] = 0.;
-          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
-            neues->v.x[2] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->v[0], repetition,optional))
+            neues->v[0] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->v[1], repetition,optional))
+            neues->v[1] = 0.;
+          if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->v[2], repetition,optional))
+            neues->v[2] = 0.;
           // here we don't care if forces are present (last in trajectories is always equal to current position)
           neues->type = elementhash[i]; // find element type
@@ -1289 +1289 @@
         istringstream input2(zeile);
         atom *neues = World::getInstance().createAtom();
-        input2 >> neues->x.x[0]; // x
-        input2 >> neues->x.x[1]; // y
-        input2 >> neues->x.x[2]; // z
+        input2 >> neues->x[0]; // x
+        input2 >> neues->x[1]; // y
+        input2 >> neues->x[2]; // z
         input2 >> l;
         neues->type = elementhash[No]; // find element type
@@ -1570 +1570 @@
              'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
              MolNo,         /* residue sequence number */
-             Walker->node->x[0],                 /* position X in Angstroem */
-             Walker->node->x[1],                 /* position Y in Angstroem */
-             Walker->node->x[2],                 /* position Z in Angstroem */
+             Walker->node->at(0),                 /* position X in Angstroem */
+             Walker->node->at(1),                 /* position Y in Angstroem */
+             Walker->node->at(2),                 /* position Z in Angstroem */
              (double)Walker->type->Valence,         /* occupancy */
              (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
@@ -1624 +1624 @@
            'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
            0,         /* residue sequence number */
-           Walker->node->x[0],                 /* position X in Angstroem */
-           Walker->node->x[1],                 /* position Y in Angstroem */
-           Walker->node->x[2],                 /* position Z in Angstroem */
+           Walker->node->at(0),                 /* position X in Angstroem */
+           Walker->node->at(1),                 /* position Y in Angstroem */
+           Walker->node->at(2),                 /* position Z in Angstroem */
            (double)Walker->type->Valence,         /* occupancy */
            (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
@@ -1677 +1677 @@
     *output << mol->name << "\t";
     *output << 0 << "\t";
-    *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
-    *output << (double)Walker->type->Valence << "\t";
+    *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
+    *output << static_cast<double>(Walker->type->Valence) << "\t";
     *output << Walker->type->symbol << "\t";
     for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
@@ -1752 +1752 @@
         *output << (*MolWalker)->name << "\t";
         *output << MolCounter << "\t";
-        *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
+        *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
         *output << (double)Walker->type->Valence << "\t";
         *output << Walker->type->symbol << "\t";

src/ellipsoid.cpp

@@ -116 +116 @@
   // put parameters into suitable ellipsoid form
   for (int i=0;i<3;i++) {
-    Center.x[i] = gsl_vector_get(x, i+0);
+    Center[i] = gsl_vector_get(x, i+0);
     EllipsoidLength[i] = gsl_vector_get(x, i+3);
     EllipsoidAngle[i] = gsl_vector_get(x, i+6);
@@ -160 +160 @@
     x = gsl_vector_alloc (9);
     for (int i=0;i<3;i++) {
-      gsl_vector_set (x, i+0, EllipsoidCenter->x[i]);
+      gsl_vector_set (x, i+0, EllipsoidCenter->at(i));
       gsl_vector_set (x, i+3, EllipsoidLength[i]);
       gsl_vector_set (x, i+6, EllipsoidAngle[i]);
@@ -195 +195 @@
       if (status == GSL_SUCCESS) {
         for (int i=0;i<3;i++) {
-          EllipsoidCenter->x[i] = gsl_vector_get (s->x,i+0);
+          EllipsoidCenter->at(i) = gsl_vector_get (s->x,i+0);
           EllipsoidLength[i] = gsl_vector_get (s->x, i+3);
           EllipsoidAngle[i] = gsl_vector_get (s->x, i+6);
@@ -427 +427 @@
       output << number << "\t";
       for (int i=0;i<3;i++)
-        output << setprecision(9) << EllipsoidCenter.x[i] << "\t";
+        output << setprecision(9) << EllipsoidCenter[i] << "\t";
       for (int i=0;i<3;i++)
         output << setprecision(9) << EllipsoidLength[i] << "\t";

src/gslmatrix.cpp

@@ -10 +10 @@
 #include "gslmatrix.hpp"
 #include "helpers.hpp"
+#include "Helpers/fast_functions.hpp"
 
 #include <cassert>

src/helpers.cpp

@@ -6 +6 @@
 
 #include "helpers.hpp"
+#include "Helpers/fast_functions.hpp"
 #include "log.hpp"
 #include "memoryusageobserver.hpp"
@@ -49 +50 @@
   while (*b < lower_bound)
     *b += step;
-};
-
-/** Returns the power of \a n with respect to \a base.
- * \param base basis
- * \param n power
- * \return \f$base^n\f$
- */
-int pot(int base, int n)
-{
-  int res = 1;
-  int j;
-  for (j=n;j--;)
-    res *= base;
-  return res;
 };
 
@@ -175 +162 @@
 };
 
-/** hard-coded determinant of a 3x3 matrix.
- * \param a[9] matrix
- * \return \f$det(a)\f$
- */
-double RDET3(const double a[NDIM*NDIM])
-{
-  return ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]);
-};
-
-/** hard-coded determinant of a 2x2 matrix.
- * \param a[4] matrix
- * \return \f$det(a)\f$
- */
-double RDET2(const double a[4])
-{
-  return ((a[0])*(a[3])-(a[1])*(a[2]));
-};
-
-/** hard-coded determinant of a 2x2 matrix.
- * \param a0 (0,0) entry of matrix
- * \param a1 (0,1) entry of matrix
- * \param a2 (1,0) entry of matrix
- * \param a3 (1,1) entry of matrix
- * \return \f$det(a)\f$
- */
-double RDET2(const double a0, const double a1, const double a2, const double a3)
-{
-  return ((a0)*(a3)-(a1)*(a2));
-};
+
 
 /** Comparison function for GSL heapsort on distances in two molecules.

src/helpers.hpp

@@ -24 +24 @@
 /********************************************** definitions *********************************/
 
-// some algebraic matrix stuff
-double RDET3(const double a[NDIM*NDIM]);
-double RDET2(const double a[4]);
-double RDET2(const double a0, const double a1, const double a2, const double a3);
-
 /********************************************** helpful functions *********************************/
 
@@ -53 +48 @@
 bool check_bounds(double *x, double *cell_size);
 void bound(double *b, double lower_bound, double upper_bound);
-int pot(int base, int n);
 int CountLinesinFile(ifstream &InputFile);
 char *FixedDigitNumber(const int FragmentNumber, const int digits);
@@ -64 +58 @@
 /********************************************** helpful template functions *********************************/
 
-/** Flips two values.
- * \param x first value
- * \param y second value
- */
-template <typename T> void flip(T &x, T &y)
-{
-  T tmp;
-  tmp = x;
-  x = y;
-  y = tmp;
-};
 
 /** returns greater of the two values.

src/leastsquaremin.cpp

@@ -25 +25 @@
   for (int i=num;i--;) {
     for(int j=NDIM;j--;)
-      sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
+      sum += (gsl_vector_get(x,j) - (vectors[i])->at(j))*(gsl_vector_get(x,j) - (vectors[i])->at(j));
   }
 

src/linearsystemofequations.cpp

@@ -54 +54 @@
 {
   assert ( columns == NDIM && "Vector class is always three-dimensional, unlike this LEqS!");
-  b->SetFromDoubleArray(x->x);
+  b->SetFromDoubleArray(x->get());
 };
 
@@ -100 +100 @@
   // copy solution
   for (size_t i=0;i<x->dimension;i++)
-    v.x[i] = x->Get(i);
+    v[i] = x->Get(i);
   return status;
 };

src/linkedcell.cpp

@@ -54 +54 @@
   Walker = set->GetPoint();
   for (int i=0;i<NDIM;i++) {
-    max.x[i] = Walker->node->x[i];
-    min.x[i] = Walker->node->x[i];
+    max[i] = Walker->node->at(i);
+    min[i] = Walker->node->at(i);
   }
   set->GoToFirst();
@@ -61 +61 @@
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      if (max.x[i] < Walker->node->x[i])
-        max.x[i] = Walker->node->x[i];
-      if (min.x[i] > Walker->node->x[i])
-        min.x[i] = Walker->node->x[i];
+      if (max[i] < Walker->node->at(i))
+        max[i] = Walker->node->at(i);
+      if (min[i] > Walker->node->at(i))
+        min[i] = Walker->node->at(i);
     }
     set->GoToNext();
@@ -72 +72 @@
   // 2. find then number of cells per axis
   for (int i=0;i<NDIM;i++) {
-    N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
+    N[i] = static_cast<int>(floor((max[i] - min[i])/RADIUS)+1);
   }
   Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
@@ -94 +94 @@
     Walker = set->GetPoint();
     for (int i=0;i<NDIM;i++) {
-      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
+      n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
     }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
@@ -128 +128 @@
   LinkedNodes::iterator Runner = set->begin();
   for (int i=0;i<NDIM;i++) {
-    max.x[i] = (*Runner)->node->x[i];
-    min.x[i] = (*Runner)->node->x[i];
+    max[i] = (*Runner)->node->at(i);
+    min[i] = (*Runner)->node->at(i);
   }
   for (LinkedNodes::iterator Runner = set->begin(); Runner != set->end(); Runner++) {
     Walker = *Runner;
     for (int i=0;i<NDIM;i++) {
-      if (max.x[i] < Walker->node->x[i])
-        max.x[i] = Walker->node->x[i];
-      if (min.x[i] > Walker->node->x[i])
-        min.x[i] = Walker->node->x[i];
+      if (max[i] < Walker->node->at(i))
+        max[i] = Walker->node->at(i);
+      if (min[i] > Walker->node->at(i))
+        min[i] = Walker->node->at(i);
     }
   }
@@ -144 +144 @@
   // 2. find then number of cells per axis
   for (int i=0;i<NDIM;i++) {
-    N[i] = (int)floor((max.x[i] - min.x[i])/RADIUS)+1;
+    N[i] = static_cast<int>(floor((max[i] - min[i])/RADIUS)+1);
   }
   Log() << Verbose(2) << "Number of cells per axis are " << N[0] << ", " << N[1] << " and " << N[2] << "." << endl;
@@ -165 +165 @@
     Walker = *Runner;
     for (int i=0;i<NDIM;i++) {
-      n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
+      n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
     }
     index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
@@ -252 +252 @@
   bool status = false;
   for (int i=0;i<NDIM;i++) {
-    n[i] = (int)floor((Walker->node->x[i] - min.x[i])/RADIUS);
+    n[i] = static_cast<int>(floor((Walker->node->at(i) - min[i])/RADIUS));
   }
   index = n[0] * N[1] * N[2] + n[1] * N[2] + n[2];
@@ -293 +293 @@
   bool status = true;
   for (int i=0;i<NDIM;i++) {
-    n[i] = (int)floor((x->x[i] - min.x[i])/RADIUS);
-    if (max.x[i] < x->x[i])
+    n[i] = static_cast<int>(floor((x->at(i) - min[i])/RADIUS));
+    if (max[i] < x->at(i))
       status = false;
-    if (min.x[i] > x->x[i])
+    if (min[i] > x->at(i))
       status = false;
   }

src/molecule.cpp

@@ -25 +25 @@
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -231 +233 @@
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
-      l = TopReplacement->x.x[i] - TopOrigin->x.x[i];
+      l = TopReplacement->x[i] - TopOrigin->x[i];
       if (fabs(l) > BondDistance) { // is component greater than bond distance
-        Orthovector1.x[i] = (l < 0) ? -1. : +1.;
+        Orthovector1[i] = (l < 0) ? -1. : +1.;
       } // (signs are correct, was tested!)
     }
@@ -305 +307 @@
 
         // determine the plane of these two with the *origin
-        AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
+        try {
+          Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
+        }
+        catch(LinearDependenceException &excp){
+          Log() << Verbose(0) << excp;
+          // TODO: figure out what to do with the Orthovector in this case
+          AllWentWell = false;
+        }
       } else {
         Orthovector1.GetOneNormalVector(&InBondvector);
@@ -313 +322 @@
       //Log() << Verbose(0) << endl;
       // orthogonal vector and bond vector between origin and replacement form the new plane
-      Orthovector1.MakeNormalVector(&InBondvector);
+      Orthovector1.MakeNormalTo(InBondvector);
       Orthovector1.Normalize();
       //Log() << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
@@ -345 +354 @@
       SecondOtherAtom->x.Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
-        FirstOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (sin(bondangle));
-        SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
+        FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
+        SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
       }
       FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
@@ -352 +361 @@
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
       for(int i=NDIM;i--;) { // and make relative to origin atom
-        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
-        SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
+        FirstOtherAtom->x[i] += TopOrigin->x[i];
+        SecondOtherAtom->x[i] += TopOrigin->x[i];
       }
       // ... and add to molecule
@@ -394 +403 @@
 //      Orthovector1.Output(out);
 //      Log() << Verbose(0) << endl;
-      AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
+      try{
+        Orthovector2 = Plane(InBondvector, Orthovector1,0).getNormal();
+      }
+      catch(LinearDependenceException &excp) {
+        Log() << Verbose(0) << excp;
+        AllWentWell = false;
+      }
 //      Log() << Verbose(3) << "Orthovector2: ";
 //      Orthovector2.Output(out);
@@ -514 +529 @@
     }
     for(j=NDIM;j--;) {
-      Walker->x.x[j] = x[j];
-      Walker->Trajectory.R.at(MDSteps-1).x[j] = x[j];
-      Walker->Trajectory.U.at(MDSteps-1).x[j] = 0;
-      Walker->Trajectory.F.at(MDSteps-1).x[j] = 0;
+      Walker->x[j] = x[j];
+      Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
+      Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
+      Walker->Trajectory.F.at(MDSteps-1)[j] = 0;
     }
     AddAtom(Walker);  // add to molecule
@@ -661 +676 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
-  cell_size[0] = dim->x[0];
+  cell_size[0] = dim->at(0);
   cell_size[1] = 0.;
-  cell_size[2] = dim->x[1];
+  cell_size[2] = dim->at(1);
   cell_size[3] = 0.;
   cell_size[4] = 0.;
-  cell_size[5] = dim->x[2];
+  cell_size[5] = dim->at(2);
 };
 
@@ -755 +770 @@
   for (int i=0;i<NDIM;i++) {
     j += i+1;
-    result = result && ((x->x[i] >= 0) && (x->x[i] < cell_size[j]));
+    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
   }
   //return result;
@@ -1005 +1020 @@
     DeterminePeriodicCenter(CenterOfGravity);
     OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
-    Log() << Verbose(5) << "Center of Gravity: ";
-    CenterOfGravity.Output();
-    Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";
-    OtherCenterOfGravity.Output();
-    Log() << Verbose(0) << endl;
+    Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl;
+    Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl;
     if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
       Log() << Verbose(4) << "Centers of gravity don't match." << endl;

src/molecule_dynamics.cpp

@@ -14 +14 @@
 #include "molecule.hpp"
 #include "parser.hpp"
+#include "Plane.hpp"
 
 /************************************* Functions for class molecule *********************************/
@@ -79 +80 @@
   //        Log() << Verbose(0) << trajectory2 << endl;
       // determine normal vector for both
-      normal.MakeNormalVector(&trajectory1, &trajectory2);
+      normal = Plane(trajectory1, trajectory2,0).getNormal();
       // print all vectors for debugging
   //        Log() << Verbose(0) << "Normal vector in between: ";
@@ -85 +86 @@
       // setup matrix
       for (int i=NDIM;i--;) {
-        gsl_matrix_set(A, 0, i, trajectory1.x[i]);
-        gsl_matrix_set(A, 1, i, trajectory2.x[i]);
-        gsl_matrix_set(A, 2, i, normal.x[i]);
-        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep).x[i] - Runner->Trajectory.R.at(Params.startstep).x[i]));
+        gsl_matrix_set(A, 0, i, trajectory1[i]);
+        gsl_matrix_set(A, 1, i, trajectory2[i]);
+        gsl_matrix_set(A, 2, i, normal[i]);
+        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
       }
       // solve the linear system by Householder transformations
@@ -514 +515 @@
       Sprinter = mol->AddCopyAtom(Walker);
       for (int n=NDIM;n--;) {
-        Sprinter->x.x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+        Sprinter->x[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
         // add to Trajectories
         //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
         if (step < MaxSteps) {
-          Walker->Trajectory.R.at(step).x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
-          Walker->Trajectory.U.at(step).x[n] = 0.;
-          Walker->Trajectory.F.at(step).x[n] = 0.;
+          Walker->Trajectory.R.at(step)[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
+          Walker->Trajectory.U.at(step)[n] = 0.;
+          Walker->Trajectory.F.at(step)[n] = 0.;
         }
       }
@@ -582 +583 @@
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Velocity.x[d] += Force.Matrix[0][i][d+5];
+        Velocity[d] += Force.Matrix[0][i][d+5];
       }
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+5] -= Velocity.x[d]/(double)AtomCount;
+        Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
       }
     // solve a constrained potential if we are meant to

src/molecule_fragmentation.cpp

@@ -1669 +1669 @@
     // remove bonds that are beyond bonddistance
     for(int i=NDIM;i--;)
-      Translationvector.x[i] = 0.;
+      Translationvector[i] = 0.;
     // scan all bonds
     Binder = first;
@@ -1676 +1676 @@
       Binder = Binder->next;
       for (int i=NDIM;i--;) {
-        tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
+        tmp = fabs(Binder->leftatom->x[i] - Binder->rightatom->x[i]);
         //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
         if (tmp > BondDistance) {
@@ -1690 +1690 @@
       // create translation vector from their periodically modified distance
       for (int i=NDIM;i--;) {
-        tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
+        tmp = Binder->leftatom->x[i] - Binder->rightatom->x[i];
         if (fabs(tmp) > BondDistance)
-          Translationvector.x[i] = (tmp < 0) ? +1. : -1.;
+          Translationvector[i] = (tmp < 0) ? +1. : -1.;
       }
       Translationvector.MatrixMultiplication(matrix);
       //Log() << Verbose(3) << "Translation vector is ";
-      Translationvector.Output();
-      Log() << Verbose(0) << endl;
+      Log() << Verbose(0) << Translationvector <<  endl;
       // apply to all atoms of first component via BFS
       for (int i=AtomCount;i--;)

src/molecule_geometry.cpp

@@ -69 +69 @@
   if (ptr != end) {   //list not empty?
     for (int i=NDIM;i--;) {
-      max->x[i] = ptr->x.x[i];
-      min->x[i] = ptr->x.x[i];
+      max->at(i) = ptr->x[i];
+      min->at(i) = ptr->x[i];
     }
     while (ptr->next != end) {  // continue with second if present
@@ -76 +76 @@
       //ptr->Output(1,1,out);
       for (int i=NDIM;i--;) {
-        max->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
-        min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[i];
+        max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
+        min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(i);
       }
     }
@@ -272 +272 @@
          if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
+              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
                 Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl;
                 if (tmp > 0)
-                  Translationvector.x[j] -= 1.;
+                  Translationvector[j] -= 1.;
                 else
-                  Translationvector.x[j] += 1.;
+                  Translationvector[j] += 1.;
               }
             }
@@ -286 +286 @@
         Testvector.MatrixMultiplication(matrix);
         Center.AddVector(&Testvector);
-        Log() << Verbose(1) << "vector is: ";
-        Testvector.Output();
-        Log() << Verbose(0) << endl;
+        Log() << Verbose(1) << "vector is: " << Testvector << endl;
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
@@ -298 +296 @@
             Testvector.MatrixMultiplication(matrix);
             Center.AddVector(&Testvector);
-            Log() << Verbose(1) << "Hydrogen vector is: ";
-            Testvector.Output();
-            Log() << Verbose(0) << endl;
+            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
           }
         }
@@ -336 +332 @@
     x.CopyVector(&ptr->x);
     //x.SubtractVector(CenterOfGravity);
-    InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-    InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-    InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-    InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-    InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-    InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-    InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-    InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-    InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
+    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
+    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
+    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
+    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
+    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
+    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
+    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
+    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
   }
   // print InertiaTensor for debugging
@@ -388 +384 @@
       x.CopyVector(&ptr->x);
       //x.SubtractVector(CenterOfGravity);
-      InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-      InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-      InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-      InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-      InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-      InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-      InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-      InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-      InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
+      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
+      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
+      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
+      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
+      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
+      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
     }
     // print InertiaTensor for debugging
@@ -423 +419 @@
   double alpha, tmp;
   Vector z_axis;
-  z_axis.x[0] = 0.;
-  z_axis.x[1] = 0.;
-  z_axis.x[2] = 1.;
+  z_axis[0] = 0.;
+  z_axis[1] = 0.;
+  z_axis[2] = 1.;
 
   // rotate on z-x plane
   Log() << Verbose(0) << "Begin of Aligning all atoms." << endl;
-  alpha = atan(-n->x[0]/n->x[2]);
+  alpha = atan(-n->at(0)/n->at(2));
   Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ";
   while (ptr->next != end) {
     ptr = ptr->next;
-    tmp = ptr->x.x[0];
-    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+    tmp = ptr->x[0];
+    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
+    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[0];
-      ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = ptr->Trajectory.R.at(j)[0];
+      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
+      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
     }
   }
   // rotate n vector
-  tmp = n->x[0];
-  n->x[0] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
-  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-  Log() << Verbose(1) << "alignment vector after first rotation: ";
-  n->Output();
-  Log() << Verbose(0) << endl;
+  tmp = n->at(0);
+  n->at(0) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
+  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
+  Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl;
 
   // rotate on z-y plane
   ptr = start;
-  alpha = atan(-n->x[1]/n->x[2]);
+  alpha = atan(-n->at(1)/n->at(2));
   Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ";
   while (ptr->next != end) {
     ptr = ptr->next;
-    tmp = ptr->x.x[1];
-    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+    tmp = ptr->x[1];
+    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
+    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[1];
-      ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = ptr->Trajectory.R.at(j)[1];
+      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
+      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
     }
   }
   // rotate n vector (for consistency check)
-  tmp = n->x[1];
-  n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
-  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-
-  Log() << Verbose(1) << "alignment vector after second rotation: ";
-  n->Output();
-  Log() << Verbose(1) << endl;
+  tmp = n->at(1);
+  n->at(1) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
+  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
+
+  Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl;
   Log() << Verbose(0) << "End of Aligning all atoms." << endl;
 };
@@ -490 +482 @@
 
   // initialize vectors
-  a.x[0] = gsl_vector_get(x,0);
-  a.x[1] = gsl_vector_get(x,1);
-  a.x[2] = gsl_vector_get(x,2);
-  b.x[0] = gsl_vector_get(x,3);
-  b.x[1] = gsl_vector_get(x,4);
-  b.x[2] = gsl_vector_get(x,5);
+  a[0] = gsl_vector_get(x,0);
+  a[1] = gsl_vector_get(x,1);
+  a[2] = gsl_vector_get(x,2);
+  b[0] = gsl_vector_get(x,3);
+  b[1] = gsl_vector_get(x,4);
+  b[2] = gsl_vector_get(x,5);
   // go through all atoms
   while (ptr != par->mol->end) {

src/molecule_graph.cpp

@@ -17 +17 @@
 #include "memoryallocator.hpp"
 #include "molecule.hpp"
+#include "Helpers/fast_functions.hpp"
 
 struct BFSAccounting

src/moleculelist.cpp

@@ -681 +681 @@
     for (int k = 0; k < NDIM; k++) {
       j += k + 1;
-      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
-      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
+      BoxDimension[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
+      (*ListRunner)->cell_size[j] += BoxDimension[k] * 2.;
     }
     (*ListRunner)->Translate(&BoxDimension);
@@ -908 +908 @@
   // center at set box dimensions
   mol->CenterEdge(&center);
-  mol->cell_size[0] = center.x[0];
+  mol->cell_size[0] = center[0];
   mol->cell_size[1] = 0;
-  mol->cell_size[2] = center.x[1];
+  mol->cell_size[2] = center[1];
   mol->cell_size[3] = 0;
   mol->cell_size[4] = 0;
-  mol->cell_size[5] = center.x[2];
+  mol->cell_size[5] = center[2];
   insert(mol);
 }

src/tesselation.cpp

@@ -15 +15 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "vector_ops.hpp"
 #include "verbose.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 class molecule;
@@ -258 +261 @@
       helper[i].CopyVector(node->node->node);
       helper[i].SubtractVector(&BaseLineCenter);
-      helper[i].MakeNormalVector(&BaseLine);  // we want to compare the triangle's heights' angles!
+      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;
       i++;
@@ -402 +405 @@
         Info FunctionInfo(__func__);
   // get normal vector
-  NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node);
+  NormalVector = Plane(*(endpoints[0]->node->node),
+                       *(endpoints[1]->node->node),
+                       *(endpoints[2]->node->node)).getNormal();
 
   // make it always point inward (any offset vector onto plane projected onto normal vector suffices)
@@ -428 +433 @@
   Vector helper;
 
-  if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) {
+  try {
+    *Intersection = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*MolCenter, *x);
+  }
+  catch (LinearDependenceException &excp) {
+    Log() << Verbose(1) << excp;
     eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;
     return false;
@@ -450 +459 @@
   int i=0;
   do {
-    if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) {
+    try {
+      CrossPoint = GetIntersectionOfTwoLinesOnPlane(*(endpoints[i%3]->node->node),
+                                                    *(endpoints[(i+1)%3]->node->node),
+                                                    *(endpoints[(i+2)%3]->node->node),
+                                                    *Intersection);
       helper.CopyVector(endpoints[(i+1)%3]->node->node);
       helper.SubtractVector(endpoints[i%3]->node->node);
@@ -462 +475 @@
       }
       i++;
-    } else 
+    } catch (LinearDependenceException &excp){
       break;
+    }
   } while (i<3);
   if (i==3) {
@@ -492 +506 @@
   Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl;
   GetCenter(&Direction);
-  if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) {
+  try {
+    *ClosestPoint = Plane(NormalVector, *(endpoints[0]->node->node)).GetIntersection(*x, Direction);
+  }
+  catch (LinearDependenceException &excp) {
     ClosestPoint->CopyVector(x);
   }
@@ -518 +535 @@
     Direction.SubtractVector(endpoints[i%3]->node->node);
     // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal);
-    CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node);
+
+    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);
@@ -731 +750 @@
   int counter=0;
   for (; Runner[2] != endpoints.end(); ) {
-    TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);
+    TemporaryNormal = Plane(*((*Runner[0])->node->node),
+                            *((*Runner[1])->node->node),
+                            *((*Runner[2])->node->node)).getNormal();
     for (int i=0;i<3;i++) // increase each of them
       Runner[i]++;
@@ -1220 +1241 @@
       // 2. next, we have to check whether all points reside on only one side of the triangle
       // 3. construct plane vector
-      PlaneVector.MakeNormalVector(A->second->node->node,
-          baseline->second.first->second->node->node,
-          baseline->second.second->second->node->node);
+      PlaneVector = Plane(*(A->second->node->node),
+                          *(baseline->second.first->second->node->node),
+                          *(baseline->second.second->second->node->node)).getNormal();
       Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl;
       // 4. loop over all points
@@ -1391 +1412 @@
         // vector in propagation direction (out of triangle)
         // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
-        PropagationVector.MakeNormalVector(&BaseLine, &NormalVector);
+        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!
@@ -1448 +1469 @@
             // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
             flag = true;
-            VirtualNormalVector.MakeNormalVector(baseline->second->endpoints[0]->node->node, baseline->second->endpoints[1]->node->node, target->second->node->node);
+            VirtualNormalVector = Plane(*(baseline->second->endpoints[0]->node->node),
+                                        *(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);
@@ -2063 +2086 @@
         if (List != NULL) {
           for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) {
-            if ((*Runner)->node->x[i] > maxCoordinate[i]) {
+            if ((*Runner)->node->at(i) > maxCoordinate[i]) {
               Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;
-              maxCoordinate[i] = (*Runner)->node->x[i];
+              maxCoordinate[i] = (*Runner)->node->at(i);
               MaxPoint[i] = (*Runner);
             }
@@ -2083 +2106 @@
   for (int k=0;k<NDIM;k++) {
     NormalVector.Zero();
-    NormalVector.x[k] = 1.;
+    NormalVector[k] = 1.;
     BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
     Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;
@@ -2117 +2140 @@
 
     // look in one direction of baseline for initial candidate
-    SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector);  // whether we look "left" first or "right" first is not important ...
+    SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal();  // whether we look "left" first or "right" first is not important ...
 
     // adding point 1 and point 2 and add the line between them
@@ -2336 +2359 @@
 
     // construct SearchDirection and an "outward pointer"
-    SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal);
+    SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
     helper.CopyVector(&CircleCenter);
     helper.SubtractVector(ThirdNode->node);
@@ -2979 +3002 @@
                   Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;
 
-                  if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)
-                  && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
-                  ) {
+                  try {
+                    NewNormalVector = Plane(*(CandidateLine.BaseLine->endpoints[0]->node->node),
+                                              *(CandidateLine.BaseLine->endpoints[1]->node->node),
+                                              *(Candidate->node)).getNormal();
+
+                    if(!fabs(NewNormalVector.NormSquared()) > HULLEPSILON)
+                      throw LinearDependenceException(__FILE__,__LINE__);
+
                     Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;
                     radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter);
@@ -3043 +3071 @@
                       Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;
                     }
-                  } else {
+                  }
+                  catch (LinearDependenceException &excp){
+                    Log() << Verbose(1) << excp;
                     Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
                   }
@@ -3565 +3595 @@
   Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
 
@@ -3633 +3663 @@
   int counter = 0;
   while (TesselC != SetOfNeighbours->end()) {
-    helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node);
+    helper = Plane(*((*TesselA)->node),
+                   *((*TesselB)->node),
+                   *((*TesselC)->node)).getNormal();
     Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl;
     counter++;
@@ -3670 +3702 @@
   Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;
   if (AngleZero.NormSquared() > MYEPSILON)
-    OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero);
+    OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
   else
-    OrthogonalVector.MakeNormalVector(&PlaneNormal);
+    OrthogonalVector.MakeNormalTo(PlaneNormal);
   Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;
 
@@ -4002 +4034 @@
           OrthogonalVector.CopyVector((*MiddleNode)->node);
           OrthogonalVector.SubtractVector(&OldPoint);
-          OrthogonalVector.MakeNormalVector(&Reference);
+          OrthogonalVector.MakeNormalTo(Reference);
           angle = GetAngle(Point, Reference, OrthogonalVector);
           //if (angle < M_PI)  // no wrong-sided triangles, please?

src/tesselationhelpers.cpp

@@ -14 +14 @@
 #include "tesselationhelpers.hpp"
 #include "vector.hpp"
+#include "vector_ops.hpp"
 #include "verbose.hpp"
 
@@ -53 +54 @@
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+    gsl_matrix_set(A, i, 0, a[i]);
+    gsl_matrix_set(A, i, 1, b[i]);
+    gsl_matrix_set(A, i, 2, c[i]);
   }
   m11 = DetGet(A, 1);
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, b.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    gsl_matrix_set(A, i, 1, b[i]);
+    gsl_matrix_set(A, i, 2, c[i]);
   }
   m12 = DetGet(A, 1);
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, c.x[i]);
+    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    gsl_matrix_set(A, i, 1, a[i]);
+    gsl_matrix_set(A, i, 2, c[i]);
   }
   m13 = DetGet(A, 1);
 
   for(int i=0;i<3;i++) {
-    gsl_matrix_set(A, i, 0, a.x[i]*a.x[i] + b.x[i]*b.x[i] + c.x[i]*c.x[i]);
-    gsl_matrix_set(A, i, 1, a.x[i]);
-    gsl_matrix_set(A, i, 2, b.x[i]);
+    gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    gsl_matrix_set(A, i, 1, a[i]);
+    gsl_matrix_set(A, i, 2, b[i]);
   }
   m14 = DetGet(A, 1);
@@ -83 +84 @@
     eLog() << Verbose(1) << "three points are colinear." << endl;
 
-  center->x[0] =  0.5 * m12/ m11;
-  center->x[1] = -0.5 * m13/ m11;
-  center->x[2] =  0.5 * m14/ m11;
+  center->at(0) =  0.5 * m12/ m11;
+  center->at(1) = -0.5 * m13/ m11;
+  center->at(2) =  0.5 * m14/ m11;
 
   if (fabs(a.Distance(center) - RADIUS) > MYEPSILON)
@@ -281 +282 @@
   Vector SideA,SideB,HeightA, HeightB;
   for (int i=0;i<NDIM;i++)
-    intersection.x[i] = gsl_vector_get(x, i);
+    intersection[i] = gsl_vector_get(x, i);
 
   SideA.CopyVector(&(I->x1));
@@ -336 +337 @@
     /* Starting point */
     x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1.x[0]);
-    gsl_vector_set(x, 1, point1.x[1]);
-    gsl_vector_set(x, 2, point1.x[2]);
+    gsl_vector_set(x, 0, point1[0]);
+    gsl_vector_set(x, 1, point1[1]);
+    gsl_vector_set(x, 2, point1[2]);
 
     /* Set initial step sizes to 1 */
@@ -372 +373 @@
   double t1, t2;
   for (int i = 0; i < NDIM; i++) {
-    intersection.x[i] = gsl_vector_get(s->x, i);
+    intersection[i] = gsl_vector_get(s->x, i);
   }
 
@@ -700 +701 @@
   // calculate the intersection between this projected baseline and Base
   Vector *Intersection = new Vector;
-  Intersection->GetIntersectionOfTwoLinesOnPlane(Base->endpoints[0]->node->node, Base->endpoints[1]->node->node, &NewOffset, &NewDirection, &Normal);
+  *Intersection = GetIntersectionOfTwoLinesOnPlane(*(Base->endpoints[0]->node->node),
+                                                   *(Base->endpoints[1]->node->node),
+                                                     NewOffset, NewDirection);
   Normal.CopyVector(Intersection);
   Normal.SubtractVector(Base->endpoints[0]->node->node);
@@ -747 +750 @@
       *vrmlfile << "Sphere {" << endl << "  "; // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *vrmlfile << Walker->node->x[i]-center->x[i] << " ";
+        *vrmlfile << Walker->node->at(i)-center->at(i) << " ";
       *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
       cloud->GoToNext();
@@ -757 +760 @@
       for (i=0;i<3;i++) { // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+          *vrmlfile << TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j) << " ";
         *vrmlfile << "\t";
       }
@@ -792 +795 @@
     *rasterfile << "# current virtual sphere\n";
     *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
-    *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
+    *rasterfile << "2\n  " << helper[0] << " " << helper[1] << " " << helper[2] << "\t" << 5. << "\t1 0 0\n";
     *rasterfile << "9\n  terminating special property\n";
     delete(center);
@@ -820 +823 @@
       *rasterfile << "2" << endl << "  ";  // 2 is sphere type
       for (i=0;i<NDIM;i++)
-        *rasterfile << Walker->node->x[i]-center->x[i] << " ";
+        *rasterfile << Walker->node->at(i)-center->at(i) << " ";
       *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
       cloud->GoToNext();
@@ -831 +834 @@
       for (i=0;i<3;i++) {  // print each node
         for (int j=0;j<NDIM;j++)  // and for each node all NDIM coordinates
-          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->x[j]-center->x[j] << " ";
+          *rasterfile << TriangleRunner->second->endpoints[i]->node->node->at(j)-center->at(j) << " ";
         *rasterfile << "\t";
       }
@@ -877 +880 @@
       Walker = target->second->node;
       LookupList[Walker->nr] = Counter++;
-      *tecplot << Walker->node->x[0] << " " << Walker->node->x[1] << " " << Walker->node->x[2] << " " << target->second->value << endl;
+      *tecplot << Walker->node->at(0) << " " << Walker->node->at(1) << " " << Walker->node->at(2) << " " << target->second->value << endl;
     }
     *tecplot << endl;

src/unittests/ActOnAllUnitTest.cpp

@@ -99 +99 @@
 };
 
+#if 0
+
+// Unitttest deactivated for now...
+// Reasons:
+//    Vector::MakeNormalVector has been removed and is now part of Plane class
+//    VectorList::ActOnAll is in the test directory and seems to be used only for unittests
+//    MakeNormal never depended on the input Vector except in the case of linear dependence (i.e. failure)
+//
+// TODO: Rethink what exactely is being tested at this point and introduce a new unittest for the
+//       tested behaviour. Most likely this should result in a thourough unittest of the plane class
+
 /** UnitTest for VectorList::ActOnAll() and Vector::MakeNormalVector()
  */
@@ -113 +124 @@
   // check that x and y components are zero
   for (ListOfVectors::iterator Runner = VL.Vectors.begin(); Runner != VL.Vectors.end(); Runner++) {
-    CPPUNIT_ASSERT_EQUAL( (*Runner)->x[0] , 0. );
-    CPPUNIT_ASSERT_EQUAL( (*Runner)->x[1] , 0. );
+    CPPUNIT_ASSERT_EQUAL( (*Runner)->at(0) , 0. );
+    CPPUNIT_ASSERT_EQUAL( (*Runner)->at(1) , 0. );
   }
 };
+
+#endif

src/unittests/ActOnAllUnitTest.hpp

@@ -20 +20 @@
     CPPUNIT_TEST ( AddSubtractTest );
     CPPUNIT_TEST ( ScaleTest );
+#if 0
+    // test deactivated. See .cpp for reasons and workarounds
     CPPUNIT_TEST ( NormalizeTest );
+#endif
     CPPUNIT_TEST_SUITE_END();
 

src/unittests/vectorunittest.cpp

@@ -17 +17 @@
 #include "memoryusageobserver.hpp"
 #include "vector.hpp"
+#include "vector_ops.hpp"
 #include "vectorunittest.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 #ifdef HAVE_TESTRUNNER
@@ -195 +198 @@
 {
   // plane at (0,0,0) normal to (1,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionWithPlane(&unit, &zero, &zero, &two) );
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(unit, zero).GetIntersection(zero, two) );
   CPPUNIT_ASSERT_EQUAL( zero, fixture );
 
   // plane at (2,1,0) normal to (0,1,0) cuts line from (1,0,0) to (0,1,1) at ???
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionWithPlane(&otherunit, &two, &unit, &notunit) );
+  CPPUNIT_ASSERT_NO_THROW(fixture = Plane(otherunit, two).GetIntersection( unit, notunit) );
   CPPUNIT_ASSERT_EQUAL( Vector(0., 1., 1.), fixture );
 
   // four vectors equal to zero
-  CPPUNIT_ASSERT_EQUAL( false, fixture.GetIntersectionOfTwoLinesOnPlane(&zero, &zero, &zero, &zero, NULL) );
+  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(zero, zero, zero, zero), LinearDependenceException);
+  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // four vectors equal to unit
+  CPPUNIT_ASSERT_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, unit, unit, unit), LinearDependenceException);
+  //CPPUNIT_ASSERT_EQUAL( zero, fixture );
+
+  // two equal lines
+  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, two));
+  CPPUNIT_ASSERT_EQUAL( unit, fixture );
+
+  // line from (1,0,0) to (2,1,0) cuts line from (1,0,0) to (0,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, unit, otherunit) );
+  CPPUNIT_ASSERT_EQUAL( unit, fixture );
+
+  // line from (1,0,0) to (0,0,0) cuts line from (0,0,0) to (2,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW( fixture = GetIntersectionOfTwoLinesOnPlane(unit, zero, zero, two) );
   CPPUNIT_ASSERT_EQUAL( zero, fixture );
 
-  // four vectors equal to unit
-  CPPUNIT_ASSERT_EQUAL( false, fixture.GetIntersectionOfTwoLinesOnPlane(&unit, &unit, &unit, &unit, NULL) );
+  // line from (1,0,0) to (2,1,0) cuts line from (0,0,0) to (0,1,0) at ???
+  CPPUNIT_ASSERT_NO_THROW(fixture = GetIntersectionOfTwoLinesOnPlane(unit, two, zero, otherunit) );
+  CPPUNIT_ASSERT_EQUAL( Vector(0., -1., 0.), fixture );
+};
+
+/** UnitTest for vector rotations.
+ */
+void VectorTest::VectorRotationTest()
+{
+  fixture.Init(-1.,0.,0.);
+
+  // zero vector does not change
+  fixture = RotateVector(zero,unit, 1.);
   CPPUNIT_ASSERT_EQUAL( zero, fixture );
 
-  // two equal lines
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionOfTwoLinesOnPlane(&unit, &two, &unit, &two, NULL) );
+  fixture = RotateVector(zero, two, 1.);
+  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
+
+  // vector on axis does not change
+  fixture = RotateVector(unit,unit, 1.);
   CPPUNIT_ASSERT_EQUAL( unit, fixture );
 
-  // line from (1,0,0) to (2,1,0) cuts line from (1,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionOfTwoLinesOnPlane(&unit, &two, &unit, &otherunit, NULL) );
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  // line from (1,0,0) to (0,0,0) cuts line from (0,0,0) to (2,1,0) at ???
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionOfTwoLinesOnPlane(&unit, &zero, &zero, &two, NULL) );
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  // line from (1,0,0) to (2,1,0) cuts line from (0,0,0) to (0,1,0) at ???
-  CPPUNIT_ASSERT_EQUAL( true, fixture.GetIntersectionOfTwoLinesOnPlane(&unit, &two, &zero, &otherunit, NULL) );
-  CPPUNIT_ASSERT_EQUAL( Vector(0., -1., 0.), fixture );
-};
-
-/** UnitTest for vector rotations.
- */
-void VectorTest::VectorRotationTest()
-{
-  fixture.Init(-1.,0.,0.);
-
-  // zero vector does not change
-  fixture.CopyVector(&zero);
-  fixture.RotateVector(&unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero, fixture );
-
-  fixture.RotateVector(&two, 1.);
-  CPPUNIT_ASSERT_EQUAL( zero,  fixture);
-
-  // vector on axis does not change
-  fixture.CopyVector(&unit);
-  fixture.RotateVector(&unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
-  fixture.RotateVector(&unit, 1.);
-  CPPUNIT_ASSERT_EQUAL( unit, fixture );
-
   // rotations
-  fixture.CopyVector(&otherunit);
-  fixture.RotateVector(&unit, M_PI);
+  fixture = RotateVector(otherunit, unit, M_PI);
   CPPUNIT_ASSERT_EQUAL( Vector(0.,-1.,0.), fixture );
 
-  fixture.CopyVector(&otherunit);
-  fixture.RotateVector(&unit, 2. * M_PI);
+  fixture = RotateVector(otherunit, unit, 2. * M_PI);
   CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
 
-  fixture.CopyVector(&otherunit);
-  fixture.RotateVector(&unit, 0);
+  fixture = RotateVector(otherunit,unit, 0);
   CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
 
-  fixture.Init(0.,0.,1.);
-  fixture.RotateVector(&notunit, M_PI);
+  fixture = RotateVector(Vector(0.,0.,1.), notunit, M_PI);
   CPPUNIT_ASSERT_EQUAL( otherunit, fixture );
 }

src/vector.cpp

@@ -7 +7 @@
 
 #include "defs.hpp"
-#include "helpers.hpp"
-#include "info.hpp"
 #include "gslmatrix.hpp"
 #include "leastsquaremin.hpp"
-#include "log.hpp"
 #include "memoryallocator.hpp"
 #include "vector.hpp"
-#include "verbose.hpp"
+#include "Helpers/fast_functions.hpp"
+#include "Helpers/Assert.hpp"
+#include "Plane.hpp"
+#include "Exceptions/LinearDependenceException.hpp"
 
 #include <gsl/gsl_linalg.h>
@@ -25 +25 @@
 /** Constructor of class vector.
  */
-Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
+Vector::Vector()
+{
+  x[0] = x[1] = x[2] = 0.;
+};
 
 /** Constructor of class vector.
  */
-Vector::Vector(const double x1, const double x2, const double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
+Vector::Vector(const double x1, const double x2, const double x3)
+{
+  x[0] = x1;
+  x[1] = x2;
+  x[2] = x3;
+};
+
+/**
+ * Copy constructor
+ */
+Vector::Vector(const Vector& src)
+{
+  x[0] = src[0];
+  x[1] = src[1];
+  x[2] = src[2];
+}
+
+/**
+ * Assignment operator
+ */
+Vector& Vector::operator=(const Vector& src){
+  // check for self assignment
+  if(&src!=this){
+    x[0] = src[0];
+    x[1] = src[1];
+    x[2] = src[2];
+  }
+  return *this;
+}
 
 /** Desctructor of class vector.
@@ -208 +239 @@
 };
 
-/** Calculates the intersection point between a line defined by \a *LineVector and \a *LineVector2 and a plane defined by \a *Normal and \a *PlaneOffset.
- * According to [Bronstein] the vectorial plane equation is:
- *   -# \f$\stackrel{r}{\rightarrow} \cdot \stackrel{N}{\rightarrow} + D = 0\f$,
- * where \f$\stackrel{r}{\rightarrow}\f$ is the vector to be testet, \f$\stackrel{N}{\rightarrow}\f$ is the plane's normal vector and
- * \f$D = - \stackrel{a}{\rightarrow} \stackrel{N}{\rightarrow}\f$, the offset with respect to origin, if \f$\stackrel{a}{\rightarrow}\f$,
- * is an offset vector onto the plane. The line is parametrized by \f$\stackrel{x}{\rightarrow} + k \stackrel{t}{\rightarrow}\f$, where
- * \f$\stackrel{x}{\rightarrow}\f$ is the offset and \f$\stackrel{t}{\rightarrow}\f$ the directional vector (NOTE: No need to normalize
- * the latter). Inserting the parametrized form into the plane equation and solving for \f$k\f$, which we insert then into the parametrization
- * of the line yields the intersection point on the plane.
- * \param *out output stream for debugging
- * \param *PlaneNormal Plane's normal vector
- * \param *PlaneOffset Plane's offset vector
- * \param *Origin first vector of line
- * \param *LineVector second vector of line
- * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
- */
-bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
-{
-  Info FunctionInfo(__func__);
-  double factor;
-  Vector Direction, helper;
-
-  // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
-  Direction.CopyVector(LineVector);
-  Direction.SubtractVector(Origin);
-  Direction.Normalize();
-  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
-  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
-  factor = Direction.ScalarProduct(PlaneNormal);
-  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
-    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
-    return false;
-  }
-  helper.CopyVector(PlaneOffset);
-  helper.SubtractVector(Origin);
-  factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
-    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
-    CopyVector(Origin);
-    return true;
-  }
-  //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
-  Direction.Scale(factor);
-  CopyVector(Origin);
-  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
-  AddVector(&Direction);
-
-  // test whether resulting vector really is on plane
-  helper.CopyVector(this);
-  helper.SubtractVector(PlaneOffset);
-  if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
-    return true;
-  } else {
-    eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl;
-    return false;
-  }
-};
-
 /** Calculates the minimum distance of this vector to the plane.
  * \param *out output stream for debugging
@@ -280 +252 @@
   temp.CopyVector(this);
   temp.SubtractVector(PlaneOffset);
-  temp.MakeNormalVector(PlaneNormal);
+  temp.MakeNormalTo(*PlaneNormal);
   temp.Scale(-1.);
   // then add connecting vector from plane to point
@@ -292 +264 @@
 
   return (temp.Norm()*sign);
-};
-
-/** Calculates the intersection of the two lines that are both on the same plane.
- * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html 
- * \param *out output stream for debugging
- * \param *Line1a first vector of first line
- * \param *Line1b second vector of first line
- * \param *Line2a first vector of second line
- * \param *Line2b second vector of second line
- * \param *PlaneNormal normal of plane, is supplemental/arbitrary
- * \return true - \a this will contain the intersection on return, false - lines are parallel
- */
-bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
-{
-  Info FunctionInfo(__func__);
-
-  GSLMatrix *M = new GSLMatrix(4,4);
-
-  M->SetAll(1.);
-  for (int i=0;i<3;i++) {
-    M->Set(0, i, Line1a->x[i]);
-    M->Set(1, i, Line1b->x[i]);
-    M->Set(2, i, Line2a->x[i]);
-    M->Set(3, i, Line2b->x[i]);
-  }
-  
-  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    cout << "\t" << M->Get(i,j);
-  //  cout << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
-    return false;
-  }
-  delete(M);
-  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
-
-
-  // constuct a,b,c
-  Vector a;
-  Vector b;
-  Vector c;
-  Vector d;
-  a.CopyVector(Line1b);
-  a.SubtractVector(Line1a);
-  b.CopyVector(Line2b);
-  b.SubtractVector(Line2a);
-  c.CopyVector(Line2a);
-  c.SubtractVector(Line1a);
-  d.CopyVector(Line2b);
-  d.SubtractVector(Line1b);
-  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   Zero();
-   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
-   return false;
-  }
-
-  // check for parallelity
-  Vector parallel;
-  double factor = 0.;
-  if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
-    parallel.CopyVector(Line1a);
-    parallel.SubtractVector(Line2a);
-    factor = parallel.ScalarProduct(&a)/a.Norm();
-    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-      CopyVector(Line2a);
-      Log() << Verbose(1) << "Lines conincide." << endl;
-      return true;
-    } else {
-      parallel.CopyVector(Line1a);
-      parallel.SubtractVector(Line2b);
-      factor = parallel.ScalarProduct(&a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        CopyVector(Line2b);
-        Log() << Verbose(1) << "Lines conincide." << endl;
-        return true;
-      }
-    }
-    Log() << Verbose(1) << "Lines are parallel." << endl;
-    Zero();
-    return false;
-  }
-
-  // obtain s
-  double s;
-  Vector temp1, temp2;
-  temp1.CopyVector(&c);
-  temp1.VectorProduct(&b);
-  temp2.CopyVector(&a);
-  temp2.VectorProduct(&b);
-  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
-
-  // construct intersection
-  CopyVector(&a);
-  Scale(s);
-  AddVector(Line1a);
-  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
-
-  return true;
 };
 
@@ -538 +403 @@
 };
 
-/** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
- * \param *axis rotation axis
- * \param alpha rotation angle in radian
- */
-void Vector::RotateVector(const Vector * const axis, const double alpha)
-{
-  Vector a,y;
-  // normalise this vector with respect to axis
-  a.CopyVector(this);
-  a.ProjectOntoPlane(axis);
-  // construct normal vector
-  bool rotatable = y.MakeNormalVector(axis,&a);
-  // The normal vector cannot be created if there is linar dependency.
-  // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
-  if (!rotatable) {
-    return;
-  }
-  y.Scale(Norm());
-  // scale normal vector by sine and this vector by cosine
-  y.Scale(sin(alpha));
-  a.Scale(cos(alpha));
-  CopyVector(Projection(axis));
-  // add scaled normal vector onto this vector
-  AddVector(&y);
-  // add part in axis direction
-  AddVector(&a);
-};
+
+double& Vector::operator[](size_t i){
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
+}
+
+const double& Vector::operator[](size_t i) const{
+  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
+  return x[i];
+}
+
+double& Vector::at(size_t i){
+  return (*this)[i];
+}
+
+const double& Vector::at(size_t i) const{
+  return (*this)[i];
+}
+
+double* Vector::get(){
+  return x;
+}
 
 /** Compares vector \a to vector \a b component-wise.
@@ -575 +435 @@
   bool status = true;
   for (int i=0;i<NDIM;i++)
-    status = status && (fabs(a.x[i] - b.x[i]) < MYEPSILON);
+    status = status && (fabs(a[i] - b[i]) < MYEPSILON);
   return status;
 };
@@ -660 +520 @@
 };
 
-/** Prints a 3dim vector.
- * prints no end of line.
- */
-void Vector::Output() const
-{
-  Log() << Verbose(0) << "(";
-  for (int i=0;i<NDIM;i++) {
-    Log() << Verbose(0) << x[i];
-    if (i != 2)
-      Log() << Verbose(0) << ",";
-  }
-  Log() << Verbose(0) << ")";
-};
-
 ostream& operator<<(ostream& ost, const Vector& m)
 {
   ost << "(";
   for (int i=0;i<NDIM;i++) {
-    ost << m.x[i];
+    ost << m[i];
     if (i != 2)
       ost << ",";
@@ -752 +598 @@
  * \param *matrix NDIM_NDIM array
  */
-void Vector::InverseMatrixMultiplication(const double * const A)
+bool Vector::InverseMatrixMultiplication(const double * const A)
 {
   Vector C;
@@ -779 +625 @@
     for (int i=NDIM;i--;)
       x[i] = C.x[i];
+    return true;
   } else {
-    eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl;
+    return false;
   }
 };
@@ -806 +653 @@
   projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
-  Log() << Verbose(1) << "Vector: ";
-  Output();
-  Log() << Verbose(0) << "\t";
   for (int i=NDIM;i--;)
     x[i] -= 2.*projection*n->x[i];
-  Log() << Verbose(0) << "Projected vector: ";
-  Output();
-  Log() << Verbose(0) << endl;
-};
-
-/** Calculates normal vector for three given vectors (being three points in space).
- * Makes this vector orthonormal to the three given points, making up a place in 3d space.
- * \param *y1 first vector
- * \param *y2 second vector
- * \param *y3 third vector
- * \return true - success, vectors are linear independent, false - failure due to linear dependency
- */
-bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2, const Vector * const y3)
-{
-  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)) {
-    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
-    return false;
-  }
-//  Log() << Verbose(4) << "relative, first plane coordinates:";
-//  x1.Output((ofstream *)&cout);
-//  Log() << Verbose(0) << endl;
-//  Log() << Verbose(4) << "second plane coordinates:";
-//  x2.Output((ofstream *)&cout);
-//  Log() << Verbose(0) << endl;
-
-  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
-  this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
-  this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
-  Normalize();
-
-  return true;
-};
-
-
-/** Calculates orthonormal vector to two given vectors.
- * Makes this vector orthonormal to two given vectors. This is very similar to the other
- * vector::MakeNormalVector(), only there three points whereas here two difference
- * vectors are given.
- * \param *x1 first vector
- * \param *x2 second vector
- * \return true - success, vectors are linear independent, false - failure due to linear dependency
- */
-bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2)
-{
-  Vector x1,x2;
-  x1.CopyVector(y1);
-  x2.CopyVector(y2);
-  Zero();
-  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
-    return false;
-  }
-//  Log() << Verbose(4) << "relative, first plane coordinates:";
-//  x1.Output((ofstream *)&cout);
-//  Log() << Verbose(0) << endl;
-//  Log() << Verbose(4) << "second plane coordinates:";
-//  x2.Output((ofstream *)&cout);
-//  Log() << Verbose(0) << endl;
-
-  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
-  this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
-  this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
-  Normalize();
-
-  return true;
-};
-
-/** Calculates orthonormal vector to one given vectors.
+};
+
+
+/** Calculates orthonormal vector to one given vector.
  * Just subtracts the projection onto the given vector from this vector.
  * The removed part of the vector is Vector::Projection()
@@ -890 +664 @@
  * \return true - success, false - vector is zero
  */
-bool Vector::MakeNormalVector(const Vector * const y1)
+bool Vector::MakeNormalTo(const Vector &y1)
 {
   bool result = false;
-  double factor = y1->ScalarProduct(this)/y1->NormSquared();
+  double factor = y1.ScalarProduct(this)/y1.NormSquared();
   Vector x1;
-  x1.CopyVector(y1);
+  x1.CopyVector(&y1);
   x1.Scale(factor);
   SubtractVector(&x1);
@@ -917 +691 @@
   double norm;
 
-  Log() << Verbose(4);
-  GivenVector->Output();
-  Log() << Verbose(0) << endl;
   for (j=NDIM;j--;)
     Components[j] = -1;
@@ -926 +697 @@
     if (fabs(GivenVector->x[j]) > MYEPSILON)
       Components[Last++] = j;
-  Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
 
   switch(Last) {
@@ -966 +736 @@
 };
 
-/** Creates a new vector as the one with least square distance to a given set of \a vectors.
- * \param *vectors set of vectors
- * \param num number of vectors
- * \return true if success, false if failed due to linear dependency
- */
-bool Vector::LSQdistance(const Vector **vectors, int num)
-{
-  int j;
-
-  for (j=0;j<num;j++) {
-    Log() << Verbose(1) << j << "th atom's vector: ";
-    (vectors[j])->Output();
-    Log() << Verbose(0) << endl;
-  }
-
-  int np = 3;
-  struct LSQ_params par;
-
-   const gsl_multimin_fminimizer_type *T =
-     gsl_multimin_fminimizer_nmsimplex;
-   gsl_multimin_fminimizer *s = NULL;
-   gsl_vector *ss, *y;
-   gsl_multimin_function minex_func;
-
-   size_t iter = 0, i;
-   int status;
-   double size;
-
-   /* Initial vertex size vector */
-   ss = gsl_vector_alloc (np);
-   y = gsl_vector_alloc (np);
-
-   /* Set all step sizes to 1 */
-   gsl_vector_set_all (ss, 1.0);
-
-   /* Starting point */
-   par.vectors = vectors;
-   par.num = num;
-
-   for (i=NDIM;i--;)
-    gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
-
-   /* Initialize method and iterate */
-   minex_func.f = &LSQ;
-   minex_func.n = np;
-   minex_func.params = (void *)&par;
-
-   s = gsl_multimin_fminimizer_alloc (T, np);
-   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
-
-   do
-     {
-       iter++;
-       status = gsl_multimin_fminimizer_iterate(s);
-
-       if (status)
-         break;
-
-       size = gsl_multimin_fminimizer_size (s);
-       status = gsl_multimin_test_size (size, 1e-2);
-
-       if (status == GSL_SUCCESS)
-         {
-           printf ("converged to minimum at\n");
-         }
-
-       printf ("%5d ", (int)iter);
-       for (i = 0; i < (size_t)np; i++)
-         {
-           printf ("%10.3e ", gsl_vector_get (s->x, i));
-         }
-       printf ("f() = %7.3f size = %.3f\n", s->fval, size);
-     }
-   while (status == GSL_CONTINUE && iter < 100);
-
-  for (i=(size_t)np;i--;)
-    this->x[i] = gsl_vector_get(s->x, i);
-   gsl_vector_free(y);
-   gsl_vector_free(ss);
-   gsl_multimin_fminimizer_free (s);
-
-  return true;
-};
 
 /** Adds vector \a *y componentwise.
@@ -1092 +779 @@
 }
 
-
-/** Asks for position, checks for boundary.
- * \param cell_size unitary size of cubic cell, coordinates must be within 0...cell_size
- * \param check whether bounds shall be checked (true) or not (false)
- */
-void Vector::AskPosition(const double * const cell_size, const bool check)
-{
-  char coords[3] = {'x','y','z'};
-  int j = -1;
-  for (int i=0;i<3;i++) {
-    j += i+1;
-    do {
-      Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
-      cin >> x[i];
-    } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
-  }
-};
-
+// this function is completely unused so it is deactivated until new uses arrive and a new
+// place can be found
+#if 0
 /** Solves a vectorial system consisting of two orthogonal statements and a norm statement.
  * This is linear system of equations to be solved, however of the three given (skp of this vector\
@@ -1281 +953 @@
 };
 
+#endif
+
 /**
  * Checks whether this vector is within the parallelepiped defined by the given three vectors and

src/vector.hpp

@@ -24 +24 @@
 class Vector {
   public:
-    double x[NDIM];
 
   Vector();
   Vector(const double x1, const double x2, const double x3);
+  Vector(const Vector& src);
   ~Vector();
 
@@ -48 +48 @@
   void CopyVector(const Vector * const y);
   void CopyVector(const Vector &y);
-  void RotateVector(const Vector * const y, const double alpha);
   void VectorProduct(const Vector * const y);
   void ProjectOntoPlane(const Vector * const y);
@@ -63 +62 @@
   void Scale(const double factor);
   void MatrixMultiplication(const double * const M);
-  void InverseMatrixMultiplication(const double * const M);
+  bool InverseMatrixMultiplication(const double * const M);
   void KeepPeriodic(const double * const matrix);
   void LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors);
   double CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const;
-  bool GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector);
-  bool GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *Normal = NULL);
   bool GetOneNormalVector(const Vector * const x1);
-  bool MakeNormalVector(const Vector * const y1);
-  bool MakeNormalVector(const Vector * const y1, const Vector * const y2);
-  bool MakeNormalVector(const Vector * const x1, const Vector * const x2, const Vector * const x3);
-  bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
-  bool LSQdistance(const Vector ** vectors, int dim);
-  void AskPosition(const double * const cell_size, const bool check);
-  void Output() const;
+  bool MakeNormalTo(const Vector &y1);
+  //bool SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c);
   bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
   void WrapPeriodically(const double * const M, const double * const Minv);
+
+  // Accessors ussually come in pairs... and sometimes even more than that
+  double& operator[](size_t i);
+  const double& operator[](size_t i) const;
+  double& at(size_t i);
+  const double& at(size_t i) const;
+
+  // Assignment operator
+  Vector &operator=(const Vector& src);
+
+  // Access to internal structure
+  double* get();
+
+private:
+  double x[NDIM];
 
 };