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 |
|