--- trunk/SHAPES/RigidBody.cpp	2004/06/15 20:20:36	1271
+++ trunk/SHAPES/RigidBody.cpp	2004/06/18 19:01:40	1279
@@ -15,8 +15,6 @@ void RigidBody::addAtom(VDWAtom* at) {
 void RigidBody::addAtom(VDWAtom* at) {
 
   vec3 coords;
-  vec3 euler;
-  mat3x3 Atmp;
 
   myAtoms.push_back(at);
    
@@ -175,15 +173,17 @@ void RigidBody::calcRefCoords( ) {
 
 void RigidBody::calcRefCoords( ) {
 
-  int i,j,k, 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 evects[3][3];
+  double prePos[3], rotPos[3];
   double r, r2, len;
+  double iMat[3][3];
 
   // First, find the center of mass:
   
@@ -270,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++) {
@@ -284,7 +350,7 @@ void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &m
   }
 }
 
-void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
+void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
 
   double phi, theta, psi;
   
@@ -344,12 +410,9 @@ void RigidBody::getEulerAngles(double myEuler[3]) {
   
   
   double phi,theta,psi,eps;
-  double pi;
-  double cphi,ctheta,cpsi;
-  double sphi,stheta,spsi;
-  double b[3];
-  int flip[3];
-  
+  double ctheta;
+  double stheta;
+    
   // set the tolerance for Euler angles and rotation elements
   
   eps = 1.0e-8;
@@ -400,14 +463,36 @@ void RigidBody::findCOM() {
   return (x > y) ? y : x;
 }
 
+double RigidBody::findMaxExtent(){
+	int i;
+	double refAtomPos[3];
+	double maxExtent;
+	double tempExtent;
+	
+	//zero the extent variables
+	maxExtent = 0.0;
+	tempExtent = 0.0;
+	for (i=0; i<3; i++)
+		refAtomPos[i] = 0.0;
+	
+	//loop over all atoms
+	for (i=0; i<myAtoms.size(); i++){
+		getAtomRefCoor(refAtomPos, i);
+		tempExtent = sqrt(refAtomPos[0]*refAtomPos[0] + refAtomPos[1]*refAtomPos[1] 
+						  + refAtomPos[2]*refAtomPos[2]);
+		if (tempExtent > maxExtent)
+			maxExtent = tempExtent;
+	}
+	return maxExtent;
+}
+
 void RigidBody::findCOM() {
   
   size_t i;
   int j;
   double mtmp;
   double ptmp[3];
-  double vtmp[3];
-  
+   
   for(j = 0; j < 3; j++) {
     pos[j] = 0.0;
   }