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