Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/vector.cpp

    rb84d5d ra67d19  
    1515#include "vector.hpp"
    1616#include "verbose.hpp"
     17#include "World.hpp"
    1718
    1819#include <gsl/gsl_linalg.h>
    … …  
    2627 */
    2728Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
     29
     30/** Constructor of class vector.
     31 */
     32Vector::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 */
     41Vector::Vector(const Vector &a)
     42{
     43  x[0] = a.x[0];
     44  x[1] = a.x[1];
     45  x[2] = a.x[2];
     46};
    2847
    2948/** Constructor of class vector.
    … …  
    234253  Direction.SubtractVector(Origin);
    235254  Direction.Normalize();
    236   Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl;
     255  DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);
    237256  //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl;
    238257  factor = Direction.ScalarProduct(PlaneNormal);
    239258  if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane?
    240     Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl;
     259    DoLog(1) && (Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl);
    241260    return false;
    242261  }
    … …  
    245264  factor = helper.ScalarProduct(PlaneNormal)/factor;
    246265  if (fabs(factor) < MYEPSILON) { // Origin is in-plane
    247     Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl;
     266    DoLog(1) && (Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl);
    248267    CopyVector(Origin);
    249268    return true;
    … …  
    252271  Direction.Scale(factor);
    253272  CopyVector(Origin);
    254   Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl;
     273  DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);
    255274  AddVector(&Direction);
    256275
    … …  
    259278  helper.SubtractVector(PlaneOffset);
    260279  if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) {
    261     Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl;
     280    DoLog(1) && (Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl);
    262281    return true;
    263282  } else {
    264     eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl;
     283    DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);
    265284    return false;
    266285  }
    267286};
    268287
    269 /** Calculates the minimum distance of this vector to the plane.
     288/** Calculates the minimum distance vector of this vector to the plane.
    270289 * \param *out output stream for debugging
    271290 * \param *PlaneNormal normal of plane
    272291 * \param *PlaneOffset offset of plane
    273  * \return distance to plane
    274  */
    275 double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     292 * \return distance vector onto to plane
     293 */
     294Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    276295{
    277296  Vector temp;
    … …  
    291310    sign = 0.;
    292311
    293   return (temp.Norm()*sign);
     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 */
     324double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     325{
     326  return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
    294327};
    295328
    … …  
    319352 
    320353  //Log() << Verbose(1) << "Coefficent matrix is:" << endl;
     354  //ostream &output = Log() << Verbose(1);
    321355  //for (int i=0;i<4;i++) {
    322356  //  for (int j=0;j<4;j++)
    323   //    cout << "\t" << M->Get(i,j);
    324   //  cout << endl;
     357  //    output << "\t" << M->Get(i,j);
     358  //  output << endl;
    325359  //}
    326360  if (fabs(M->Determinant()) > MYEPSILON) {
    327     Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl;
     361    DoLog(1) && (Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl);
    328362    return false;
    329363  }
    330   delete(M);
    331   Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl;
     364  DoLog(1) && (Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl);
    332365
    333366
    … …  
    345378  d.CopyVector(Line2b);
    346379  d.SubtractVector(Line1b);
    347   Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl;
     380  DoLog(1) && (Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl);
    348381  if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) {
    349382   Zero();
    350    Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl;
     383   DoLog(1) && (Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl);
    351384   return false;
    352385  }
    … …  
    361394    if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    362395      CopyVector(Line2a);
    363       Log() << Verbose(1) << "Lines conincide." << endl;
     396      DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
    364397      return true;
    365398    } else {
    … …  
    369402      if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) {
    370403        CopyVector(Line2b);
    371         Log() << Verbose(1) << "Lines conincide." << endl;
     404        DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);
    372405        return true;
    373406      }
    374407    }
    375     Log() << Verbose(1) << "Lines are parallel." << endl;
     408    DoLog(1) && (Log() << Verbose(1) << "Lines are parallel." << endl);
    376409    Zero();
    377410    return false;
    … …  
    385418  temp2.CopyVector(&a);
    386419  temp2.VectorProduct(&b);
    387   Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl;
     420  DoLog(1) && (Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl);
    388421  if (fabs(temp2.NormSquared()) > MYEPSILON)
    389422    s = temp1.ScalarProduct(&temp2)/temp2.NormSquared();
    390423  else
    391424    s = 0.;
    392   Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl;
     425  DoLog(1) && (Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl);
    393426
    394427  // construct intersection
    … …  
    396429  Scale(s);
    397430  AddVector(Line1a);
    398   Log() << Verbose(1) << "Intersection is at " << *this << "." << endl;
     431  DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);
    399432
    400433  return true;
    … …  
    584617 * \return lhs + a
    585618 */
    586 const Vector& operator+=(Vector& a, const Vector& b)
     619Vector& operator+=(Vector& a, const Vector& b)
    587620{
    588621  a.AddVector(&b);
    … …  
    595628 * \return lhs - a
    596629 */
    597 const Vector& operator-=(Vector& a, const Vector& b)
     630Vector& operator-=(Vector& a, const Vector& b)
    598631{
    599632  a.SubtractVector(&b);
    … …  
    606639 * \return lhs.x[i] * m
    607640 */
    608 const Vector& operator*=(Vector& a, const double m)
     641Vector& operator*=(Vector& a, const double m)
    609642{
    610643  a.Scale(m);
    … …  
    617650 * \return a + b
    618651 */
    619 Vector const operator+(const Vector& a, const Vector& b)
    620 {
    621   Vector x(a);
    622   x.AddVector(&b);
    623   return x;
     652Vector& operator+(const Vector& a, const Vector& b)
     653{
     654  Vector *x = new Vector;
     655  x->CopyVector(&a);
     656  x->AddVector(&b);
     657  return *x;
    624658};
    625659
    … …  
    629663 * \return a - b
    630664 */
    631 Vector const operator-(const Vector& a, const Vector& b)
    632 {
    633   Vector x(a);
    634   x.SubtractVector(&b);
    635   return x;
     665Vector& operator-(const Vector& a, const Vector& b)
     666{
     667  Vector *x = new Vector;
     668  x->CopyVector(&a);
     669  x->SubtractVector(&b);
     670  return *x;
    636671};
    637672
    … …  
    641676 * \return m * a
    642677 */
    643 Vector const operator*(const Vector& a, const double m)
    644 {
    645   Vector x(a);
    646   x.Scale(m);
    647   return x;
     678Vector& operator*(const Vector& a, const double m)
     679{
     680  Vector *x = new Vector;
     681  x->CopyVector(&a);
     682  x->Scale(m);
     683  return *x;
    648684};
    649685
    … …  
    653689 * \return m * a
    654690 */
    655 Vector const operator*(const double m, const Vector& a )
    656 {
    657   Vector x(a);
    658   x.Scale(m);
    659   return x;
     691Vector& operator*(const double m, const Vector& a )
     692{
     693  Vector *x = new Vector;
     694  x->CopyVector(&a);
     695  x->Scale(m);
     696  return *x;
    660697};
    661698
    … …  
    665702void Vector::Output() const
    666703{
    667   Log() << Verbose(0) << "(";
     704  DoLog(0) && (Log() << Verbose(0) << "(");
    668705  for (int i=0;i<NDIM;i++) {
    669     Log() << Verbose(0) << x[i];
     706    DoLog(0) && (Log() << Verbose(0) << x[i]);
    670707    if (i != 2)
    671       Log() << Verbose(0) << ",";
    672   }
    673   Log() << Verbose(0) << ")";
     708      DoLog(0) && (Log() << Verbose(0) << ",");
     709  }
     710  DoLog(0) && (Log() << Verbose(0) << ")");
    674711};
    675712
    … …  
    780817      x[i] = C.x[i];
    781818  } else {
    782     eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl;
     819    DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);
    783820  }
    784821};
    … …  
    806843  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
    807844  // withdraw projected vector twice from original one
    808   Log() << Verbose(1) << "Vector: ";
     845  DoLog(1) && (Log() << Verbose(1) << "Vector: ");
    809846  Output();
    810   Log() << Verbose(0) << "\t";
     847  DoLog(0) && (Log() << Verbose(0) << "\t");
    811848  for (int i=NDIM;i--;)
    812849    x[i] -= 2.*projection*n->x[i];
    813   Log() << Verbose(0) << "Projected vector: ";
     850  DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");
    814851  Output();
    815   Log() << Verbose(0) << endl;
     852  DoLog(0) && (Log() << Verbose(0) << endl);
    816853};
    817854
    … …  
    832869  x2.SubtractVector(y2);
    833870  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    834     eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
     871    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
    835872    return false;
    836873  }
    … …  
    866903  Zero();
    867904  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    868     eLog() << Verbose(2) << "Given vectors are linear dependent." << endl;
     905    DoeLog(2) && (eLog()<< Verbose(2) << "Given vectors are linear dependent." << endl);
    869906    return false;
    870907  }
    … …  
    917954  double norm;
    918955
    919   Log() << Verbose(4);
     956  DoLog(4) && (Log() << Verbose(4));
    920957  GivenVector->Output();
    921   Log() << Verbose(0) << endl;
     958  DoLog(0) && (Log() << Verbose(0) << endl);
    922959  for (j=NDIM;j--;)
    923960    Components[j] = -1;
    … …  
    926963    if (fabs(GivenVector->x[j]) > MYEPSILON)
    927964      Components[Last++] = j;
    928   Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
     965  DoLog(4) && (Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl);
    929966
    930967  switch(Last) {
    … …  
    9761013
    9771014  for (j=0;j<num;j++) {
    978     Log() << Verbose(1) << j << "th atom's vector: ";
     1015    DoLog(1) && (Log() << Verbose(1) << j << "th atom's vector: ");
    9791016    (vectors[j])->Output();
    980     Log() << Verbose(0) << endl;
     1017    DoLog(0) && (Log() << Verbose(0) << endl);
    9811018  }
    9821019
    … …  
    10731110void Vector::CopyVector(const Vector * const y)
    10741111{
    1075   // check for self assignment
    1076   if(y!=this){
    1077     for (int i=NDIM;i--;)
    1078       this->x[i] = y->x[i];
    1079   }
     1112  for (int i=NDIM;i--;)
     1113    this->x[i] = y->x[i];
    10801114}
    10811115
    … …  
    10851119void Vector::CopyVector(const Vector &y)
    10861120{
    1087   // check for self assignment
    1088   if(&y!=this) {
    1089     for (int i=NDIM;i--;)
    1090       this->x[i] = y.x[i];
    1091   }
     1121  for (int i=NDIM;i--;)
     1122    this->x[i] = y.x[i];
    10921123}
    10931124
    … …  
    11041135    j += i+1;
    11051136    do {
    1106       Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
     1137      DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");
    11071138      cin >> x[i];
    11081139    } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
    … …  
    11351166  B2 = cos(beta) * x2->Norm() * c;
    11361167  C = c * c;
    1137   Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
     1168  DoLog(2) && (Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl);
    11381169  int flag = 0;
    11391170  if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    … …  
    11741205  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    11751206  D3 = y->x[0]/x1->x[0]*A-B1;
    1176   Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
     1207  DoLog(2) && (Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n");
    11771208  if (fabs(D1) < MYEPSILON) {
    1178     Log() << Verbose(2) << "D1 == 0!\n";
     1209    DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");
    11791210    if (fabs(D2) > MYEPSILON) {
    1180       Log() << Verbose(3) << "D2 != 0!\n";
     1211      DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");
    11811212      x[2] = -D3/D2;
    11821213      E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    11831214      E2 = -x1->x[1]/x1->x[0];
    1184       Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     1215      DoLog(3) && (Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n");
    11851216      F1 = E1*E1 + 1.;
    11861217      F2 = -E1*E2;
    11871218      F3 = E1*E1 + D3*D3/(D2*D2) - C;
    1188       Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     1219      DoLog(3) && (Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
    11891220      if (fabs(F1) < MYEPSILON) {
    1190         Log() << Verbose(4) << "F1 == 0!\n";
    1191         Log() << Verbose(4) << "Gleichungssystem linear\n";
     1221        DoLog(4) && (Log() << Verbose(4) << "F1 == 0!\n");
     1222        DoLog(4) && (Log() << Verbose(4) << "Gleichungssystem linear\n");
    11921223        x[1] = F3/(2.*F2);
    11931224      } else {
    11941225        p = F2/F1;
    11951226        q = p*p - F3/F1;
    1196         Log() << Verbose(4) << "p " << p << "\tq " << q << endl;
     1227        DoLog(4) && (Log() << Verbose(4) << "p " << p << "\tq " << q << endl);
    11971228        if (q < 0) {
    1198           Log() << Verbose(4) << "q < 0" << endl;
     1229          DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);
    11991230          return false;
    12001231        }
    … …  
    12031234      x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    12041235    } else {
    1205       Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n";
     1236      DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");
    12061237      return false;
    12071238    }
    … …  
    12091240    E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    12101241    E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    1211     Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     1242    DoLog(2) && (Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n");
    12121243    F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    12131244    F2 = -(E1*E2 + D2*D3/(D1*D1));
    12141245    F3 = E1*E1 + D3*D3/(D1*D1) - C;
    1215     Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     1246    DoLog(2) && (Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");
    12161247    if (fabs(F1) < MYEPSILON) {
    1217       Log() << Verbose(3) << "F1 == 0!\n";
    1218       Log() << Verbose(3) << "Gleichungssystem linear\n";
     1248      DoLog(3) && (Log() << Verbose(3) << "F1 == 0!\n");
     1249      DoLog(3) && (Log() << Verbose(3) << "Gleichungssystem linear\n");
    12191250      x[2] = F3/(2.*F2);
    12201251    } else {
    12211252      p = F2/F1;
    12221253      q = p*p - F3/F1;
    1223       Log() << Verbose(3) << "p " << p << "\tq " << q << endl;
     1254      DoLog(3) && (Log() << Verbose(3) << "p " << p << "\tq " << q << endl);
    12241255      if (q < 0) {
    1225         Log() << Verbose(3) << "q < 0" << endl;
     1256        DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);
    12261257        return false;
    12271258      }
    … …  
    12611292    for (j=2;j>=0;j--) {
    12621293      k = (i & pot(2,j)) << j;
    1263       Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
     1294      DoLog(2) && (Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl);
    12641295      sign[j] = (k == 0) ? 1. : -1.;
    12651296    }
    1266     Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
     1297    DoLog(2) && (Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n");
    12671298    // apply sign matrix
    12681299    for (j=NDIM;j--;)
    … …  
    12701301    // calculate angle and check
    12711302    ang = x2->Angle (this);
    1272     Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
     1303    DoLog(1) && (Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t");
    12731304    if (fabs(ang - cos(beta)) < MYEPSILON) {
    12741305      break;
Note: See TracChangeset for help on using the changeset viewer.