| 1 |
vc = Pc / Mc; |
| 2 |
ac = -momentumTarget_ / Mc + vc; |
| 3 |
acrec = -momentumTarget_ / Mc; |
| 4 |
|
| 5 |
// We now need the inverse of the inertia tensor to calculate the angular velocity of the cold slab; |
| 6 |
|
| 7 |
Ici = Ic.inverse(); |
| 8 |
omegac = Ici * Lc; |
| 9 |
bc = -(Ici * angularMomentumTarget_) + omegac; |
| 10 |
bcrec = bc - omegac; |
| 11 |
|
| 12 |
cNumerator = Kc - kineticTarget_; |
| 13 |
if (doLinearPart) |
| 14 |
cNumerator -= 0.5 * Mc * ac.lengthSquare(); |
| 15 |
|
| 16 |
if (doAngularPart) |
| 17 |
cNumerator -= 0.5 * ( dot(bc, Ic * bc)); |
| 18 |
|
| 19 |
cDenominator = Kc; |
| 20 |
|
| 21 |
if (doLinearPart) |
| 22 |
cDenominator -= 0.5 * Mc * vc.lengthSquare(); |
| 23 |
|
| 24 |
if (doAngularPart) |
| 25 |
cDenominator -= 0.5*(dot(omegac, Ic * omegac)); |
| 26 |
|
| 27 |
c = sqrt(cNumerator / cDenominator); |
| 28 |
|
| 29 |
if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients |
| 30 |
|
| 31 |
vh = Ph / Mh; |
| 32 |
ah = momentumTarget_ / Mh + vh; |
| 33 |
ahrec = momentumTarget_ / Mh; |
| 34 |
|
| 35 |
// We now need the inverse of the inertia tensor to calculate the angular velocity of the hot slab; |
| 36 |
|
| 37 |
Ihi = Ih.inverse(); |
| 38 |
omegah = Ihi * Lh; |
| 39 |
bh = (Ihi * angularMomentumTarget_) + omegah; |
| 40 |
bhrec = bh - omegah; |
| 41 |
|
| 42 |
hNumerator = Kh + kineticTarget_; |
| 43 |
if (doLinearPart) |
| 44 |
hNumerator -= 0.5 * Mh * ah.lengthSquare(); |
| 45 |
|
| 46 |
if (doAngularPart) |
| 47 |
hNumerator -= 0.5 * ( dot(bh, Ih * bh)); |
| 48 |
|
| 49 |
hDenominator = Kh; |
| 50 |
if (doLinearPart) |
| 51 |
hDenominator -= 0.5 * Mh * vh.lengthSquare(); |
| 52 |
if (doAngularPart) |
| 53 |
hDenominator -= 0.5*(dot(omegah, Ih * omegah)); |
| 54 |
|
| 55 |
h = sqrt(hNumerator / hDenominator); |
| 56 |
if ((h > 0.9) && (h < 1.1)) {//restrict scaling coefficients |
| 57 |
|
| 58 |
rPos = (*sdi)->getPos() - coordinateOrigin_; |
| 59 |
if (doLinearPart) |
| 60 |
vel = ((*sdi)->getVel() - vc) * c + ac; |
| 61 |
if (doAngularPart) |
| 62 |
vel = ((*sdi)->getVel() - cross(omegac, rPos)) * c + cross(bc, rPos); |
| 63 |
|
