Changes in src/vector.cpp [a67d19:b34306]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/vector.cpp
ra67d19 rb34306 253 253 Direction.SubtractVector(Origin); 254 254 Direction.Normalize(); 255 DoLog(1) && (Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl);255 Log() << Verbose(1) << "INFO: Direction is " << Direction << "." << endl; 256 256 //Log() << Verbose(1) << "INFO: PlaneNormal is " << *PlaneNormal << " and PlaneOffset is " << *PlaneOffset << "." << endl; 257 257 factor = Direction.ScalarProduct(PlaneNormal); 258 258 if (fabs(factor) < MYEPSILON) { // Uniqueness: line parallel to plane? 259 DoLog(1) && (Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl);259 Log() << Verbose(1) << "BAD: Line is parallel to plane, no intersection." << endl; 260 260 return false; 261 261 } … … 264 264 factor = helper.ScalarProduct(PlaneNormal)/factor; 265 265 if (fabs(factor) < MYEPSILON) { // Origin is in-plane 266 DoLog(1) && (Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl);266 Log() << Verbose(1) << "GOOD: Origin of line is in-plane." << endl; 267 267 CopyVector(Origin); 268 268 return true; … … 271 271 Direction.Scale(factor); 272 272 CopyVector(Origin); 273 DoLog(1) && (Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl);273 Log() << Verbose(1) << "INFO: Scaled direction is " << Direction << "." << endl; 274 274 AddVector(&Direction); 275 275 … … 278 278 helper.SubtractVector(PlaneOffset); 279 279 if (helper.ScalarProduct(PlaneNormal) < MYEPSILON) { 280 DoLog(1) && (Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl);280 Log() << Verbose(1) << "GOOD: Intersection is " << *this << "." << endl; 281 281 return true; 282 282 } else { 283 DoeLog(2) && (eLog()<< Verbose(2) << "Intersection point " << *this << " is not on plane." << endl);283 eLog() << Verbose(2) << "Intersection point " << *this << " is not on plane." << endl; 284 284 return false; 285 285 } … … 352 352 353 353 //Log() << Verbose(1) << "Coefficent matrix is:" << endl; 354 //ostream &output = Log() << Verbose(1);355 354 //for (int i=0;i<4;i++) { 356 355 // for (int j=0;j<4;j++) 357 // output << "\t" << M->Get(i,j);358 // output << endl;356 // cout << "\t" << M->Get(i,j); 357 // cout << endl; 359 358 //} 360 359 if (fabs(M->Determinant()) > MYEPSILON) { 361 DoLog(1) && (Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl);360 Log() << Verbose(1) << "Determinant of coefficient matrix is NOT zero." << endl; 362 361 return false; 363 362 } 364 DoLog(1) && (Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl);363 Log() << Verbose(1) << "INFO: Line1a = " << *Line1a << ", Line1b = " << *Line1b << ", Line2a = " << *Line2a << ", Line2b = " << *Line2b << "." << endl; 365 364 366 365 … … 378 377 d.CopyVector(Line2b); 379 378 d.SubtractVector(Line1b); 380 DoLog(1) && (Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl);379 Log() << Verbose(1) << "INFO: a = " << a << ", b = " << b << ", c = " << c << "." << endl; 381 380 if ((a.NormSquared() < MYEPSILON) || (b.NormSquared() < MYEPSILON)) { 382 381 Zero(); 383 DoLog(1) && (Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl);382 Log() << Verbose(1) << "At least one of the lines is ill-defined, i.e. offset equals second vector." << endl; 384 383 return false; 385 384 } … … 394 393 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { 395 394 CopyVector(Line2a); 396 DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);395 Log() << Verbose(1) << "Lines conincide." << endl; 397 396 return true; 398 397 } else { … … 402 401 if ((factor >= -MYEPSILON) && (factor - 1. < MYEPSILON)) { 403 402 CopyVector(Line2b); 404 DoLog(1) && (Log() << Verbose(1) << "Lines conincide." << endl);403 Log() << Verbose(1) << "Lines conincide." << endl; 405 404 return true; 406 405 } 407 406 } 408 DoLog(1) && (Log() << Verbose(1) << "Lines are parallel." << endl);407 Log() << Verbose(1) << "Lines are parallel." << endl; 409 408 Zero(); 410 409 return false; … … 418 417 temp2.CopyVector(&a); 419 418 temp2.VectorProduct(&b); 420 DoLog(1) && (Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl);419 Log() << Verbose(1) << "INFO: temp1 = " << temp1 << ", temp2 = " << temp2 << "." << endl; 421 420 if (fabs(temp2.NormSquared()) > MYEPSILON) 422 421 s = temp1.ScalarProduct(&temp2)/temp2.NormSquared(); 423 422 else 424 423 s = 0.; 425 DoLog(1) && (Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl);424 Log() << Verbose(1) << "Factor s is " << temp1.ScalarProduct(&temp2) << "/" << temp2.NormSquared() << " = " << s << "." << endl; 426 425 427 426 // construct intersection … … 429 428 Scale(s); 430 429 AddVector(Line1a); 431 DoLog(1) && (Log() << Verbose(1) << "Intersection is at " << *this << "." << endl);430 Log() << Verbose(1) << "Intersection is at " << *this << "." << endl; 432 431 433 432 return true; … … 702 701 void Vector::Output() const 703 702 { 704 DoLog(0) && (Log() << Verbose(0) << "(");703 Log() << Verbose(0) << "("; 705 704 for (int i=0;i<NDIM;i++) { 706 DoLog(0) && (Log() << Verbose(0) << x[i]);705 Log() << Verbose(0) << x[i]; 707 706 if (i != 2) 708 DoLog(0) && (Log() << Verbose(0) << ",");709 } 710 DoLog(0) && (Log() << Verbose(0) << ")");707 Log() << Verbose(0) << ","; 708 } 709 Log() << Verbose(0) << ")"; 711 710 }; 712 711 … … 817 816 x[i] = C.x[i]; 818 817 } else { 819 DoeLog(1) && (eLog()<< Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl);818 eLog() << Verbose(1) << "inverse of matrix does not exists: det A = " << detA << "." << endl; 820 819 } 821 820 }; … … 843 842 projection = ScalarProduct(n)/n->ScalarProduct(n); // remove constancy from n (keep as logical one) 844 843 // withdraw projected vector twice from original one 845 DoLog(1) && (Log() << Verbose(1) << "Vector: ");844 Log() << Verbose(1) << "Vector: "; 846 845 Output(); 847 DoLog(0) && (Log() << Verbose(0) << "\t");846 Log() << Verbose(0) << "\t"; 848 847 for (int i=NDIM;i--;) 849 848 x[i] -= 2.*projection*n->x[i]; 850 DoLog(0) && (Log() << Verbose(0) << "Projected vector: ");849 Log() << Verbose(0) << "Projected vector: "; 851 850 Output(); 852 DoLog(0) && (Log() << Verbose(0) << endl);851 Log() << Verbose(0) << endl; 853 852 }; 854 853 … … 869 868 x2.SubtractVector(y2); 870 869 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);870 eLog() << Verbose(2) << "Given vectors are linear dependent." << endl; 872 871 return false; 873 872 } … … 903 902 Zero(); 904 903 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);904 eLog() << Verbose(2) << "Given vectors are linear dependent." << endl; 906 905 return false; 907 906 } … … 954 953 double norm; 955 954 956 DoLog(4) && (Log() << Verbose(4));955 Log() << Verbose(4); 957 956 GivenVector->Output(); 958 DoLog(0) && (Log() << Verbose(0) << endl);957 Log() << Verbose(0) << endl; 959 958 for (j=NDIM;j--;) 960 959 Components[j] = -1; … … 963 962 if (fabs(GivenVector->x[j]) > MYEPSILON) 964 963 Components[Last++] = j; 965 DoLog(4) && (Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl);964 Log() << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl; 966 965 967 966 switch(Last) { … … 1013 1012 1014 1013 for (j=0;j<num;j++) { 1015 DoLog(1) && (Log() << Verbose(1) << j << "th atom's vector: ");1014 Log() << Verbose(1) << j << "th atom's vector: "; 1016 1015 (vectors[j])->Output(); 1017 DoLog(0) && (Log() << Verbose(0) << endl);1016 Log() << Verbose(0) << endl; 1018 1017 } 1019 1018 … … 1135 1134 j += i+1; 1136 1135 do { 1137 DoLog(0) && (Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ");1136 Log() << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: "; 1138 1137 cin >> x[i]; 1139 1138 } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check)); … … 1166 1165 B2 = cos(beta) * x2->Norm() * c; 1167 1166 C = c * c; 1168 DoLog(2) && (Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl);1167 Log() << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl; 1169 1168 int flag = 0; 1170 1169 if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping … … 1205 1204 D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2]; 1206 1205 D3 = y->x[0]/x1->x[0]*A-B1; 1207 DoLog(2) && (Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n");1206 Log() << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n"; 1208 1207 if (fabs(D1) < MYEPSILON) { 1209 DoLog(2) && (Log() << Verbose(2) << "D1 == 0!\n");1208 Log() << Verbose(2) << "D1 == 0!\n"; 1210 1209 if (fabs(D2) > MYEPSILON) { 1211 DoLog(3) && (Log() << Verbose(3) << "D2 != 0!\n");1210 Log() << Verbose(3) << "D2 != 0!\n"; 1212 1211 x[2] = -D3/D2; 1213 1212 E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2; 1214 1213 E2 = -x1->x[1]/x1->x[0]; 1215 DoLog(3) && (Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n");1214 Log() << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n"; 1216 1215 F1 = E1*E1 + 1.; 1217 1216 F2 = -E1*E2; 1218 1217 F3 = E1*E1 + D3*D3/(D2*D2) - C; 1219 DoLog(3) && (Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");1218 Log() << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n"; 1220 1219 if (fabs(F1) < MYEPSILON) { 1221 DoLog(4) && (Log() << Verbose(4) << "F1 == 0!\n");1222 DoLog(4) && (Log() << Verbose(4) << "Gleichungssystem linear\n");1220 Log() << Verbose(4) << "F1 == 0!\n"; 1221 Log() << Verbose(4) << "Gleichungssystem linear\n"; 1223 1222 x[1] = F3/(2.*F2); 1224 1223 } else { 1225 1224 p = F2/F1; 1226 1225 q = p*p - F3/F1; 1227 DoLog(4) && (Log() << Verbose(4) << "p " << p << "\tq " << q << endl);1226 Log() << Verbose(4) << "p " << p << "\tq " << q << endl; 1228 1227 if (q < 0) { 1229 DoLog(4) && (Log() << Verbose(4) << "q < 0" << endl);1228 Log() << Verbose(4) << "q < 0" << endl; 1230 1229 return false; 1231 1230 } … … 1234 1233 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2]; 1235 1234 } else { 1236 DoLog(2) && (Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n");1235 Log() << Verbose(2) << "Gleichungssystem unterbestimmt\n"; 1237 1236 return false; 1238 1237 } … … 1240 1239 E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1; 1241 1240 E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2]; 1242 DoLog(2) && (Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n");1241 Log() << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n"; 1243 1242 F1 = E2*E2 + D2*D2/(D1*D1) + 1.; 1244 1243 F2 = -(E1*E2 + D2*D3/(D1*D1)); 1245 1244 F3 = E1*E1 + D3*D3/(D1*D1) - C; 1246 DoLog(2) && (Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n");1245 Log() << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n"; 1247 1246 if (fabs(F1) < MYEPSILON) { 1248 DoLog(3) && (Log() << Verbose(3) << "F1 == 0!\n");1249 DoLog(3) && (Log() << Verbose(3) << "Gleichungssystem linear\n");1247 Log() << Verbose(3) << "F1 == 0!\n"; 1248 Log() << Verbose(3) << "Gleichungssystem linear\n"; 1250 1249 x[2] = F3/(2.*F2); 1251 1250 } else { 1252 1251 p = F2/F1; 1253 1252 q = p*p - F3/F1; 1254 DoLog(3) && (Log() << Verbose(3) << "p " << p << "\tq " << q << endl);1253 Log() << Verbose(3) << "p " << p << "\tq " << q << endl; 1255 1254 if (q < 0) { 1256 DoLog(3) && (Log() << Verbose(3) << "q < 0" << endl);1255 Log() << Verbose(3) << "q < 0" << endl; 1257 1256 return false; 1258 1257 } … … 1292 1291 for (j=2;j>=0;j--) { 1293 1292 k = (i & pot(2,j)) << j; 1294 DoLog(2) && (Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl);1293 Log() << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl; 1295 1294 sign[j] = (k == 0) ? 1. : -1.; 1296 1295 } 1297 DoLog(2) && (Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n");1296 Log() << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n"; 1298 1297 // apply sign matrix 1299 1298 for (j=NDIM;j--;) … … 1301 1300 // calculate angle and check 1302 1301 ang = x2->Angle (this); 1303 DoLog(1) && (Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t");1302 Log() << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t"; 1304 1303 if (fabs(ang - cos(beta)) < MYEPSILON) { 1305 1304 break;
Note:
See TracChangeset
for help on using the changeset viewer.