Changes in src/vector.cpp [a67d19:717e0c]

File:

• src/vector.cpp (modified) (25 diffs)

src/vector.cpp

@@ -8 +8 @@
 #include "defs.hpp"
 #include "helpers.hpp"
-#include "info.hpp"
-#include "gslmatrix.hpp"
+#include "memoryallocator.hpp"
 #include "leastsquaremin.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "vector.hpp"
 #include "verbose.hpp"
-#include "World.hpp"
-
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_permutation.h>
-#include <gsl/gsl_vector.h>
 
 /************************************ Functions for class vector ************************************/
@@ -27 +19 @@
  */
 Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
-
-/** Constructor of class vector.
- */
-Vector::Vector(const Vector * const a)
-{
-  x[0] = a->x[0];
-  x[1] = a->x[1];
-  x[2] = a->x[2];
-};
-
-/** Constructor of class vector.
- */
-Vector::Vector(const Vector &a)
-{
-  x[0] = a.x[0];
-  x[1] = a.x[1];
-  x[2] = a.x[2];
-};
 
 /** Constructor of class vector.
@@ -241 +215 @@
  * \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)
+ * \return true -  \a this contains intersection point on return, false - line is parallel to 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;
@@ -253 +226 @@
   Direction.SubtractVector(Origin);
   Direction.Normalize();
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);
-  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
+  //Log() << Verbose(4) << "INFO: Direction is " << Direction << "." << endl;
   factor = Direction.ScalarProduct(PlaneNormal);
-  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
-    DoLog(1) && (Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl);
+  if (factor < MYEPSILON) { // Uniqueness: line parallel to plane?
+    eLog() << Verbose(2) << "Line is parallel to plane, no intersection." << endl;
     return false;
   }
@@ -263 +235 @@
   helper.SubtractVector(Origin);
   factor = helper.ScalarProduct(PlaneNormal)/factor;
-  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
-    DoLog(1) && (Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl);
+  if (factor < MYEPSILON) { // Origin is in-plane
+    //Log() << Verbose(2) << "Origin of line is in-plane, simple." << endl;
     CopyVector(Origin);
     return true;
@@ -271 +243 @@
   Direction.Scale(factor);
   CopyVector(Origin);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);
+  //Log() << Verbose(4) << "INFO: Scaled direction is " << Direction << "." << endl;
   AddVector(&Direction);
 
@@ -278 +250 @@
   helper.SubtractVector(PlaneOffset);
   if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl);
+    //Log() << Verbose(2) << "INFO: Intersection at " << *this << " is good." << endl;
     return true;
   } else {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);
+    eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl;
     return false;
   }
 };
 
-/** Calculates the minimum distance vector of this vector to the plane.
+/** Calculates the minimum distance of this vector to the plane.
  * \param *out output stream for debugging
  * \param *PlaneNormal normal of plane
  * \param *PlaneOffset offset of plane
- * \return distance vector onto to plane
- */
-Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
+ * \return distance to plane
+ */
+double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
 {
   Vector temp;
@@ -310 +282 @@
     sign = 0.;
 
-  temp.Normalize();
-  temp.Scale(sign);
-  return temp;
-};
-
-/** Calculates the minimum distance of this vector to the plane.
- * \sa Vector::GetDistanceVectorToPlane()
- * \param *out output stream for debugging
- * \param *PlaneNormal normal of plane
- * \param *PlaneOffset offset of plane
- * \return distance to plane
- */
-double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
-{
-  return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
+  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 
+ * We construct auxiliary plane with its vector normal to one line direction and the PlaneNormal, then a vector
+ * from the first line's offset onto the plane. Finally, scale by factor is 1/cos(angle(line1,line2..)) = 1/SP(...), and
+ * project onto the first line's direction and add its offset.
  * \param *out output stream for debugging
  * \param *Line1a first vector of first line
@@ -339 +299 @@
 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;
-  //ostream &output = Log() << Verbose(1);
-  //for (int i=0;i<4;i++) {
-  //  for (int j=0;j<4;j++)
-  //    output << "\t" << M->Get(i,j);
-  //  output << endl;
-  //}
-  if (fabs(M->Determinant()) > MYEPSILON) {
-    DoLog(1) && (Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl);
+  bool result = true;
+  Vector Direction, OtherDirection;
+  Vector AuxiliaryNormal;
+  Vector Distance;
+  const Vector *Normal = NULL;
+  Vector *ConstructedNormal = NULL;
+  bool FreeNormal = false;
+
+  // construct both direction vectors
+  Zero();
+  Direction.CopyVector(Line1b);
+  Direction.SubtractVector(Line1a);
+  if (Direction.IsZero())
     return false;
-  }
-  DoLog(1) && (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);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl);
-  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
-   Zero();
-   DoLog(1) && (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);
-      DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
-      return true;
+  OtherDirection.CopyVector(Line2b);
+  OtherDirection.SubtractVector(Line2a);
+  if (OtherDirection.IsZero())
+    return false;
+
+  Direction.Normalize();
+  OtherDirection.Normalize();
+
+  //Log() << Verbose(4) << "INFO: Normalized Direction " << Direction << " and OtherDirection " << OtherDirection << "." << endl;
+
+  if (fabs(OtherDirection.ScalarProduct(&Direction) - 1.) < MYEPSILON) { // lines are parallel
+    if ((Line1a == Line2a) || (Line1a == Line2b))
+      CopyVector(Line1a);
+    else if ((Line1b == Line2b) || (Line1b == Line2b))
+        CopyVector(Line1b);
+    else
+      return false;
+    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
+    return true;
+  } else {
+    // check whether we have a plane normal vector
+    if (PlaneNormal == NULL) {
+      ConstructedNormal = new Vector;
+      ConstructedNormal->MakeNormalVector(&Direction, &OtherDirection);
+      Normal = ConstructedNormal;
+      FreeNormal = true;
+    } else
+      Normal = PlaneNormal;
+
+    AuxiliaryNormal.MakeNormalVector(&OtherDirection, Normal);
+    //Log() << Verbose(4) << "INFO: PlaneNormal is " << *Normal << " and AuxiliaryNormal " << AuxiliaryNormal << "." << endl;
+
+    Distance.CopyVector(Line2a);
+    Distance.SubtractVector(Line1a);
+    //Log() << Verbose(4) << "INFO: Distance is " << Distance << "." << endl;
+    if (Distance.IsZero()) {
+      // offsets are equal, match found
+      CopyVector(Line1a);
+      result = true;
     } else {
-      parallel.CopyVector(Line1a);
-      parallel.SubtractVector(Line2b);
-      factor = parallel.ScalarProduct(&a)/a.Norm();
-      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
-        CopyVector(Line2b);
-        DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
-        return true;
-      }
+      CopyVector(Distance.Projection(&AuxiliaryNormal));
+      //Log() << Verbose(4) << "INFO: Projected Distance is " << *this << "." << endl;
+      double factor = Direction.ScalarProduct(&AuxiliaryNormal);
+      //Log() << Verbose(4) << "INFO: Scaling factor is " << factor << "." << endl;
+      Scale(1./(factor*factor));
+      //Log() << Verbose(4) << "INFO: Scaled Distance is " << *this << "." << endl;
+      CopyVector(Projection(&Direction));
+      //Log() << Verbose(4) << "INFO: Distance, projected into Direction, is " << *this << "." << endl;
+      if (this->IsZero())
+        result = false;
+      else
+        result = true;
+      AddVector(Line1a);
     }
