| 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 |
< |
} |
| 173 |
> |
IAtom(0, 0) += mtmp * r2; |
| 174 |
> |
IAtom(1, 1) += mtmp * r2; |
| 175 |
> |
IAtom(2, 2) += mtmp * r2; |
| 176 |
|
|
| 177 |
< |
//project the inertial moment of directional atoms into this rigid body |
| 179 |
< |
for (std::size_t i = 0; i < atoms_.size(); i++) { |
| 177 |
> |
//project the inertial moment of directional atoms into this rigid body |
| 178 |
|
if (atoms_[i]->isDirectional()) { |
| 179 |
< |
RectMatrix<double, 3, 3> Iproject = refOrients_[i].transpose() * atoms_[i]->getI(); |
| 180 |
< |
Itmp(0, 0) += Iproject(0, 0); |
| 181 |
< |
Itmp(1, 1) += Iproject(1, 1); |
| 182 |
< |
Itmp(2, 2) += Iproject(2, 2); |
| 179 |
> |
//IAtom += atoms_[i]->getI(); |
| 180 |
> |
Itmp += IAtom; |
| 181 |
> |
Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; |
| 182 |
> |
} else { |
| 183 |
> |
Itmp += IAtom; |
| 184 |
|
} |
| 185 |
|
} |
| 186 |
|
|
| 187 |
+ |
std::cout << Itmp <<std::endl; |
| 188 |
|
//diagonalize |
| 189 |
|
Vector3d evals; |
| 190 |
|
Mat3x3d::diagonalize(Itmp, evals, sU_); |
| 273 |
|
if (atoms_[i]->isDirectional()) { |
| 274 |
|
|
| 275 |
|
dAtom = (DirectionalAtom *) atoms_[i]; |
| 276 |
< |
dAtom->setA(a * refOrients_[i]); |
| 277 |
< |
//dAtom->rotateBy( A ); |
| 276 |
> |
dAtom->setA(refOrients_[i] * a); |
| 277 |
|
} |
| 278 |
|
|
| 279 |
|
} |
| 300 |
|
if (atoms_[i]->isDirectional()) { |
| 301 |
|
|
| 302 |
|
dAtom = (DirectionalAtom *) atoms_[i]; |
| 303 |
< |
dAtom->setA(a * refOrients_[i], frame); |
| 303 |
> |
dAtom->setA(refOrients_[i] * a, frame); |
| 304 |
|
} |
| 305 |
|
|
| 306 |
|
} |
