Changes in src/vector.cpp [a67d19:7b36fe]
- File:
-
- 1 edited
-
src/vector.cpp (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/vector.cpp
ra67d19 r7b36fe 15 15 #include "vector.hpp" 16 16 #include "verbose.hpp" 17 #include "World.hpp"18 17 19 18 #include <gsl/gsl_linalg.h> … … 27 26 */ 28 27 Vector::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 };47 28 48 29 /** Constructor of class vector. … … 253 234 Direction.SubtractVector(Origin); 254 235 Direction.Normalize(); 255 DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);236 Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl; 256 237 //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl; 257 238 factor = Direction.ScalarProduct(PlaneNormal); 258 239 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; 260 241 return false; 261 242 } … … 264 245 factor = helper.ScalarProduct(PlaneNormal)/factor; 265 246 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; 267 248 CopyVector(Origin); 268 249 return true; … … 271 252 Direction.Scale(factor); 272 253 CopyVector(Origin); 273 DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);254 Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl; 274 255 AddVector(&Direction); 275 256 … … 278 259 helper.SubtractVector(PlaneOffset); 279 260 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; 281 262 return true; 282 263 } 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; 284 265 return false; 285 266 } 286 267 }; 287 268 288 /** Calculates the minimum distance vectorof this vector to the plane.269 /** Calculates the minimum distance of this vector to the plane. 289 270 * \param *out output stream for debugging 290 271 * \param *PlaneNormal normal of plane 291 272 * \param *PlaneOffset offset of plane 292 * \return distance vector ontoto plane293 */ 294 Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const273 * \return distance to plane 274 */ 275 double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const 295 276 { 296 277 Vector temp; … … 310 291 sign = 0.; 311 292 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); 327 294 }; 328 295 … … 352 319 353 320 //Log() << Verbose(1) << "Coefficent matrix is:" << endl; 354 //ostream &output = Log() << Verbose(1);355 321 //for (int i=0;i<4;i++) { 356 322 // 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; 359 325 //} 360 326 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; 362 328 return false; 363 329 } 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; 365 331 366 332 … … 378 344 d.CopyVector(Line2b); 379 345 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; 381 347 if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) { 382 348 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; 384 350 return false; 385 351 } … … 394 360 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { 395 361 CopyVector(Line2a); 396 DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);362 Log() << Verbose(1) << "Lines conincide." << endl; 397 363 return true; 398 364 } else { … … 402 368 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { 403 369 CopyVector(Line2b); 404 DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);370 Log() << Verbose(1) << "Lines conincide." << endl; 405 371 return true; 406 372 } 407 373 } 408 DoLog(1) && (Log() << Verbose(1) << "Lines are parallel." << endl);374 Log() << Verbose(1) << "Lines are parallel." << endl; 409 375 Zero(); 410 376 return false; … … 418 384 temp2.CopyVector(&a); 419 385 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; 421 387 if (fabs(temp2.NormSquared()) > MYEPSILON) 422 388 s = temp1.ScalarProduct(&temp2)/temp2.NormSquared(); 423 389 else 424 390 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; 426 392 427 393 // construct intersection … … 429 395 Scale(s); 430 396 AddVector(Line1a); 431 DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);397 Log() << Verbose(1) << "Intersection is at " << *this << "." << endl; 432 398 433 399 return true; … … 702 668 void Vector::Output() const 703 669 { 704 DoLog(0) && (Log() << Verbose(0) << "(");670 Log() << Verbose(0) << "("; 705 671 for (int i=0;i<NDIM;i++) { 706 DoLog(0) && (Log() << Verbose(0) << x[i]);672 Log() << Verbose(0) << x[i]; 707 673 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) << ")"; 711 677 }; 712 678 … … 817 783 x[i] = C.x[i]; 818 784 } 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; 820 786 } 821 787 }; … … 843 809 projection = ScalarProduct(n)/n->ScalarProduct(n); // remove constancy from n (keep as logical one) 844 810 // withdraw projected vector twice from original one 845 DoLog(1) && (Log() << Verbose(1) << "Vector: ");811 Log() << Verbose(1) << "Vector: "; 846 812 Output(); 847 DoLog(0) && (Log() << Verbose(0) << "\t");813 Log() << Verbose(0) << "\t"; 848 814 for (int i=NDIM;i--;) 849 815 x[i] -= 2.*projection*n->x[i]; 850 DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");816 Log() << Verbose(0) << "Projected vector: "; 851 817 Output(); 852 DoLog(0) && (Log() << Verbose(0) << endl);818 Log() << Verbose(0) << endl; 853 819 }; 854 820 … … 869 835 x2.SubtractVector(y2); 870 836 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; 872 838 return false; 873 839 } … … 903 869 Zero(); 904 870 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; 906 872 return false; 907 873 } … … 954 920 double norm; 955 921 956 DoLog(4) && (Log() << Verbose(4));922 Log() << Verbose(4); 957 923 GivenVector->Output(); 958 DoLog(0) && (Log() << Verbose(0) << endl);924 Log() << Verbose(0) << endl; 959 925 for (j=NDIM;j--;) 960 926 Components[j] = -1; … … 963 929 if (fabs(GivenVector->x[j]) > MYEPSILON) 964 930 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; 966 932 967 933 switch(Last) { … … 1013 979 1014 980 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: "; 1016 982 (vectors[j])->Output(); 1017 DoLog(0) && (Log() << Verbose(0) << endl);983 Log() << Verbose(0) << endl; 1018 984 } 1019 985 … … 1135 1101 j += i+1; 1136 1102 do { 1137 DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");1103 Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: "; 1138 1104 cin >> x[i]; 1139 1105 } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check)); … … 1166 1132 B2 = cos(beta) * x2->Norm() * c; 1167 1133 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; 1169 1135 int flag = 0; 1170 1136 if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping … … 1205 1171 D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 1206 1172 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"; 1208 1174 if (fabs(D1) < MYEPSILON) { 1209 DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");1175 Log() << Verbose(2) << "D1 == 0!\n"; 1210 1176 if (fabs(D2) > MYEPSILON) { 1211 DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");1177 Log() << Verbose(3) << "D2 != 0!\n"; 1212 1178 x[2] = -D3/D2; 1213 1179 E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2; 1214 1180 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"; 1216 1182 F1 = E1*E1 + 1.; 1217 1183 F2 = -E1*E2; 1218 1184 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"; 1220 1186 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"; 1223 1189 x[1] = F3/(2.*F2); 1224 1190 } else { 1225 1191 p = F2/F1; 1226 1192 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; 1228 1194 if (q < 0) { 1229 DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);1195 Log() << Verbose(4) << "q < 0" << endl; 1230 1196 return false; 1231 1197 } … … 1234 1200 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2]; 1235 1201 } else { 1236 DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");1202 Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n"; 1237 1203 return false; 1238 1204 } … … 1240 1206 E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1; 1241 1207 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"; 1243 1209 F1 = E2*E2 + D2*D2/(D1*D1) + 1.; 1244 1210 F2 = -(E1*E2 + D2*D3/(D1*D1)); 1245 1211 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"; 1247 1213 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"; 1250 1216 x[2] = F3/(2.*F2); 1251 1217 } else { 1252 1218 p = F2/F1; 1253 1219 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; 1255 1221 if (q < 0) { 1256 DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);1222 Log() << Verbose(3) << "q < 0" << endl; 1257 1223 return false; 1258 1224 } … … 1292 1258 for (j=2;j>=0;j--) { 1293 1259 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; 1295 1261 sign[j] = (k == 0) ? 1. : -1.; 1296 1262 } 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"; 1298 1264 // apply sign matrix 1299 1265 for (j=NDIM;j--;) … … 1301 1267 // calculate angle and check 1302 1268 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"; 1304 1270 if (fabs(ang - cos(beta)) < MYEPSILON) { 1305 1271 break;
Note:
See TracChangeset
for help on using the changeset viewer.
