--- trunk/OOPSE/libmdtools/DirectionalAtom.cpp	2003/12/12 15:42:13	878
+++ trunk/OOPSE/libmdtools/DirectionalAtom.cpp	2004/02/10 21:33:44	1046
@@ -429,16 +429,22 @@ void DirectionalAtom::getGrad( double grad[6] ) {
   ephi[0] = 0.0;
   ephi[1] = 0.0;
   ephi[2] = 1.0;
-  etheta[0] = -sphi;
-  etheta[1] = cphi;
+
+  etheta[0] = cphi;
+  etheta[1] = sphi;
   etheta[2] = 0.0;
-  epsi[0] = ctheta * cphi;
-  epsi[1] = ctheta * sphi;
-  epsi[2] = -stheta;
+  
+  epsi[0] = stheta * cphi;
+  epsi[1] = stheta * sphi;
+  epsi[2] = ctheta;
   
   for (int j = 0 ; j<3; j++)
     grad[j] = frc[j];
 
+  grad[3] = 0;
+  grad[4] = 0;
+  grad[5] = 0;
+
   for (int j = 0; j < 3; j++ ) {
     
     grad[3] += trq[j]*ephi[j];
@@ -449,16 +455,23 @@ void DirectionalAtom::getGrad( double grad[6] ) {
 
 }
 
-
+/**
+  * getEulerAngles computes a set of Euler angle values consistent
+  *  with an input rotation matrix.  They are returned in the following
+  * order:
+  *  myEuler[0] = phi;
+  *  myEuler[1] = theta;
+  *  myEuler[2] = psi;
+*/
 void DirectionalAtom::getEulerAngles(double myEuler[3]) {
 
-  // getEulerAngles computes a set of Euler angle values consistent
-  // with an input rotation matrix.  They are returned in the following
-  // order:
-  //  myEuler[0] = phi;
-  //  myEuler[1] = theta;
-  //  myEuler[2] = psi;
+  // We use so-called "x-convention", which is the most common definition. 
+  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first 
+  // rotation is by an angle phi about the z-axis, the second is by an angle  
+  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
+  //z-axis (again). 
   
+  
   double phi,theta,psi,eps;
   double pi;
   double cphi,ctheta,cpsi;
@@ -469,80 +482,40 @@ void DirectionalAtom::getEulerAngles(double myEuler[3]
   // set the tolerance for Euler angles and rotation elements
   
   eps = 1.0e-8;
-    
-  // get a trial value of theta from a single rotation element
-  
-  theta = asin(min(1.0,max(-1.0,-Amat[Axz])));
-  ctheta = cos(theta);
-  stheta = -Amat[Axz];
-  
-  // set the phi/psi difference when theta is either 90 or -90
-  
-  if (fabs(ctheta) <= eps) {
-    phi = 0.0;
-    if (fabs(Amat[Azx]) < eps) {
-      psi = asin(min(1.0,max(-1.0,-Amat[Ayx]/Amat[Axz])));
-    } else {
-      if (fabs(Amat[Ayx]) < eps) {
-        psi = acos(min(1.0,max(-1.0,-Amat[Azx]/Amat[Axz])));
-      } else {
-        psi = atan(Amat[Ayx]/Amat[Azx]);
-      }    
-    }
-  } 
 
-  // set the phi and psi values for all other theta values
+  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
+  ctheta = Amat[Azz]; 
+  stheta = sqrt(1.0 - ctheta * ctheta);
+
+  // when sin(theta) is close to 0, we need to consider singularity
+  // In this case, we can assign an arbitary value to phi (or psi), and then determine 
+  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
+  // in cases of singularity.  
+  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi. 
+  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
+  // change the sign of both of the parameters passed to atan2.
   
-  else {
-    if (fabs(Amat[Axx]) < eps) {
-      phi = asin(min(1.0,max(-1.0,Amat[Axy]/ctheta)));
-    } else {
-      if (fabs(Amat[Axy]) < eps) {
-        phi = acos(min(1.0,max(-1.0,Amat[Axx]/ctheta)));
-      } else {
-        phi = atan(Amat[Axy]/Amat[Axx]);
-      }
-    }
-    if (fabs(Amat[Azz]) < eps) {
-      psi = asin(min(1.0,max(-1.0,Amat[Ayz]/ctheta)));
-    } else {
-      if (fabs(Amat[Ayz]) < eps) {
-        psi = acos(min(1.0,max(-1.0,Amat[Azz]/ctheta)));
-      }
-      psi = atan(Amat[Ayz]/Amat[Azz]);
-    }
+  if (fabs(stheta) <= eps){
+    psi = 0.0;
+    phi = atan2(-Amat[Ayx], Amat[Axx]);  
   }
-
-  // find sine and cosine of the trial phi and psi values
-
-  cphi = cos(phi);
-  sphi = sin(phi);
-  cpsi = cos(psi);
-  spsi = sin(psi);
-
-  // reconstruct the diagonal of the rotation matrix
-
-  b[0] = ctheta * cphi;
-  b[1] = spsi*stheta*sphi + cpsi*cphi;
-  b[2] = ctheta * cpsi;
-
-  // compare the correct matrix diagonal to rebuilt diagonal
-
-  for (int i = 0; i < 3; i++) {
-    flip[i] = 0;
-    if (fabs(Amat[3*i + i] - b[i]) > eps)  flip[i] = 1;
+  // we only have one unique solution
+  else{    
+      phi = atan2(Amat[Azx], -Amat[Azy]);
+      psi = atan2(Amat[Axz], Amat[Ayz]);
   }
 
-  // alter Euler angles to get correct rotation matrix values
-  
-  if (flip[0] && flip[1]) phi = phi - copysign(M_PI,phi);
-  if (flip[0] && flip[2]) theta = -theta + copysign(M_PI, theta);
-  if (flip[1] && flip[2]) psi = psi - copysign(M_PI, psi);
+  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
+  //if (phi < 0)
+  //  phi += M_PI;
 
+  //if (psi < 0)
+  //  psi += M_PI;
+
   myEuler[0] = phi;
   myEuler[1] = theta;
   myEuler[2] = psi;
-
+  
   return;
 }