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

File:

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

src/vector.cpp

@@ -1 +1 @@
 /** \file vector.cpp
- *
+ * 
  * Function implementations for the class vector.
- *
- */
-
+ * 
+ */
+ 
 #include "molecules.hpp"
-
+ 
 
 /************************************ Functions for class vector ************************************/
@@ -22 +22 @@
 Vector::~Vector() {};
 
-/** Calculates square of distance between this and another vector.
+/** Calculates 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 Vector::Distance(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 |\f$
- */
-double Vector::Distance(const Vector *y) const
-{
-  double res = 0.;
-  for (int i=NDIM;i--;)
-    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
-  return (sqrt(res));
+  return (res);  
 };
 
@@ -81 +69 @@
         if (tmp < res) res = tmp;
       }
-  return (res);
+  return (res);  
 };
 
@@ -124 +112 @@
   for (int i=NDIM;i--;)
     res += x[i]*y->x[i];
-  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);
-
-};
-
+  return (res);  
+};
 
 /** projects this vector onto plane defined by \a *y.
- * \param *y normal vector of plane
+ * \param *y array to normal vector of plane
  * \return \f$\langle x, y \rangle\f$
  */
@@ -153 +123 @@
   Vector tmp;
   tmp.CopyVector(y);
-  tmp.Normalize();
-  tmp.Scale(ScalarProduct(&tmp));
+  tmp.Scale(Projection(y));
   this->SubtractVector(&tmp);
 };
@@ -175 +144 @@
   for (int i=NDIM;i--;)
     res += this->x[i]*this->x[i];
-  return (sqrt(res));
+  return (sqrt(res));  
 };
 
@@ -209 +178 @@
  */
 void Vector::Init(double x1, double x2, double x3)
-{
+{ 
   x[0] = x1;
   x[1] = x2;
@@ -219 +188 @@
  * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
  */
-double Vector::Angle(const Vector *y) const
-{
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
+double Vector::Angle(Vector *y) const
+{
+  return acos(this->ScalarProduct(y)/Norm()/y->Norm()); 
 };
 
@@ -269 +238 @@
 
 /** Sums two vectors \a  and \b component-wise.
- * \param a first vector
+ * \param a first vector 
  * \param b second vector
  * \return a + b
@@ -282 +251 @@
 
 /** Factors given vector \a a times \a m.
- * \param a vector
+ * \param a vector 
  * \param m factor
  * \return a + b
@@ -313 +282 @@
 };
 
-/** 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 << ")";
+ofstream& operator<<(ofstream& ost,Vector& m)
+{
+        m.Output(&ost);
         return ost;
 };
@@ -351 +309 @@
 };
 
-/** Translate atom by given vector.
+/** Translate atom by given vector. 
  * \param trans[] translation vector.
  */
@@ -358 +316 @@
   for (int i=NDIM;i--;)
     x[i] += trans->x[i];
-};
+}; 
 
 /** Do a matrix multiplication.
@@ -397 +355 @@
     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];
@@ -421 +379 @@
 {
   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.
  */
@@ -440 +398 @@
   Output((ofstream *)&cout);
   cout << endl;
-};
+}; 
 
 /** Calculates normal vector for three given vectors (being three points in space).
@@ -472 +430 @@
   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
   Normalize();
-
+  
   return true;
 };
@@ -530 +488 @@
 /** 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)
@@ -551 +509 @@
       Components[Last++] = j;
   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
-
+        
   switch(Last) {
     case 3:  // threecomponent system
@@ -564 +522 @@
     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;
@@ -583 +541 @@
 {
 //  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);
@@ -595 +553 @@
  * \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: ";
@@ -607 +565 @@
   int np = 3;
         struct LSQ_params par;
-
+    
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -613 +571 @@
    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++)
@@ -664 +622 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
-
+ 
   for (i=(size_t)np;i--;)
     this->x[i] = gsl_vector_get(s->x, i);
@@ -730 +688 @@
  * \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)
@@ -748 +706 @@
   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;
@@ -781 +739 @@
   // 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;
@@ -797 +755 @@
       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;
@@ -824 +782 @@
     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;
@@ -874 +832 @@
     }
     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];
@@ -880 +838 @@
     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;
     }