--- trunk/src/primitives/RigidBody.cpp	2005/01/13 19:40:37	253
+++ trunk/src/primitives/RigidBody.cpp	2005/02/13 06:57:48	318
@@ -172,6 +172,16 @@ void  RigidBody::calcRefCoords() {
         Itmp(0, 0) += mtmp * r2;
         Itmp(1, 1) += mtmp * r2;
         Itmp(2, 2) += mtmp * r2;
+    }
+
+    //project the inertial moment of directional atoms into this rigid body
+    for (std::size_t i = 0; i < atoms_.size(); i++) {
+        if (atoms_[i]->isDirectional()) {
+            RectMatrix<double, 3, 3> 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);
+        }
     }
 
     //diagonalize 
@@ -264,13 +274,99 @@ void  RigidBody::updateAtoms() {
           dAtom = (DirectionalAtom *) atoms_[i];
           dAtom->setA(a * refOrients_[i]);
           //dAtom->rotateBy( A );      
+        }
+
+    }
+  
+}
+
+
+void  RigidBody::updateAtoms(int frame) {
+    unsigned int i;
+    Vector3d ref;
+    Vector3d apos;
+    DirectionalAtom* dAtom;
+    Vector3d pos = getPos(frame);
+    RotMat3x3d a = getA(frame);
+    
+    for (i = 0; i < atoms_.size(); i++) {
+     
+        ref = body2Lab(refCoords_[i]);
+
+        apos = pos + ref;
+
+        atoms_[i]->setPos(apos, frame);
+
+        if (atoms_[i]->isDirectional()) {
+          
+          dAtom = (DirectionalAtom *) atoms_[i];
+          dAtom->setA(a * refOrients_[i], frame);
         }
 
     }
   
+}
+
+void RigidBody::updateAtomVel() {
+    Mat3x3d skewMat;;
+
+    Vector3d ji = getJ();
+    Mat3x3d I =  getI();
+
+    skewMat(0, 0) =0;
+    skewMat(0, 1) = ji[2] /I(2, 2);
+    skewMat(0, 2) = -ji[1] /I(1, 1);
+
+    skewMat(1, 0) = -ji[2] /I(2, 2);
+    skewMat(1, 1) = 0;
+    skewMat(1, 2) = ji[0]/I(0, 0);
+
+    skewMat(2, 0) =ji[1] /I(1, 1);
+    skewMat(2, 1) = -ji[0]/I(0, 0);
+    skewMat(2, 2) = 0;
+
+    Mat3x3d mat = (getA() * skewMat).transpose();
+    Vector3d rbVel = getVel();
+
+
+    Vector3d velRot;        
+    for (int i =0 ; i < refCoords_.size(); ++i) {
+        atoms_[i]->setVel(rbVel + mat * refCoords_[i]);
+    }
+
 }
 
+void RigidBody::updateAtomVel(int frame) {
+    Mat3x3d skewMat;;
 
+    Vector3d ji = getJ(frame);
+    Mat3x3d I =  getI();
+
+    skewMat(0, 0) =0;
+    skewMat(0, 1) = ji[2] /I(2, 2);
+    skewMat(0, 2) = -ji[1] /I(1, 1);
+
+    skewMat(1, 0) = -ji[2] /I(2, 2);
+    skewMat(1, 1) = 0;
+    skewMat(1, 2) = ji[0]/I(0, 0);
+
+    skewMat(2, 0) =ji[1] /I(1, 1);
+    skewMat(2, 1) = -ji[0]/I(0, 0);
+    skewMat(2, 2) = 0;
+
+    Mat3x3d mat = (getA(frame) * skewMat).transpose();
+    Vector3d rbVel = getVel(frame);
+
+
+    Vector3d velRot;        
+    for (int i =0 ; i < refCoords_.size(); ++i) {
+        atoms_[i]->setVel(rbVel + mat * refCoords_[i], frame);
+    }
+
+}
+
+        
+
 bool RigidBody::getAtomPos(Vector3d& pos, unsigned int index) {
     if (index < atoms_.size()) {