--- trunk/OOPSE/libmdtools/RigidBody.cpp	2004/06/21 15:54:27	1283
+++ trunk/OOPSE/libmdtools/RigidBody.cpp	2004/06/21 18:52:21	1284
@@ -645,8 +645,18 @@ void RigidBody::accept(BaseVisitor* v){
 
   //for(atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter)
   //  (*atomIter)->accept(v);
+}
+void RigidBody::getAtomRefCoor(double pos[3], int index){
+  vec3 ref;
+
+  ref = refCoords[index];
+  pos[0] = ref[0];
+  pos[1] = ref[1];
+  pos[2] = ref[2];
+  
 }
 
+
 void RigidBody::getAtomPos(double theP[3], int index){
   vec3 ref;
 
@@ -662,12 +672,40 @@ void RigidBody::getAtomRefCoor(double pos[3], int inde
 }
 
 
-void RigidBody::getAtomRefCoor(double pos[3], int index){
+void RigidBody::getAtomVel(double theV[3], int index){
   vec3 ref;
+  double velRot[3];
+  double skewMat[3][3];
+  double aSkewMat[3][3];
+  double aSkewTransMat[3][3];
+  
+  //velRot = $(A\cdot skew(I^{-1}j))^{T}refCoor$
 
+  if (index >= myAtoms.size())
+    cerr << index << " is an invalid index, current rigid body contains " << myAtoms.size() << "atoms" << endl;
+
   ref = refCoords[index];
-  pos[0] = ref[0];
-  pos[1] = ref[1];
-  pos[2] = ref[2];
+
+  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;
   
+  matMul3(A, skewMat, aSkewMat);
+
+  transposeMat3(aSkewMat, aSkewTransMat);
+
+  matVecMul3(aSkewTransMat, ref.vec, velRot);
+  theV[0] = vel[0] + velRot[0];
+  theV[1] = vel[1] + velRot[1];
+  theV[2] = vel[2] + velRot[2];
 }
+
+