--- trunk/OOPSE/libmdtools/DirectionalAtom.cpp	2003/08/12 19:56:49	689
+++ trunk/OOPSE/libmdtools/DirectionalAtom.cpp	2004/02/10 21:33:44	1046
@@ -1,8 +1,10 @@
-#include <cmath>
+#include <math.h>
 
 #include "Atom.hpp"
 #include "simError.h"
 
+
+
 void DirectionalAtom::zeroForces() {
   if( hasCoords ){
     frc[offsetX] = 0.0; 
@@ -39,32 +41,31 @@ void DirectionalAtom::setCoords(void){
   else{
     sprintf( painCave.errMsg,
 	     "Attempted to set Atom %d  coordinates with an unallocated "
-	     "SimState object.\n" );
+	     "SimState object.\n", index );
     painCave.isFatal = 1;
     simError();
   }
 
   hasCoords = true;
 
-  mu[index] = myMu;
+  *mu = myMu;
 
 }
 
 double DirectionalAtom::getMu( void ) {
 
   if( hasCoords ){
-    return mu[index];
+    return *mu;
   }
   else{
     return myMu;
   }
-  return 0;
 }
 
 void DirectionalAtom::setMu( double the_mu ) { 
 
   if( hasCoords ){
-    mu[index] = the_mu;
+    *mu = the_mu;
     myMu = the_mu;
   }
   else{
@@ -401,4 +402,127 @@ void DirectionalAtom::getI( double the_I[3][3] ){
   the_I[2][0] = Izx;
   the_I[2][1] = Izy;
   the_I[2][2] = Izz;
+}
+
+void DirectionalAtom::getGrad( double grad[6] ) {
+
+  double myEuler[3];
+  double phi, theta, psi;
+  double cphi, sphi, ctheta, stheta;
+  double ephi[3];
+  double etheta[3];
+  double epsi[3];
+
+  this->getEulerAngles(myEuler);
+
+  phi = myEuler[0];
+  theta = myEuler[1];
+  psi = myEuler[2];
+
+  cphi = cos(phi);
+  sphi = sin(phi);
+  ctheta = cos(theta);
+  stheta = sin(theta);
+
+  // get unit vectors along the phi, theta and psi rotation axes
+
+  ephi[0] = 0.0;
+  ephi[1] = 0.0;
+  ephi[2] = 1.0;
+
+  etheta[0] = cphi;
+  etheta[1] = sphi;
+  etheta[2] = 0.0;
+  
+  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];
+    grad[4] += trq[j]*etheta[j];
+    grad[5] += trq[j]*epsi[j];
+    
+  }
+
+}
+
+/**
+  * 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]) {
+
+  // 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;
+  double sphi,stheta,spsi;
+  double b[3];
+  int flip[3];
+
+  // set the tolerance for Euler angles and rotation elements
+  
+  eps = 1.0e-8;
+
+  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.
+  
+  if (fabs(stheta) <= eps){
+    psi = 0.0;
+    phi = atan2(-Amat[Ayx], Amat[Axx]);  
+  }
+  // we only have one unique solution
+  else{    
+      phi = atan2(Amat[Azx], -Amat[Azy]);
+      psi = atan2(Amat[Axz], Amat[Ayz]);
+  }
+
+  //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;
 }
+
+double DirectionalAtom::max(double x, double  y) {  
+  return (x > y) ? x : y;
+}
+
+double DirectionalAtom::min(double x, double  y) {  
+  return (x > y) ? y : x;
+}