(* by Nasser M. Abbasi, version: July 26,2013 *)
Manipulate[
tick;
Module[{sol, t, \[Theta]1sol, \[Theta]2sol, v1sol, v2sol, eq1, eq2,
ic, \[Theta]1, \[Theta]2, result},
If[state == "running" || state == "step" || state == "initial",
eq1 = 3/2 r^2 m1 \[Theta]1''[t] +
r^2 ((k1 + k2) \[Theta]1[t] - k2 \[Theta]2[t]) == 0;
eq2 = 3/2 r^2 m2 \[Theta]2''[t] +
r^2 (-k2 \[Theta]1[t] + (k2 + k3) \[Theta]2[t]) == 0;
ic = {\[Theta]1[0] == \[Theta]1current, \[Theta]2[
0] == \[Theta]2current, \[Theta]1'[0] ==
v1current, \[Theta]2'[0] == v2current};
sol = First@
NDSolve[{eq1, eq2,
ic}, {\[Theta]1[t], \[Theta]2[t], \[Theta]1'[t], \[Theta]2'[
t]}, {t, 0, stepSize}]
];
If[state == "running" || state == "step",
\[Theta]1current = \[Theta]1[t] /. sol /. t -> stepSize;
\[Theta]2current = \[Theta]2[t] /. sol /. t -> stepSize;
v1current = \[Theta]1'[t] /. sol /. t -> stepSize;
v2current = \[Theta]2'[t] /. sol /. t -> stepSize;
x1 = springRelaxedL + r \[Theta]1[t] /. sol /. t -> 0;
x2 = x1 + springRelaxedL + r \[Theta]2[t] /. sol /. t -> 0;
currentTime = currentTime + stepSize;
If[state == "running",
tick += del
];
];
Text@drawSystem[currentTime, r, x1,
x2, \[Theta]1current, \[Theta]2current, v1current, v2current, m1,
m2, k1, k2, k3, showGrid]
],
Text@Grid[{
{
Grid[{
{
Button[Text[Style["run", 11]], state = "running"; tick += del,
ImageSize -> {50, 35}],
Button[Text[Style["pause", 11]], state = "paused";
tick += del, ImageSize -> {50, 35}]
},
{Button[Text[Style["step", 11]], state = "step"; tick += del,
ImageSize -> {50, 35}],
Button[Text[Style["reset", 11]], state = "initial";
tick = 0; \[Theta]10 = 0; \[Theta]1current = 0; \[Theta]20 =
0; \[Theta]2current = 0; v1current = 40.0*(2 Pi)/60;
v10 = 40.0;
v2current = 0; v20 = 0; currentTime = 0; k1 = 100; k2 = 200;
k3 = 300; m1 = 6; m2 = 1; x1 = springRelaxedL;
x2 = 2 springRelaxedL; tick += del,
ImageSize -> {50, 35}
]
}
}, Spacings -> {0.4, .2}, Alignment -> Center
]
},
{
Grid[{
{
"step size",
Manipulator[
Dynamic[stepSize, {stepSize = #} &], {0.001, 0.02, 0.001},
ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[stepSize, {4, 3}]], "sec"
},
{Row[{Subscript["\[Theta]", 1]', "(0)"}],
Manipulator[
Dynamic[v10, {v10 = #; v1current = v10*2 Pi/60} &], {0,
50, .1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[v10, {2, 1}]], "rpm"
},
{Row[{Subscript["\[Theta]", 2]', "(0)"}],
Manipulator[
Dynamic[v20 , {v20 = #; v2current = v20*2 Pi/60} &], {0,
50, .1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[v20, {2, 1}]], "rpm"
},
{Text@TraditionalForm@Style[Subscript[k, 1], 11],
Manipulator[Dynamic[k1 , {k1 = #; tick += del} &],
{100, 500, 1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[k1, 3]], "N/m"
},
{TraditionalForm@Style[Subscript[k, 2], 11],
Manipulator[Dynamic[k2 , {k2 = #; tick += del} &],
{100, 500, 1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[k2, 3]], "N/m"
},
{TraditionalForm@Style[Subscript[k, 3], 11],
Manipulator[Dynamic[k3 , {k3 = #; tick += del} &],
{100, 500, 1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[k3, 3]], "N/m"
},
{TraditionalForm@Style[Subscript[m, 1], 11],
Manipulator[Dynamic[m1 , {m1 = #; tick += del} &],
{.1, 10, .1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[m1, {3, 1}]], "kg"
},
{TraditionalForm@Style[Subscript[m, 2], 11],
Manipulator[Dynamic[m2 , {m2 = #; tick += del} &],
{.1, 10, .1}, ImageSize -> Tiny, ContinuousAction -> True],
Dynamic[padIt2[m2, {3, 1}]], "kg"
}
}, Frame -> True,
FrameStyle -> Directive[Thickness[.001], Gray],
Alignment -> Center, Spacings -> {.4, .1}
]
}
}],
Row[{"show grid", Spacer[4],
Checkbox[Dynamic[showGrid, {showGrid = #; tick += del} &]]}],
{{showGrid, True}, None},
{{springRelaxedL, 2}, None},
{{x1, 2}, None},(*keep same as springRelaxedL*)
{{x2, 4}, None},(*keep same as 2*springRelaxedL*)
{{r, .5}, None},
{{state, "initial"}, None},
{{stepSize, 0.015}, None},
{{m1, 6}, None},
{{m2, 1}, None},
{{k1, 100}, None},
{{k3, 300}, None},
{{k2, 200}, None},
{{\[Theta]10, 0}, None}, (*initial radial position m1*)
{{\[Theta]20, 0}, None}, (*initial radial position m2*)
{{v10, 40}, None}, (*initial angular position RPM, m2*)
{{v20, 0}, None}, (*initial angula position RPM, m2*)
{{\[Theta]1current, 0}, None}, (*current radial position m1*)
{{\[Theta]2current, 0}, None}, (*current radial position m2*)
{{v1current, 40 (2 Pi/60)},
None}, (*current angular position, rad/sec m2*)
{{v2current, 0 (2 Pi/60)},
None}, (*current angula position rad/sec m2*)
{{currentTime, 0}, None}, (*simulation time*)
{{del, $MachineEpsilon}, None},
{{tick, 0}, None},
ControlPlacement -> Left,
SynchronousUpdating -> False,
SynchronousInitialization -> False,
ContinuousAction -> False,
Alignment -> Center,
ImageMargins -> 0,
FrameMargins -> 0,
Paneled -> True,
Frame -> False,
AutorunSequencing -> {1},
TrackedSymbols :> {tick},
Initialization :>
{
integerStrictPositive = (IntegerQ[#] && # > 0 &);
integerPositive = (IntegerQ[#] && # >= 0 &);
numericStrictPositive = (Element[#, Reals] && # > 0 &);
numericPositive = (Element[#, Reals] && # >= 0 &);
numericStrictNegative = (Element[#, Reals] && # < 0 &);
numericNegative = (Element[#, Reals] && # <= 0 &);
bool = (Element[#, Booleans] &);
numeric = (Element[#, Reals] &);
integer = (Element[#, Integers] &);
(*--------------------------------------------*)
(* helper function for formatting *)
(*--------------------------------------------*)
padIt1[v_?numeric, f_List] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"-", "+"},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*--------------------------------------------*)
(* helper function for formatting *)
(*--------------------------------------------*)
padIt2[v_?numeric, f_List] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"", ""},
NumberPadding -> {"0", "0"}, SignPadding -> True];
padIt2[v_?numeric, f_Integer] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"", ""},
NumberPadding -> {"0", "0"}, SignPadding -> True];
(*--------------------------------------------*)
getKE[m1_, m2_, v1_, v2_] :=
Module[{i1 = 1/2 m1 r^2, i2 = 1/2 m2 r^2},
1/2 i1 (v1)^2 + 1/2 m1 (r v1)^2 + 1/2 i2 (v2)^2 +
1/2 m2 (r v2)^2];
getPE[\[Theta]1_, \[Theta]2_, k1_, k2_, k3_] :=
1/2 k1 (r \[Theta]1)^2 + 1/2 k2 (r \[Theta]2 - r \[Theta]1)^2 +
1/2 k3 (r \[Theta]2)^2;
(*--------------------------------------------*)
drawSystem[currentTime_,
r_(*radius of sphere*),
x1_(*distance of center of 1st cylinder from left wall*),
x2_(*distance of center of 2nd cylinder from left wall*),
\[Theta]1sol_(*angle of rotation of first cyliner*),
\[Theta]2sol_(*angle of rotation of second cyliner*),
v1sol_(*anglular velocity of first cyliner*),
v2sol_(*anglular velocity of first cyliner*),
m1_, m2_, k1_, k2_, k3_, showGrid_] := Module[{ke, pe, h1},
ke = getKE[m1, m2, v1sol, v2sol];
pe = getPE[\[Theta]1sol, \[Theta]2sol, k1, k2, k3];
h1 = Text@Style[Grid[{
{"time (sec)",
"P.E. (J)",
"K.E. (J)",
"energy (J)",
Row[{Subscript["\[Theta]", 1]', " (rpm)"}],
Row[{Subscript["\[Theta]", 2]', " (rpm)"}]
},
{padIt2[currentTime, {6, 2}],
padIt2[pe, {6, 3}],
padIt2[ke, {6, 3}],
padIt2[ke + pe, {6, 3}],
padIt2[v1sol*60/(2 Pi), {6, 3}],
padIt2[v2sol*60/(2 Pi), {6, 3}]}
},
Frame -> All,
FrameStyle -> Directive[Thickness[.001], Gray],
Spacings -> {1, 1.2},
Alignment -> Center]
, 12];
Grid[{
{
h1
},
{
Graphics[{
{AbsoluteThickness[k1/250],
Line@makeSpring[0, r, x1 - r, r, r/2]},
{EdgeForm[Black], Opacity[m1/10], Red, Disk[{x1, r}, r]},
{Black,
Line[{{x1, r}, {x1 + r Sin[\[Theta]1sol],
r + r Cos[\[Theta]1sol]}}]},
{AbsoluteThickness[k2/250],
Line@makeSpring[x1 + r, r, x2 - r, r, r/2]},
{EdgeForm[Black], Opacity[m2/10], Red, Disk[{x2, r}, r]},
{Black,
Line[{{x2, r}, {x2 + r Sin[\[Theta]2sol],
r + r Cos[\[Theta]2sol]}}]},
{AbsoluteThickness[k3/250],
Line@makeSpring[x2 + r, r, 3 springRelaxedL, r, r/2]}
},
PlotRange -> {{0, 3 springRelaxedL}, {-0.5 r, 3 r}},
ImageSize -> 400,
Axes -> True,
Frame -> True,
AxesOrigin -> {0, 0},
AspectRatio -> Automatic,
GridLines ->
If[showGrid, {Range[0, 3 springRelaxedL, .5], Automatic},
None],
GridLinesStyle -> LightGray
]
}
}, Alignment -> Center, Spacings -> {.1, .2}
]
];
(*--------------------------------------------*)
makeSpring[xFirst_?numeric, yFirst_?numeric, xEnd_?numeric,
yEnd_?numeric, szel_?numeric] :=
Module[{hx, veghossz, hossz, hy, dh, tbl},
hx = xEnd - xFirst;
If[Abs[hx] <= $MachineEpsilon, hx = 10^-6];
hy = yEnd - yFirst;
If[Abs[hy] <= $MachineEpsilon, hy = 10^-6];
veghossz = 0.03;
hossz = Sqrt[hx^2 + hy^2];
dh = (hossz - 2*veghossz)/20;
tbl =
Table[If[
OddQ[i], {xFirst + hx*(i*dh + veghossz)/hossz + hy*szel/hossz,
yFirst + hy*(i*dh + veghossz)/hossz -
hx*szel/hossz}, {xFirst + hx*(i*dh + veghossz)/hossz -
hy*szel/hossz,
yFirst + hy*(i*dh + veghossz)/hossz + hx*szel/hossz}], {i, 2,
18}];
{{xFirst, yFirst}}~
Join~{{xFirst + hx*(dh + veghossz)/hossz,
yFirst + hy*(dh + veghossz)/hossz}}~Join~tbl~
Join~{{xFirst + hx*(19*dh + veghossz)/hossz,
yFirst + hy*(19*dh + veghossz)/hossz}}~Join~{{xEnd, yEnd}}
];
}
]