| 297 |
|
vel =integrableObject->getVel(); |
| 298 |
|
mass = integrableObject->getMass(); |
| 299 |
|
if (integrableObject->isDirectional()){ |
| 300 |
– |
//calculate angular velocity in lab frame |
| 300 |
|
Mat3x3d I = integrableObject->getI(); |
| 301 |
|
Vector3d angMom = integrableObject->getJ(); |
| 302 |
< |
Vector3d omega; |
| 302 |
> |
A = integrableObject->getA(); |
| 303 |
> |
Atrans = A.transpose(); |
| 304 |
> |
|
| 305 |
> |
Vector3d omegaBody; |
| 306 |
|
|
| 307 |
|
if (integrableObject->isLinear()) { |
| 308 |
|
int linearAxis = integrableObject->linearAxis(); |
| 309 |
|
int l = (linearAxis +1 )%3; |
| 310 |
|
int m = (linearAxis +2 )%3; |
| 311 |
< |
omega[l] = angMom[l] /I(l, l); |
| 312 |
< |
omega[m] = angMom[m] /I(m, m); |
| 311 |
> |
omegaBody[l] = angMom[l] /I(l, l); |
| 312 |
> |
omegaBody[m] = angMom[m] /I(m, m); |
| 313 |
|
|
| 314 |
|
} else { |
| 315 |
< |
omega[0] = angMom[0] /I(0, 0); |
| 316 |
< |
omega[1] = angMom[1] /I(1, 1); |
| 317 |
< |
omega[2] = angMom[2] /I(2, 2); |
| 315 |
> |
omegaBody[0] = angMom[0] /I(0, 0); |
| 316 |
> |
omegaBody[1] = angMom[1] /I(1, 1); |
| 317 |
> |
omegaBody[2] = angMom[2] /I(2, 2); |
| 318 |
|
} |
| 319 |
|
|
| 320 |
< |
//std::cerr << "I = " << I(0,0) << "\t" << I(1,1) << "\t" << I(2,2) << "\n\n"; |
| 320 |
> |
Vector3d omegaLab = Atrans * omegaBody; |
| 321 |
|
|
| 322 |
< |
//apply friction force and torque at center of resistance |
| 321 |
< |
A = integrableObject->getA(); |
| 322 |
< |
Atrans = A.transpose(); |
| 323 |
< |
//std::cerr << "A = " << integrableObject->getA() << "\n"; |
| 324 |
< |
//std::cerr << "Atrans = " << A.transpose() << "\n\n"; |
| 325 |
< |
Vector3d rcr = Atrans * hydroProps_[index]->getCOR(); |
| 326 |
< |
//std::cerr << "cor = " << hydroProps_[index]->getCOR() << "\n\n\n\n"; |
| 327 |
< |
//std::cerr << "rcr = " << rcr << "\n\n"; |
| 328 |
< |
Vector3d vcdLab = vel + cross(omega, rcr); |
| 329 |
< |
|
| 330 |
< |
//std::cerr << "velL = " << vel << "\n\n"; |
| 331 |
< |
//std::cerr << "vcdL = " << vcdLab << "\n\n"; |
| 332 |
< |
Vector3d vcdBody = A* vcdLab; |
| 333 |
< |
//std::cerr << "vcdB = " << vcdBody << "\n\n"; |
| 334 |
< |
Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega); |
| 322 |
> |
// apply friction force and torque at center of resistance |
| 323 |
|
|
| 324 |
< |
//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n"; |
| 325 |
< |
//std::cerr << "ffB = " << frictionForceBody << "\n\n"; |
| 326 |
< |
Vector3d frictionForceLab = Atrans*frictionForceBody; |
| 327 |
< |
//std::cerr << "ffL = " << frictionForceLab << "\n\n"; |
| 328 |
< |
//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; |
| 324 |
> |
Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR(); |
| 325 |
> |
Vector3d vcdLab = vel + cross(omegaLab, rcrLab); |
| 326 |
> |
|
| 327 |
> |
Vector3d vcdBody = A * vcdLab; |
| 328 |
> |
Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody); |
| 329 |
> |
|
| 330 |
> |
Vector3d frictionForceLab = Atrans * frictionForceBody; |
| 331 |
|
integrableObject->addFrc(frictionForceLab); |
| 332 |
< |
//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; |
| 333 |
< |
//std::cerr << "ome = " << omega << "\n\n"; |
| 334 |
< |
Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega); |
| 345 |
< |
//std::cerr << "ftB = " << frictionTorqueBody << "\n\n"; |
| 346 |
< |
Vector3d frictionTorqueLab = Atrans*frictionTorqueBody; |
| 347 |
< |
//std::cerr << "ftL = " << frictionTorqueLab << "\n\n"; |
| 348 |
< |
//std::cerr << "ftL2 = " << frictionTorqueLab+cross(rcr,frictionForceLab) << "\n\n"; |
| 349 |
< |
//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; |
| 350 |
< |
integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab)); |
| 351 |
< |
//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; |
| 332 |
> |
Vector3d frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody); |
| 333 |
> |
Vector3d frictionTorqueLab = Atrans * frictionTorqueBody; |
| 334 |
> |
integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab)); |
| 335 |
|
|
| 336 |
|
//apply random force and torque at center of resistance |
| 337 |
|
Vector3d randomForceBody; |
| 338 |
|
Vector3d randomTorqueBody; |
| 339 |
|
genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_); |
| 340 |
< |
//std::cerr << "rfB = " << randomForceBody << "\n\n"; |
| 341 |
< |
//std::cerr << "rtB = " << randomTorqueBody << "\n\n"; |
| 359 |
< |
Vector3d randomForceLab = Atrans*randomForceBody; |
| 360 |
< |
Vector3d randomTorqueLab = Atrans* randomTorqueBody; |
| 340 |
> |
Vector3d randomForceLab = Atrans * randomForceBody; |
| 341 |
> |
Vector3d randomTorqueLab = Atrans * randomTorqueBody; |
| 342 |
|
integrableObject->addFrc(randomForceLab); |
| 343 |
< |
//std::cerr << "rfL = " << randomForceLab << "\n\n"; |
| 363 |
< |
//std::cerr << "rtL = " << randomTorqueLab << "\n\n"; |
| 364 |
< |
//std::cerr << "rtL2 = " << randomTorqueLab + cross(rcr, randomForceLab) << "\n\n"; |
| 365 |
< |
integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab )); |
| 343 |
> |
integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab )); |
| 344 |
|
|
| 345 |
|
} else { |
| 346 |
|
//spherical atom |
| 347 |
|
Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel); |
| 370 |
– |
//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n"; |
| 348 |
|
Vector3d randomForce; |
| 349 |
|
Vector3d randomTorque; |
| 350 |
|
genRandomForceAndTorque(randomForce, randomTorque, index, variance_); |
| 376 |
|
Z[0] = randNumGen_.randNorm(0, variance); |
| 377 |
|
Z[1] = randNumGen_.randNorm(0, variance); |
| 378 |
|
Z[2] = randNumGen_.randNorm(0, variance); |
| 402 |
– |
//Z[3] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
| 403 |
– |
//Z[4] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
| 404 |
– |
//Z[5] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
| 379 |
|
Z[3] = randNumGen_.randNorm(0, variance); |
| 380 |
|
Z[4] = randNumGen_.randNorm(0, variance); |
| 381 |
|
Z[5] = randNumGen_.randNorm(0, variance); |
| 382 |
|
|
| 409 |
– |
|
| 383 |
|
generalForce = hydroProps_[index]->getS()*Z; |
| 384 |
|
|
| 385 |
|
force[0] = generalForce[0]; |
