Changeset 1f591b for molecuilder/src/vector.cpp

Timestamp:

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

Author:

Tillmann Crueger <crueger@…>

Children:

e7ea64

Parents:

0f55b2

Message:

Prepared interface of Vector Class for transition to VectorComposites

File:

• molecuilder/src/vector.cpp (modified) (33 diffs)

molecuilder/src/vector.cpp

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