ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DirectionalAtom.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/DirectionalAtom.cpp (file contents):
Revision 670 by mmeineke, Thu Aug 7 21:47:18 2003 UTC vs.
Revision 1046 by tim, Tue Feb 10 21:33:44 2004 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2  
3   #include "Atom.hpp"
4   #include "simError.h"
5  
6 +
7 +
8   void DirectionalAtom::zeroForces() {
9    if( hasCoords ){
10      frc[offsetX] = 0.0;
# Line 23 | Line 25 | double DirectionalAtom::getMu( void ) {
25    }
26   }
27  
28 < double DirectionalAtom::getMu( void ) {
28 > void DirectionalAtom::setCoords(void){
29  
30 <  if( hasCoords ){
31 <    return mu[index];
30 >  if( myConfig->isAllocated() ){
31 >
32 >    myConfig->getAtomPointers( index,
33 >                     &pos,
34 >                     &vel,
35 >                     &frc,
36 >                     &trq,
37 >                     &Amat,
38 >                     &mu,  
39 >                     &ul );
40    }
41    else{
32 –    
42      sprintf( painCave.errMsg,
43 <             "Attempt to get Mu for atom %d before coords set.\n",
44 <             index );
43 >             "Attempted to set Atom %d  coordinates with an unallocated "
44 >             "SimState object.\n", index );
45      painCave.isFatal = 1;
46      simError();
47    }
48 <  return 0;
48 >
49 >  hasCoords = true;
50 >
51 >  *mu = myMu;
52 >
53   }
54  
55 + double DirectionalAtom::getMu( void ) {
56 +
57 +  if( hasCoords ){
58 +    return *mu;
59 +  }
60 +  else{
61 +    return myMu;
62 +  }
63 + }
64 +
65   void DirectionalAtom::setMu( double the_mu ) {
66  
67    if( hasCoords ){
68 <    mu[index] = the_mu;
68 >    *mu = the_mu;
69 >    myMu = the_mu;
70    }
71    else{
72 <    
49 <    sprintf( painCave.errMsg,
50 <             "Attempt to set Mu for atom %d before coords set.\n",
51 <             index );
52 <    painCave.isFatal = 1;
53 <    simError();
72 >    myMu = the_mu;
73    }
74   }
75  
# Line 383 | Line 402 | void DirectionalAtom::getI( double the_I[3][3] ){
402    the_I[2][0] = Izx;
403    the_I[2][1] = Izy;
404    the_I[2][2] = Izz;
405 + }
406 +
407 + void DirectionalAtom::getGrad( double grad[6] ) {
408 +
409 +  double myEuler[3];
410 +  double phi, theta, psi;
411 +  double cphi, sphi, ctheta, stheta;
412 +  double ephi[3];
413 +  double etheta[3];
414 +  double epsi[3];
415 +
416 +  this->getEulerAngles(myEuler);
417 +
418 +  phi = myEuler[0];
419 +  theta = myEuler[1];
420 +  psi = myEuler[2];
421 +
422 +  cphi = cos(phi);
423 +  sphi = sin(phi);
424 +  ctheta = cos(theta);
425 +  stheta = sin(theta);
426 +
427 +  // get unit vectors along the phi, theta and psi rotation axes
428 +
429 +  ephi[0] = 0.0;
430 +  ephi[1] = 0.0;
431 +  ephi[2] = 1.0;
432 +
433 +  etheta[0] = cphi;
434 +  etheta[1] = sphi;
435 +  etheta[2] = 0.0;
436 +  
437 +  epsi[0] = stheta * cphi;
438 +  epsi[1] = stheta * sphi;
439 +  epsi[2] = ctheta;
440 +  
441 +  for (int j = 0 ; j<3; j++)
442 +    grad[j] = frc[j];
443 +
444 +  grad[3] = 0;
445 +  grad[4] = 0;
446 +  grad[5] = 0;
447 +
448 +  for (int j = 0; j < 3; j++ ) {
449 +    
450 +    grad[3] += trq[j]*ephi[j];
451 +    grad[4] += trq[j]*etheta[j];
452 +    grad[5] += trq[j]*epsi[j];
453 +    
454 +  }
455 +
456   }
457 +
458 + /**
459 +  * getEulerAngles computes a set of Euler angle values consistent
460 +  *  with an input rotation matrix.  They are returned in the following
461 +  * order:
462 +  *  myEuler[0] = phi;
463 +  *  myEuler[1] = theta;
464 +  *  myEuler[2] = psi;
465 + */
466 + void DirectionalAtom::getEulerAngles(double myEuler[3]) {
467 +
468 +  // We use so-called "x-convention", which is the most common definition.
469 +  // In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
470 +  // rotation is by an angle phi about the z-axis, the second is by an angle  
471 +  // theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
472 +  //z-axis (again).
473 +  
474 +  
475 +  double phi,theta,psi,eps;
476 +  double pi;
477 +  double cphi,ctheta,cpsi;
478 +  double sphi,stheta,spsi;
479 +  double b[3];
480 +  int flip[3];
481 +
482 +  // set the tolerance for Euler angles and rotation elements
483 +  
484 +  eps = 1.0e-8;
485 +
486 +  theta = acos(min(1.0,max(-1.0,Amat[Azz])));
487 +  ctheta = Amat[Azz];
488 +  stheta = sqrt(1.0 - ctheta * ctheta);
489 +
490 +  // when sin(theta) is close to 0, we need to consider singularity
491 +  // In this case, we can assign an arbitary value to phi (or psi), and then determine
492 +  // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
493 +  // in cases of singularity.  
494 +  // we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
495 +  // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
496 +  // change the sign of both of the parameters passed to atan2.
497 +  
498 +  if (fabs(stheta) <= eps){
499 +    psi = 0.0;
500 +    phi = atan2(-Amat[Ayx], Amat[Axx]);  
501 +  }
502 +  // we only have one unique solution
503 +  else{    
504 +      phi = atan2(Amat[Azx], -Amat[Azy]);
505 +      psi = atan2(Amat[Axz], Amat[Ayz]);
506 +  }
507 +
508 +  //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
509 +  //if (phi < 0)
510 +  //  phi += M_PI;
511 +
512 +  //if (psi < 0)
513 +  //  psi += M_PI;
514 +
515 +  myEuler[0] = phi;
516 +  myEuler[1] = theta;
517 +  myEuler[2] = psi;
518 +  
519 +  return;
520 + }
521 +
522 + double DirectionalAtom::max(double x, double  y) {  
523 +  return (x > y) ? x : y;
524 + }
525 +
526 + double DirectionalAtom::min(double x, double  y) {  
527 +  return (x > y) ? y : x;
528 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines