Changes in src/vector.cpp [2985c8:cc2ee5]

File:

• src/vector.cpp (modified) (33 diffs, 1 prop)

src/vector.cpp

@@ -1 +1 @@
 /** \file vector.cpp
- * 
+ *
  * Function implementations for the class vector.
- * 
- */
- 
+ *
+ */
+
 #include "molecules.hpp"
- 
+
 
 /************************************ Functions for class vector ************************************/
@@ -21 +21 @@
  */
 Vector::~Vector() {};
+
+/** Calculates square of distance between this and another vector.
+ * \param *y array to second vector
+ * \return \f$| x - y |^2\f$
+ */
+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]);
+  return (res);
+};
 
 /** Calculates distance between this and another vector.
  * \param *y array to second vector
- * \return \f$| x - y |^2\f$
+ * \return \f$| x - y |\f$
  */
 double Vector::Distance(const Vector *y) const
@@ -31 +43 @@
   for (int i=NDIM;i--;)
     res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (res);  
+  return (sqrt(res));
 };
 
@@ -69 +81 @@
         if (tmp < res) res = tmp;
       }
-  return (res);  
+  return (res);
 };
 
@@ -112 +124 @@
   for (int i=NDIM;i--;)
     res += x[i]*y->x[i];
-  return (res);  
-};
+  return (res);
+};
+
+
+/** Calculates VectorProduct between this and another vector.
+ *  -# returns the Product in place of vector from which it was initiated
+ *  -# ATTENTION: Only three dim.
+ *  \param *y array to vector with which to calculate crossproduct
+ *  \return \f$ x \times y \f&
+ */
+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);
+
+};
+
 
 /** projects this vector onto plane defined by \a *y.
- * \param *y array to normal vector of plane
+ * \param *y normal vector of plane
  * \return \f$\langle x, y \rangle\f$
  */
@@ -123 +153 @@
   Vector tmp;
   tmp.CopyVector(y);
-  tmp.Scale(Projection(y));
+  tmp.Normalize();
+  tmp.Scale(ScalarProduct(&tmp));
   this->SubtractVector(&tmp);
 };
@@ -144 +175 @@
   for (int i=NDIM;i--;)
     res += this->x[i]*this->x[i];
-  return (sqrt(res));  
+  return (sqrt(res));
 };
 
@@ -178 +209 @@
  */
 void Vector::Init(double x1, double x2, double x3)
-{ 
+{
   x[0] = x1;
   x[1] = x2;
@@ -188 +219 @@
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
-double Vector::Angle(Vector *y) const
-{
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
+double Vector::Angle(const Vector *y) const
+{
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
 };
 
@@ -238 +269 @@
 
 /** Sums two vectors \a  and \b component-wise.
- * \param a first vector 
+ * \param a first vector
  * \param b second vector
  * \return a + b
@@ -251 +282 @@
 
 /** Factors given vector \a a times \a m.
- * \param a vector 
+ * \param a vector
  * \param m factor
  * \return a + b
@@ -282 +313 @@
 };
 
-ofstream& operator<<(ofstream& ost,Vector& m)
-{
-        m.Output(&ost);
+/** Prints a 3dim vector to a stream.
+ * \param ost output stream
+ * \param v Vector to be printed
+ * \return output stream
+ */
+ostream& operator<<(ostream& ost,Vector& m)
+{
+  ost << "(";
+  for (int i=0;i<NDIM;i++) {
+    ost << m.x[i];
+    if (i != 2)
+      ost << ",";
+  }
+  ost << ")";
         return ost;
 };
@@ -309 +351 @@
 };
 
-/** Translate atom by given vector. 
+/** Translate atom by given vector.
  * \param trans[] translation vector.
  */
@@ -316 +358 @@
   for (int i=NDIM;i--;)
     x[i] += trans->x[i];
-}; 
+};
 
 /** Do a matrix multiplication.
@@ -355 +397 @@
     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-  
+
     // do the matrix multiplication
     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
@@ -379 +421 @@
 {
   for(int i=NDIM;i--;)
-    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i]; 
-};
-
-/** Mirrors atom against a given plane. 
+    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
+};
+
+/** Mirrors atom against a given plane.
  * \param n[] normal vector of mirror plane.
  */
@@ -398 +440 @@
   Output((ofstream *)&cout);
   cout << endl;
