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

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 1113 by tim, Thu Apr 15 16:18:26 2004 UTC vs.
Revision 1125 by gezelter, Mon Apr 19 22:13:01 2004 UTC

# Line 687 | Line 687 | template<typename T> void Integrator<T>::rotationPropa
687  
688    double angle;
689    double A[3][3], I[3][3];
690 +  int i, j, k;
691  
692    // use the angular velocities to propagate the rotation matrix a
693    // full time step
694  
695    sd->getA(A);
696    sd->getI(I);
696
697  // rotate about the x-axis
698  angle = dt2 * ji[0] / I[0][0];
699  this->rotate( 1, 2, angle, ji, A );
697  
698 <  // rotate about the y-axis
699 <  angle = dt2 * ji[1] / I[1][1];
700 <  this->rotate( 2, 0, angle, ji, A );
701 <
702 <  // rotate about the z-axis
703 <  angle = dt * ji[2] / I[2][2];
704 <  this->rotate( 0, 1, angle, ji, A);
698 >  if (sd->isLinear()) {
699 >    i = sd->linearAxis();
700 >    j = (i+1)%3;
701 >    k = (i+2)%3;
702 >    
703 >    angle = dt2 * ji[j] / I[j][j];
704 >    this->rotate( k, i, angle, ji, A );
705  
706 <  // rotate about the y-axis
707 <  angle = dt2 * ji[1] / I[1][1];
711 <  this->rotate( 2, 0, angle, ji, A );
706 >    angle = dt * ji[k] / I[k][k];
707 >    this->rotate( i, j, angle, ji, A);
708  
709 <  // rotate about the x-axis
710 <  angle = dt2 * ji[0] / I[0][0];
711 <  this->rotate( 1, 2, angle, ji, A );
712 <
709 >    angle = dt2 * ji[j] / I[j][j];
710 >    this->rotate( k, i, angle, ji, A );
711 >
712 >  } else {
713 >    // rotate about the x-axis
714 >    angle = dt2 * ji[0] / I[0][0];
715 >    this->rotate( 1, 2, angle, ji, A );
716 >    
717 >    // rotate about the y-axis
718 >    angle = dt2 * ji[1] / I[1][1];
719 >    this->rotate( 2, 0, angle, ji, A );
720 >    
721 >    // rotate about the z-axis
722 >    angle = dt * ji[2] / I[2][2];
723 >    this->rotate( 0, 1, angle, ji, A);
724 >    
725 >    // rotate about the y-axis
726 >    angle = dt2 * ji[1] / I[1][1];
727 >    this->rotate( 2, 0, angle, ji, A );
728 >    
729 >    // rotate about the x-axis
730 >    angle = dt2 * ji[0] / I[0][0];
731 >    this->rotate( 1, 2, angle, ji, A );
732 >    
733 >  }
734    sd->setA( A  );
735   }
736  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines