--- trunk/src/primitives/RigidBody.cpp 2005/04/15 22:04:00 507 +++ trunk/src/primitives/RigidBody.cpp 2005/10/23 21:08:08 695 @@ -164,27 +164,25 @@ namespace oopse { } // Moment of Inertia calculation - Mat3x3d Itmp(0.0); - + Mat3x3d Itmp(0.0); for (std::size_t i = 0; i < atoms_.size(); i++) { + Mat3x3d IAtom(0.0); mtmp = atoms_[i]->getMass(); - Itmp -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; + IAtom -= outProduct(refCoords_[i], refCoords_[i]) * mtmp; double r2 = refCoords_[i].lengthSquare(); - Itmp(0, 0) += mtmp * r2; - Itmp(1, 1) += mtmp * r2; - Itmp(2, 2) += mtmp * r2; - } + IAtom(0, 0) += mtmp * r2; + IAtom(1, 1) += mtmp * r2; + IAtom(2, 2) += mtmp * r2; + Itmp += IAtom; - //project the inertial moment of directional atoms into this rigid body - for (std::size_t i = 0; i < atoms_.size(); i++) { + //project the inertial moment of directional atoms into this rigid body if (atoms_[i]->isDirectional()) { - RectMatrix Iproject = refOrients_[i].transpose() * atoms_[i]->getI(); - Itmp(0, 0) += Iproject(0, 0); - Itmp(1, 1) += Iproject(1, 1); - Itmp(2, 2) += Iproject(2, 2); - } + Itmp += refOrients_[i].transpose() * atoms_[i]->getI() * refOrients_[i]; + } } + // std::cout << Itmp << std::endl; + //diagonalize Vector3d evals; Mat3x3d::diagonalize(Itmp, evals, sU_); @@ -273,8 +271,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(a * refOrients_[i]); - //dAtom->rotateBy( A ); + dAtom->setA(refOrients_[i] * a); } } @@ -301,7 +298,7 @@ namespace oopse { if (atoms_[i]->isDirectional()) { dAtom = (DirectionalAtom *) atoms_[i]; - dAtom->setA(a * refOrients_[i], frame); + dAtom->setA(refOrients_[i] * a, frame); } }