-}; 
+};
 
 /** Calculates normal vector for three given vectors (being three points in space).
@@ -430 +472 @@
   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
   Normalize();
-  
+
   return true;
 };
@@ -488 +530 @@
 /** Creates this vector as one of the possible orthonormal ones to the given one.
  * Just scan how many components of given *vector are unequal to zero and
- * try to get the skp of both to be zero accordingly.  
+ * try to get the skp of both to be zero accordingly.
  * \param *vector given vector
  * \return true - success, false - failure (null vector given)
@@ -509 +551 @@
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
-        
+
   switch(Last) {
     case 3:  // threecomponent system
@@ -522 +564 @@
     case 1: // one component system
       // set sole non-zero component to 0, and one of the other zero component pendants to 1
-      x[(Components[0]+2)%NDIM] = 0.;   
-      x[(Components[0]+1)%NDIM] = 1.;   
-      x[Components[0]] = 0.;   
+      x[(Components[0]+2)%NDIM] = 0.;
+      x[(Components[0]+1)%NDIM] = 1.;
+      x[Components[0]] = 0.;
       return true;
       break;
@@ -541 +583 @@
 {
 //  cout << Verbose(3) << "For comparison: ";
-//  cout << "A " << A->Projection(this) << "\t"; 
-//  cout << "B " << B->Projection(this) << "\t"; 
-//  cout << "C " << C->Projection(this) << "\t"; 
+//  cout << "A " << A->Projection(this) << "\t";
+//  cout << "B " << B->Projection(this) << "\t";
+//  cout << "C " << C->Projection(this) << "\t";
 //  cout << endl;
   return A->Projection(this);
@@ -553 +595 @@
  * \return true if success, false if failed due to linear dependency
  */
-bool Vector::LSQdistance(Vector **vectors, int num) 
+bool Vector::LSQdistance(Vector **vectors, int num)
 {
         int j;
-                                
+
         for (j=0;j<num;j++) {
                 cout << Verbose(1) << j << "th atom's vector: ";
@@ -565 +607 @@
   int np = 3;
         struct LSQ_params par;
-    
+
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -571 +613 @@
    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.); 
-         
+                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++)
@@ -622 +664 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
- 
+
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
@@ -688 +730 @@
  * \param alpha first angle
  * \param beta second angle
- * \param c norm of final vector 
+ * \param c norm of final vector
  * \return a vector with \f$\langle x1,x2 \rangle=A\f$, \f$\langle x1,y \rangle = B\f$ and with norm \a c.
- * \bug this is not yet working properly 
+ * \bug this is not yet working properly
  */
 bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
@@ -706 +748 @@
   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     if (fabs(x1->x[1]) > MYEPSILON) {
-      flag = 1;    
+      flag = 1;
     } else if (fabs(x1->x[2]) > MYEPSILON) {
        flag = 2;
@@ -739 +781 @@
   // now comes the case system
   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
-  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 
+  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
   D3 = y->x[0]/x1->x[0]*A-B1;
   cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
   if (fabs(D1) < MYEPSILON) {
-    cout << Verbose(2) << "D1 == 0!\n"; 
+    cout << Verbose(2) << "D1 == 0!\n";
     if (fabs(D2) > MYEPSILON) {
-      cout << Verbose(3) << "D2 != 0!\n"; 
+      cout << Verbose(3) << "D2 != 0!\n";
       x[2] = -D3/D2;
       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
@@ -755 +797 @@
       cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
       if (fabs(F1) < MYEPSILON) {
-        cout << Verbose(4) << "F1 == 0!\n"; 
+        cout << Verbose(4) << "F1 == 0!\n";
         cout << Verbose(4) << "Gleichungssystem linear\n";
-        x[1] = F3/(2.*F2); 
+        x[1] = F3/(2.*F2);
       } else {
         p = F2/F1;
         q = p*p - F3/F1;
-        cout << Verbose(4) << "p " << p << "\tq " << q << endl;  
+        cout << Verbose(4) << "p " << p << "\tq " << q << endl;
         if (q < 0) {
           cout << Verbose(4) << "q < 0" << endl;
@@ -782 +824 @@
     cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     if (fabs(F1) < MYEPSILON) {
-      cout << Verbose(3) << "F1 == 0!\n"; 
+      cout << Verbose(3) << "F1 == 0!\n";
       cout << Verbose(3) << "Gleichungssystem linear\n";
-      x[2] = F3/(2.*F2);     
+      x[2] = F3/(2.*F2);
     } else {
       p = F2/F1;
       q = p*p - F3/F1;
-      cout << Verbose(3) << "p " << p << "\tq " << q << endl;  
+      cout << Verbose(3) << "p " << p << "\tq " << q << endl;
       if (q < 0) {
         cout << Verbose(3) << "q < 0" << endl;
@@ -832 +874 @@
     }
     cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
-    // apply sign matrix 
+    // apply sign matrix
     for (j=NDIM;j--;)
       x[j] *= sign[j];
@@ -838 +880 @@
     ang = x2->Angle (this);
     cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
-    if (fabs(ang - cos(beta)) < MYEPSILON) {  
+    if (fabs(ang - cos(beta)) < MYEPSILON) {
       break;
     }