--- trunk/OOPSE/libmdtools/RigidBody.cpp	2004/04/19 20:54:58	1120
+++ trunk/OOPSE/libmdtools/RigidBody.cpp	2004/06/21 18:52:21	1284
@@ -8,6 +8,7 @@ RigidBody::RigidBody() : StuntDouble() {
   objType = OT_RIGIDBODY;
   is_linear = false;
   linear_axis =  -1;
+  momIntTol = 1e-6;
 }
 
 RigidBody::~RigidBody() {
@@ -261,6 +262,18 @@ void RigidBody::body2Lab( double r[3] ){
 
 }
 
+double RigidBody::getZangle( ){
+    return zAngle;
+}
+
+void RigidBody::setZangle( double zAng ){
+    zAngle = zAng;
+}
+
+void RigidBody::addZangle( double zAng ){
+    zAngle += zAng;
+}
+
 void RigidBody::calcRefCoords( ) {
 
   int i,j,k, it, n_linear_coords;
@@ -632,4 +645,67 @@ void RigidBody::accept(BaseVisitor* v){
 
   //for(atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter)
   //  (*atomIter)->accept(v);
-}
\ No newline at end of file
+}
+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;
+
+  if (index >= myAtoms.size())
+    cerr << index << " is an invalid index, current rigid body contains " << myAtoms.size() << "atoms" << endl;
+
+  ref = refCoords[index];
+  body2Lab(ref.vec);
+  
+  theP[0] = pos[0] + ref[0];
+  theP[1] = pos[1] + ref[1];
+  theP[2] = pos[2] + ref[2];
+}
+
+
+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];
+
+  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];
+}
+
+