-    DoLog(1) && (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);
-  DoLog(1) && (Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl);
-  if (fabs(temp2.NormSquared()) > MYEPSILON)
-    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
-  else
-    s = 0.;
-  DoLog(1) && (Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl);
-
-  // construct intersection
-  CopyVector(&a);
-  Scale(s);
-  AddVector(Line1a);
-  DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);
-
-  return true;
+
+    if (FreeNormal)
+      delete(ConstructedNormal);
+  }
+  if (result)
+    Log() << Verbose(4) << "INFO: Intersection is " << *this << "." << endl;
+
+  return result;
 };
 
@@ -537 +480 @@
   else
     return false;
-};
-
-/** Checks whether vector is normal to \a *normal.
- * @return true - vector is normalized, false - vector is not
- */
-bool Vector::IsEqualTo(const Vector * const a) const
-{
-  bool status = true;
-  for (int i=0;i<NDIM;i++) {
-    if (fabs(x[i] - a->x[i]) > MYEPSILON)
-      status = false;
-  }
-  return status;
 };
 
@@ -702 +632 @@
 void Vector::Output() const
 {
-  DoLog(0) && (Log() << Verbose(0) << "(");
+  Log() << Verbose(0) << "(";
   for (int i=0;i<NDIM;i++) {
-    DoLog(0) && (Log() << Verbose(0) << x[i]);
+    Log() << Verbose(0) << x[i];
     if (i != 2)
-      DoLog(0) && (Log() << Verbose(0) << ",");
-  }
-  DoLog(0) && (Log() << Verbose(0) << ")");
+      Log() << Verbose(0) << ",";
+  }
+  Log() << Verbose(0) << ")";
 };
 
@@ -817 +747 @@
       x[i] = C.x[i];
   } else {
-    DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);
+    eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl;
   }
 };
@@ -843 +773 @@
   projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
   // withdraw projected vector twice from original one
-  DoLog(1) && (Log() << Verbose(1) << "Vector: ");
+  Log() << Verbose(1) << "Vector: ";
   Output();
-  DoLog(0) && (Log() << Verbose(0) << "\t");
+  Log() << Verbose(0) << "\t";
   for (int i=NDIM;i--;)
     x[i] -= 2.*projection*n->x[i];
-  DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");
+  Log() << Verbose(0) << "Projected vector: ";
   Output();
