Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/vector.cpp

    • Property mode changed from 100755 to 100644
    rc4d4df r7ea9e6  
    2121/** Constructor of class vector.
    2222 */
    23 Vector::Vector(double x1, double x2, double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
     23Vector::Vector(const double x1, const double x2, const double x3) { x[0] = x1; x[1] = x2; x[2] = x3; };
    2424
    2525/** Desctructor of class vector.
    … …  
    3131 * \return \f$| x - y |^2\f$
    3232 */
    33 double Vector::DistanceSquared(const Vector *y) const
     33double Vector::DistanceSquared(const Vector * const y) const
    3434{
    3535  double res = 0.;
    … …  
    4343 * \return \f$| x - y |\f$
    4444 */
    45 double Vector::Distance(const Vector *y) const
     45double Vector::Distance(const Vector * const y) const
    4646{
    4747  double res = 0.;
    … …  
    5656 * \return \f$| x - y |\f$
    5757 */
    58 double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
     58double Vector::PeriodicDistance(const Vector * const y, const double * const cell_size) const
    5959{
    6060  double res = Distance(y), tmp, matrix[NDIM*NDIM];
    … …  
    9494 * \return \f$| x - y |^2\f$
    9595 */
    96 double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
     96double Vector::PeriodicDistanceSquared(const Vector * const y, const double * const cell_size) const
    9797{
    9898  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    … …  
    131131 * Tries to translate a vector into each adjacent neighbouring cell.
    132132 */
    133 void Vector::KeepPeriodic(ofstream *out, double *matrix)
     133void Vector::KeepPeriodic(ofstream *out, const double * const matrix)
    134134{
    135135//  int N[NDIM];
    … …  
    162162 * \return \f$\langle x, y \rangle\f$
    163163 */
    164 double Vector::ScalarProduct(const Vector *y) const
     164double Vector::ScalarProduct(const Vector * const y) const
    165165{
    166166  double res = 0.;
    … …  
    177177 *  \return \f$ x \times y \f&
    178178 */
    179 void Vector::VectorProduct(const Vector *y)
     179void Vector::VectorProduct(const Vector * const y)
    180180{
    181181  Vector tmp;
    … …  
    184184  tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
    185185  this->CopyVector(&tmp);
    186 
    187186};
    188187
    … …  
    192191 * \return \f$\langle x, y \rangle\f$
    193192 */
    194 void Vector::ProjectOntoPlane(const Vector *y)
     193void Vector::ProjectOntoPlane(const Vector * const y)
    195194{
    196195  Vector tmp;
    … …  
    217216 * \return true -  \a this contains intersection point on return, false - line is parallel to plane
    218217 */
    219 bool Vector::GetIntersectionWithPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset, Vector *Origin, Vector *LineVector)
     218bool Vector::GetIntersectionWithPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector)
    220219{
    221220  double factor;
    … …  
    264263 * \return distance to plane
    265264 */
    266 double Vector::DistanceToPlane(ofstream *out, Vector *PlaneNormal, Vector *PlaneOffset)
     265double Vector::DistanceToPlane(ofstream *out, const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    267266{
    268267  Vector temp;
    … …  
    276275  temp.AddVector(this);
    277276  temp.SubtractVector(PlaneOffset);
    278 
    279   return temp.Norm();
     277  double sign = temp.ScalarProduct(PlaneNormal);
     278  if (fabs(sign) > MYEPSILON)
     279    sign /= fabs(sign);
     280  else
     281    sign = 0.;
     282
     283  return (temp.Norm()*sign);
    280284};
    281285
    … …  
    292296 * \return true - \a this will contain the intersection on return, false - lines are parallel
    293297 */
    294 bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, Vector *Line1a, Vector *Line1b, Vector *Line2a, Vector *Line2b, const Vector *PlaneNormal)
     298bool Vector::GetIntersectionOfTwoLinesOnPlane(ofstream *out, const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *PlaneNormal)
    295299{
    296300  bool result = true;
    … …  
    375379 * \param *y array to second vector
    376380 */
    377 void Vector::ProjectIt(const Vector *y)
     381void Vector::ProjectIt(const Vector * const y)
    378382{
    379383  Vector helper(*y);
    … …  
    386390 * \return Vector
    387391 */
    388 Vector Vector::Projection(const Vector *y) const
     392Vector Vector::Projection(const Vector * const y) const
    389393{
    390394  Vector helper(*y);
    … …  
    435439/** Zeros all components of this vector.
    436440 */
    437 void Vector::One(double one)
     441void Vector::One(const double one)
    438442{
    439443  for (int i=NDIM;i--;)
    … …  
    443447/** Initialises all components of this vector.
    444448 */
    445 void Vector::Init(double x1, double x2, double x3)
     449void Vector::Init(const double x1, const double x2, const double x3)
    446450{
    447451  x[0] = x1;
    … …  
    469473 * @return true - vector is normalized, false - vector is not
    470474 */
    471 bool Vector::IsNormalTo(const Vector *normal) const
     475bool Vector::IsNormalTo(const Vector * const normal) const
    472476{
    473477  if (ScalarProduct(normal) < MYEPSILON)
    … …  
    481485 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
    482486 */
    483 double Vector::Angle(const Vector *y) const
     487double Vector::Angle(const Vector * const y) const
    484488{
    485489  double norm1 = Norm(), norm2 = y->Norm();
    … …  
    500504 * \param alpha rotation angle in radian
    501505 */
    502 void Vector::RotateVector(const Vector *axis, const double alpha)
     506void Vector::RotateVector(const Vector * const axis, const double alpha)
    503507{
    504508  Vector a,y;
    … …  
    656660 * \param *factor pointer to scaling factor
    657661 */
    658 void Vector::Scale(double **factor)
     662void Vector::Scale(const double ** const factor)
    659663{
    660664  for (int i=NDIM;i--;)
    … …  
    662666};
    663667
    664 void Vector::Scale(double *factor)
     668void Vector::Scale(const double * const factor)
    665669{
    666670  for (int i=NDIM;i--;)
    … …  
    668672};
    669673
    670 void Vector::Scale(double factor)
     674void Vector::Scale(const double factor)
    671675{
    672676  for (int i=NDIM;i--;)
    … …  
    677681 * \param trans[] translation vector.
    678682 */
    679 void Vector::Translate(const Vector *trans)
     683void Vector::Translate(const Vector * const trans)
    680684{
    681685  for (int i=NDIM;i--;)
    … …  
    687691 * \param *Minv inverse matrix
    688692 */
    689 void Vector::WrapPeriodically(const double *M, const double *Minv)
     693void Vector::WrapPeriodically(const double * const M, const double * const Minv)
    690694{
    691695  MatrixMultiplication(Minv);
    … …  
    704708 * \param *matrix NDIM_NDIM array
    705709 */
    706 void Vector::MatrixMultiplication(const double *M)
     710void Vector::MatrixMultiplication(const double * const M)
    707711{
    708712  Vector C;
    … …  
    716720};
    717721
    718 /** Calculate the inverse of a 3x3 matrix.
     722/** Do a matrix multiplication with the \a *A' inverse.
    719723 * \param *matrix NDIM_NDIM array
    720724 */
    721 double * Vector::InverseMatrix(double *A)
    722 {
    723   double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
     725void Vector::InverseMatrixMultiplication(const double * const A)
     726{
     727  Vector C;
     728  double B[NDIM*NDIM];
    724729  double detA = RDET3(A);
    725730  double detAReci;
    726731
    727   for (int i=0;i<NDIM*NDIM;++i)
    728     B[i] = 0.;
    729732  // calculate the inverse B
    730733  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    … …  
    739742    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    740743    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
    741   }
    742   return B;
    743 };
    744 
    745 /** Do a matrix multiplication with the \a *A' inverse.
    746  * \param *matrix NDIM_NDIM array
    747  */
    748 void Vector::InverseMatrixMultiplication(const double *A)
    749 {
    750   Vector C;
    751   double B[NDIM*NDIM];
    752   double detA = RDET3(A);
    753   double detAReci;
    754 
    755   // calculate the inverse B
    756   if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    757     detAReci = 1./detA;
    758     B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    759     B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    760     B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    761     B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    762     B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    763     B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    764     B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    765     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    766     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
    767744
    768745    // do the matrix multiplication
    … …  
    786763 * \param *factors three-component vector with the factor for each given vector
    787764 */
    788 void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
     765void Vector::LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors)
    789766{
    790767  for(int i=NDIM;i--;)
    … …  
    795772 * \param n[] normal vector of mirror plane.
    796773 */
    797 void Vector::Mirror(const Vector *n)
     774void Vector::Mirror(const Vector * const n)
    798775{
    799776  double projection;
    … …  
    817794 * \return true - success, vectors are linear independent, false - failure due to linear dependency
    818795 */
    819 bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
     796bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2, const Vector * const y3)
    820797{
    821798  Vector x1, x2;
    … …  
    853830 * \return true - success, vectors are linear independent, false - failure due to linear dependency
    854831 */
    855 bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
     832bool Vector::MakeNormalVector(const Vector * const y1, const Vector * const y2)
    856833{
    857834  Vector x1,x2;
    … …  
    884861 * \return true - success, false - vector is zero
    885862 */
    886 bool Vector::MakeNormalVector(const Vector *y1)
     863bool Vector::MakeNormalVector(const Vector * const y1)
    887864{
    888865  bool result = false;
    … …  
    904881 * \return true - success, false - failure (null vector given)
    905882 */
    906 bool Vector::GetOneNormalVector(const Vector *GivenVector)
     883bool Vector::GetOneNormalVector(const Vector * const GivenVector)
    907884{
    908885  int Components[NDIM]; // contains indices of non-zero components
    … …  
    950927 * \return scaling parameter for this vector
    951928 */
    952 double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
     929double Vector::CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const
    953930{
    954931//  cout << Verbose(3) << "For comparison: ";
    … …  
    965942 * \return true if success, false if failed due to linear dependency
    966943 */
    967 bool Vector::LSQdistance(Vector **vectors, int num)
     944bool Vector::LSQdistance(const Vector **vectors, int num)
    968945{
    969946  int j;
    … …  
    10471024 * \param *y vector
    10481025 */
    1049 void Vector::AddVector(const Vector *y)
     1026void Vector::AddVector(const Vector * const y)
    10501027{
    10511028  for (int i=NDIM;i--;)
    … …  
    10561033 * \param *y vector
    10571034 */
    1058 void Vector::SubtractVector(const Vector *y)
     1035void Vector::SubtractVector(const Vector * const y)
    10591036{
    10601037  for (int i=NDIM;i--;)
    … …  
    10651042 * \param *y vector
    10661043 */
    1067 void Vector::CopyVector(const Vector *y)
     1044void Vector::CopyVector(const Vector * const y)
    10681045{
    10691046  for (int i=NDIM;i--;)
    … …  
    10741051 * \param y vector
    10751052 */
    1076 void Vector::CopyVector(const Vector y)
     1053void Vector::CopyVector(const Vector &y)
    10771054{
    10781055  for (int i=NDIM;i--;)
    … …  
    10851062 * \param check whether bounds shall be checked (true) or not (false)
    10861063 */
    1087 void Vector::AskPosition(double *cell_size, bool check)
     1064void Vector::AskPosition(const double * const cell_size, const bool check)
    10881065{
    10891066  char coords[3] = {'x','y','z'};
    … …  
    11131090 * \bug this is not yet working properly
    11141091 */
    1115 bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
     1092bool Vector::SolveSystem(Vector * x1, Vector * x2, Vector * y, const double alpha, const double beta, const double c)
    11161093{
    11171094  double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    … …  
    11391116      break;
    11401117    case 2:
    1141       flip(&x1->x[0],&x1->x[1]);
    1142       flip(&x2->x[0],&x2->x[1]);
    1143       flip(&y->x[0],&y->x[1]);
    1144       //flip(&x[0],&x[1]);
    1145       flip(&x1->x[1],&x1->x[2]);
    1146       flip(&x2->x[1],&x2->x[2]);
    1147       flip(&y->x[1],&y->x[2]);
    1148       //flip(&x[1],&x[2]);
     1118      flip(x1->x[0],x1->x[1]);
     1119      flip(x2->x[0],x2->x[1]);
     1120      flip(y->x[0],y->x[1]);
     1121      //flip(x[0],x[1]);
     1122      flip(x1->x[1],x1->x[2]);
     1123      flip(x2->x[1],x2->x[2]);
     1124      flip(y->x[1],y->x[2]);
     1125      //flip(x[1],x[2]);
    11491126    case 1:
    1150       flip(&x1->x[0],&x1->x[1]);
    1151       flip(&x2->x[0],&x2->x[1]);
    1152       flip(&y->x[0],&y->x[1]);
    1153       //flip(&x[0],&x[1]);
    1154       flip(&x1->x[1],&x1->x[2]);
    1155       flip(&x2->x[1],&x2->x[2]);
    1156       flip(&y->x[1],&y->x[2]);
    1157       //flip(&x[1],&x[2]);
     1127      flip(x1->x[0],x1->x[1]);
     1128      flip(x2->x[0],x2->x[1]);
     1129      flip(y->x[0],y->x[1]);
     1130      //flip(x[0],x[1]);
     1131      flip(x1->x[1],x1->x[2]);
     1132      flip(x2->x[1],x2->x[2]);
     1133      flip(y->x[1],y->x[2]);
     1134      //flip(x[1],x[2]);
    11581135      break;
    11591136  }
    … …  
    12241201      break;
    12251202    case 2:
    1226       flip(&x1->x[0],&x1->x[1]);
    1227       flip(&x2->x[0],&x2->x[1]);
    1228       flip(&y->x[0],&y->x[1]);
    1229       flip(&x[0],&x[1]);
    1230       flip(&x1->x[1],&x1->x[2]);
    1231       flip(&x2->x[1],&x2->x[2]);
    1232       flip(&y->x[1],&y->x[2]);
    1233       flip(&x[1],&x[2]);
     1203      flip(x1->x[0],x1->x[1]);
     1204      flip(x2->x[0],x2->x[1]);
     1205      flip(y->x[0],y->x[1]);
     1206      flip(x[0],x[1]);
     1207      flip(x1->x[1],x1->x[2]);
     1208      flip(x2->x[1],x2->x[2]);
     1209      flip(y->x[1],y->x[2]);
     1210      flip(x[1],x[2]);
    12341211    case 1:
    1235       flip(&x1->x[0],&x1->x[1]);
    1236       flip(&x2->x[0],&x2->x[1]);
    1237       flip(&y->x[0],&y->x[1]);
    1238       //flip(&x[0],&x[1]);
    1239       flip(&x1->x[1],&x1->x[2]);
    1240       flip(&x2->x[1],&x2->x[2]);
    1241       flip(&y->x[1],&y->x[2]);
    1242       flip(&x[1],&x[2]);
     1212      flip(x1->x[0],x1->x[1]);
     1213      flip(x2->x[0],x2->x[1]);
     1214      flip(y->x[0],y->x[1]);
     1215      //flip(x[0],x[1]);
     1216      flip(x1->x[1],x1->x[2]);
     1217      flip(x2->x[1],x2->x[2]);
     1218      flip(y->x[1],y->x[2]);
     1219      flip(x[1],x[2]);
    12431220      break;
    12441221  }
    … …  
    12761253 * @param three vectors forming the matrix that defines the shape of the parallelpiped
    12771254 */
    1278 bool Vector::IsInParallelepiped(Vector offset, double *parallelepiped)
     1255bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
    12791256{
    12801257  Vector a;
Note: See TracChangeset for help on using the changeset viewer.