--- trunk/SHAPES/RigidBody.cpp	2004/06/17 21:27:38	1276
+++ trunk/SHAPES/RigidBody.cpp	2004/06/18 19:01:40	1279
@@ -173,14 +173,17 @@ void RigidBody::calcRefCoords( ) {
 
 void RigidBody::calcRefCoords( ) {
 
-  int i, j, it, n_linear_coords;
+  int i, j, it, n_linear_coords, pAxis, maxAxis, midAxis;
   double mtmp;
   vec3 apos;
   double refCOM[3];
   vec3 ptmp;
   double Itmp[3][3];
+  double pAxisMat[3][3], pAxisRotMat[3][3];
   double evals[3];
+  double prePos[3], rotPos[3];
   double r, r2, len;
+  double iMat[3][3];
 
   // First, find the center of mass:
   
@@ -267,6 +270,72 @@ void RigidBody::calcRefCoords( ) {
                exit(-1);     
   }
   
+  //sort and reorder the moment axes
+  if (evals[0] < evals[1] && evals[0] < evals[2])
+    pAxis = 0;
+  else if (evals[1] < evals[2])
+    pAxis = 1;
+  else 
+    pAxis = 2;
+  
+  if (evals[0] > evals[1] && evals[0] > evals[2])
+    maxAxis = 0;
+  else if (evals[1] > evals[2])
+    maxAxis = 1;
+  else 
+    maxAxis = 2;
+  
+  midAxis = 0;
+  if (midAxis == pAxis || midAxis == pAxis)
+    midAxis = 1;
+  if (midAxis == pAxis || midAxis == pAxis)
+    midAxis = 2;
+  
+  if (pAxis != maxAxis){
+    //zero out our matrices
+    for (i=0; i<3; i++){
+      for (j=0; j<3; j++) {
+        pAxisMat[i][j] = 0.0;
+        pAxisRotMat[i][j] = 0.0;
+      }
+    }
+    
+    //let z be the smallest and x be the largest eigenvalue axes
+    for (i=0; i<3; i++){
+      pAxisMat[i][2] = I[i][pAxis];
+      pAxisMat[i][1] = I[i][midAxis];
+      pAxisMat[i][0] = I[i][maxAxis];
+    }
+    
+    //calculate the proper rotation matrix
+    transposeMat3(pAxisMat, pAxisRotMat);
+    
+    //rotate the rigid body to the principle axis frame
+    for (i = 0; i < myAtoms.size(); i++) {
+      apos = refCoords[i];
+      for (j=0; j<3; j++)
+        prePos[j] = apos[j];
+      
+      matVecMul3(pAxisRotMat, prePos, rotPos);
+      
+      for (j=0; j < 3; j++) 
+        apos[j] = rotPos[j];
+      
+      refCoords[i] = apos;
+    }
+    
+    //the lab and the body frame match up at this point, so A = Identity Matrix
+    for (i=0; i<3; i++){
+      for (j=0; j<3; j++){
+        if (i == j)
+          iMat[i][j] = 1.0;
+        else
+          iMat[i][j] = 0.0;
+      }
+    }
+    setA(iMat);
+  }
+           
   // renormalize column vectors:
   
   for (i=0; i < 3; i++) {