-  DoLog(0) && (Log() << Verbose(0) << endl);
+  Log() << Verbose(0) << endl;
 };
 
@@ -869 +799 @@
   x2.SubtractVector(y2);
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
+    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
     return false;
   }
@@ -903 +833 @@
   Zero();
   if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
-    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
+    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
     return false;
   }
@@ -954 +884 @@
   double norm;
 
-  DoLog(4) && (Log() << Verbose(4));
+  Log() << Verbose(4);
   GivenVector->Output();
-  DoLog(0) && (Log() << Verbose(0) << endl);
+  Log() << Verbose(0) << endl;
   for (j=NDIM;j--;)
     Components[j] = -1;
@@ -963 +893 @@
     if (fabs(GivenVector->x[j]) > MYEPSILON)
       Components[Last++] = j;
-  DoLog(4) && (Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl);
+  Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
 
   switch(Last) {
@@ -1013 +943 @@
 
   for (j=0;j<num;j++) {
-    DoLog(1) && (Log() << Verbose(1) << j << "th atom's vector: ");
+    Log() << Verbose(1) << j << "th atom's vector: ";
     (vectors[j])->Output();
-    DoLog(0) && (Log() << Verbose(0) << endl);
+    Log() << Verbose(0) << endl;
   }
 
@@ -1135 +1065 @@
     j += i+1;
     do {
-      DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");
+      Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
       cin >> x[i];
     } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
@@ -1166 +1096 @@
   B2 = cos(beta) * x2->Norm() * c;
   C = c * c;
-  DoLog(2) && (Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl);
+  Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
   int flag = 0;
   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
@@ -1205 +1135 @@
   D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
   D3 = y->x[0]/x1->x[0]*A-B1;
-  DoLog(2) && (Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n");
+  Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
   if (fabs(D1) < MYEPSILON) {
-    DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");
+    Log() << Verbose(2) << "D1 == 0!\n";
     if (fabs(D2) > MYEPSILON) {
-      DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");
+      Log() << Verbose(3) << "D2 != 0!\n";
       x[2] = -D3/D2;
       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
       E2 = -x1->x[1]/x1->x[0];
-      DoLog(3) && (Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n");
+      Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
       F1 = E1*E1 + 1.;
       F2 = -E1*E2;
       F3 = E1*E1 + D3*D3/(D2*D2) - C;
-      DoLog(3) && (Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
+      Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
       if (fabs(F1) < MYEPSILON) {
-        DoLog(4) && (Log() << Verbose(4) << "F1 == 0!\n");
-        DoLog(4) && (Log() << Verbose(4) << "Gleichungssystem linear\n");
+        Log() << Verbose(4) << "F1 == 0!\n";
+        Log() << Verbose(4) << "Gleichungssystem linear\n";
         x[1] = F3/(2.*F2);
       } else {
         p = F2/F1;
         q = p*p - F3/F1;
-        DoLog(4) && (Log() << Verbose(4) << "p " << p << "\tq " << q << endl);
+        Log() << Verbose(4) << "p " << p << "\tq " << q << endl;
         if (q < 0) {
-          DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);
+          Log() << Verbose(4) << "q < 0" << endl;
           return false;
         }
@@ -1234 +1164 @@
       x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
     } else {
-      DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");
+      Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n";
       return false;
     }
@@ -1240 +1170 @@
     E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
     E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
-    DoLog(2) && (Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n");
+    Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
     F2 = -(E1*E2 + D2*D3/(D1*D1));
     F3 = E1*E1 + D3*D3/(D1*D1) - C;
-    DoLog(2) && (Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
+    Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     if (fabs(F1) < MYEPSILON) {
-      DoLog(3) && (Log() << Verbose(3) << "F1 == 0!\n");
-      DoLog(3) && (Log() << Verbose(3) << "Gleichungssystem linear\n");
+      Log() << Verbose(3) << "F1 == 0!\n";
+      Log() << Verbose(3) << "Gleichungssystem linear\n";
       x[2] = F3/(2.*F2);
     } else {
       p = F2/F1;
       q = p*p - F3/F1;
-      DoLog(3) && (Log() << Verbose(3) << "p " << p << "\tq " << q << endl);
+      Log() << Verbose(3) << "p " << p << "\tq " << q << endl;
       if (q < 0) {
-        DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);
+        Log() << Verbose(3) << "q < 0" << endl;
         return false;
       }
@@ -1292 +1222 @@
     for (j=2;j>=0;j--) {
       k = (i & pot(2,j)) << j;
-      DoLog(2) && (Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl);
+      Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
       sign[j] = (k == 0) ? 1. : -1.;
     }
-    DoLog(2) && (Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n");
+    Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
     // apply sign matrix
     for (j=NDIM;j--;)
@@ -1301 +1231 @@
     // calculate angle and check
     ang = x2->Angle (this);
-    DoLog(1) && (Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t");
+    Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
     if (fabs(ang - cos(beta)) < MYEPSILON) {
       break;