| 1 |
kstocke1 |
3926 |
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 |
|
|
|
