| 6 |
|
|
| 7 |
|
RigidBody::RigidBody() : StuntDouble() { |
| 8 |
|
objType = OT_RIGIDBODY; |
| 9 |
< |
com_good = false; |
| 10 |
< |
precalc_done = false; |
| 9 |
> |
is_linear = false; |
| 10 |
> |
linear_axis = -1; |
| 11 |
> |
momIntTol = 1e-6; |
| 12 |
|
} |
| 13 |
|
|
| 14 |
|
RigidBody::~RigidBody() { |
| 98 |
|
trq[i] = 0.0; |
| 99 |
|
} |
| 100 |
|
|
| 100 |
– |
forces_good = false; |
| 101 |
– |
|
| 101 |
|
} |
| 102 |
|
|
| 103 |
< |
void RigidBody::setEulerAngles( double phi, double theta, double psi ){ |
| 103 |
> |
void RigidBody::setEuler( double phi, double theta, double psi ){ |
| 104 |
|
|
| 105 |
|
A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
| 106 |
|
A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
| 190 |
|
|
| 191 |
|
for (int i = 0; i < 3; i++) |
| 192 |
|
for (int j = 0; j < 3; j++) |
| 193 |
< |
the_A[i][j] = the_A[i][j]; |
| 193 |
> |
the_A[i][j] = A[i][j]; |
| 194 |
|
|
| 195 |
|
} |
| 196 |
|
|
| 228 |
|
|
| 229 |
|
void RigidBody::getI( double the_I[3][3] ){ |
| 230 |
|
|
| 232 |
– |
if (precalc_done) { |
| 233 |
– |
|
| 231 |
|
for (int i = 0; i < 3; i++) |
| 232 |
|
for (int j = 0; j < 3; j++) |
| 233 |
|
the_I[i][j] = I[i][j]; |
| 234 |
|
|
| 238 |
– |
} else { |
| 239 |
– |
|
| 240 |
– |
} |
| 235 |
|
} |
| 236 |
|
|
| 237 |
|
void RigidBody::lab2Body( double r[3] ){ |
| 264 |
|
|
| 265 |
|
void RigidBody::calcRefCoords( ) { |
| 266 |
|
|
| 267 |
< |
int i,j,k; |
| 267 |
> |
int i,j,k, it, n_linear_coords; |
| 268 |
|
double mtmp; |
| 269 |
|
vec3 apos; |
| 270 |
|
double refCOM[3]; |
| 271 |
+ |
vec3 ptmp; |
| 272 |
+ |
double Itmp[3][3]; |
| 273 |
+ |
double evals[3]; |
| 274 |
+ |
double evects[3][3]; |
| 275 |
+ |
double r, r2, len; |
| 276 |
|
|
| 277 |
+ |
// First, find the center of mass: |
| 278 |
+ |
|
| 279 |
|
mass = 0.0; |
| 280 |
|
for (j=0; j<3; j++) |
| 281 |
|
refCOM[j] = 0.0; |
| 294 |
|
for(j = 0; j < 3; j++) |
| 295 |
|
refCOM[j] /= mass; |
| 296 |
|
|
| 297 |
+ |
// Next, move the origin of the reference coordinate system to the COM: |
| 298 |
+ |
|
| 299 |
|
for (i = 0; i < myAtoms.size(); i++) { |
| 300 |
|
apos = refCoords[i]; |
| 301 |
|
for (j=0; j < 3; j++) { |
| 304 |
|
refCoords[i] = apos; |
| 305 |
|
} |
| 306 |
|
|
| 307 |
+ |
// Moment of Inertia calculation |
| 308 |
+ |
|
| 309 |
+ |
for (i = 0; i < 3; i++) |
| 310 |
+ |
for (j = 0; j < 3; j++) |
| 311 |
+ |
Itmp[i][j] = 0.0; |
| 312 |
+ |
|
| 313 |
+ |
for (it = 0; it < myAtoms.size(); it++) { |
| 314 |
+ |
|
| 315 |
+ |
mtmp = myAtoms[it]->getMass(); |
| 316 |
+ |
ptmp = refCoords[it]; |
| 317 |
+ |
r= norm3(ptmp.vec); |
| 318 |
+ |
r2 = r*r; |
| 319 |
+ |
|
| 320 |
+ |
for (i = 0; i < 3; i++) { |
| 321 |
+ |
for (j = 0; j < 3; j++) { |
| 322 |
+ |
|
| 323 |
+ |
if (i==j) Itmp[i][j] += mtmp * r2; |
| 324 |
+ |
|
| 325 |
+ |
Itmp[i][j] -= mtmp * ptmp.vec[i]*ptmp.vec[j]; |
| 326 |
+ |
} |
| 327 |
+ |
} |
| 328 |
+ |
} |
| 329 |
+ |
|
| 330 |
+ |
diagonalize3x3(Itmp, evals, sU); |
| 331 |
+ |
|
| 332 |
+ |
// zero out I and then fill the diagonals with the moments of inertia: |
| 333 |
+ |
|
| 334 |
+ |
n_linear_coords = 0; |
| 335 |
+ |
|
| 336 |
+ |
for (i = 0; i < 3; i++) { |
| 337 |
+ |
for (j = 0; j < 3; j++) { |
| 338 |
+ |
I[i][j] = 0.0; |
| 339 |
+ |
} |
| 340 |
+ |
I[i][i] = evals[i]; |
| 341 |
+ |
|
| 342 |
+ |
if (fabs(evals[i]) < momIntTol) { |
| 343 |
+ |
is_linear = true; |
| 344 |
+ |
n_linear_coords++; |
| 345 |
+ |
linear_axis = i; |
| 346 |
+ |
} |
| 347 |
+ |
} |
| 348 |
+ |
|
| 349 |
+ |
if (n_linear_coords > 1) { |
| 350 |
+ |
sprintf( painCave.errMsg, |
| 351 |
+ |
"RigidBody error.\n" |
| 352 |
+ |
"\tOOPSE found more than one axis in this rigid body with a vanishing \n" |
| 353 |
+ |
"\tmoment of inertia. This can happen in one of three ways:\n" |
| 354 |
+ |
"\t 1) Only one atom was specified, or \n" |
| 355 |
+ |
"\t 2) All atoms were specified at the same location, or\n" |
| 356 |
+ |
"\t 3) The programmers did something stupid.\n" |
| 357 |
+ |
"\tIt is silly to use a rigid body to describe this situation. Be smarter.\n" |
| 358 |
+ |
); |
| 359 |
+ |
painCave.isFatal = 1; |
| 360 |
+ |
simError(); |
| 361 |
+ |
} |
| 362 |
+ |
|
| 363 |
+ |
// renormalize column vectors: |
| 364 |
+ |
|
| 365 |
+ |
for (i=0; i < 3; i++) { |
| 366 |
+ |
len = 0.0; |
| 367 |
+ |
for (j = 0; j < 3; j++) { |
| 368 |
+ |
len += sU[i][j]*sU[i][j]; |
| 369 |
+ |
} |
| 370 |
+ |
len = sqrt(len); |
| 371 |
+ |
for (j = 0; j < 3; j++) { |
| 372 |
+ |
sU[i][j] /= len; |
| 373 |
+ |
} |
| 374 |
+ |
} |
| 375 |
|
} |
| 376 |
|
|
| 377 |
|
void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){ |
| 437 |
|
// (Actually, on second thought, don't. Integrator does this now.) |
| 438 |
|
// lab2Body(trq); |
| 439 |
|
|
| 369 |
– |
forces_good = true; |
| 370 |
– |
|
| 440 |
|
} |
| 441 |
|
|
| 442 |
|
void RigidBody::updateAtoms() { |
| 625 |
|
vel[j] /= mass; |
| 626 |
|
} |
| 627 |
|
|
| 559 |
– |
com_good = true; |
| 628 |
|
} |
| 561 |
– |
|
| 562 |
– |
void RigidBody::findOrient() { |
| 563 |
– |
|
| 564 |
– |
size_t it; |
| 565 |
– |
int i, j; |
| 566 |
– |
double ptmp[3]; |
| 567 |
– |
double Itmp[3][3]; |
| 568 |
– |
double evals[3]; |
| 569 |
– |
double evects[3][3]; |
| 570 |
– |
double r2, mtmp, len; |
| 629 |
|
|
| 630 |
< |
if (!com_good) findCOM(); |
| 630 |
> |
void RigidBody::accept(BaseVisitor* v){ |
| 631 |
> |
vector<Atom*>::iterator atomIter; |
| 632 |
> |
v->visit(this); |
| 633 |
|
|
| 634 |
< |
// Calculate inertial tensor matrix elements: |
| 635 |
< |
|
| 576 |
< |
for (i = 0; i < 3; i++) |
| 577 |
< |
for (j = 0; j < 3; j++) |
| 578 |
< |
Itmp[i][j] = 0.0; |
| 579 |
< |
|
| 580 |
< |
for (it = 0; it < myAtoms.size(); it++) { |
| 581 |
< |
|
| 582 |
< |
mtmp = myAtoms[it]->getMass(); |
| 583 |
< |
myAtoms[it]->getPos(ptmp); |
| 584 |
< |
|
| 585 |
< |
for (j = 0; j < 3; j++) |
| 586 |
< |
ptmp[j] = pos[j] - ptmp[j]; |
| 587 |
< |
|
| 588 |
< |
r2 = norm3(ptmp); |
| 589 |
< |
|
| 590 |
< |
for (i = 0; i < 3; i++) { |
| 591 |
< |
for (j = 0; j < 3; j++) { |
| 592 |
< |
|
| 593 |
< |
if (i==j) Itmp[i][j] = mtmp * r2; |
| 594 |
< |
|
| 595 |
< |
Itmp[i][j] -= mtmp * ptmp[i]*ptmp[j]; |
| 596 |
< |
} |
| 597 |
< |
} |
| 598 |
< |
} |
| 599 |
< |
|
| 600 |
< |
diagonalize3x3(Itmp, evals, sU); |
| 601 |
< |
|
| 602 |
< |
// zero out I and then fill the diagonals with the moments of inertia: |
| 603 |
< |
|
| 604 |
< |
for (i = 0; i < 3; i++) { |
| 605 |
< |
for (j = 0; j < 3; j++) { |
| 606 |
< |
I[i][j] = 0.0; |
| 607 |
< |
} |
| 608 |
< |
I[i][i] = evals[i]; |
| 609 |
< |
} |
| 610 |
< |
|
| 611 |
< |
// renormalize column vectors: |
| 612 |
< |
|
| 613 |
< |
for (i=0; i < 3; i++) { |
| 614 |
< |
len = 0.0; |
| 615 |
< |
for (j = 0; j < 3; j++) { |
| 616 |
< |
len += sU[i][j]*sU[i][j]; |
| 617 |
< |
} |
| 618 |
< |
len = sqrt(len); |
| 619 |
< |
for (j = 0; j < 3; j++) { |
| 620 |
< |
sU[i][j] /= len; |
| 621 |
< |
} |
| 622 |
< |
} |
| 623 |
< |
|
| 624 |
< |
// sU now contains the coordinates of the 'special' frame; |
| 625 |
< |
|
| 626 |
< |
orient_good = true; |
| 627 |
< |
|
| 634 |
> |
//for(atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter) |
| 635 |
> |
// (*atomIter)->accept(v); |
| 636 |
|
} |
