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