--- trunk/SHAPES/RigidBody.cpp	2004/06/17 21:27:38	1276
+++ trunk/SHAPES/RigidBody.cpp	2004/06/23 23:19:43	1292
@@ -1,4 +1,5 @@
 #include <math.h>
+#include <iostream>
 #include "RigidBody.hpp"
 #include "VDWAtom.hpp"
 #include "MatVec3.h"
@@ -173,15 +174,18 @@ 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 r, r2, len;
-
+  double iMat[3][3];
+  double test[3];
+  
   // First, find the center of mass:
   
   mass = 0.0;
@@ -193,7 +197,6 @@ void RigidBody::calcRefCoords( ) {
     mass += mtmp;
 
     apos = refCoords[i];
-    
     for(j = 0; j < 3; j++) {
       refCOM[j] += apos[j]*mtmp;     
     }    
@@ -257,7 +260,7 @@ void RigidBody::calcRefCoords( ) {
   if (n_linear_coords > 1) {
           printf(
                "RigidBody error.\n"
-               "\tOOPSE found more than one axis in this rigid body with a vanishing \n"
+               "\tSHAPES found more than one axis in this rigid body with a vanishing \n"
                "\tmoment of inertia.  This can happen in one of three ways:\n"
                "\t 1) Only one atom was specified, or \n"
                "\t 2) All atoms were specified at the same location, or\n"
@@ -266,7 +269,7 @@ void RigidBody::calcRefCoords( ) {
                );
                exit(-1);     
   }
-  
+
   // renormalize column vectors:
   
   for (i=0; i < 3; i++) {
@@ -277,8 +280,65 @@ void RigidBody::calcRefCoords( ) {
     len = sqrt(len);
     for (j = 0; j < 3; j++) {
       sU[i][j] /= len;
+    }
+  }
+  
+  //sort and reorder the moment axes
+
+  // The only problem below is for molecules like C60 with 3 nearly identical 
+  // non-zero moments of inertia.  In this case it doesn't really matter which is 
+  // the principal axis, so they get assigned nearly randomly depending on the 
+  // floating point comparison between eigenvalues
+  if (! is_linear) {
+    pAxis = 0;
+    maxAxis = 0;
+  
+    for (i = 0; i < 3; i++) {
+      if (evals[i] < evals[pAxis]) pAxis = i;
+      if (evals[i] > evals[maxAxis]) maxAxis = i;
     }
+  
+    midAxis = 0;
+    for (i=0; i < 3; i++) {
+      if (pAxis != i && maxAxis != i) midAxis = i;
+    }
+  } else {
+    pAxis = linear_axis;
+    // linear molecules have one zero moment of inertia and two identical
+    // moments of inertia.  In this case, it doesn't matter which is chosen
+    // as mid and which is max, so just permute from the pAxis:
+    midAxis = (pAxis + 1)%3;
+    maxAxis = (pAxis + 2)%3;
+  }
+      
+  //let z be the smallest and x be the largest eigenvalue axes
+  for (i=0; i<3; i++){
+    pAxisMat[i][2] = sU[i][pAxis];
+    pAxisMat[i][1] = sU[i][midAxis];
+    pAxisMat[i][0] = sU[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++) {
+    matVecMul3(pAxisRotMat, refCoords[i].vec, refCoords[i].vec);
+    myAtoms[i]->setPos(refCoords[i].vec);
+  }
+  
+  identityMat3(iMat);
+  setA(iMat);
+  
+  //and resort the moments of intertia to match the new orientation
+  for (i=0; i<3; i++)
+    if (evals[i]<momIntTol)
+      evals[i] = 0.0;
+  I[0][0] = evals[maxAxis];
+  I[1][1] = evals[midAxis];
+  I[2][2] = evals[pAxis];
 }
 
 void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
@@ -473,3 +533,21 @@ void RigidBody::getAtomRefCoor(double pos[3], int inde
   pos[2] = ref[2];
   
 }
+
+double RigidBody::getAtomRpar(int index){
+  
+  return myAtoms[index]->getRpar();
+  
+}
+
+double RigidBody::getAtomEps(int index){
+  
+  return myAtoms[index]->getEps();
+  
+}
+
+char *RigidBody::getAtomBase(int index){
+  
+  return myAtoms[index]->getBase();
+  
+}