Ignore:
Timestamp:
Apr 29, 2010, 1:55:21 PM (16 years ago)
Author:
Tillmann Crueger <crueger@…>
Children:
070651, 5d1a94
Parents:
90c4460 (diff), 32842d8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'VectorRefactoring' into StructureRefactoring

Conflicts:

molecuilder/src/Legacy/oldmenu.cpp
molecuilder/src/Makefile.am
molecuilder/src/analysis_correlation.cpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/molecule.cpp
molecuilder/src/molecule_fragmentation.cpp
molecuilder/src/molecule_geometry.cpp
molecuilder/src/molecule_graph.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/tesselationhelpers.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/vector.cpp

    r90c4460 r0d111b  
    66
    77
    8 #include "defs.hpp"
    9 #include "helpers.hpp"
    10 #include "info.hpp"
    11 #include "gslmatrix.hpp"
    12 #include "leastsquaremin.hpp"
    13 #include "log.hpp"
    14 #include "memoryallocator.hpp"
    158#include "vector.hpp"
    169#include "verbose.hpp"
    1710#include "World.hpp"
    18 
    19 #include <gsl/gsl_linalg.h>
    20 #include <gsl/gsl_matrix.h>
    21 #include <gsl/gsl_permutation.h>
    22 #include <gsl/gsl_vector.h>
     11#include "Helpers/Assert.hpp"
     12#include "Helpers/fast_functions.hpp"
     13
     14#include <iostream>
     15
     16using namespace std;
     17
    2318
    2419/************************************ Functions for class vector ************************************/
    … …  
    2621/** Constructor of class vector.
    2722 */
    28 Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
     23Vector::Vector()
     24{
     25  x[0] = x[1] = x[2] = 0.;
     26};
     27
     28/**
     29 * Copy constructor
     30 */
     31
     32Vector::Vector(const Vector& src)
     33{
     34  x[0] = src[0];
     35  x[1] = src[1];
     36  x[2] = src[2];
     37}
    2938
    3039/** Constructor of class vector.
    3140 */
    32 Vector::Vector(const Vector * const a)
    33 {
    34   x[0] = a->x[0];
    35   x[1] = a->x[1];
    36   x[2] = a->x[2];
    37 };
    38 
    39 /** Constructor of class vector.
    40  */
    41 Vector::Vector(const Vector &a)
    42 {
    43   x[0] = a.x[0];
    44   x[1] = a.x[1];
    45   x[2] = a.x[2];
    46 };
    47 
    48 /** Constructor of class vector.
    49  */
    50 Vector::Vector(const double x1, const double x2, const double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
     41Vector::Vector(const double x1, const double x2, const double x3)
     42{
     43  x[0] = x1;
     44  x[1] = x2;
     45  x[2] = x3;
     46};
     47
     48/**
     49 * Assignment operator
     50 */
     51Vector& Vector::operator=(const Vector& src){
     52  // check for self assignment
     53  if(&src!=this){
     54    x[0] = src[0];
     55    x[1] = src[1];
     56    x[2] = src[2];
     57  }
     58  return *this;
     59}
    5160
    5261/** Desctructor of class vector.
    … …  
    5867 * \return \f$| x - y |^2\f$
    5968 */
    60 double Vector::DistanceSquared(const Vector * const y) const
     69double Vector::DistanceSquared(const Vector &y) const
    6170{
    6271  double res = 0.;
    6372  for (int i=NDIM;i--;)
    64     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
     73    res += (x[i]-y[i])*(x[i]-y[i]);
    6574  return (res);
    6675};
    … …  
    7079 * \return \f$| x - y |\f$
    7180 */
    72 double Vector::Distance(const Vector * const y) const
    73 {
    74   double res = 0.;
    75   for (int i=NDIM;i--;)
    76     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
    77   return (sqrt(res));
     81double Vector::Distance(const Vector &y) const
     82{
     83  return (sqrt(DistanceSquared(y)));
    7884};
    7985
    … …  
    8389 * \return \f$| x - y |\f$
    8490 */
    85 double Vector::PeriodicDistance(const Vector * const y, const double * const cell_size) const
     91double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
    8692{
    8793  double res = Distance(y), tmp, matrix[NDIM*NDIM];
    88   Vector Shiftedy, TranslationVector;
    89   int N[NDIM];
    90   matrix[0] = cell_size[0];
    91   matrix[1] = cell_size[1];
    92   matrix[2] = cell_size[3];
    93   matrix[3] = cell_size[1];
    94   matrix[4] = cell_size[2];
    95   matrix[5] = cell_size[4];
    96   matrix[6] = cell_size[3];
    97   matrix[7] = cell_size[4];
    98   matrix[8] = cell_size[5];
    99   // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    100   for (N[0]=-1;N[0]<=1;N[0]++)
    101     for (N[1]=-1;N[1]<=1;N[1]++)
    102       for (N[2]=-1;N[2]<=1;N[2]++) {
    103         // create the translation vector
    104         TranslationVector.Zero();
    105         for (int i=NDIM;i--;)
    106           TranslationVector.x[i] = (double)N[i];
    107         TranslationVector.MatrixMultiplication(matrix);
    108         // add onto the original vector to compare with
    109         Shiftedy.CopyVector(y);
    110         Shiftedy.AddVector(&TranslationVector);
    111         // get distance and compare with minimum so far
    112         tmp = Distance(&Shiftedy);
    113         if (tmp < res) res = tmp;
    114       }
    115   return (res);
     94    Vector Shiftedy, TranslationVector;
     95    int N[NDIM];
     96    matrix[0] = cell_size[0];
     97    matrix[1] = cell_size[1];
     98    matrix[2] = cell_size[3];
     99    matrix[3] = cell_size[1];
     100    matrix[4] = cell_size[2];
     101    matrix[5] = cell_size[4];
     102    matrix[6] = cell_size[3];
     103    matrix[7] = cell_size[4];
     104    matrix[8] = cell_size[5];
     105    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     106    for (N[0]=-1;N[0]<=1;N[0]++)
     107      for (N[1]=-1;N[1]<=1;N[1]++)
     108        for (N[2]=-1;N[2]<=1;N[2]++) {
     109          // create the translation vector
     110          TranslationVector.Zero();
     111          for (int i=NDIM;i--;)
     112            TranslationVector[i] = (double)N[i];
     113          TranslationVector.MatrixMultiplication(matrix);
     114          // add onto the original vector to compare with
     115          Shiftedy = y + TranslationVector;
     116          // get distance and compare with minimum so far
     117          tmp = Distance(Shiftedy);
     118          if (tmp < res) res = tmp;
     119        }
     120    return (res);
    116121};
    117122
    … …  
    121126 * \return \f$| x - y |^2\f$
    122127 */
    123 double Vector::PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const
     128double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
    124129{
    125130  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    126   Vector Shiftedy, TranslationVector;
    127   int N[NDIM];
    128   matrix[0] = cell_size[0];
    129   matrix[1] = cell_size[1];
    130   matrix[2] = cell_size[3];
    131   matrix[3] = cell_size[1];
    132   matrix[4] = cell_size[2];
    133   matrix[5] = cell_size[4];
    134   matrix[6] = cell_size[3];
    135   matrix[7] = cell_size[4];
    136   matrix[8] = cell_size[5];
    137   // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    138   for (N[0]=-1;N[0]<=1;N[0]++)
    139     for (N[1]=-1;N[1]<=1;N[1]++)
    140       for (N[2]=-1;N[2]<=1;N[2]++) {
    141         // create the translation vector
    142         TranslationVector.Zero();
    143         for (int i=NDIM;i--;)
    144           TranslationVector.x[i] = (double)N[i];
    145         TranslationVector.MatrixMultiplication(matrix);
    146         // add onto the original vector to compare with
    147         Shiftedy.CopyVector(y);
    148         Shiftedy.AddVector(&TranslationVector);
    149         // get distance and compare with minimum so far
    150         tmp = DistanceSquared(&Shiftedy);
    151         if (tmp < res) res = tmp;
    152       }
    153   return (res);
     131    Vector Shiftedy, TranslationVector;
     132    int N[NDIM];
     133    matrix[0] = cell_size[0];
     134    matrix[1] = cell_size[1];
     135    matrix[2] = cell_size[3];
     136    matrix[3] = cell_size[1];
     137    matrix[4] = cell_size[2];
     138    matrix[5] = cell_size[4];
     139    matrix[6] = cell_size[3];
     140    matrix[7] = cell_size[4];
     141    matrix[8] = cell_size[5];
     142    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     143    for (N[0]=-1;N[0]<=1;N[0]++)
     144      for (N[1]=-1;N[1]<=1;N[1]++)
     145        for (N[2]=-1;N[2]<=1;N[2]++) {
     146          // create the translation vector
     147          TranslationVector.Zero();
     148          for (int i=NDIM;i--;)
     149            TranslationVector[i] = (double)N[i];
     150          TranslationVector.MatrixMultiplication(matrix);
     151          // add onto the original vector to compare with
     152          Shiftedy = y + TranslationVector;
     153          // get distance and compare with minimum so far
     154          tmp = DistanceSquared(Shiftedy);
     155          if (tmp < res) res = tmp;
     156        }
     157    return (res);
    154158};
    155159
    … …  
    160164void Vector::KeepPeriodic(const double * const matrix)
    161165{
    162 //  int N[NDIM];
    163 //  bool flag = false;
    164   //vector Shifted, TranslationVector;
    165   Vector TestVector;
    166 //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
    167 //  Log() << Verbose(2) << "Vector is: ";
    168 //  Output(out);
    169 //  Log() << Verbose(0) << endl;
    170   TestVector.CopyVector(this);
    171   TestVector.InverseMatrixMultiplication(matrix);
    172   for(int i=NDIM;i--;) { // correct periodically
    173     if (TestVector.x[i] < 0) {  // get every coefficient into the interval [0,1)
    174       TestVector.x[i] += ceil(TestVector.x[i]);
    175     } else {
    176       TestVector.x[i] -= floor(TestVector.x[i]);
     166  //  int N[NDIM];
     167  //  bool flag = false;
     168    //vector Shifted, TranslationVector;
     169  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
     170  //  Log() << Verbose(2) << "Vector is: ";
     171  //  Output(out);
     172  //  Log() << Verbose(0) << endl;
     173    InverseMatrixMultiplication(matrix);
     174    for(int i=NDIM;i--;) { // correct periodically
     175      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
     176        at(i) += ceil(at(i));
     177      } else {
     178        at(i) -= floor(at(i));
     179      }
    177180    }
    178   }
    179   TestVector.MatrixMultiplication(matrix);
    180   CopyVector(&TestVector);
    181 //  Log() << Verbose(2) << "New corrected vector is: ";
    182 //  Output(out);
    183 //  Log() << Verbose(0) << endl;
    184 //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
     181    MatrixMultiplication(matrix);
     182  //  Log() << Verbose(2) << "New corrected vector is: ";
     183  //  Output(out);
     184  //  Log() << Verbose(0) << endl;
     185  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
    185186};
    186187
    … …  
    189190 * \return \f$\langle x, y \rangle\f$
    190191 */
    191 double Vector::ScalarProduct(const Vector * const y) const
     192double Vector::ScalarProduct(const Vector &y) const
    192193{
    193194  double res = 0.;
    194195  for (int i=NDIM;i--;)
    195     res += x[i]*y->x[i];
     196    res += x[i]*y[i];
    196197  return (res);
    197198};
    … …  
    204205 *  \return \f$ x \times y \f&
    205206 */
    206 void Vector::VectorProduct(const Vector * const y)
     207void Vector::VectorProduct(const Vector &y)
    207208{
    208209  Vector tmp;
    209   tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
    210   tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
    211   tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
    212   this->CopyVector(&tmp);
     210  tmp[0] = x[1]* (y[2]) - x[2]* (y[1]);
     211  tmp[1] = x[2]* (y[0]) - x[0]* (y[2]);
     212  tmp[2] = x[0]* (y[1]) - x[1]* (y[0]);
     213  (*this) = tmp;
    213214};
    214215
    … …  
    218219 * \return \f$\langle x, y \rangle\f$
    219220 */
    220 void Vector::ProjectOntoPlane(const Vector * const y)
     221void Vector::ProjectOntoPlane(const Vector &y)
    221222{
    222223  Vector tmp;
    223   tmp.CopyVector(y);
     224  tmp = y;
    224225  tmp.Normalize();
    225   tmp.Scale(ScalarProduct(&tmp));
    226   this->SubtractVector(&tmp);
    227 };
    228 
    229 /** Calculates the intersection point between a line defined by \a *LineVector and \a *LineVector2 and a plane defined by \a *Normal and \a *PlaneOffset.
    230  * According to [Bronstein] the vectorial plane equation is:
    231  *   -# \f$\stackrel{r}{\rightarrow} \cdot \stackrel{N}{\rightarrow} + D = 0\f$,
    232  * where \f$\stackrel{r}{\rightarrow}\f$ is the vector to be testet, \f$\stackrel{N}{\rightarrow}\f$ is the plane's normal vector and
    233  * \f$D = - \stackrel{a}{\rightarrow} \stackrel{N}{\rightarrow}\f$, the offset with respect to origin, if \f$\stackrel{a}{\rightarrow}\f$,
    234  * is an offset vector onto the plane. The line is parametrized by \f$\stackrel{x}{\rightarrow} + k \stackrel{t}{\rightarrow}\f$, where
    235  * \f$\stackrel{x}{\rightarrow}\f$ is the offset and \f$\stackrel{t}{\rightarrow}\f$ the directional vector (NOTE: No need to normalize
    236  * the latter). Inserting the parametrized form into the plane equation and solving for \f$k\f$, which we insert then into the parametrization
    237  * of the line yields the intersection point on the plane.
    238  * \param *out output stream for debugging
    239  * \param *PlaneNormal Plane's normal vector
    240  * \param *PlaneOffset Plane's offset vector
    241  * \param *Origin first vector of line
    242  * \param *LineVector second vector of line
    243  * \return true -  \a this contains intersection point on return, false - line is parallel to plane (even if in-plane)
    244  */
    245 bool Vector::GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
    246 {
    247   Info FunctionInfo(__func__);
    248   double factor;
    249   Vector Direction, helper;
    250 
    251   // find intersection of a line defined by Offset and Direction with a  plane defined by triangle
    252   Direction.CopyVector(LineVector);
    253   Direction.SubtractVector(Origin);
    254   Direction.Normalize();
    255   DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);
    256   //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
    257   factor = Direction.ScalarProduct(PlaneNormal);
    258   if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
    259     DoLog(1) && (Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl);
    260     return false;
    261   }
    262   helper.CopyVector(PlaneOffset);
    263   helper.SubtractVector(Origin);
    264   factor = helper.ScalarProduct(PlaneNormal)/factor;
    265   if (fabs(factor) < MYEPSILON) { // Origin is in-plane
    266     DoLog(1) && (Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl);
    267     CopyVector(Origin);
    268     return true;
    269   }
    270   //factor = Origin->ScalarProduct(PlaneNormal)*(-PlaneOffset->ScalarProduct(PlaneNormal))/(Direction.ScalarProduct(PlaneNormal));
    271   Direction.Scale(factor);
    272   CopyVector(Origin);
    273   DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);
    274   AddVector(&Direction);
    275 
    276   // test whether resulting vector really is on plane
    277   helper.CopyVector(this);
    278   helper.SubtractVector(PlaneOffset);
    279   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
    280     DoLog(1) && (Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl);
    281     return true;
    282   } else {
    283     DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);
    284     return false;
    285   }
     226  tmp.Scale(ScalarProduct(tmp));
     227  *this -= tmp;
    286228};
    287229
    … …  
    290232 * \param *PlaneNormal normal of plane
    291233 * \param *PlaneOffset offset of plane
     234 * \return distance to plane
    292235 * \return distance vector onto to plane
    293236 */
    294 Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    295 {
    296   Vector temp;
    297 
    298   // first create part that is orthonormal to PlaneNormal with withdraw
    299   temp.CopyVector(this);
    300   temp.SubtractVector(PlaneOffset);
    301   temp.MakeNormalVector(PlaneNormal);
     237Vector Vector::GetDistanceVectorToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
     238{
     239  Vector temp = (*this) - PlaneOffset;
     240  temp.MakeNormalTo(PlaneNormal);
    302241  temp.Scale(-1.);
    303242  // then add connecting vector from plane to point
    304   temp.AddVector(this);
    305   temp.SubtractVector(PlaneOffset);
     243  temp += (*this)-PlaneOffset;
    306244  double sign = temp.ScalarProduct(PlaneNormal);
    307245  if (fabs(sign) > MYEPSILON)
    … …  
    314252  return temp;
    315253};
     254
    316255
    317256/** Calculates the minimum distance of this vector to the plane.
    … …  
    322261 * \return distance to plane
    323262 */
    324 double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     263double Vector::DistanceToPlane(const Vector &PlaneNormal, const Vector &PlaneOffset) const
    325264{
    326265  return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
    327 };
    328 
    329 /** Calculates the intersection of the two lines that are both on the same plane.
    330  * This is taken from Weisstein, Eric W. "Line-Line Intersection." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Line-LineIntersection.html
    331  * \param *out output stream for debugging
    332  * \param *Line1a first vector of first line
    333  * \param *Line1b second vector of first line
    334  * \param *Line2a first vector of second line
    335  * \param *Line2b second vector of second line
    336  * \param *PlaneNormal normal of plane, is supplemental/arbitrary
    337  * \return true - \a this will contain the intersection on return, false - lines are parallel
    338  */
    339 bool Vector::GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
    340 {
    341   Info FunctionInfo(__func__);
    342 
    343   GSLMatrix *M = new GSLMatrix(4,4);
    344 
    345   M->SetAll(1.);
    346   for (int i=0;i<3;i++) {
    347     M->Set(0, i, Line1a->x[i]);
    348     M->Set(1, i, Line1b->x[i]);
    349     M->Set(2, i, Line2a->x[i]);
    350     M->Set(3, i, Line2b->x[i]);
    351   }
    352  
    353   //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
    354   //ostream &output = Log() << Verbose(1);
    355   //for (int i=0;i<4;i++) {
    356   //  for (int j=0;j<4;j++)
    357   //    output << "\t" << M->Get(i,j);
    358   //  output << endl;
    359   //}
    360   if (fabs(M->Determinant()) > MYEPSILON) {
    361     DoLog(1) && (Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl);
    362     return false;
    363   }
    364   delete(M);
    365   DoLog(1) && (Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl);
    366 
    367 
    368   // constuct a,b,c
    369   Vector a;
    370   Vector b;
    371   Vector c;
    372   Vector d;
    373   a.CopyVector(Line1b);
    374   a.SubtractVector(Line1a);
    375   b.CopyVector(Line2b);
    376   b.SubtractVector(Line2a);
    377   c.CopyVector(Line2a);
    378   c.SubtractVector(Line1a);
    379   d.CopyVector(Line2b);
    380   d.SubtractVector(Line1b);
    381   DoLog(1) && (Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl);
    382   if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
    383    Zero();
    384    DoLog(1) && (Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl);
    385    return false;
    386   }
    387 
    388   // check for parallelity
    389   Vector parallel;
    390   double factor = 0.;
    391   if (fabs(a.ScalarProduct(&b)*a.ScalarProduct(&b)/a.NormSquared()/b.NormSquared() - 1.) < MYEPSILON) {
    392     parallel.CopyVector(Line1a);
    393     parallel.SubtractVector(Line2a);
    394     factor = parallel.ScalarProduct(&a)/a.Norm();
    395     if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    396       CopyVector(Line2a);
    397       DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
    398       return true;
    399     } else {
    400       parallel.CopyVector(Line1a);
    401       parallel.SubtractVector(Line2b);
    402       factor = parallel.ScalarProduct(&a)/a.Norm();
    403       if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    404         CopyVector(Line2b);
    405         DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
    406         return true;
    407       }
    408     }
    409     DoLog(1) && (Log() << Verbose(1) << "Lines are parallel." << endl);
    410     Zero();
    411     return false;
    412   }
    413 
    414   // obtain s
    415   double s;
    416   Vector temp1, temp2;
    417   temp1.CopyVector(&c);
    418   temp1.VectorProduct(&b);
    419   temp2.CopyVector(&a);
    420   temp2.VectorProduct(&b);
    421   DoLog(1) && (Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl);
    422   if (fabs(temp2.NormSquared()) > MYEPSILON)
    423     s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
    424   else
    425     s = 0.;
    426   DoLog(1) && (Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl);
    427 
    428   // construct intersection
    429   CopyVector(&a);
    430   Scale(s);
    431   AddVector(Line1a);
    432   DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);
    433 
    434   return true;
    435266};
    436267
    … …  
    438269 * \param *y array to second vector
    439270 */
    440 void Vector::ProjectIt(const Vector * const y)
    441 {
    442   Vector helper(*y);
    443   helper.Scale(-(ScalarProduct(y)));
    444   AddVector(&helper);
     271void Vector::ProjectIt(const Vector &y)
     272{
     273  (*this) += (-ScalarProduct(y))*y;
    445274};
    446275
    … …  
    449278 * \return Vector
    450279 */
    451 Vector Vector::Projection(const Vector * const y) const
    452 {
    453   Vector helper(*y);
    454   helper.Scale((ScalarProduct(y)/y->NormSquared()));
     280Vector Vector::Projection(const Vector &y) const
     281{
     282  Vector helper = y;
     283  helper.Scale((ScalarProduct(y)/y.NormSquared()));
    455284
    456285  return helper;
    … …  
    462291double Vector::Norm() const
    463292{
    464   double res = 0.;
    465   for (int i=NDIM;i--;)
    466     res += this->x[i]*this->x[i];
    467   return (sqrt(res));
     293  return (sqrt(NormSquared()));
    468294};
    469295
    … …  
    473299double Vector::NormSquared() const
    474300{
    475   return (ScalarProduct(this));
     301  return (ScalarProduct(*this));
    476302};
    477303
    … …  
    480306void Vector::Normalize()
    481307{
    482   double res = 0.;
    483   for (int i=NDIM;i--;)
    484     res += this->x[i]*this->x[i];
    485   if (fabs(res) > MYEPSILON)
    486     res = 1./sqrt(res);
    487   Scale(&res);
     308  double factor = Norm();
     309  (*this) *= 1/factor;
    488310};
    489311
    … …  
    492314void Vector::Zero()
    493315{
    494   for (int i=NDIM;i--;)
    495     this->x[i] = 0.;
     316  at(0)=at(1)=at(2)=0;
    496317};
    497318
    … …  
    500321void Vector::One(const double one)
    501322{
    502   for (int i=NDIM;i--;)
    503     this->x[i] = one;
    504 };
    505 
    506 /** Initialises all components of this vector.
    507  */
    508 void Vector::Init(const double x1, const double x2, const double x3)
    509 {
    510   x[0] = x1;
    511   x[1] = x2;
    512   x[2] = x3;
     323  at(0)=at(1)=at(2)=one;
    513324};
    514325
    … …  
    532343 * @return true - vector is normalized, false - vector is not
    533344 */
    534 bool Vector::IsNormalTo(const Vector * const normal) const
     345bool Vector::IsNormalTo(const Vector &normal) const
    535346{
    536347  if (ScalarProduct(normal) < MYEPSILON)
    … …  
    543354 * @return true - vector is normalized, false - vector is not
    544355 */
    545 bool Vector::IsEqualTo(const Vector * const a) const
     356bool Vector::IsEqualTo(const Vector &a) const
    546357{
    547358  bool status = true;
    548359  for (int i=0;i<NDIM;i++) {
    549     if (fabs(x[i] - a->x[i]) > MYEPSILON)
     360    if (fabs(x[i] - a[i]) > MYEPSILON)
    550361      status = false;
    551362  }
    … …  
    557368 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
    558369 */
    559 double Vector::Angle(const Vector * const y) const
    560 {
    561   double norm1 = Norm(), norm2 = y->Norm();
     370double Vector::Angle(const Vector &y) const
     371{
     372  double norm1 = Norm(), norm2 = y.Norm();
    562373  double angle = -1;
    563374  if ((fabs(norm1) > MYEPSILON) && (fabs(norm2) > MYEPSILON))
    … …  
    572383};
    573384
    574 /** Rotates the vector relative to the origin around the axis given by \a *axis by an angle of \a alpha.
    575  * \param *axis rotation axis
    576  * \param alpha rotation angle in radian
    577  */
    578 void Vector::RotateVector(const Vector * const axis, const double alpha)
    579 {
    580   Vector a,y;
    581   // normalise this vector with respect to axis
    582   a.CopyVector(this);
    583   a.ProjectOntoPlane(axis);
    584   // construct normal vector
    585   bool rotatable = y.MakeNormalVector(axis,&a);
    586   // The normal vector cannot be created if there is linar dependency.
    587   // Then the vector to rotate is on the axis and any rotation leads to the vector itself.
    588   if (!rotatable) {
    589     return;
    590   }
    591   y.Scale(Norm());
    592   // scale normal vector by sine and this vector by cosine
    593   y.Scale(sin(alpha));
    594   a.Scale(cos(alpha));
    595   CopyVector(Projection(axis));
    596   // add scaled normal vector onto this vector
    597   AddVector(&y);
    598   // add part in axis direction
    599   AddVector(&a);
    600 };
     385
     386double& Vector::operator[](size_t i){
     387  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
     388  return x[i];
     389}
     390
     391const double& Vector::operator[](size_t i) const{
     392  ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
     393  return x[i];
     394}
     395
     396double& Vector::at(size_t i){
     397  return (*this)[i];
     398}
     399
     400const double& Vector::at(size_t i) const{
     401  return (*this)[i];
     402}
     403
     404double* Vector::get(){
     405  return x;
     406}
    601407
    602408/** Compares vector \a to vector \a b component-wise.
    … …  
    605411 * \return a == b
    606412 */
    607 bool operator==(const Vector& a, const Vector& b)
    608 {
    609   bool status = true;
    610   for (int i=0;i<NDIM;i++)
    611     status = status && (fabs(a.x[i] - b.x[i]) < MYEPSILON);
    612   return status;
     413bool Vector::operator==(const Vector& b) const
     414{
     415  return IsEqualTo(b);
    613416};
    614417
    … …  
    618421 * \return lhs + a
    619422 */
    620 const Vector& operator+=(Vector& a, const Vector& b)
    621 {
    622   a.AddVector(&b);
    623   return a;
     423const Vector& Vector::operator+=(const Vector& b)
     424{
     425  this->AddVector(b);
     426  return *this;
    624427};
    625428
    … …  
    629432 * \return lhs - a
    630433 */
    631 const Vector& operator-=(Vector& a, const Vector& b)
    632 {
    633   a.SubtractVector(&b);
    634   return a;
     434const Vector& Vector::operator-=(const Vector& b)
     435{
     436  this->SubtractVector(b);
     437  return *this;
    635438};
    636439
    … …  
    651454 * \return a + b
    652455 */
    653 Vector const operator+(const Vector& a, const Vector& b)
    654 {
    655   Vector x(a);
    656   x.AddVector(&b);
     456Vector const Vector::operator+(const Vector& b) const
     457{
     458  Vector x = *this;
     459  x.AddVector(b);
    657460  return x;
    658461};
    … …  
    663466 * \return a - b
    664467 */
    665 Vector const operator-(const Vector& a, const Vector& b)
    666 {
    667   Vector x(a);
    668   x.SubtractVector(&b);
     468Vector const Vector::operator-(const Vector& b) const
     469{
     470  Vector x = *this;
     471  x.SubtractVector(b);
    669472  return x;
    670473};
    … …  
    694497};
    695498
    696 /** Prints a 3dim vector.
    697  * prints no end of line.
    698  */
    699 void Vector::Output() const
    700 {
    701   DoLog(0) && (Log() << Verbose(0) << "(");
    702   for (int i=0;i<NDIM;i++) {
    703     DoLog(0) && (Log() << Verbose(0) << x[i]);
    704     if (i != 2)
    705       DoLog(0) && (Log() << Verbose(0) << ",");
    706   }
    707   DoLog(0) && (Log() << Verbose(0) << ")");
    708 };
    709 
    710499ostream& operator<<(ostream& ost, const Vector& m)
    711500{
    712501  ost << "(";
    713502  for (int i=0;i<NDIM;i++) {
    714     ost << m.x[i];
     503    ost << m[i];
    715504    if (i != 2)
    716505      ost << ",";
    … …  
    720509};
    721510
    722 /** Scales each atom coordinate by an individual \a factor.
    723  * \param *factor pointer to scaling factor
    724  */
    725 void Vector::Scale(const double ** const factor)
     511
     512void Vector::ScaleAll(const double *factor)
    726513{
    727514  for (int i=NDIM;i--;)
    728     x[i] *= (*factor)[i];
    729 };
    730 
    731 void Vector::Scale(const double * const factor)
    732 {
    733   for (int i=NDIM;i--;)
    734     x[i] *= *factor;
    735 };
     515    x[i] *= factor[i];
     516};
     517
     518
    736519
    737520void Vector::Scale(const double factor)
    … …  
    739522  for (int i=NDIM;i--;)
    740523    x[i] *= factor;
    741 };
    742 
    743 /** Translate atom by given vector.
    744  * \param trans[] translation vector.
    745  */
    746 void Vector::Translate(const Vector * const trans)
    747 {
    748   for (int i=NDIM;i--;)
    749     x[i] += trans->x[i];
    750524};
    751525
    … …  
    773547void Vector::MatrixMultiplication(const double * const M)
    774548{
    775   Vector C;
    776549  // do the matrix multiplication
    777   C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    778   C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
    779   C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    780   // transfer the result into this
    781   for (int i=NDIM;i--;)
    782     x[i] = C.x[i];
     550  at(0) = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
     551  at(1) = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
     552  at(2) = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    783553};
    784554
    … …  
    786556 * \param *matrix NDIM_NDIM array
    787557 */
    788 void Vector::InverseMatrixMultiplication(const double * const A)
    789 {
    790   Vector C;
     558bool Vector::InverseMatrixMultiplication(const double * const A)
     559{
    791560  double B[NDIM*NDIM];
    792561  double detA = RDET3(A);
    … …  
    807576
    808577    // do the matrix multiplication
    809     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
    810     C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
    811     C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    812     // transfer the result into this
    813     for (int i=NDIM;i--;)
    814       x[i] = C.x[i];
     578    at(0) = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
     579    at(1) = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
     580    at(2) = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
     581
     582    return true;
    815583  } else {
    816     DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);
     584    return false;
    817585  }
    818586};
    … …  
    826594 * \param *factors three-component vector with the factor for each given vector
    827595 */
    828 void Vector::LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors)
    829 {
    830   for(int i=NDIM;i--;)
    831     x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
     596void Vector::LinearCombinationOfVectors(const Vector &x1, const Vector &x2, const Vector &x3, const double * const factors)
     597{
     598  (*this) = (factors[0]*x1) +
     599            (factors[1]*x2) +
     600            (factors[2]*x3);
    832601};
    833602
    … …  
    835604 * \param n[] normal vector of mirror plane.
    836605 */
    837 void Vector::Mirror(const Vector * const n)
     606void Vector::Mirror(const Vector &n)
    838607{
    839608  double projection;
    840   projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
     609  projection = ScalarProduct(n)/n.NormSquared();    // remove constancy from n (keep as logical one)
    841610  // withdraw projected vector twice from original one
    842   DoLog(1) && (Log() << Verbose(1) << "Vector: ");
    843   Output();
    844   DoLog(0) && (Log() << Verbose(0) << "\t");
    845611  for (int i=NDIM;i--;)
    846     x[i] -= 2.*projection*n->x[i];
    847   DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");
    848   Output();
    849   DoLog(0) && (Log() << Verbose(0) << endl);
    850 };
    851 
    852 /** Calculates normal vector for three given vectors (being three points in space).
    853  * Makes this vector orthonormal to the three given points, making up a place in 3d space.
    854  * \param *y1 first vector
    855  * \param *y2 second vector
    856  * \param *y3 third vector
    857  * \return true - success, vectors are linear independent, false - failure due to linear dependency
    858  */
    859 bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2, const Vector * const y3)
    860 {
    861   Vector x1, x2;
    862 
    863   x1.CopyVector(y1);
    864   x1.SubtractVector(y2);
    865   x2.CopyVector(y3);
    866   x2.SubtractVector(y2);
    867   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    868     DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
    869     return false;
    870   }
    871 //  Log() << Verbose(4) << "relative, first plane coordinates:";
    872 //  x1.Output((ofstream *)&cout);
    873 //  Log() << Verbose(0) << endl;
    874 //  Log() << Verbose(4) << "second plane coordinates:";
    875 //  x2.Output((ofstream *)&cout);
    876 //  Log() << Verbose(0) << endl;
    877 
    878   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    879   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    880   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    881   Normalize();
    882 
    883   return true;
    884 };
    885 
    886 
    887 /** Calculates orthonormal vector to two given vectors.
    888  * Makes this vector orthonormal to two given vectors. This is very similar to the other
    889  * vector::MakeNormalVector(), only there three points whereas here two difference
    890  * vectors are given.
    891  * \param *x1 first vector
    892  * \param *x2 second vector
    893  * \return true - success, vectors are linear independent, false - failure due to linear dependency
    894  */
    895 bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2)
    896 {
    897   Vector x1,x2;
    898   x1.CopyVector(y1);
    899   x2.CopyVector(y2);
    900   Zero();
    901   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    902     DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
    903     return false;
    904   }
    905 //  Log() << Verbose(4) << "relative, first plane coordinates:";
    906 //  x1.Output((ofstream *)&cout);
    907 //  Log() << Verbose(0) << endl;
    908 //  Log() << Verbose(4) << "second plane coordinates:";
    909 //  x2.Output((ofstream *)&cout);
    910 //  Log() << Verbose(0) << endl;
    911 
    912   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    913   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    914   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    915   Normalize();
    916 
    917   return true;
     612    at(i) -= 2.*projection*n[i];
    918613};
    919614
    … …  
    924619 * \return true - success, false - vector is zero
    925620 */
    926 bool Vector::MakeNormalVector(const Vector * const y1)
     621bool Vector::MakeNormalTo(const Vector &y1)
    927622{
    928623  bool result = false;
    929   double factor = y1->ScalarProduct(this)/y1->NormSquared();
     624  double factor = y1.ScalarProduct(*this)/y1.NormSquared();
    930625  Vector x1;
    931   x1.CopyVector(y1);
    932   x1.Scale(factor);
    933   SubtractVector(&x1);
     626  x1 = factor * y1;
     627  SubtractVector(x1);
    934628  for (int i=NDIM;i--;)
    935629    result = result || (fabs(x[i]) > MYEPSILON);
    … …  
    944638 * \return true - success, false - failure (null vector given)
    945639 */
    946 bool Vector::GetOneNormalVector(const Vector * const GivenVector)
     640bool Vector::GetOneNormalVector(const Vector &GivenVector)
    947641{
    948642  int Components[NDIM]; // contains indices of non-zero components
    … …  
    951645  double norm;
    952646
    953   DoLog(4) && (Log() << Verbose(4));
    954   GivenVector->Output();
    955   DoLog(0) && (Log() << Verbose(0) << endl);
    956647  for (j=NDIM;j--;)
    957648    Components[j] = -1;
    958649  // find two components != 0
    959650  for (j=0;j<NDIM;j++)
    960     if (fabs(GivenVector->x[j]) > MYEPSILON)
     651    if (fabs(GivenVector[j]) > MYEPSILON)
    961652      Components[Last++] = j;
    962   DoLog(4) && (Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl);
    963653
    964654  switch(Last) {
    965655    case 3:  // threecomponent system
    966656    case 2:  // two component system
    967       norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
     657      norm = sqrt(1./(GivenVector[Components[1]]*GivenVector[Components[1]]) + 1./(GivenVector[Components[0]]*GivenVector[Components[0]]));
    968658      x[Components[2]] = 0.;
    969659      // in skp both remaining parts shall become zero but with opposite sign and third is zero
    970       x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
    971       x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
     660      x[Components[1]] = -1./GivenVector[Components[1]] / norm;
     661      x[Components[0]] = 1./GivenVector[Components[0]] / norm;
    972662      return true;
    973663      break;
    … …  
    984674};
    985675
    986 /** Determines parameter needed to multiply this vector to obtain intersection point with plane defined by \a *A, \a *B and \a *C.
    987  * \param *A first plane vector
    988  * \param *B second plane vector
    989  * \param *C third plane vector
    990  * \return scaling parameter for this vector
    991  */
    992 double Vector::CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const
    993 {
    994 //  Log() << Verbose(3) << "For comparison: ";
    995 //  Log() << Verbose(0) << "A " << A->Projection(this) << "\t";
    996 //  Log() << Verbose(0) << "B " << B->Projection(this) << "\t";
    997 //  Log() << Verbose(0) << "C " << C->Projection(this) << "\t";
    998 //  Log() << Verbose(0) << endl;
    999   return A->ScalarProduct(this);
    1000 };
    1001 
    1002 /** Creates a new vector as the one with least square distance to a given set of \a vectors.
    1003  * \param *vectors set of vectors
    1004  * \param num number of vectors
    1005  * \return true if success, false if failed due to linear dependency
    1006  */
    1007 bool Vector::LSQdistance(const Vector **vectors, int num)
    1008 {
    1009   int j;
    1010 
    1011   for (j=0;j<num;j++) {
    1012     DoLog(1) && (Log() << Verbose(1) << j << "th atom's vector: ");
    1013     (vectors[j])->Output();
    1014     DoLog(0) && (Log() << Verbose(0) << endl);
    1015   }
    1016 
    1017   int np = 3;
    1018   struct LSQ_params par;
    1019 
    1020    const gsl_multimin_fminimizer_type *T =
    1021      gsl_multimin_fminimizer_nmsimplex;
    1022    gsl_multimin_fminimizer *s = NULL;
    1023    gsl_vector *ss, *y;
    1024    gsl_multimin_function minex_func;
    1025 
    1026    size_t iter = 0, i;
    1027    int status;
    1028    double size;
    1029 
    1030    /* Initial vertex size vector */
    1031    ss = gsl_vector_alloc (np);
    1032    y = gsl_vector_alloc (np);
    1033 
    1034    /* Set all step sizes to 1 */
    1035    gsl_vector_set_all (ss, 1.0);
    1036 
    1037    /* Starting point */
    1038    par.vectors = vectors;
    1039    par.num = num;
    1040 
    1041    for (i=NDIM;i--;)
    1042     gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
    1043 
    1044    /* Initialize method and iterate */
    1045    minex_func.f = &LSQ;
    1046    minex_func.n = np;
    1047    minex_func.params = (void *)&par;
    1048 
    1049    s = gsl_multimin_fminimizer_alloc (T, np);
    1050    gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
    1051 
    1052    do
    1053      {
    1054        iter++;
    1055        status = gsl_multimin_fminimizer_iterate(s);
    1056 
    1057        if (status)
    1058          break;
    1059 
    1060        size = gsl_multimin_fminimizer_size (s);
    1061        status = gsl_multimin_test_size (size, 1e-2);
    1062 
    1063        if (status == GSL_SUCCESS)
    1064          {
    1065            printf ("converged to minimum at\n");
    1066          }
    1067 
    1068        printf ("%5d ", (int)iter);
    1069        for (i = 0; i < (size_t)np; i++)
    1070          {
    1071            printf ("%10.3e ", gsl_vector_get (s->x, i));
    1072          }
    1073        printf ("f() = %7.3f size = %.3f\n", s->fval, size);
    1074      }
    1075    while (status == GSL_CONTINUE && iter < 100);
    1076 
    1077   for (i=(size_t)np;i--;)
    1078     this->x[i] = gsl_vector_get(s->x, i);
    1079    gsl_vector_free(y);
    1080    gsl_vector_free(ss);
    1081    gsl_multimin_fminimizer_free (s);
    1082 
    1083   return true;
    1084 };
    1085 
    1086676/** Adds vector \a *y componentwise.
    1087677 * \param *y vector
    1088678 */
    1089 void Vector::AddVector(const Vector * const y)
    1090 {
    1091   for (int i=NDIM;i--;)
    1092     this->x[i] += y->x[i];
     679void Vector::AddVector(const Vector &y)
     680{
     681  for(int i=NDIM;i--;)
     682    x[i] += y[i];
    1093683}
    1094684
    … …  
    1096686 * \param *y vector
    1097687 */
    1098 void Vector::SubtractVector(const Vector * const y)
    1099 {
    1100   for (int i=NDIM;i--;)
    1101     this->x[i] -= y->x[i];
    1102 }
    1103 
    1104 /** Copy vector \a *y componentwise.
    1105  * \param *y vector
    1106  */
    1107 void Vector::CopyVector(const Vector * const y)
    1108 {
    1109   // check for self assignment
    1110   if(y!=this){
    1111     for (int i=NDIM;i--;)
    1112       this->x[i] = y->x[i];
    1113   }
    1114 }
    1115 
    1116 /** Copy vector \a y componentwise.
    1117  * \param y vector
    1118  */
    1119 void Vector::CopyVector(const Vector &y)
    1120 {
    1121   // check for self assignment
    1122   if(&y!=this) {
    1123     for (int i=NDIM;i--;)
    1124       this->x[i] = y.x[i];
    1125   }
    1126 }
    1127 
    1128 
    1129 /** Asks for position, checks for boundary.
    1130  * \param cell_size unitary size of cubic cell, coordinates must be within 0...cell_size
    1131  * \param check whether bounds shall be checked (true) or not (false)
    1132  */
    1133 void Vector::AskPosition(const double * const cell_size, const bool check)
    1134 {
    1135   char coords[3] = {'x','y','z'};
    1136   int j = -1;
    1137   for (int i=0;i<3;i++) {
    1138     j += i+1;
    1139     do {
    1140       DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");
    1141       cin >> x[i];
    1142     } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
    1143   }
    1144 };
    1145 
    1146 /** Solves a vectorial system consisting of two orthogonal statements and a norm statement.
    1147  * This is linear system of equations to be solved, however of the three given (skp of this vector\
    1148  * with either of the three hast to be zero) only two are linear independent. The third equation
    1149  * is that the vector should be of magnitude 1 (orthonormal). This all leads to a case-based solution
    1150  * where very often it has to be checked whether a certain value is zero or not and thus forked into
    1151  * another case.
    1152  * \param *x1 first vector
    1153  * \param *x2 second vector
    1154  * \param *y third vector
    1155  * \param alpha first angle
    1156  * \param beta second angle
    1157  * \param c norm of final vector
    1158  * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
    1159  * \bug this is not yet working properly
    1160  */
    1161 bool Vector::SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c)
    1162 {
    1163   double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    1164   double ang; // angle on testing
    1165   double sign[3];
    1166   int i,j,k;
    1167   A = cos(alpha) * x1->Norm() * c;
    1168   B1 = cos(beta + M_PI/2.) * y->Norm() * c;
    1169   B2 = cos(beta) * x2->Norm() * c;
    1170   C = c * c;
    1171   DoLog(2) && (Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl);
    1172   int flag = 0;
    1173   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    1174     if (fabs(x1->x[1]) > MYEPSILON) {
    1175       flag = 1;
    1176     } else if (fabs(x1->x[2]) > MYEPSILON) {
    1177        flag = 2;
    1178     } else {
    1179       return false;
    1180     }
    1181   }
    1182   switch (flag) {
    1183     default:
    1184     case 0:
    1185       break;
    1186     case 2:
    1187       flip(x1->x[0],x1->x[1]);
    1188       flip(x2->x[0],x2->x[1]);
    1189       flip(y->x[0],y->x[1]);
    1190       //flip(x[0],x[1]);
    1191       flip(x1->x[1],x1->x[2]);
    1192       flip(x2->x[1],x2->x[2]);
    1193       flip(y->x[1],y->x[2]);
    1194       //flip(x[1],x[2]);
    1195     case 1:
    1196       flip(x1->x[0],x1->x[1]);
    1197       flip(x2->x[0],x2->x[1]);
    1198       flip(y->x[0],y->x[1]);
    1199       //flip(x[0],x[1]);
    1200       flip(x1->x[1],x1->x[2]);
    1201       flip(x2->x[1],x2->x[2]);
    1202       flip(y->x[1],y->x[2]);
    1203       //flip(x[1],x[2]);
    1204       break;
    1205   }
    1206   // now comes the case system
    1207   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
    1208   D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    1209   D3 = y->x[0]/x1->x[0]*A-B1;
    1210   DoLog(2) && (Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n");
    1211   if (fabs(D1) < MYEPSILON) {
    1212     DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");
    1213     if (fabs(D2) > MYEPSILON) {
    1214       DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");
    1215       x[2] = -D3/D2;
    1216       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    1217       E2 = -x1->x[1]/x1->x[0];
    1218       DoLog(3) && (Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n");
    1219       F1 = E1*E1 + 1.;
    1220       F2 = -E1*E2;
    1221       F3 = E1*E1 + D3*D3/(D2*D2) - C;
    1222       DoLog(3) && (Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
    1223       if (fabs(F1) < MYEPSILON) {
    1224         DoLog(4) && (Log() << Verbose(4) << "F1 == 0!\n");
    1225         DoLog(4) && (Log() << Verbose(4) << "Gleichungssystem linear\n");
    1226         x[1] = F3/(2.*F2);
    1227       } else {
    1228         p = F2/F1;
    1229         q = p*p - F3/F1;
    1230         DoLog(4) && (Log() << Verbose(4) << "p " << p << "\tq " << q << endl);
    1231         if (q < 0) {
    1232           DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);
    1233           return false;
    1234         }
    1235         x[1] = p + sqrt(q);
    1236       }
    1237       x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    1238     } else {
    1239       DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");
    1240       return false;
    1241     }
    1242   } else {
    1243     E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    1244     E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    1245     DoLog(2) && (Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n");
    1246     F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    1247     F2 = -(E1*E2 + D2*D3/(D1*D1));
    1248     F3 = E1*E1 + D3*D3/(D1*D1) - C;
    1249     DoLog(2) && (Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
    1250     if (fabs(F1) < MYEPSILON) {
    1251       DoLog(3) && (Log() << Verbose(3) << "F1 == 0!\n");
    1252       DoLog(3) && (Log() << Verbose(3) << "Gleichungssystem linear\n");
    1253       x[2] = F3/(2.*F2);
    1254     } else {
    1255       p = F2/F1;
    1256       q = p*p - F3/F1;
    1257       DoLog(3) && (Log() << Verbose(3) << "p " << p << "\tq " << q << endl);
    1258       if (q < 0) {
    1259         DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);
    1260         return false;
    1261       }
    1262       x[2] = p + sqrt(q);
    1263     }
    1264     x[1] = (-D2 * x[2] - D3)/D1;
    1265     x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    1266   }
    1267   switch (flag) { // back-flipping
    1268     default:
    1269     case 0:
    1270       break;
    1271     case 2:
    1272       flip(x1->x[0],x1->x[1]);
    1273       flip(x2->x[0],x2->x[1]);
    1274       flip(y->x[0],y->x[1]);
    1275       flip(x[0],x[1]);
    1276       flip(x1->x[1],x1->x[2]);
    1277       flip(x2->x[1],x2->x[2]);
    1278       flip(y->x[1],y->x[2]);
    1279       flip(x[1],x[2]);
    1280     case 1:
    1281       flip(x1->x[0],x1->x[1]);
    1282       flip(x2->x[0],x2->x[1]);
    1283       flip(y->x[0],y->x[1]);
    1284       //flip(x[0],x[1]);
    1285       flip(x1->x[1],x1->x[2]);
    1286       flip(x2->x[1],x2->x[2]);
    1287       flip(y->x[1],y->x[2]);
    1288       flip(x[1],x[2]);
    1289       break;
    1290   }
    1291   // one z component is only determined by its radius (without sign)
    1292   // thus check eight possible sign flips and determine by checking angle with second vector
    1293   for (i=0;i<8;i++) {
    1294     // set sign vector accordingly
    1295     for (j=2;j>=0;j--) {
    1296       k = (i & pot(2,j)) << j;
    1297       DoLog(2) && (Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl);
    1298       sign[j] = (k == 0) ? 1. : -1.;
    1299     }
    1300     DoLog(2) && (Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n");
    1301     // apply sign matrix
    1302     for (j=NDIM;j--;)
    1303       x[j] *= sign[j];
    1304     // calculate angle and check
    1305     ang = x2->Angle (this);
    1306     DoLog(1) && (Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t");
    1307     if (fabs(ang - cos(beta)) < MYEPSILON) {
    1308       break;
    1309     }
    1310     // unapply sign matrix (is its own inverse)
    1311     for (j=NDIM;j--;)
    1312       x[j] *= sign[j];
    1313   }
    1314   return true;
    1315 };
     688void Vector::SubtractVector(const Vector &y)
     689{
     690  for(int i=NDIM;i--;)
     691    x[i] -= y[i];
     692}
    1316693
    1317694/**
    … …  
    1324701bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
    1325702{
    1326   Vector a;
    1327   a.CopyVector(this);
    1328   a.SubtractVector(&offset);
     703  Vector a = (*this)-offset;
    1329704  a.InverseMatrixMultiplication(parallelepiped);
    1330705  bool isInside = true;
    1331706
    1332707  for (int i=NDIM;i--;)
    1333     isInside = isInside && ((a.x[i] <= 1) && (a.x[i] >= 0));
     708    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
    1334709
    1335710  return isInside;
Note: See TracChangeset for help on using the changeset viewer.