Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/vector.cpp

    ra67d19 r7b36fe  
    1515#include "vector.hpp"
    1616#include "verbose.hpp"
    17 #include "World.hpp"
    1817
    1918#include <gsl/gsl_linalg.h>
     
    2726 */
    2827Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
    29 
    30 /** Constructor of class vector.
    31  */
    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 };
    4728
    4829/** Constructor of class vector.
     
    253234  Direction.SubtractVector(Origin);
    254235  Direction.Normalize();
    255   DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);
     236  Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
    256237  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
    257238  factor = Direction.ScalarProduct(PlaneNormal);
    258239  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
    259     DoLog(1) && (Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl);
     240    Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
    260241    return false;
    261242  }
     
    264245  factor = helper.ScalarProduct(PlaneNormal)/factor;
    265246  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
    266     DoLog(1) && (Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl);
     247    Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
    267248    CopyVector(Origin);
    268249    return true;
     
    271252  Direction.Scale(factor);
    272253  CopyVector(Origin);
    273   DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);
     254  Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
    274255  AddVector(&Direction);
    275256
     
    278259  helper.SubtractVector(PlaneOffset);
    279260  if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
    280     DoLog(1) && (Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl);
     261    Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
    281262    return true;
    282263  } else {
    283     DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);
     264    eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl;
    284265    return false;
    285266  }
    286267};
    287268
    288 /** Calculates the minimum distance vector of this vector to the plane.
     269/** Calculates the minimum distance of this vector to the plane.
    289270 * \param *out output stream for debugging
    290271 * \param *PlaneNormal normal of plane
    291272 * \param *PlaneOffset offset of plane
    292  * \return distance vector onto to plane
    293  */
    294 Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     273 * \return distance to plane
     274 */
     275double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    295276{
    296277  Vector temp;
     
    310291    sign = 0.;
    311292
    312   temp.Normalize();
    313   temp.Scale(sign);
    314   return temp;
    315 };
    316 
    317 /** Calculates the minimum distance of this vector to the plane.
    318  * \sa Vector::GetDistanceVectorToPlane()
    319  * \param *out output stream for debugging
    320  * \param *PlaneNormal normal of plane
    321  * \param *PlaneOffset offset of plane
    322  * \return distance to plane
    323  */
    324 double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    325 {
    326   return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
     293  return (temp.Norm()*sign);
    327294};
    328295
     
    352319 
    353320  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
    354   //ostream &output = Log() << Verbose(1);
    355321  //for (int i=0;i<4;i++) {
    356322  //  for (int j=0;j<4;j++)
    357   //    output << "\t" << M->Get(i,j);
    358   //  output << endl;
     323  //    cout << "\t" << M->Get(i,j);
     324  //  cout << endl;
    359325  //}
    360326  if (fabs(M->Determinant()) > MYEPSILON) {
    361     DoLog(1) && (Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl);
     327    Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
    362328    return false;
    363329  }
    364   DoLog(1) && (Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl);
     330  Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
    365331
    366332
     
    378344  d.CopyVector(Line2b);
    379345  d.SubtractVector(Line1b);
    380   DoLog(1) && (Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl);
     346  Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
    381347  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
    382348   Zero();
    383    DoLog(1) && (Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl);
     349   Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
    384350   return false;
    385351  }
     
    394360    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    395361      CopyVector(Line2a);
    396       DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
     362      Log() << Verbose(1) << "Lines conincide." << endl;
    397363      return true;
    398364    } else {
     
    402368      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    403369        CopyVector(Line2b);
    404         DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
     370        Log() << Verbose(1) << "Lines conincide." << endl;
    405371        return true;
    406372      }
    407373    }
    408     DoLog(1) && (Log() << Verbose(1) << "Lines are parallel." << endl);
     374    Log() << Verbose(1) << "Lines are parallel." << endl;
    409375    Zero();
    410376    return false;
     
    418384  temp2.CopyVector(&a);
    419385  temp2.VectorProduct(&b);
    420   DoLog(1) && (Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl);
     386  Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
    421387  if (fabs(temp2.NormSquared()) > MYEPSILON)
    422388    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
    423389  else
    424390    s = 0.;
    425   DoLog(1) && (Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl);
     391  Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
    426392
    427393  // construct intersection
     
    429395  Scale(s);
    430396  AddVector(Line1a);
    431   DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);
     397  Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
    432398
    433399  return true;
     
    702668void Vector::Output() const
    703669{
    704   DoLog(0) && (Log() << Verbose(0) << "(");
     670  Log() << Verbose(0) << "(";
    705671  for (int i=0;i<NDIM;i++) {
    706     DoLog(0) && (Log() << Verbose(0) << x[i]);
     672    Log() << Verbose(0) << x[i];
    707673    if (i != 2)
    708       DoLog(0) && (Log() << Verbose(0) << ",");
    709   }
    710   DoLog(0) && (Log() << Verbose(0) << ")");
     674      Log() << Verbose(0) << ",";
     675  }
     676  Log() << Verbose(0) << ")";
    711677};
    712678
     
    817783      x[i] = C.x[i];
    818784  } else {
    819     DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);
     785    eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl;
    820786  }
    821787};
     
    843809  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
    844810  // withdraw projected vector twice from original one
    845   DoLog(1) && (Log() << Verbose(1) << "Vector: ");
     811  Log() << Verbose(1) << "Vector: ";
    846812  Output();
    847   DoLog(0) && (Log() << Verbose(0) << "\t");
     813  Log() << Verbose(0) << "\t";
    848814  for (int i=NDIM;i--;)
    849815    x[i] -= 2.*projection*n->x[i];
    850   DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");
     816  Log() << Verbose(0) << "Projected vector: ";
    851817  Output();
    852   DoLog(0) && (Log() << Verbose(0) << endl);
     818  Log() << Verbose(0) << endl;
    853819};
    854820
     
    869835  x2.SubtractVector(y2);
    870836  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    871     DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
     837    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
    872838    return false;
    873839  }
     
    903869  Zero();
    904870  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    905     DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
     871    eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
    906872    return false;
    907873  }
     
    954920  double norm;
    955921
    956   DoLog(4) && (Log() << Verbose(4));
     922  Log() << Verbose(4);
    957923  GivenVector->Output();
    958   DoLog(0) && (Log() << Verbose(0) << endl);
     924  Log() << Verbose(0) << endl;
    959925  for (j=NDIM;j--;)
    960926    Components[j] = -1;
     
    963929    if (fabs(GivenVector->x[j]) > MYEPSILON)
    964930      Components[Last++] = j;
    965   DoLog(4) && (Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl);
     931  Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
    966932
    967933  switch(Last) {
     
    1013979
    1014980  for (j=0;j<num;j++) {
    1015     DoLog(1) && (Log() << Verbose(1) << j << "th atom's vector: ");
     981    Log() << Verbose(1) << j << "th atom's vector: ";
    1016982    (vectors[j])->Output();
    1017     DoLog(0) && (Log() << Verbose(0) << endl);
     983    Log() << Verbose(0) << endl;
    1018984  }
    1019985
     
    11351101    j += i+1;
    11361102    do {
    1137       DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");
     1103      Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
    11381104      cin >> x[i];
    11391105    } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
     
    11661132  B2 = cos(beta) * x2->Norm() * c;
    11671133  C = c * c;
    1168   DoLog(2) && (Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl);
     1134  Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
    11691135  int flag = 0;
    11701136  if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     
    12051171  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    12061172  D3 = y->x[0]/x1->x[0]*A-B1;
    1207   DoLog(2) && (Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n");
     1173  Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
    12081174  if (fabs(D1) < MYEPSILON) {
    1209     DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");
     1175    Log() << Verbose(2) << "D1 == 0!\n";
    12101176    if (fabs(D2) > MYEPSILON) {
    1211       DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");
     1177      Log() << Verbose(3) << "D2 != 0!\n";
    12121178      x[2] = -D3/D2;
    12131179      E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    12141180      E2 = -x1->x[1]/x1->x[0];
    1215       DoLog(3) && (Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n");
     1181      Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    12161182      F1 = E1*E1 + 1.;
    12171183      F2 = -E1*E2;
    12181184      F3 = E1*E1 + D3*D3/(D2*D2) - C;
    1219       DoLog(3) && (Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
     1185      Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    12201186      if (fabs(F1) < MYEPSILON) {
    1221         DoLog(4) && (Log() << Verbose(4) << "F1 == 0!\n");
    1222         DoLog(4) && (Log() << Verbose(4) << "Gleichungssystem linear\n");
     1187        Log() << Verbose(4) << "F1 == 0!\n";
     1188        Log() << Verbose(4) << "Gleichungssystem linear\n";
    12231189        x[1] = F3/(2.*F2);
    12241190      } else {
    12251191        p = F2/F1;
    12261192        q = p*p - F3/F1;
    1227         DoLog(4) && (Log() << Verbose(4) << "p " << p << "\tq " << q << endl);
     1193        Log() << Verbose(4) << "p " << p << "\tq " << q << endl;
    12281194        if (q < 0) {
    1229           DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);
     1195          Log() << Verbose(4) << "q < 0" << endl;
    12301196          return false;
    12311197        }
     
    12341200      x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    12351201    } else {
    1236       DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");
     1202      Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n";
    12371203      return false;
    12381204    }
     
    12401206    E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    12411207    E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    1242     DoLog(2) && (Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n");
     1208    Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    12431209    F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    12441210    F2 = -(E1*E2 + D2*D3/(D1*D1));
    12451211    F3 = E1*E1 + D3*D3/(D1*D1) - C;
    1246     DoLog(2) && (Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
     1212    Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    12471213    if (fabs(F1) < MYEPSILON) {
    1248       DoLog(3) && (Log() << Verbose(3) << "F1 == 0!\n");
    1249       DoLog(3) && (Log() << Verbose(3) << "Gleichungssystem linear\n");
     1214      Log() << Verbose(3) << "F1 == 0!\n";
     1215      Log() << Verbose(3) << "Gleichungssystem linear\n";
    12501216      x[2] = F3/(2.*F2);
    12511217    } else {
    12521218      p = F2/F1;
    12531219      q = p*p - F3/F1;
    1254       DoLog(3) && (Log() << Verbose(3) << "p " << p << "\tq " << q << endl);
     1220      Log() << Verbose(3) << "p " << p << "\tq " << q << endl;
    12551221      if (q < 0) {
    1256         DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);
     1222        Log() << Verbose(3) << "q < 0" << endl;
    12571223        return false;
    12581224      }
     
    12921258    for (j=2;j>=0;j--) {
    12931259      k = (i & pot(2,j)) << j;
    1294       DoLog(2) && (Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl);
     1260      Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
    12951261      sign[j] = (k == 0) ? 1. : -1.;
    12961262    }
    1297     DoLog(2) && (Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n");
     1263    Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
    12981264    // apply sign matrix
    12991265    for (j=NDIM;j--;)
     
    13011267    // calculate angle and check
    13021268    ang = x2->Angle (this);
    1303     DoLog(1) && (Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t");
     1269    Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
    13041270    if (fabs(ang - cos(beta)) < MYEPSILON) {
    13051271      break;
Note: See TracChangeset for help on using the changeset viewer.