--- trunk/SHAPES/RigidBody.cpp	2004/06/18 19:01:40	1279
+++ trunk/SHAPES/RigidBody.cpp	2004/06/21 15:54:27	1283
@@ -1,4 +1,5 @@
 #include <math.h>
+#include <iostream>
 #include "RigidBody.hpp"
 #include "VDWAtom.hpp"
 #include "MatVec3.h"
@@ -47,6 +48,8 @@ void RigidBody::setEuler( double phi, double theta, do
     A[2][1] = -cos(phi) * sin(theta);
     A[2][2] = cos(theta);
 
+    printf("A[2][x] = %lf\t%lf\t%lf\n", A[2][0], A[2][1], A[2][2]);
+
 }
 
 void RigidBody::getQ( double q[4] ){
@@ -181,10 +184,10 @@ void RigidBody::calcRefCoords( ) {
   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];
-
+  double test[3];
+  
   // First, find the center of mass:
   
   mass = 0.0;
@@ -196,7 +199,6 @@ void RigidBody::calcRefCoords( ) {
     mass += mtmp;
 
     apos = refCoords[i];
-    
     for(j = 0; j < 3; j++) {
       refCOM[j] += apos[j]*mtmp;     
     }    
@@ -260,7 +262,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"
@@ -269,73 +271,7 @@ 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++) {
@@ -348,6 +284,64 @@ void RigidBody::calcRefCoords( ) {
       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);
+  
+  
+  for (i=0; i<myAtoms.size(); i++){
+    apos = refCoords[i];
+    printf("%f\t%f\t%f\n",apos[0],apos[1],apos[2]);
+  }
+    
+  //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);
+  }
+   
+  for (i=0; i<myAtoms.size(); i++){
+    apos = refCoords[i];
+    printf("%f\t%f\t%f\n",apos[0],apos[1],apos[2]);
+  }
+  
+  identityMat3(iMat);
+  setA(iMat);
 }
 
 void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
@@ -542,3 +536,15 @@ 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();
+  
+}