| 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 ){ |
| 60 |
|
else{ |
| 61 |
|
return myMu; |
| 62 |
|
} |
| 61 |
– |
return 0; |
| 63 |
|
} |
| 64 |
|
|
| 65 |
|
void DirectionalAtom::setMu( double the_mu ) { |
| 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 |
+ |
} |
