--- trunk/SHAPES/RigidBody.cpp	2004/06/15 20:20:36	1271
+++ trunk/SHAPES/RigidBody.cpp	2004/06/22 18:04:58	1285
@@ -1,4 +1,5 @@
 #include <math.h>
+#include <iostream>
 #include "RigidBody.hpp"
 #include "VDWAtom.hpp"
 #include "MatVec3.h"
@@ -15,8 +16,6 @@ void RigidBody::addAtom(VDWAtom* at) {
 void RigidBody::addAtom(VDWAtom* at) {
 
   vec3 coords;
-  vec3 euler;
-  mat3x3 Atmp;
 
   myAtoms.push_back(at);
    
@@ -175,16 +174,18 @@ 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 r, r2, len;
-
+  double iMat[3][3];
+  double test[3];
+  
   // First, find the center of mass:
   
   mass = 0.0;
@@ -196,7 +197,6 @@ void RigidBody::calcRefCoords( ) {
     mass += mtmp;
 
     apos = refCoords[i];
-    
     for(j = 0; j < 3; j++) {
       refCOM[j] += apos[j]*mtmp;     
     }    
@@ -260,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"
@@ -269,7 +269,7 @@ void RigidBody::calcRefCoords( ) {
                );
                exit(-1);     
   }
-  
+
   // renormalize column vectors:
   
   for (i=0; i < 3; i++) {
@@ -280,11 +280,59 @@ 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);
 }
 
-void RigidBody::doEulerToRotMat(vec3 &euler, mat3x3 &myA ){
+void RigidBody::doEulerToRotMat(double euler[3], double myA[3][3] ){
 
   double phi, theta, psi;
   
@@ -344,12 +392,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 +445,36 @@ double RigidBody::min(double x, double  y) {  
   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;
   }
@@ -457,3 +524,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();
+  
+}