--- trunk/OOPSE/libmdtools/RigidBody.cpp	2004/04/15 02:45:10	1112
+++ trunk/OOPSE/libmdtools/RigidBody.cpp	2004/04/15 16:18:26	1113
@@ -6,8 +6,6 @@ RigidBody::RigidBody() : StuntDouble() {
 
 RigidBody::RigidBody() : StuntDouble() {
   objType = OT_RIGIDBODY;
-  com_good = false;
-  precalc_done = false;
 }
 
 RigidBody::~RigidBody() {
@@ -97,11 +95,9 @@ void RigidBody::zeroForces() {
     trq[i] = 0.0;
   }
 
-  forces_good = false;
-
 }
 
-void RigidBody::setEulerAngles( double phi, double theta, double psi ){
+void RigidBody::setEuler( double phi, double theta, double psi ){
   
     A[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
     A[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
@@ -191,7 +187,7 @@ void RigidBody::getA( double the_A[3][3] ){
   
   for (int i = 0; i < 3; i++) 
     for (int j = 0; j < 3; j++) 
-      the_A[i][j] = the_A[i][j];
+      the_A[i][j] = A[i][j];
 
 }
 
@@ -229,15 +225,10 @@ void RigidBody::getI( double the_I[3][3] ){  
 
 void RigidBody::getI( double the_I[3][3] ){  
 
-  if (precalc_done) {
-
     for (int i = 0; i < 3; i++) 
       for (int j = 0; j < 3; j++) 
         the_I[i][j] = I[i][j];
 
-  } else {
-
-  }
 }
 
 void RigidBody::lab2Body( double r[3] ){
@@ -270,11 +261,18 @@ void RigidBody::calcRefCoords( ) {
 
 void RigidBody::calcRefCoords( ) {
 
-  int i,j,k;
+  int i,j,k, it;
   double mtmp;
   vec3 apos;
   double refCOM[3];
+  vec3 ptmp;
+  double Itmp[3][3];
+  double evals[3];
+  double evects[3][3];
+  double r, r2, len;
 
+  // First, find the center of mass:
+  
   mass = 0.0;
   for (j=0; j<3; j++)
     refCOM[j] = 0.0;
@@ -293,14 +291,62 @@ void RigidBody::calcRefCoords( ) {
   for(j = 0; j < 3; j++) 
     refCOM[j] /= mass;
 
+// Next, move the origin of the reference coordinate system to the COM:
+
   for (i = 0; i < myAtoms.size(); i++) {
     apos = refCoords[i];
     for (j=0; j < 3; j++) {
       apos[j] = apos[j] - refCOM[j];
     }
     refCoords[i] = apos;
+  }
+
+// Moment of Inertia calculation
+
+  for (i = 0; i < 3; i++) 
+    for (j = 0; j < 3; j++)
+      Itmp[i][j] = 0.0;  
+  
+  for (it = 0; it < myAtoms.size(); it++) {
+
+    mtmp = myAtoms[it]->getMass();
+    ptmp = refCoords[it];
+    r= norm3(ptmp.vec);
+    r2 = r*r;
+    
+    for (i = 0; i < 3; i++) {
+      for (j = 0; j < 3; j++) {
+        
+        if (i==j) Itmp[i][j] += mtmp * r2;
+
+        Itmp[i][j] -= mtmp * ptmp.vec[i]*ptmp.vec[j];
+      }
+    }
   }
+  
+  diagonalize3x3(Itmp, evals, sU);
+  
+  // zero out I and then fill the diagonals with the moments of inertia:
 
+  for (i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      I[i][j] = 0.0;  
+    }
+    I[i][i] = evals[i];
+  }
+  
+  // renormalize column vectors:
+  
+  for (i=0; i < 3; i++) {
+    len = 0.0;
+    for (j = 0; j < 3; j++) {
+      len += sU[i][j]*sU[i][j];
+    }
+    len = sqrt(len);
+    for (j = 0; j < 3; j++) {
+      sU[i][j] /= len;
+    }
+  }
 }
 
 void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
@@ -366,8 +412,6 @@ void RigidBody::calcForcesAndTorques() {
   // (Actually, on second thought, don't.  Integrator does this now.)
   // lab2Body(trq);
 
-  forces_good = true;
-
 }
 
 void RigidBody::updateAtoms() {
@@ -556,73 +600,4 @@ void RigidBody::findCOM() {
     vel[j] /= mass;
   }
 
-  com_good = true;
 }
-  
-void RigidBody::findOrient() {
-  
-  size_t it;
-  int i, j;
-  double ptmp[3];
-  double Itmp[3][3];
-  double evals[3];
-  double evects[3][3];
-  double r2, mtmp, len;
-
-  if (!com_good) findCOM();
-
-  // Calculate inertial tensor matrix elements:
-
-  for (i = 0; i < 3; i++) 
-    for (j = 0; j < 3; j++)
-      Itmp[i][j] = 0.0;  
-  
-  for (it = 0; it < myAtoms.size(); it++) {
-    
-    mtmp = myAtoms[it]->getMass();    
-    myAtoms[it]->getPos(ptmp);
-
-    for (j = 0; j < 3; j++) 
-      ptmp[j] = pos[j] - ptmp[j];
-
-    r2 = norm3(ptmp);
-    
-    for (i = 0; i < 3; i++) {
-      for (j = 0; j < 3; j++) {
-        
-        if (i==j) Itmp[i][j] = mtmp * r2;
-
-        Itmp[i][j] -= mtmp * ptmp[i]*ptmp[j];
-      }
-    }
-  }
-  
-  diagonalize3x3(Itmp, evals, sU);
-  
-  // zero out I and then fill the diagonals with the moments of inertia:
-
-  for (i = 0; i < 3; i++) {
-    for (j = 0; j < 3; j++) {
-      I[i][j] = 0.0;  
-    }
-    I[i][i] = evals[i];
-  }
-  
-  // renormalize column vectors:
-  
-  for (i=0; i < 3; i++) {
-    len = 0.0;
-    for (j = 0; j < 3; j++) {
-      len += sU[i][j]*sU[i][j];
-    }
-    len = sqrt(len);
-    for (j = 0; j < 3; j++) {
-      sU[i][j] /= len;
-    }
-  }
-  
-  // sU now contains the coordinates of the 'special' frame;
-  
-  orient_good = true;
-
-}