164 |
|
} |
165 |
|
|
166 |
|
// Moment of Inertia calculation |
167 |
< |
Mat3x3d Itmp(0.0); |
168 |
< |
|
167 |
> |
Mat3x3d Itmp(0.0); |
168 |
|
for (std::size_t i = 0; i < atoms_.size(); i++) { |
169 |
+ |
Mat3x3d IAtom(0.0); |
170 |
|
mtmp = atoms_[i]->getMass(); |
171 |
< |
Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; |
171 |
> |
IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; |
172 |
|
double r2 = refCoords_[i].lengthSquare(); |
173 |
< |
Itmp(0, 0) += mtmp * r2; |
174 |
< |
Itmp(1, 1) += mtmp * r2; |
175 |
< |
Itmp(2, 2) += mtmp * r2; |
176 |
< |
} |
177 |
< |
|
178 |
< |
//project the inertial moment of directional atoms into this rigid body |
179 |
< |
for (std::size_t i = 0; i < atoms_.size(); i++) { |
173 |
> |
IAtom(0, 0) += mtmp * r2; |
174 |
> |
IAtom(1, 1) += mtmp * r2; |
175 |
> |
IAtom(2, 2) += mtmp * r2; |
176 |
> |
Itmp += IAtom; |
177 |
> |
|
178 |
> |
//project the inertial moment of directional atoms into this rigid body |
179 |
|
if (atoms_[i]->isDirectional()) { |
180 |
< |
Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; |
181 |
< |
} |
180 |
> |
Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; |
181 |
> |
} |
182 |
|
} |
183 |
|
|
184 |
+ |
// std::cout << Itmp << std::endl; |
185 |
+ |
|
186 |
|
//diagonalize |
187 |
|
Vector3d evals; |
188 |
|
Mat3x3d::diagonalize(Itmp, evals, sU_); |
483 |
|
"RigidBody error.\n" |
484 |
|
"\tAtom %s does not have a position specified.\n" |
485 |
|
"\tThis means RigidBody cannot set up reference coordinates.\n", |
486 |
< |
ats->getType() ); |
486 |
> |
ats->getType().c_str() ); |
487 |
|
painCave.isFatal = 1; |
488 |
|
simError(); |
489 |
|
} |
503 |
|
"RigidBody error.\n" |
504 |
|
"\tAtom %s does not have an orientation specified.\n" |
505 |
|
"\tThis means RigidBody cannot set up reference orientations.\n", |
506 |
< |
ats->getType() ); |
506 |
> |
ats->getType().c_str() ); |
507 |
|
painCave.isFatal = 1; |
508 |
|
simError(); |
509 |
|
} |