(*By Nasser M. Abbasi, version: January 19, 2015 *)
Manipulate[
gTick;
Module[{g1, $x, $xv, $y, $yv, $\[Theta], $\[Omega], $trace, cm, initialYEffective, g = 9.81, L0},
L0 = initialRelaxedLength[len0, m0, M0, k1, g];
initialYEffective = initialY + L0;
If[setIC,
setIC = False;
cm = calculateCenterOfMass[m0, M0, initialYEffective, initialTheta Degree, initialBobX ];
list@setIC[{initialBobX, initialBobSpeed, initialY, initialYSpeed, initialTheta Degree, initialThetaDot, cm, r0, L0}];
isRelaxedSpring = If[Abs[initialBobX] <= $MachineEpsilon, True, False];
isRelaxedTopSpring = If[Abs[initialY - L0] <= $MachineEpsilon, True, False]
];
If[runningState == "RUNNING" || runningState == "STEP",
solver@step[m0, len0, k1, M0, delT];
If[runningState == "STEP", runningState = "STOP", gTick += del]
];
g1 = display@makeDisplay[plotType, m0, M0, len0, r0, k1, showPlots, showCounters, traceCM, traceBob, w];
{simTime, $x, $xv, $y, $yv, $\[Theta], $\[Omega], cm, $trace} = list@getCurrent[];
FinishDynamic[];
Text@g1
],
(*---------- control layout ------------*)
Item[Grid[{
{Button[Style["run", 12], {If[Not[plotType == "velocity/acceleration"], runningState = "RUNNING"]; gTick += del},
ImageSize -> {55, 35}]},
{Button[Style["stop", 12], {If[Not[plotType == "velocity/acceleration"], runningState = "STOP"]; gTick += del},
ImageSize -> {55, 35}]},
{Button[Style["step", 12], {If[Not[plotType == "velocity/acceleration"], runningState = "STEP"]; gTick += del},
ImageSize -> {55, 35}]},
{Button[Style["reset", 12],(*bring simulation back to initial conditions*)
{
setIC = True;
runningState = "STOP";
gTick += del
}, ImageSize -> {55, 35}]
},
{Grid[{
{Style["time (sec)", 11]},
{Style[Dynamic@padIt2[Mod[simTime, 1000], {7, 4}], 11]}
}, Spacings -> {0, 0.5}]
}
}, Spacings -> {0.1, .7}
], ControlPlacement -> Right],
Grid[{
{Grid[
{
{Row[{"speed", Spacer[15], "(slow)"}],
Manipulator[Dynamic[delT, {delT = #} &], {0.001, 0.05, 0.001}, ImageSize -> Tiny, ContinuousAction -> False],
"(fast)"
},
{Row[{"display zoom", Spacer[10], "(in)"}],
Manipulator[Dynamic[w,
{w = Chop@#; gTick += del} &], {0.05, 3, 0.01}, ImageSize -> Tiny, ContinuousAction -> False],
"(out)"
}
}, Spacings -> {.6, .5}, Alignment -> Left, Frame -> True, FrameStyle -> Directive[Thickness[.005], Gray]
]
},
{Style["top spring initial conditions", 12]},
{Grid[{
{"spring mass",
Manipulator[Dynamic[m0,
{m0 = Chop@#; setIC = True; gTick += del} &], {1, 5, 0.1}, ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt2[m0, {2, 2}], 11]
},
{"spring stiffness",
Manipulator[Dynamic[k1,
{k1 = Chop@#; setIC = True; gTick += del} &], {300, 10000, 10}, ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt2[k1, 4], 11]
},
{"angular velocity",
Manipulator[Dynamic[initialThetaDot,
{initialThetaDot = #; setIC = True; gTick += del} &], {-0.5, 0.5, .1}, ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt1[initialThetaDot, {2, 2}], 11]
},
{"angle (degree)",
Manipulator[Dynamic[initialTheta,
{initialTheta = #; setIC = True; gTick += del} &], {-45, 45, 1}, ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt1[initialTheta, 3], 11]
},
{"initial extension",
Manipulator[Dynamic[initialY,
{(initialY = #) &, (initialY = Chop[#]; setIC = True; gTick += del) &}], {-0.25, 0.25, 0.01}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt1[initialY, {2, 2}], 11]
}
}, Spacings -> {.6, .3}, Alignment -> Left, Frame -> True, FrameStyle -> Directive[Thickness[.005], Gray]
]
},
{Style["bottom spring initial conditions", 12]},
{Grid[{
{"spring stiffness",
Manipulator[Dynamic[k, {k = Chop@#; setIC = True; gTick += del} &], {300, 3000, 10}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[k, 4], 11]
},
{"bob mass",
Manipulator[Dynamic[M0, {M0 = Chop@#; setIC = True; gTick += del} &], {0, 50, 1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[M0, 2], 11]
},
{"bob initial position",
Manipulator[Dynamic[initialBobX, {initialBobX = #; setIC = True; gTick += del} &], {-0.5, 0.5, 0.1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt1[initialBobX, {1, 1}], 11]
},
{
Grid[{{
Row[{"relax", Spacer[2],
Checkbox[Dynamic[isRelaxedSpring, {isRelaxedSpring = #; initialBobX = 0; setIC = True; gTick += del} &]]}], Spacer[9],
Row[{"trace bob", Spacer[2],
Checkbox[Dynamic[traceBob, {traceBob = #; gTick += del} &], Enabled -> Dynamic[M0 > 0]]}], Spacer[9],
Row[{"show CM", Spacer[4],
Checkbox[Dynamic[traceCM, {traceCM = #; gTick += del} &]]}]
}}, Spacings -> {0.4, 0}], SpanFromLeft
},
{"bob initial speed",
Manipulator[Dynamic[initialBobSpeed, {initialBobSpeed = Chop@#; setIC = True; gTick += del} &], {-0.1, 0.1, 0.01},
ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt1[initialBobSpeed, {2, 2}], 11]
}
}, Spacings -> {.6, .3}, Alignment -> Left, Frame -> True, FrameStyle -> Directive[Thickness[.005], Gray]
]
},
{Grid[{
{Style["plot type", 12], Style["test case", 12]},
{
PopupMenu[Dynamic[plotType, {plotType = #;
If[plotType == "velocity/acceleration",
If[runningState == "RUNNING",
runningState = "SUSPEND_RUNNING"
],
If[runningState == "SUSPEND_RUNNING", runningState = "RUNNING"]
]; gTick += del} &],
{ "disk angular velocity" -> Style["disk angular velocity", 12],
"disk angle" -> Style["disk angle", 12],
"bob position" -> Style["bob position", 12],
"bob velocity" -> Style["bob velocity", 12],
"vertical displacement" -> Style["vertical displacement", 12]
}, ImageSize -> All, ContinuousAction -> False, Enabled -> Dynamic[showPlots]]
,
PopupMenu[Dynamic[testCase, {testCase = #;
runningState = "STOP";
Which[testCase == 1,
(isRelaxedSpring = False; delT = 0.006; m0 = 1; k = 3000; k1 = 10000; w = 1.2;
initialThetaDot = 0; initialTheta = 0; initialBobX = -0.5; initialBobSpeed = 0; M0 = 50; traceCM = False;
traceBob = False; initialY = 0.25; setIC = True; runningState = "STOP"; plotType = "bob position"
),
testCase == 2,
(isRelaxedSpring = False; delT = 0.05; m0 = 5; k = 500; w = 0.9; k1 = 10000;
initialThetaDot = 0.0; initialTheta = 15; initialBobX = -.2; initialBobSpeed = 0; M0 = 0; traceCM = False;
traceBob = False; setIC = True; runningState = "STOP"; plotType = "disk angular velocity"
),
testCase == 3,
(isRelaxedSpring = True; delT = 0.03; m0 = 5; k = 3000; k1 = 10000; M0 = 5; traceCM = False; traceBob = False; w = 1.1;
initialThetaDot = 0; initialTheta = 0; initialBobX = 0; initialBobSpeed = 0; setIC = True; runningState = "STOP";
plotType = "disk angular velocity"
)
];
gTick += del} &],
{
1 -> Style["1", 12], 2 -> Style["2", 12], 3 -> Style["3", 12]
}, ImageSize -> All, ContinuousAction -> False]
},
{Grid[{
{
Row[{"show plot ", Spacer[4], Checkbox[Dynamic[showPlots, {showPlots = #; gTick += del} &]]}], Spacer[12],
Row[{"show counters ", Spacer[4], Checkbox[Dynamic[showCounters, {showCounters = #; gTick += del} &]]}], SpanFromLeft
}
}
], SpanFromLeft
},
{Grid[{
{"plot window size",
Manipulator[Dynamic[xScale,
{xScale = #; list@setSize[xScale]; setIC = True; gTick += del} &], {10, 1000, 1}, ImageSize -> Tiny,
ContinuousAction -> False, Enabled -> Dynamic[showPlots]],
Style[Dynamic@padIt2[xScale, 4], 11]
}
}], SpanFromLeft
}
}, Frame -> True, FrameStyle -> Directive[Thickness[.005], Gray]
]
}
}, Spacings -> {0, {2 -> 0.7, 4 -> 0.7, 6 -> 0.5, 8 -> 0}}, Alignment -> Center, Frame -> None],
(*----control variables--*)
{{m0, 1}, None},(*mass of top spring*)
{{M0, 10}, None},
{{len0, 0.6}, None},(*relaxed length, fixed does not change*)
{{delT, 0.01}, None},
{{k, 500}, None},
{{k1, 500}, None},
{{initialThetaDot, 0.0}, None},
{{initialTheta, 15}, None},
{{initialBobX, -0.35}, None},
{{initialY, .1}, None},
{{initialYDot, 0}, None},
{{initialBobSpeed, 0}, None},
{{initialYSpeed, 0}, None},
{{isRelaxedSpring, False}, None},
{{isRelaxedTopSpring, False}, None},
{{plotType, "bob position"}, None},
{{testCase, 1}, None},
{{showPlots, False}, None},
{{showCounters, True}, None},
{{xScale, 400}, None},
{{traceBob, False}, None},
{{traceCM, False}, None},
{{setIC, False}, None},
(*----- refresh control ---*)
{{gTick, 0}, ControlType -> None},
{{del, $MachineEpsilon}, None},
(*----- state variables for sim---*)
{{runningState, "STOP"}, None},(*"RUNNING","STEP","SUSPEND","STOP"*)
{{wasHitLastTime, False}, None},
{{r0, .05}, None},
{{w, 1.3}, None},
{{simTime, 0}, None},
TrackedSymbols :> {gTick},
ControlPlacement -> Left,
SynchronousUpdating -> False,
SynchronousInitialization -> False,
ContinuousAction -> False,
Alignment -> Center,
ImageMargins -> 0,
FrameMargins -> 0,
Paneled -> True,
Frame -> False,
AutorunSequencing -> {1},
Initialization :>
{
(*definitions used for parameter checking*)
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] &);
(* This list object is used to store and the manage the time series generated during running of the simulation so that display \
object can use it to make plots and the animation *)
listClass[$size_?integer(*maximum size of the list*),
$r0_?numericStrictPositive,
$L0_?numericStrictPositive] :=
Module[{size, lists, r0, L0, self},
(*lists has the structure {idx,{ {t,x,vx,y,vy,\[Theta],\[Omega],{cmX,cmY},{xTrace,yTrace}}},... ,....}*)
self@getSize[] := size;
self@setSize[n_] := ( size = n ; lists = {0, Table[Table[0, {9}], {n}] });
self@add[{time_?numericPositive,(* current time*)
x_?numeric,(*bob position*)
vx_?numeric,(*bob speed*)
y_?numeric,(*arm spring position*)
vy_?numeric,(*arm spring speed*)
\[Theta]_?numeric,(*disk angle*)
\[Omega]_?numeric,(*disk angular speed*)
cm_List(*center of mass coordinates*)
}
] := Module[{},
lists[[1]]++;
If[lists[[1]] > size, lists[[1]] = 1];
lists[[ 2, lists[[1]], 1 ]] = time;
lists[[ 2, lists[[1]], 2 ]] = x;
lists[[ 2, lists[[1]], 3 ]] = vx;
lists[[ 2, lists[[1]], 4 ]] = y;
lists[[ 2, lists[[1]], 5 ]] = vy;
lists[[ 2, lists[[1]], 6 ]] = \[Theta];
lists[[ 2, lists[[1]], 7 ]] = \[Omega];
lists[[ 2, lists[[1]], 8 ]] = cm;
lists[[ 2, lists[[1]], 9 ]] = {(L0 + y + r0) Sin[\[Theta]] + x Cos[\[Theta]], -(L0 + y + r0) Cos[\[Theta]] + x Sin[\[Theta]]}
];
self@setIC[{x_?numeric, vx_?numeric, y_?numeric, vy_?numeric, \[Theta]_?numeric, \[Omega]_?numeric, cm_List,
rr0_?numericStrictPositive, LL0_?numericStrictPositive}] := (
list@clear[];
r0 = rr0;
L0 = LL0;
list@add[{0, x, vx, y, vy, \[Theta], \[Omega], cm}]
);
self@getCurrent[] := lists[[2, lists[[1]] ]];
self@clear[] := lists[[1]] = 0;
self@getPositionList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 2}]] ] ;
self@getSpeedList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 3}]] ] ;
self@getYList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 4}]] ] ;
self@getYSpeedList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 5}]] ] ;
self@getThetaList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 6}]] ] ;
self@getOmegaList[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 7}]] ] ;
self@getBobTraceData[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, 9]] ] ;
self@getCMTraceData[] := Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, 8]] ] ;
(*--- constructor---*)
self@setSize[$size];
r0 = $r0;
L0 = $L0;
self
];
(*-----------------------------------------------------------*)
solverClass[] :=
Module[{solve, self},
(*----------------------------------------------*)
solve[M0_?numericStrictPositive,
m0_?numericPositive,
len0_?numericStrictPositive,
eq1_,
eq2_,
eq3_,
t_,
x_,
y_,
\[Theta]_,
delT_?numericStrictPositive(*time length to integrate over for NDSolve*),
k1_?numericStrictPositive] :=
Module[{currentSolution, currentX, currentXV, currentY, currentYV, current\[Theta], currentOmega,
ic, $x, $xv, $y, $yv, $\[Theta], $\[Omega], simTime, trace, g = 9.81, cm, L0},
L0 = initialRelaxedLength[len0, m0, M0, k1, g];
{simTime, $x, $xv, $y, $yv, $\[Theta], $\[Omega], cm, trace} = list@getCurrent[];
ic = {x[0] == $x, x'[0] == $xv, y[0] == $y, y'[0] == $yv, \[Theta][0] == Mod[$\[Theta], 2 Pi], \[Theta]'[0] == $\[Omega]};
currentSolution =
Quiet@NDSolve[Flatten@{eq1, eq2, eq3, ic}, {x, x', y, y', \[Theta], \[Theta]'}, {t, 0, delT}, Method -> {"BDF"},
MaxSteps -> Infinity];
currentSolution = First@currentSolution;
currentX = x[delT] /. currentSolution;
currentXV = x'[delT] /. currentSolution;
currentY = y[delT] /. currentSolution;
currentYV = y'[delT] /. currentSolution;
current\[Theta] = \[Theta][delT] /. currentSolution;
currentOmega = \[Theta]'[delT] /. currentSolution;
cm = calculateCenterOfMass[m0, M0, L0 + currentY, current\[Theta], currentX];
list@add[{simTime + delT, currentX, currentXV, currentY, currentYV, current\[Theta], currentOmega, cm}]
];
(*----------------------------------------------*)
solve[
m0_?numericPositive,
len0_?numericStrictPositive,
eq2_,
eq3_,
t_,
y_,
\[Theta]_,
delT_?numericStrictPositive(*time length to integrate over for NDSolve*),
k1_?numericStrictPositive] :=
Module[{currentSolution, currentY, currentYV, current\[Theta], currentOmega, ic, $x, $xv, $y, $yv, $\[Theta], $\[Omega],
simTime, trace, g = 9.81, cm, L0},
L0 = initialRelaxedLength[len0, m0, M0, k1, g];
{simTime, $x, $xv, $y, $yv, $\[Theta], $\[Omega], cm, trace} = list@getCurrent[];
ic = {y[0] == $y, y'[0] == $yv, \[Theta][0] == Mod[$\[Theta], 2 Pi], \[Theta]'[0] == $\[Omega]};
currentSolution =
Quiet@NDSolve[Flatten@{eq2, eq3, ic}, {y, y', \[Theta], \[Theta]'}, {t, 0, delT}, Method -> {"BDF"}, MaxSteps -> Infinity];
currentSolution = First@currentSolution;
currentY = y[delT] /. currentSolution;
currentYV = y'[delT] /. currentSolution;
current\[Theta] = \[Theta][delT] /. currentSolution;
currentOmega = \[Theta]'[delT] /. currentSolution;
cm = { (L0 + currentY)/2 Sin[current\[Theta]], -((currentY + L0)/2) Cos[current\[Theta]]};
list@add[{simTime + delT, $x, $xv, currentY, currentYV, current\[Theta], currentOmega, cm}]
];
(*----------------------------------------------*)
self@step[
m0_?numericPositive,(*pendulum mass*)
len0_?numericStrictPositive,
k1_?numericPositive,(*pendulum arm k*)
M0_?numericPositive,(*bob mass*)
delT_?numericStrictPositive(*time length to integrate over for NDSolve*)
] := Module[{eq1, eq2, eq3, t, y, x, \[Theta], current\[Theta], g = 9.81},
(*these are the equations of motion for the case of bob having mass and not*)
(*these were derived using Lagrnagian method. See my report for how this was done, link is below*)
If[M0 > $MachineEpsilon,
eq1 = g M0 Sin[\[Theta][t]] + k x[t] -
M0 Derivative[1][\[Theta]][t] (-Derivative[1][y][t] + x[t] Derivative[1][\[Theta]][t]) +
M0 (Derivative[1][y][t] Derivative[1][\[Theta]][t] + (x^\[Prime]\[Prime])[
t] + (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) (\[Theta]^\[Prime]\[Prime])[t]) == 0;
eq2 = 1/2 g m0 Sin[\[Theta][t]] (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) +
g M0 (Cos[\[Theta][t]] x[t] + Sin[\[Theta][t]] (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t])) +
2/3 m0 (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) Derivative[1][y][t] Derivative[1][\[Theta]][t] +
1/3 m0 (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t])^2 (\[Theta]^\[Prime]\[Prime])[t] +
1/2 M0 (2 Derivative[1][x][t] (-Derivative[1][y][t] + x[t] Derivative[1][\[Theta]][t]) +
2 Derivative[1][y][t] (Derivative[1][x][t] + (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) Derivative[1][\[Theta]][t]) +
2 x[t] (Derivative[1][x][t] Derivative[1][\[Theta]][t] - (y^\[Prime]\[Prime])[t] +
x[t] (\[Theta]^\[Prime]\[Prime])[t]) +
2 (len0 + (g m0)/(3 k1) + (g M0)/k1 +
y[t]) (Derivative[1][y][t] Derivative[1][\[Theta]][t] + (x^\[Prime]\[Prime])[
t] + (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) (\[Theta]^\[Prime]\[Prime])[t])) == 0;
eq3 = 1/2 g m0 (1 - Cos[\[Theta][t]]) + g M0 (1 - Cos[\[Theta][t]]) + k1 y[t] -
1/3 m0 (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) Derivative[1][\[Theta]][t]^2 -
M0 Derivative[1][\[Theta]][
t] (Derivative[1][x][t] + (len0 + (g m0)/(3 k1) + (g M0)/k1 + y[t]) Derivative[1][\[Theta]][t]) +
1/3 m0 (y^\[Prime]\[Prime])[t] -
M0 (Derivative[1][x][t] Derivative[1][\[Theta]][t] - (y^\[Prime]\[Prime])[t] + x[t] (\[Theta]^\[Prime]\[Prime])[t]) == 0;
solve[M0, m0, len0, eq1, eq2, eq3, t, x, y, \[Theta], delT, k1]
,
eq2 = 1/2 g m0 Sin[\[Theta][t]] (len0 + (g m0)/(3 k1) + y[t]) +
2/3 m0 (len0 + (g m0)/(3 k1) + y[t]) Derivative[1][y][t] Derivative[1][\[Theta]][t] +
1/3 m0 (len0 + (g m0)/(3 k1) + y[t])^2 (\[Theta]^\[Prime]\[Prime])[t] == 0;
eq3 = 1/2 g m0 (1 - Cos[\[Theta][t]]) + k1 y[t] - 1/3 m0 (len0 + (g m0)/(3 k1) + y[t]) Derivative[1][\[Theta]][t]^2 +
1/3 m0 (y^\[Prime]\[Prime])[t] == 0;
solve[m0, len0, eq2, eq3, t, y, \[Theta], delT, k1]
]
];
self
];
(*--------- displayClass ---------------------------------------*)
displayClass[] := Module[{makeCounters, makeDiskDiagram, makePlot, self},
(*--------- private ---------------------------------------*)
makeCounters[
m0_?numericPositive,(*pendulum arm mass*)
M0_?numericPositive,(*bob mass*)
len0_?numericStrictPositive,(*pendulum length mass*)
k_?numericPositive,
k1_?numericPositive] :=
Module[{h1, h2, currentMomentOfInertia, x, xv, y, yv, \[Theta], \[Omega], PE, KE, trace, g = 9.81, simTime, angularMomentum,
cm, L0},
{simTime, x, xv, y, yv, \[Theta], \[Omega], cm, trace} = list@getCurrent[];
(*these are the KE and PE for the system. See my report for how this was done, link is below*)
If[M0 > $MachineEpsilon,
L0 = len0 + (m0 g)/(3 k1) + (M0 g)/k1;
KE = 1/6 m0 yv^2 + 1/6 m0 (len0 + (g m0)/(3 k1) + (g M0)/k1 + y)^2 \[Omega]^2 +
1/2 M0 ((-yv + x \[Omega])^2 + (xv + (len0 + (g m0)/(3 k1) + (g M0)/k1 + y) \[Omega])^2);
PE = 1/2 k x^2 + 1/2 k1 y^2 + 1/2 g m0 (1 - Cos[\[Theta]]) (len0 + (g m0)/(3 k1) + (g M0)/k1 + y) +
g M0 (Sin[\[Theta]] x + (1 - Cos[\[Theta]]) (len0 + (g m0)/(3 k1) + (g M0)/k1 + y))
,
L0 = len0 + (m0 g)/(3 k1);
KE = 1/6 m0 yv^2 + 1/6 m0 (len0 + (g m0)/(3 k1) + y)^2 \[Omega]^2;
PE = 1/2 k1 y^2 + 1/2 g m0 (1 - Cos[\[Theta]]) (len0 + (g m0)/(3 k1) + y)
];
currentMomentOfInertia = m0 (L0 + y)^2/3;
angularMomentum = m0 (L0 + y)^2/3*\[Omega] + M0 xv (L0 + y) + M0 (L0 + y)^2 \[Omega] + M0 (x \[Omega] - yv) x;
h1 = Style[Grid[{
{"\[Theta]", "\[Omega]", "bob position", "bob velocity", Row[{Style["y", Italic], "(", Style["t", Italic], ")"}]},
{"(degree)", "(rad/sec)", "(meter)", "(meter/sec)", "(meter)"},
{ padIt2[180./Pi*\[Theta], {5, 2}], padIt1[\[Omega], {6, 3}], padIt1[x, {5, 3}], padIt1[xv, {5, 3}], padIt1[y, {4, 3}]}
},
Frame -> {All, None, { {{1, 2}, {1, 5}} -> True , {{3, 3}, {1, 5}} -> True}},
FrameStyle -> Gray,
Spacings -> 1,
ItemSize -> {{All, 2 ;; -1} -> 5},
Alignment -> Center], 12];
h2 = Style[Grid[{
{Row[{Style["I", Italic], "\[InvisibleSpace] \[InvisibleSpace]\[Omega] (J sec)"}],
Row[{"P.E. (J)"}],
Row[{"K.E. (J)"}],
Row[{"energy (J)"}]
},
{ padIt2[angularMomentum, {7, 4}],
padIt1[PE, {7, 3}],
padIt1[KE, {6, 3}],
padIt1[PE + KE, {8, 4}]
}},
Frame -> All,
FrameStyle -> Gray,
Spacings -> 1,
ItemSize -> {{All, 2 ;; -1} -> 6},
Alignment -> Center], 12];
Grid[{{h1}, {h2}}, Spacings -> {0, .1}]
];
(*---------------private----------------------------------*)
makeDiskDiagram[m0_?numericStrictPositive,
M0_?numericPositive,
r0_?numericStrictPositive,
{width_?numericStrictPositive, height_?numericStrictPositive},
traceBob_?bool,
traceCM_?bool,
w_?numericStrictPositive,
k1_?numericStrictPositive,
len0_?numericStrictPositive
] := Module[{a = 2 r0, b = 0.5 r0, x, y, xx, xv, yy, yv, \[Theta], \[Omega], traceData, traceOn = False, spring, pin, pin0,
bob, simTime, frame0, frame1, springArm, bobHasMass, cm, cg, g = 9.81, L0, options},
bobHasMass = M0 > $MachineEpsilon;
If[bobHasMass, L0 = len0 + (m0 g)/(3 k1) + (M0 g)/k1, L0 = len0 + (m0 g)/(3 k1)];
{simTime, xx, xv, yy, yv, \[Theta], \[Omega], cm, currentTrace} = list@getCurrent[];
cg = Text[Style["\[CirclePlus]", 14], cm];
If[traceBob,
traceData = list@getBobTraceData[];
traceOn = True
];
If[bobHasMass,
x = xx*Cos[\[Theta]];
y = xx*Sin[\[Theta]];
bob = {Red, Opacity[M0/0.1], EdgeForm[Thin], Disk[{L0 + yy + r0, xx}, r0]};
(*bob spring*)
spring = {Gray, Line@makeSpring[L0 + yy + r0, 0, L0 + yy + r0, xx, r0, 14]};
(*tube frame*)
frame0 = {Blue, Line[{ {L0 + yy , -w/2}, {L0 + yy, w/2}}]};
frame1 = {Blue, Line[{ {L0 + yy + a , -w/2}, {L0 + yy + a, w/2}}]};
pin0 = {Blue, EdgeForm[Thin], Disk[{L0 + yy + r0, 0}, b/2]};
];
pin = {Blue, EdgeForm[Thin], Disk[{0, 0}, b/2]};
(*pendulum arm*)
springArm = {Thickness[0.004], Line@makeSpring[0, 0, L0 + yy, 0, r0, 24]};
options = {Axes -> False,
PlotRange -> {{-w, w}, {-w , w/2}},
ImageSize -> {width, height},
ImagePadding -> All,
ImageMargins -> 0,
AxesOrigin -> {0, 0},
PlotRangeClipping -> False,
PlotRangePadding -> None,
Frame -> True,
AspectRatio -> 1.15,
GridLines -> Automatic,
GridLinesStyle -> Directive[Gray, Dashed]};
If[bobHasMass,
Graphics[{
Rotate[frame0, -(Pi/2) + \[Theta], {0, 0}],
Rotate[frame1, -(Pi/2) + \[Theta], {0, 0}],
Rotate[{Black, spring}, -(Pi/2) + \[Theta], {0, 0}],
Rotate[{Black, springArm}, -(Pi/2) + \[Theta], {0, 0}],
Rotate[bob, -(Pi/2) + \[Theta], {0, 0}],
Rotate[pin, -(Pi/2) + \[Theta], {0, 0}],
Rotate[pin0, -(Pi/2) + \[Theta], {0, 0}],
If[traceCM, cg, Sequence @@ {}],
If[traceOn, {Thin, Magenta, Line[traceData]}, Sequence @@ {}]
}, options
],
Graphics[{
Rotate[{Black, springArm}, -(Pi/2) + \[Theta], {0, 0}],
Rotate[pin, -(Pi/2) + \[Theta], {0, 0}],
If[traceCM, cg, Sequence @@ {}]
}, options
]
]
];
(*-------------------- private -----------------------------*)
makePlot[plotType_String,
w_?numericStrictPositive] := Module[{plotTitle, data},
Which[plotType == "disk angular velocity",
data = list@getOmegaList[];
If[Length[data] == 0, data = {{0, 0}}];
plotTitle = {{Style["\[Omega] (rad/sec)", 12], None}, {Style["time (sec)", 12],
Style["disk angular velocity vs. time", 12]}};
plotRange = {{data[[1, 1]], data[[-1, 1]]}, All}
,
plotType == "disk angle",
data = list@getThetaList[];
data[[All, 2]] = Map[180.0/Pi*# &, data[[All, 2]]];
If[Length[data] == 0, data = {{0, 0}}];
plotTitle = {{Style["\[Theta] (deg)", 12], None}, {Style["time (sec)", 12],
Style[Row[{"disk angle \[Theta](", Style["t", Italic], ") vs. time"}], 12]}};
plotRange = {{data[[1, 1]], data[[-1, 1]]}, {0, 360}}
,
plotType == "bob position",
data = list@getPositionList[];
If[Length[data] == 0, data = {{0, 0}}];
plotTitle = {{Style["position (meter)", 12], None}, {Style["time (sec)", 12], Style["bob position vs. time", 12]}};
plotRange = {{data[[1, 1]], data[[-1, 1]]}, {-w/2, w/2}}
,
plotType == "bob velocity",
data = list@getSpeedList[];
If[Length[data] == 0, data = {{0, 0}}];
plotTitle = {{Style[Row[{"velocity (", Style["m/s", Italic], ")"}], 12], None}, {Style["time (sec)", 12],
Style["bob velocity vs. time", 12]}};
plotRange = {{data[[1, 1]], data[[-1, 1]]}, All}
,
plotType == "vertical displacement",
data = list@getYList[];
If[Length[data] == 0, data = {{0, 0}}];
plotTitle = {{Style[Row[{Style["y", Italic], " (meter)"}], 12], None}, {Style["time (sec)", 12],
Style[Row[{"vertical displacement ", Style["y", Italic], "(", Style["t", Italic], ") vs. time"}], 12]}};
plotRange = {{data[[1, 1]], data[[-1, 1]]}, All}
];
ListPlot[ data,
PlotRange -> plotRange,
AxesOrigin -> {0, 0},
Joined -> True,
Frame -> True,
GridLines -> Automatic,
ImageSize -> {300, 160},
ImagePadding -> {{50, 10}, {35, 20}},
ImageMargins -> {{2, 1}, {1, 1}},
Axes -> False,
FrameLabel -> plotTitle,
PlotStyle -> Red,
AspectRatio -> 0.38]
];
(*------- public ----------------------------------*)
self@makeDisplay[plotType_String,
m0_?numericPositive,
M0_?numericPositive,
len0_?numericStrictPositive,
r0_?numericStrictPositive,
k1_?numericPositive,
showPlots_?bool,
showCounters_?bool,
traceCM_?bool,
traceBob_?bool,
w_] :=
If[showPlots,
If[showCounters,
Grid[{
{makeCounters[m0, M0, len0, k, k1]},
{makePlot[plotType, w]},
{makeDiskDiagram[m0, M0, r0, {contentSizeW, 0.695 contentSizeW}, traceBob, traceCM, w, k1, len0]}
}, Spacings -> 0, Frame -> None]
,
Grid[{
{makePlot[plotType, w]},
{makeDiskDiagram[m0, M0, r0, {contentSizeW, 1.02 contentSizeW}, traceBob, traceCM, w, k1, len0]}
}, Spacings -> 0, Frame -> None]
],
If[showCounters,
Grid[{
{makeCounters[m0, M0, len0, k, k1]},
{makeDiskDiagram[m0, M0, r0, {contentSizeW, 1.25 contentSizeW}, traceBob, traceCM, w, k1, len0]}
}, Spacings -> 0, Frame -> None]
,
Grid[{
{makeDiskDiagram[m0, M0, r0, {1.01 contentSizeW, 1.575 contentSizeW}, traceBob, traceCM, w, k1, len0]}
}, Spacings -> 0, Frame -> None]
]
];
self
];
(*-------------------------*)
initialRelaxedLength[len0_?numericStrictPositive, m0_?numericStrictPositive, M0_?numericPositive, k_?numericStrictPositive,
g_?numericStrictPositive] := (len0 + (m0 g)/(3 k) + (M0 g)/k);
(*-------------------------*)
calculateCenterOfMass[m0_?numericStrictPositive, M0_?numericPositive, y_?numeric, \[Theta]_?numeric, x_?numeric ] :=
{1/(m0 + M0) (m0 y/2 Sin[\[Theta]] + M0 (y Sin[\[Theta]] + x Cos[\[Theta]])),
1/(m0 + M0) (-m0 y/2 Cos[\[Theta]] + M0 (-y Cos[\[Theta]] + x Sin[\[Theta]]))};
(*--------------------------------------------*)
(* helper function for formatting *)
(*--------------------------------------------*)
padIt1[v_?numeric, f_List] :=
AccountingForm[Chop[v] , f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, SignPadding -> True];
padIt1[v_?numeric, f_Integer] :=
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];
(*This function is called to make spring,
based on code by Arpad Kosa from WRI demo*)(*at Wolfram web site modified by me.This returns a Line which is the spring*)
(*szel, larger number means bigger spring width*)
makeSpring[xFirst_?numeric, yFirst_?numeric, xEnd_?numeric, yEnd_?numeric, szel_?numeric, nTurns_?integerStrictPositive] :=
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)/nTurns;
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, nTurns - 2}
];
{{xFirst, yFirst}}~Join~{{xFirst + hx*(dh + veghossz)/hossz, yFirst + hy*(dh + veghossz)/hossz}}~Join~tbl~
Join~{{xFirst + hx*((nTurns - 1)*dh + veghossz)/hossz, yFirst + hy*((nTurns - 1)*dh + veghossz)/hossz}}~Join~{{xEnd, yEnd}}
];
(*--- constant parameters size and width of display ---*)
contentSizeW = 300;
contentSizeH = 405;
(* objects used by the simulation. These must be here in the initialization section *)
(*listClass[size,r0,L0], where L0=len0+(m0 g)/(3 k1)+(M0 g)/k1*)
list = listClass[500, 0.05, 0.6 + (1 9.81)/(3 500) + (10 9.81)/500];
(*list@add[{time,x,vx,y,vy,\[Theta],\[Omega],cm*)
list@add[{0, -0.35, 0.0, 0.1, 0.0, 15 Degree, 0.0, {0, 0}}];
solver = solverClass[];
display = displayClass[];
}
]