--- trunk/OOPSE-1.0/utils/sysbuilder/MoLocator.cpp	2004/07/29 03:31:50	1432
+++ trunk/OOPSE-1.0/utils/sysbuilder/MoLocator.cpp	2004/07/29 18:16:16	1435
@@ -5,6 +5,7 @@
 
 #include "simError.h"
 #include "MoLocator.hpp"
+#include "MatVec3.h"
 
 MoLocator::MoLocator( MoleculeStamp* theStamp, ForceFields* theFF){
 
@@ -20,12 +21,15 @@ void MoLocator::placeMol( const Vector3d& offset, cons
   Vector3d angMomentum(0.0, 0.0, 0.0);
   double quaternion[4];
   vector<StuntDouble*> myIntegrableObjects;
-
+  double rotMat[3][3];
+  
   quaternion[0] = 1.0;
   quaternion[1] = 0.0;
   quaternion[2] = 0.0;
   quaternion[3] = 0.0;
 
+  latVec2RotMat(ort, rotMat);
+  
   myIntegrableObjects = mol->getIntegrableObjects();
 
   if(myIntegrableObjects.size() != nIntegrableObjects){
@@ -39,14 +43,17 @@ void MoLocator::placeMol( const Vector3d& offset, cons
     
   for(int i=0; i<nIntegrableObjects; i++) {
 
-      newCoor = refCoords[i] + offset;
-      myIntegrableObjects[i]->setPos( newCoor.vec);
-      myIntegrableObjects[i]->setVel(velocity.vec);
+    //calculate the reference coordinate for integrable objects after rotation
+    matVecMul3(rotMat, refCoords[i].vec, newCoor.vec);    
+    newCoor +=  offset;
 
-      if(myIntegrableObjects[i]->isDirectional()){
-        myIntegrableObjects[i]->setQ(quaternion);
-        myIntegrableObjects[i]->setJ(angMomentum.vec);  
-      }
+    myIntegrableObjects[i]->setPos( newCoor.vec);
+    myIntegrableObjects[i]->setVel(velocity.vec);
+
+    if(myIntegrableObjects[i]->isDirectional()){
+      myIntegrableObjects[i]->setA(rotMat);
+      myIntegrableObjects[i]->setJ(angMomentum.vec);  
+    }
   }
 
 }
@@ -63,8 +70,6 @@ void MoLocator::calcRefCoords( void ){
   double currAtomMass;
   double molMass;
   
-  mass.resize(nIntegrableObjects);
-  
   nAtoms= myStamp->getNAtoms();
   nRigidBodies = myStamp->getNRigidBodies();
 
@@ -140,3 +145,25 @@ void MoLocator::calcRefCoords( void ){
     refCoords[i] -= refMolCom;
 }
 
+
+void latVec2RotMat(const Vector3d& lv, double rotMat[3][3]){
+
+  double theta, phi, psi;
+  
+  theta =acos(lv.z);
+  phi = atan2(lv.y, lv.x);
+  psi = 0;
+
+  rotMat[0][0] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
+  rotMat[0][1] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
+  rotMat[0][2] = sin(theta) * sin(psi);
+  
+  rotMat[1][0] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
+  rotMat[1][1] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
+  rotMat[1][2] = sin(theta) * cos(psi);
+  
+  rotMat[2][0] = sin(phi) * sin(theta);
+  rotMat[2][1] = -cos(phi) * sin(theta);
+  rotMat[2][2] = cos(theta);
+}
+