OOPSE/libmdtools/DirectionalAtom.cpp

#include <math.h>

#include "Atom.hpp"
#include "simError.h"

void DirectionalAtom::zeroForces() {
  if( hasCoords ){
    frc[offsetX] = 0.0; 
    frc[offsetY] = 0.0; 
    frc[offsetZ] = 0.0;
    
    trq[offsetX] = 0.0; 
    trq[offsetY] = 0.0; 
    trq[offsetZ] = 0.0;
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to zero frc and trq for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}

void DirectionalAtom::setCoords(void){

  if( myConfig->isAllocated() ){

    myConfig->getAtomPointers( index,
                     &pos, 
                     &vel, 
                     &frc, 
                     &trq, 
                     &Amat,
                     &mu,  
                     &ul );
  }
  else{
    sprintf( painCave.errMsg,
             "Attempted to set Atom %d  coordinates with an unallocated "
             "SimState object.\n", index );
    painCave.isFatal = 1;
    simError();
  }

  hasCoords = true;

  *mu = myMu;

}

double DirectionalAtom::getMu( void ) {

  if( hasCoords ){
    return *mu;
  }
  else{
    return myMu;
  }
}

void DirectionalAtom::setMu( double the_mu ) { 

  if( hasCoords ){
    *mu = the_mu;
    myMu = the_mu;
  }
  else{
    myMu = the_mu;
  }
}

void DirectionalAtom::setA( double the_A[3][3] ){

  if( hasCoords ){
    Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
    Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
    Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
    
    this->updateU();  
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to set Amat for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}

void DirectionalAtom::setI( double the_I[3][3] ){
  
  Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
  Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
  Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
}

void DirectionalAtom::setQ( double the_q[4] ){

  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;

  if( hasCoords ){
    q0Sqr = the_q[0] * the_q[0];
    q1Sqr = the_q[1] * the_q[1];
    q2Sqr = the_q[2] * the_q[2];
    q3Sqr = the_q[3] * the_q[3];
    
    
    Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
    Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
    Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
    
    Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
    Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
    Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
    
    Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
    Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
    Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
    
    this->updateU();
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to set Q for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }

}

void DirectionalAtom::getA( double the_A[3][3] ){
  
  if( hasCoords ){
    the_A[0][0] = Amat[Axx];
    the_A[0][1] = Amat[Axy];
    the_A[0][2] = Amat[Axz];
    
    the_A[1][0] = Amat[Ayx];
    the_A[1][1] = Amat[Ayy];
    the_A[1][2] = Amat[Ayz];
    
    the_A[2][0] = Amat[Azx];
    the_A[2][1] = Amat[Azy];
    the_A[2][2] = Amat[Azz];
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to get Amat for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }

}

void DirectionalAtom::printAmatIndex( void ){

  if( hasCoords ){
    std::cerr << "Atom[" << index << "] index =>\n" 
              << "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
              << "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
              << "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to print Amat indices for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}


void DirectionalAtom::getU( double the_u[3] ){
  
  the_u[0] = sux;
  the_u[1] = suy;
  the_u[2] = suz;

  this->body2Lab( the_u );
}

void DirectionalAtom::getQ( double q[4] ){
  
  double t, s;
  double ad1, ad2, ad3;

  if( hasCoords ){
    
    t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
    if( t > 0.0 ){
      
      s = 0.5 / sqrt( t );
      q[0] = 0.25 / s;
      q[1] = (Amat[Ayz] - Amat[Azy]) * s;
      q[2] = (Amat[Azx] - Amat[Axz]) * s;
      q[3] = (Amat[Axy] - Amat[Ayx]) * s;
    }
    else{
      
      ad1 = fabs( Amat[Axx] );
      ad2 = fabs( Amat[Ayy] );
      ad3 = fabs( Amat[Azz] );
      
      if( ad1 >= ad2 && ad1 >= ad3 ){
        
        s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
        q[0] = (Amat[Ayz] + Amat[Azy]) / s;
        q[1] = 0.5 / s;
        q[2] = (Amat[Axy] + Amat[Ayx]) / s;
        q[3] = (Amat[Axz] + Amat[Azx]) / s;
      }
      else if( ad2 >= ad1 && ad2 >= ad3 ){
        
        s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
        q[0] = (Amat[Axz] + Amat[Azx]) / s;
        q[1] = (Amat[Axy] + Amat[Ayx]) / s;
        q[2] = 0.5 / s;
        q[3] = (Amat[Ayz] + Amat[Azy]) / s;
      }
      else{
        
        s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
        q[0] = (Amat[Axy] + Amat[Ayx]) / s;
        q[1] = (Amat[Axz] + Amat[Azx]) / s;
        q[2] = (Amat[Ayz] + Amat[Azy]) / s;
        q[3] = 0.5 / s;
      }
    }
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to get Q for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}


void DirectionalAtom::setEuler( double phi, double theta, double psi ){
  
  if( hasCoords ){
    Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
    Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
    Amat[Axz] = sin(theta) * sin(psi);
    
    Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
    Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
    Amat[Ayz] = sin(theta) * cos(psi);
    
    Amat[Azx] = sin(phi) * sin(theta);
    Amat[Azy] = -cos(phi) * sin(theta);
    Amat[Azz] = cos(theta);
    
    this->updateU();
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to set Euler angles for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}


void DirectionalAtom::lab2Body( double r[3] ){

  double rl[3]; // the lab frame vector 
  
  if( hasCoords ){
    rl[0] = r[0];
    rl[1] = r[1];
    rl[2] = r[2];
    
    r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
    r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
    r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to convert lab2body for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }

}

void DirectionalAtom::body2Lab( double r[3] ){

  double rb[3]; // the body frame vector 
  
  if( hasCoords ){
    rb[0] = r[0];
    rb[1] = r[1];
    rb[2] = r[2];
    
    r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
    r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
    r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to convert body2lab for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}

void DirectionalAtom::updateU( void ){

  if( hasCoords ){
    ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
    ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
    ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to updateU for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}

void DirectionalAtom::getJ( double theJ[3] ){
  
  theJ[0] = jx;
  theJ[1] = jy;
  theJ[2] = jz;
}

void DirectionalAtom::setJ( double theJ[3] ){
  
  jx = theJ[0];
  jy = theJ[1];
  jz = theJ[2];
}

void DirectionalAtom::getTrq( double theT[3] ){
  
  if( hasCoords ){
    theT[0] = trq[offsetX];
    theT[1] = trq[offsetY];
    theT[2] = trq[offsetZ];
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to get Trq for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}

void DirectionalAtom::addTrq( double theT[3] ){
  
  if( hasCoords ){
    trq[offsetX] += theT[0];
    trq[offsetY] += theT[1];
    trq[offsetZ] += theT[2];
  }
  else{
    
    sprintf( painCave.errMsg,
             "Attempt to add Trq for atom %d before coords set.\n",
             index );
    painCave.isFatal = 1;
    simError();
  }
}


void DirectionalAtom::getI( double the_I[3][3] ){
  
  the_I[0][0] = Ixx;
  the_I[0][1] = Ixy;
  the_I[0][2] = Ixz;

  the_I[1][0] = Iyx;
  the_I[1][1] = Iyy;
  the_I[1][2] = Iyz;

  the_I[2][0] = Izx;
  the_I[2][1] = Izy;
  the_I[2][2] = Izz;
}
Revision:	829
Committed:	Tue Oct 28 16:03:37 2003 UTC (20 years, 8 months ago) by gezelter
File size:	8914 byte(s)
Log Message:	replace c++ header stuff with more portable c header stuff Also, mod file fixes and portability changes Some fortran changes will need to be reversed.
#	User	Rev	Content
1	gezelter	829	#include <math.h>
2	mmeineke	377
3			#include "Atom.hpp"
4	mmeineke	670	#include "simError.h"
5	mmeineke	377
6	mmeineke	670	void DirectionalAtom::zeroForces() {
7			if( hasCoords ){
8			frc[offsetX] = 0.0;
9			frc[offsetY] = 0.0;
10			frc[offsetZ] = 0.0;
11
12			trq[offsetX] = 0.0;
13			trq[offsetY] = 0.0;
14			trq[offsetZ] = 0.0;
15			}
16			else{
17
18			sprintf( painCave.errMsg,
19			"Attempt to zero frc and trq for atom %d before coords set.\n",
20			index );
21			painCave.isFatal = 1;
22			simError();
23			}
24			}
25	mmeineke	377
26	tim	689	void DirectionalAtom::setCoords(void){
27	mmeineke	414
28	tim	689	if( myConfig->isAllocated() ){
29
30			myConfig->getAtomPointers( index,
31			&pos,
32			&vel,
33			&frc,
34			&trq,
35			&Amat,
36			&mu,
37			&ul );
38	mmeineke	670	}
39			else{
40			sprintf( painCave.errMsg,
41	tim	689	"Attempted to set Atom %d coordinates with an unallocated "
42	mmeineke	787	"SimState object.\n", index );
43	mmeineke	670	painCave.isFatal = 1;
44			simError();
45			}
46	tim	689
47			hasCoords = true;
48
49	mmeineke	690	*mu = myMu;
50	tim	689
51			}
52
53			double DirectionalAtom::getMu( void ) {
54
55			if( hasCoords ){
56	mmeineke	690	return *mu;
57	tim	689	}
58			else{
59			return myMu;
60			}
61	mmeineke	670	}
62
63			void DirectionalAtom::setMu( double the_mu ) {
64
65			if( hasCoords ){
66	mmeineke	690	*mu = the_mu;
67	tim	689	myMu = the_mu;
68	mmeineke	670	}
69			else{
70	tim	689	myMu = the_mu;
71	mmeineke	670	}
72			}
73
74	mmeineke	377	void DirectionalAtom::setA( double the_A[3][3] ){
75
76	mmeineke	670	if( hasCoords ){
77			Amat[Axx] = the_A[0][0]; Amat[Axy] = the_A[0][1]; Amat[Axz] = the_A[0][2];
78			Amat[Ayx] = the_A[1][0]; Amat[Ayy] = the_A[1][1]; Amat[Ayz] = the_A[1][2];
79			Amat[Azx] = the_A[2][0]; Amat[Azy] = the_A[2][1]; Amat[Azz] = the_A[2][2];
80
81			this->updateU();
82			}
83			else{
84
85			sprintf( painCave.errMsg,
86			"Attempt to set Amat for atom %d before coords set.\n",
87			index );
88			painCave.isFatal = 1;
89			simError();
90			}
91	mmeineke	377	}
92
93			void DirectionalAtom::setI( double the_I[3][3] ){
94
95			Ixx = the_I[0][0]; Ixy = the_I[0][1]; Ixz = the_I[0][2];
96			Iyx = the_I[1][0]; Iyy = the_I[1][1]; Iyz = the_I[1][2];
97			Izx = the_I[2][0]; Izy = the_I[2][1]; Izz = the_I[2][2];
98			}
99
100			void DirectionalAtom::setQ( double the_q[4] ){
101
102			double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
103
104	mmeineke	670	if( hasCoords ){
105			q0Sqr = the_q[0] * the_q[0];
106			q1Sqr = the_q[1] * the_q[1];
107			q2Sqr = the_q[2] * the_q[2];
108			q3Sqr = the_q[3] * the_q[3];
109
110
111			Amat[Axx] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
112			Amat[Axy] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
113			Amat[Axz] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
114
115			Amat[Ayx] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
116			Amat[Ayy] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
117			Amat[Ayz] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
118
119			Amat[Azx] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
120			Amat[Azy] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
121			Amat[Azz] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
122
123			this->updateU();
124			}
125			else{
126
127			sprintf( painCave.errMsg,
128			"Attempt to set Q for atom %d before coords set.\n",
129			index );
130			painCave.isFatal = 1;
131			simError();
132			}
133	mmeineke	377
134			}
135
136			void DirectionalAtom::getA( double the_A[3][3] ){
137
138	mmeineke	670	if( hasCoords ){
139			the_A[0][0] = Amat[Axx];
140			the_A[0][1] = Amat[Axy];
141			the_A[0][2] = Amat[Axz];
142
143			the_A[1][0] = Amat[Ayx];
144			the_A[1][1] = Amat[Ayy];
145			the_A[1][2] = Amat[Ayz];
146
147			the_A[2][0] = Amat[Azx];
148			the_A[2][1] = Amat[Azy];
149			the_A[2][2] = Amat[Azz];
150			}
151			else{
152
153			sprintf( painCave.errMsg,
154			"Attempt to get Amat for atom %d before coords set.\n",
155			index );
156			painCave.isFatal = 1;
157			simError();
158			}
159	mmeineke	377
160			}
161
162	mmeineke	597	void DirectionalAtom::printAmatIndex( void ){
163	mmeineke	377
164	mmeineke	670	if( hasCoords ){
165			std::cerr << "Atom[" << index << "] index =>\n"
166			<< "[ " << Axx << ", " << Axy << ", " << Axz << " ]\n"
167			<< "[ " << Ayx << ", " << Ayy << ", " << Ayz << " ]\n"
168			<< "[ " << Azx << ", " << Azy << ", " << Azz << " ]\n";
169			}
170			else{
171
172			sprintf( painCave.errMsg,
173			"Attempt to print Amat indices for atom %d before coords set.\n",
174			index );
175			painCave.isFatal = 1;
176			simError();
177			}
178	mmeineke	597	}
179
180
181	mmeineke	377	void DirectionalAtom::getU( double the_u[3] ){
182
183			the_u[0] = sux;
184			the_u[1] = suy;
185			the_u[2] = suz;
186
187			this->body2Lab( the_u );
188			}
189
190			void DirectionalAtom::getQ( double q[4] ){
191
192			double t, s;
193			double ad1, ad2, ad3;
194
195	mmeineke	670	if( hasCoords ){
196	mmeineke	377
197	mmeineke	670	t = Amat[Axx] + Amat[Ayy] + Amat[Azz] + 1.0;
198			if( t > 0.0 ){
199	mmeineke	377
200	mmeineke	670	s = 0.5 / sqrt( t );
201			q[0] = 0.25 / s;
202			q[1] = (Amat[Ayz] - Amat[Azy]) * s;
203			q[2] = (Amat[Azx] - Amat[Axz]) * s;
204			q[3] = (Amat[Axy] - Amat[Ayx]) * s;
205	mmeineke	377	}
206			else{
207
208	mmeineke	670	ad1 = fabs( Amat[Axx] );
209			ad2 = fabs( Amat[Ayy] );
210			ad3 = fabs( Amat[Azz] );
211
212			if( ad1 >= ad2 && ad1 >= ad3 ){
213
214			s = 2.0 * sqrt( 1.0 + Amat[Axx] - Amat[Ayy] - Amat[Azz] );
215			q[0] = (Amat[Ayz] + Amat[Azy]) / s;
216			q[1] = 0.5 / s;
217			q[2] = (Amat[Axy] + Amat[Ayx]) / s;
218			q[3] = (Amat[Axz] + Amat[Azx]) / s;
219			}
220			else if( ad2 >= ad1 && ad2 >= ad3 ){
221
222			s = sqrt( 1.0 + Amat[Ayy] - Amat[Axx] - Amat[Azz] ) * 2.0;
223			q[0] = (Amat[Axz] + Amat[Azx]) / s;
224			q[1] = (Amat[Axy] + Amat[Ayx]) / s;
225			q[2] = 0.5 / s;
226			q[3] = (Amat[Ayz] + Amat[Azy]) / s;
227			}
228			else{
229
230			s = sqrt( 1.0 + Amat[Azz] - Amat[Axx] - Amat[Ayy] ) * 2.0;
231			q[0] = (Amat[Axy] + Amat[Ayx]) / s;
232			q[1] = (Amat[Axz] + Amat[Azx]) / s;
233			q[2] = (Amat[Ayz] + Amat[Azy]) / s;
234			q[3] = 0.5 / s;
235			}
236	mmeineke	377	}
237			}
238	mmeineke	670	else{
239
240			sprintf( painCave.errMsg,
241			"Attempt to get Q for atom %d before coords set.\n",
242			index );
243			painCave.isFatal = 1;
244			simError();
245			}
246	mmeineke	377	}
247
248
249			void DirectionalAtom::setEuler( double phi, double theta, double psi ){
250
251	mmeineke	670	if( hasCoords ){
252			Amat[Axx] = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
253			Amat[Axy] = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
254			Amat[Axz] = sin(theta) * sin(psi);
255
256			Amat[Ayx] = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
257			Amat[Ayy] = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
258			Amat[Ayz] = sin(theta) * cos(psi);
259
260			Amat[Azx] = sin(phi) * sin(theta);
261			Amat[Azy] = -cos(phi) * sin(theta);
262			Amat[Azz] = cos(theta);
263
264			this->updateU();
265			}
266			else{
267
268			sprintf( painCave.errMsg,
269			"Attempt to set Euler angles for atom %d before coords set.\n",
270			index );
271			painCave.isFatal = 1;
272			simError();
273			}
274	mmeineke	377	}
275
276
277			void DirectionalAtom::lab2Body( double r[3] ){
278
279			double rl[3]; // the lab frame vector
280
281	mmeineke	670	if( hasCoords ){
282			rl[0] = r[0];
283			rl[1] = r[1];
284			rl[2] = r[2];
285
286			r[0] = (Amat[Axx] * rl[0]) + (Amat[Axy] * rl[1]) + (Amat[Axz] * rl[2]);
287			r[1] = (Amat[Ayx] * rl[0]) + (Amat[Ayy] * rl[1]) + (Amat[Ayz] * rl[2]);
288			r[2] = (Amat[Azx] * rl[0]) + (Amat[Azy] * rl[1]) + (Amat[Azz] * rl[2]);
289			}
290			else{
291
292			sprintf( painCave.errMsg,
293			"Attempt to convert lab2body for atom %d before coords set.\n",
294			index );
295			painCave.isFatal = 1;
296			simError();
297			}
298	mmeineke	377
299			}
300
301			void DirectionalAtom::body2Lab( double r[3] ){
302
303			double rb[3]; // the body frame vector
304
305	mmeineke	670	if( hasCoords ){
306			rb[0] = r[0];
307			rb[1] = r[1];
308			rb[2] = r[2];
309
310			r[0] = (Amat[Axx] * rb[0]) + (Amat[Ayx] * rb[1]) + (Amat[Azx] * rb[2]);
311			r[1] = (Amat[Axy] * rb[0]) + (Amat[Ayy] * rb[1]) + (Amat[Azy] * rb[2]);
312			r[2] = (Amat[Axz] * rb[0]) + (Amat[Ayz] * rb[1]) + (Amat[Azz] * rb[2]);
313			}
314			else{
315
316			sprintf( painCave.errMsg,
317			"Attempt to convert body2lab for atom %d before coords set.\n",
318			index );
319			painCave.isFatal = 1;
320			simError();
321			}
322	mmeineke	377	}
323
324			void DirectionalAtom::updateU( void ){
325
326	mmeineke	670	if( hasCoords ){
327			ul[offsetX] = (Amat[Axx] * sux) + (Amat[Ayx] * suy) + (Amat[Azx] * suz);
328			ul[offsetY] = (Amat[Axy] * sux) + (Amat[Ayy] * suy) + (Amat[Azy] * suz);
329			ul[offsetZ] = (Amat[Axz] * sux) + (Amat[Ayz] * suy) + (Amat[Azz] * suz);
330			}
331			else{
332
333			sprintf( painCave.errMsg,
334			"Attempt to updateU for atom %d before coords set.\n",
335			index );
336			painCave.isFatal = 1;
337			simError();
338			}
339	mmeineke	377	}
340
341	mmeineke	599	void DirectionalAtom::getJ( double theJ[3] ){
342
343			theJ[0] = jx;
344			theJ[1] = jy;
345			theJ[2] = jz;
346			}
347
348			void DirectionalAtom::setJ( double theJ[3] ){
349
350			jx = theJ[0];
351			jy = theJ[1];
352			jz = theJ[2];
353			}
354
355			void DirectionalAtom::getTrq( double theT[3] ){
356
357	mmeineke	670	if( hasCoords ){
358			theT[0] = trq[offsetX];
359			theT[1] = trq[offsetY];
360			theT[2] = trq[offsetZ];
361			}
362			else{
363
364			sprintf( painCave.errMsg,
365			"Attempt to get Trq for atom %d before coords set.\n",
366			index );
367			painCave.isFatal = 1;
368			simError();
369			}
370	mmeineke	599	}
371
372			void DirectionalAtom::addTrq( double theT[3] ){
373
374	mmeineke	670	if( hasCoords ){
375			trq[offsetX] += theT[0];
376			trq[offsetY] += theT[1];
377			trq[offsetZ] += theT[2];
378			}
379			else{
380
381			sprintf( painCave.errMsg,
382			"Attempt to add Trq for atom %d before coords set.\n",
383			index );
384			painCave.isFatal = 1;
385			simError();
386			}
387	mmeineke	599	}
388
389
390			void DirectionalAtom::getI( double the_I[3][3] ){
391
392			the_I[0][0] = Ixx;
393			the_I[0][1] = Ixy;
394			the_I[0][2] = Ixz;
395
396			the_I[1][0] = Iyx;
397			the_I[1][1] = Iyy;
398			the_I[1][2] = Iyz;
399
400			the_I[2][0] = Izx;
401			the_I[2][1] = Izy;
402			the_I[2][2] = Izz;
403			}