(* by Nasser M. Abbasi, version: Nov 10, 2012 *)
Manipulate[
gTick;
Text@Module[{g1, wasHit = False},
If[setIC,
setIC = False;
list@setIC[{initialBobX, initialBobSpeed, initialTheta,
initialThetaDot, r0, L0}];
isRelaxedSpring =
If[Abs[initialBobX] <= $MachineEpsilon, True, False]
];
If[runningState == "RUNNING" || runningState == "STEP",
wasHit = solver@step[m0, L0, k, M0, delT, w, r0];
If[wasHit,(*animate a hit on the edge unless last time it was no \
hit*)
If[wasHitLastTime == True,
wasHit = False,
wasHitLastTime = True
],
wasHitLastTime = False
];
If[runningState == "STEP", runningState = "STOP", gTick += del]
];
g1 = display@
makeDisplay[plotType, m0, M0, L0, r0, wasHit, k, showPlots,
showCounters, trace, w];
FinishDynamic[];
g1
],
(*---------- control layout ------------*)
Text@Grid[{
{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}]
}}, Spacings -> {0.5, .2}
]},
{Style["simulation parameters", 12]},
{Grid[
{
{"simulation speed",
Manipulator[Dynamic[delT, {delT = #} &], {0.001, 0.05, 0.001},
ImageSize -> Tiny, ContinuousAction -> False],
Style[Dynamic@padIt2[delT, {3, 3}], 11]
},
{"pendulum mass",
Manipulator[Dynamic[m0,
{m0 = Chop@#;
setIC = True;
gTick += del
} &], {1, 100, 1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[m0, 3], 11]
},
{"pendulum length",
Manipulator[Dynamic[L0,
{L0 = Chop@#;
setIC = True;
gTick += del
} &], {0.1, 2, 0.1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[L0, {2, 1}], 11]
},
{"tube length",
Manipulator[Dynamic[w,
{(w = #) &,
(
w = Chop@#;
If[initialBobX > (w/2 - r0),
initialBobX = w/2 - r0
,
If[initialBobX < (-w/2 + r0),
initialBobX = -w/2 + r0
]
];
setIC = True;
gTick += del
) &
}], {0.1, 1, 0.1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[w, {1, 1}], 11]
},
{"bob mass",
Manipulator[Dynamic[M0,
{M0 = Chop@#;
setIC = True;
gTick += del
} &], {1, 100, 1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[M0, 3], 11]
},
{"spring stiffness",
Manipulator[Dynamic[k,
{k = Chop@#;
setIC = True;
gTick += del
} &], {0, 500, .1}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[k, {4, 1}], 11]
}
}, Spacings -> {.6, .5}, Alignment -> Left, Frame -> True,
FrameStyle -> Directive[Thickness[.005], Gray]
]
},
{Style["initial conditions", 12]},
{Grid[{
{"angular velocity",
Manipulator[Dynamic[initialThetaDot,
{initialThetaDot = Chop@#;
setIC = True;
gTick += del
} &], {-0.1, 0.1, .01}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt1[initialThetaDot, {2, 2}], 11]
},
{"angle (degree)",
Manipulator[Dynamic[initialTheta,
{initialTheta = Chop@#;
setIC = True;
gTick += del
} &], {0, 2 Pi, 2 Pi/100.}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt2[initialTheta*180.0/Pi, {4, 1}], 11]
},
{"bob position",
Manipulator[Dynamic[initialBobX,
{(initialBobX = #) &,
(
initialBobX = Chop[#];
If[initialBobX > (w/2 - r0),
initialBobX = w/2 - r0
,
If[initialBobX < (-w/2 + r0),
initialBobX = -w/2 + r0
]
];
setIC = True;
gTick += del
) &
}], {-0.47, 0.47, 0.01}, ImageSize -> Tiny,
ContinuousAction -> False],
Style[Dynamic@padIt1[initialBobX, {2, 2}], 11]
},
{
Grid[{{
Row[{"relax", Spacer[2],
Checkbox[Dynamic[isRelaxedSpring,
{isRelaxedSpring = #;
initialBobX = 0;
setIC = True;
gTick += del
} &]]}],
Row[{"trace", Spacer[2],
Checkbox[Dynamic[trace,
{trace = #;
If[#, list@clearTrace[]];
gTick += del
} &]]}]
}}, Spacings -> {0.4, 0}], SpanFromLeft
},
{"bob 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]
}, ImageSize -> All, ContinuousAction -> False,
Enabled -> Dynamic[showPlots]]
,
PopupMenu[Dynamic[testCase, {testCase = #;
runningState = "STOP";
Which[testCase == 1,
(
isRelaxedSpring = False; delT = 0.01; m0 = 40; k = 1;
w = .7;
initialThetaDot = 0.01; initialTheta = 45 Degree;
initialBobX = 0.3; initialBobSpeed = 0; M0 = 10;
trace = False;
setIC = True;
runningState = "STOP";
plotType = "bob position"
),
testCase == 2,
(
isRelaxedSpring = False; delT = 0.02; m0 = 16; k = 500;
w = 1;
initialThetaDot = 0.01; initialTheta = 45 Degree;
initialBobX = -.2; initialBobSpeed = 0; M0 = 30;
trace = True;
setIC = True;
runningState = "STOP";
plotType = "disk angular velocity"
),
testCase == 3,
(
isRelaxedSpring = True; delT = 0.03; m0 = 50; k = 0.5;
M0 = 100; trace = False; w = .3;
initialThetaDot = 0; initialTheta = 270 Degree;
initialBobX = .10; initialBobSpeed = 0.1;
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];
list@add[{0, initialBobX, initialBobSpeed,
initialTheta, initialThetaDot}];
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, 25}, None},
{{M0, 40}, None},
{{L0, 0.5}, None},
{{delT, 0.01}, None},
{{k, 1}, None},
{{w, .5}, None},
{{initialThetaDot, 0.09}, None},
{{initialTheta, 45 Degree}, None},
{{initialBobX, 0.05}, None},
{{initialBobSpeed, 0}, None},
{{isRelaxedSpring, False}, None},
{{plotType, "bob position"}, None},
{{testCase, 1}, None},
{{showPlots, False}, None},
{{showCounters, True}, None},
{{xScale, 400}, None},
{{trace, 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, .03}, None},
TrackedSymbols :> {gTick},
ControlPlacement -> Left,
SynchronousUpdating -> False,
SynchronousInitialization -> True,
ContinuousAction -> False,
Alignment -> Center,
ImageMargins -> 0,
FrameMargins -> 0,
Paneled -> True,
Frame -> False,
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,
v,\[Theta],\[Omega],{xTrace,yTrace}}},... ,....}*)
self@getSize[] := size;
self@setSize[n_] := ( size = n ;
lists = {0, Table[Table[0, {6}], {n}] });
self@add[{time_?numericPositive,(* current time*)
x_?numeric,(*bob position*)
v_?numeric,(*bob speed*)
\[Theta]_?numeric,(*disk angle*)
\[Omega]_?numeric(*disk angular speed*)
}
] := 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 ]] = v;
lists[[ 2, lists[[1]], 4 ]] = \[Theta];
lists[[ 2, lists[[1]], 5 ]] = \[Omega];
lists[[ 2, lists[[1]],
6 ]] = {(L0 + r0) Sin[\[Theta]] +
x Cos[\[Theta]], -(L0 + r0) Cos[\[Theta]] + x Sin[\[Theta]]}
];
self@
setIC[{x_?numeric,
v_?numeric, \[Theta]_?numeric, \[Omega]_?numeric,
rr0_?numericStrictPositive, LL0_?numericStrictPositive}] := (
list@clear[];
r0 = rr0;
L0 = LL0;
list@add[{0, x, v, \[Theta], \[Omega]}]
);
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@getThetaList[] :=
Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 4}]] ] ;
self@getOmegaList[] :=
Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, {1, 5}]] ] ;
self@getTraceData[] :=
Module[{len = lists[[1]]}, lists[[ 2, 1 ;; len, 6]] ] ;
(*--- constructor---*)
self@setSize[$size];
r0 = $r0;
L0 = $L0;
self
];
(*-----------------------------------------------------------*)
solverClass[] :=
Module[{solve, self},
(*----------------------------------------------*)
solve[M0_?numericStrictPositive,
m0_?numericPositive,
L0_?numericStrictPositive,
eq1_,
eq2_,
t_,
x_,
\[Theta]_,
w_?numericStrictPositive,
delT_?
numericStrictPositive(*time length to integrate over for \
NDSolve*),
r0_?numericStrictPositive] :=
Module[{wasHit = False, currentSolution, currentX, currentV,
current\[Theta], currentOmega, sol,
ic, $x, $v, $\[Theta], $\[Omega], simTime, trace, newV,
newOmega, eqq1, eqq2},
{simTime, $x, $v, $\[Theta], $\[Omega], trace} =
list@getCurrent[];
ic = {x[0] == $x,
x'[0] == $v, \[Theta][0] ==
Mod[$\[Theta], 2 Pi], \[Theta]'[0] == $\[Omega]};
currentSolution =
Quiet@NDSolve[
Flatten@{eq1, eq2, ic}, {x, x', \[Theta], \[Theta]'}, {t, 0,
delT}, Method -> {"BDF"}, MaxSteps -> Infinity];
currentSolution = First@currentSolution;
currentX = x[delT] /. currentSolution;
currentV = x'[delT] /. currentSolution;
current\[Theta] = \[Theta][delT] /. currentSolution;
currentOmega = \[Theta]'[delT] /. currentSolution;
(* if bob hits the edge of the disk or the base of the spring,
fix speed and reset things. Assume perfect elastic edge.
Solve such that the total moment of momuntum of the system \
before and after collision remain the same*)
If[currentX > (w/2 - r0) || currentX < (-w/2 + r0),
wasHit = True;
eqq1 = newV == -currentV;
eqq2 = (m0 L0^2/3) currentOmega + M0 (currentV) L0 +
M0 (L0^2 currentOmega) +
M0 currentX^2 currentOmega == (m0 L0^2/3) newOmega +
M0 (newV) L0 + M0 (L0^2 newOmega) + M0 currentX^2 newOmega;
sol = {newV, newOmega} /.
NSolve[{eqq1, eqq2}, {newV, newOmega}];
If[Length[sol] > 1,
If[currentV > 0,
sol = Select[sol, First[#] <= 0 &],
sol = Select[sol, First[#] > 0 &]
]
];
sol = First[sol];
currentV = sol[[1]];
currentOmega = sol[[2]]
];
list@
add[{simTime + delT, currentX, currentV, current\[Theta],
currentOmega}];
wasHit
];
(*----------------------------------------------*)
self@step[
m0_?numericPositive,(*pendulum mass*)
L0_?numericStrictPositive,(*pendulum length*)
k_?numericPositive,(*spring k*)
M0_?numericStrictPositive,(*bob mass*)
delT_?
numericStrictPositive,(*time length to integrate over for \
NDSolve*)
w_?numericPositive,
r0_?numericStrictPositive
] :=
Module[{eq1, eq2, t, x, \[Theta], current\[Theta],
wasHit = False, $x, $v, $\[Theta], $\[Omega], simTime, trace,
g = 9.81, solidPendulum, springAttached},
{simTime, $x, $v, $\[Theta], $\[Omega], trace} =
list@getCurrent[];
solidPendulum = Abs[m0] > $MachineEpsilon;
springAttached = Abs[k] > $MachineEpsilon;
If[solidPendulum,
If[springAttached,
eq1 =
g M0 Sin[\[Theta][t]] + k x[t] -
M0 x[t] Derivative[1][\[Theta]][t]^2 +
M0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t]) == 0;
eq2 =
1/2 g L0 m0 Sin[\[Theta][t]] -
g M0 (-L0 Sin[\[Theta][t]] - Cos[\[Theta][t]] x[t]) +
1/3 L0^2 m0 (\[Theta]^\[Prime]\[Prime])[t] +
1/2 M0 (4 x[t] Derivative[1][x][t] Derivative[
1][\[Theta]][t] +
2 x[t]^2 (\[Theta]^\[Prime]\[Prime])[t] +
2 L0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t])) == 0;
wasHit =
solve[M0, m0, L0, eq1, eq2, t, x, \[Theta], w, delT, r0]
,
eq1 =
g M0 Sin[\[Theta][t]] -
M0 x[t] Derivative[1][\[Theta]][t]^2 +
M0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t]) == 0;
eq2 =
1/2 g L0 m0 Sin[\[Theta][t]] -
g M0 (-L0 Sin[\[Theta][t]] - Cos[\[Theta][t]] x[t]) +
1/3 L0^2 m0 (\[Theta]^\[Prime]\[Prime])[t] +
1/2 M0 (4 x[t] Derivative[1][x][t] Derivative[
1][\[Theta]][t] +
2 x[t]^2 (\[Theta]^\[Prime]\[Prime])[t] +
2 L0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t])) == 0;
wasHit =
solve[M0, m0, L0, eq1, eq2, t, x, \[Theta], w, delT, r0]
]
,
If[springAttached,
eq1 =
g M0 Sin[\[Theta][t]] + k x[t] -
M0 x[t] Derivative[1][\[Theta]][t]^2 +
M0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t]) == 0;
eq2 = -g M0 (-L0 Sin[\[Theta][t]] - Cos[\[Theta][t]] x[t]) +
1/2 M0 (4 x[t] Derivative[1][x][t] Derivative[
1][\[Theta]][t] +
2 x[t]^2 (\[Theta]^\[Prime]\[Prime])[t] +
2 L0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t])) == 0;
wasHit = solve[M0, m0, L0, eq1, eq2, t, x, \[Theta], w, delT, r0]
,
eq1 =
g M0 Sin[\[Theta][t]] -
M0 x[t] Derivative[1][\[Theta]][t]^2 +
M0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t]) == 0;
eq2 = -g M0 (-L0 Sin[\[Theta][t]] - Cos[\[Theta][t]] x[t]) +
1/2 M0 (4 x[t] Derivative[1][x][t] Derivative[
1][\[Theta]][t] +
2 x[t]^2 (\[Theta]^\[Prime]\[Prime])[t] +
2 L0 ((x^\[Prime]\[Prime])[t] +
L0 (\[Theta]^\[Prime]\[Prime])[t])) == 0;
wasHit =
solve[M0, m0, L0, eq1, eq2, t, x, \[Theta], w, delT, r0]
]
];
wasHit
];
self
];
(*--------- displayClass ---------------------------------------*)
displayClass[] :=
Module[{makeCounters, makeDiskDiagram, makePlot, self},
(*--------- private ---------------------------------------*)
makeCounters[
m0_?numericPositive,(*pendulum arm mass*)
M0_?numericStrictPositive,(*bob mass*)
L0_?numericStrictPositive,(*pendulum length mass*)
k_?numericPositive] :=
Module[{h1, h2, currentMomentOfInertia, x,
v, \[Theta], \[Omega], PE, KE, trace, g = 9.81, simTime,
angularMomentum, solidPendulum, springAttached},
solidPendulum = Abs[m0] > $MachineEpsilon;
springAttached = Abs[k] > $MachineEpsilon;
{simTime, x, v, \[Theta], \[Omega], trace} = list@getCurrent[];
If[solidPendulum,
KE =
1/6 L0^2 m0 \[Omega]^2 +
1/2 M0 (x^2 \[Omega]^2 + (v + L0 \[Omega])^2);
currentMomentOfInertia = m0*L0^2/3;
angularMomentum =
m0 L0^2/3*\[Omega] + M0 (L0 v + (L0^2 + x^2) \[Omega]);
If[springAttached,
PE = -(1/2) g L0 m0 Cos[\[Theta]] + 1/2 k x^2 -
g M0 (L0 Cos[\[Theta]] - Sin[\[Theta]] x),
PE = -(1/2) g L0 m0 Cos[\[Theta]] -
g M0 (L0 Cos[\[Theta]] - Sin[\[Theta]] x)
],
KE = 1/2 M0 (x^2 \[Omega]^2 + (v + L0 \[Omega])^2);
currentMomentOfInertia = 0;
angularMomentum = M0 (L0 v + (L0^2 + x^2) \[Omega]);
If[springAttached,
PE = 1/2 k x^2 - g M0 (L0 Cos[\[Theta]] - Sin[\[Theta]] x),
PE = -g M0 (L0 Cos[\[Theta]] - Sin[\[Theta]] x)
]
];
h1 = Style[Grid[{
{"\[Theta]", "\[Omega]", "bob position", "bob velocity",
Style["\[CapitalIota]", Italic]
},
{"(degree)", "(rad/sec)", "(meter)", "(meter/sec)",
"(kg \!\(\*SuperscriptBox[\(meter\), \(2\)]\))"
},
{ padIt2[180./Pi*\[Theta], {5, 2}],
padIt1[\[Omega], {7, 5}], padIt1[x, {5, 4}],
padIt1[v, {6, 3}], padIt2[currentMomentOfInertia, {7, 4}]
}
},
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[{
{"time (sec)",
Row[{Style["I", Italic],
"\[InvisibleSpace] \[InvisibleSpace]\[Omega] (joule second)"}],
Row[{"P.E. (", Style["J", Italic], ")"}],
Row[{"K.E. (", Style["J", Italic], ")"}],
Row[{"energy (", Style["J", Italic], ")"}]
},
{ padIt2[Mod[simTime, 1000], {7, 4}],
padIt2[angularMomentum, {6, 4}], padIt1[PE, {7, 3}],
padIt1[KE, {6, 3}], padIt1[PE + KE, {6, 2}]
}},
Frame -> All,
FrameStyle -> Gray,
Spacings -> 1,
ItemSize -> {{All, 2 ;; -1} -> 6},
Alignment -> Center], 12];
Grid[{{h1}, {h2}}, Spacings -> {0, .1}]
];
(*---------------private----------------------------------*)
makeDiskDiagram[k_?numericPositive,
m0_?numericStrictPositive,
M0_?numericStrictPositive,
L0_?numericStrictPositive,
r0_?numericStrictPositive,
wasHit_?bool,
{ContentSizeW_?numericStrictPositive,
ContentSizeH_?numericStrictPositive},
trace_?bool,
w_] :=
Module[{a = 2 r0, b = 1.1 r0, x, y, splash, xx,
v, \[Theta], \[Omega], traceData, currentTrace, spring, pin,
bob, ext, simTime, frame0, frame1, frame2, frame3, frame4},
{simTime, xx, v, \[Theta], \[Omega], currentTrace} =
list@getCurrent[];
traceData = list@getTraceData[];
If[wasHit,
splash = {
{ Dotted, Red,
Rotate[
Line[{{L0 + r0, Sign[xx] w/2}, {L0 + r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]},
{ Dotted, Red,
Rotate[Line[{{L0 + r0, Sign[xx] w/2}, {L0 + 2 r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]},
{ Dotted, Red,
Rotate[Line[{{L0 + r0, Sign[xx] w/2}, {L0 + 3 r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]},
{ Dotted, Red,
Rotate[Line[{{L0 + r0, Sign[xx] w/2}, {L0 - r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]},
{ Dotted, Red,
Rotate[Line[{{L0 + r0, Sign[xx] w/2}, {L0 - 2 r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]},
{ Dotted, Red,
Rotate[Line[{{L0 + r0, Sign[xx] w/2}, {L0 - 3 r0,
Sign[xx] (w/2 + 2 r0)}}], -(Pi/2) + \[Theta], {0, 0}]}
};
If[xx < 0, xx = -w/2 + r0, xx = w/2 - r0]
,
splash = Sequence @@ {};
];
x = xx*Cos[\[Theta]];
y = xx*Sin[\[Theta]];
traceData = list@getTraceData[];
bob = {Red, Opacity[M0/100], EdgeForm[Thin],
Disk[{L0 + r0, xx}, r0]};
pin = {LightGray, EdgeForm[Thin], Disk[{0, 0}, b]};
ext = 1.1 Sqrt[(L0 + 2 r0)^2 + (w/2)^2];
If[k == 0,
spring = Sequence @@ {},
spring = Line@makeSpring[L0 + r0, -w/2, L0 + r0, xx, r0]
];
frame0 = {Thin, LightGray,
Line[{ {L0, 0}, {L0 , -w/2}, {L0 + a, -w/2}, {L0 + a,
w/2}, {L0 + a, w/2}, {L0, w/2}, {L0, 0}}]};
frame1 = {Blue, Opacity[m0/100], EdgeForm[Thin],
Polygon[{ {0, -b/2}, {L0 , -b/2}, {L0, b/2}, {0, b/2}}]};
frame2 = {Black, Line[{ {L0 , -w/2}, {L0 + a, -w/2}}]};
frame3 = {Black, Line[{ {L0 + a , w/2}, {L0, w/2}}]};
frame4 = {Thin, Red, Line[{ {L0 + r0 , -w/2}, {L0 + r0, w/2}}]};
Framed@Graphics[{
Rotate[frame0, -(Pi/2) + \[Theta], {0, 0}],
Rotate[{White, frame1}, -(Pi/2) + \[Theta], {0, 0}],
Rotate[frame2, -(Pi/2) + \[Theta], {0, 0}],
Rotate[frame3, -(Pi/2) + \[Theta], {0, 0}],
Rotate[frame4, -(Pi/2) + \[Theta], {0, 0}],
Rotate[{Black, spring}, -(Pi/2) + \[Theta], {0, 0}],
Rotate[bob, -(Pi/2) + \[Theta], {0, 0}],
Rotate[pin, -(Pi/2) + \[Theta], {0, 0}],
splash,
If[trace, {Thin, Darker@Red, Line[traceData]},
Sequence @@ {}]
},
Axes -> False,
PlotRange -> {{-ext, ext}, {-ext , ext}},
ImageSize -> {ContentSizeW, 1.01 ContentSizeH},
ImagePadding -> {{5, 5}, {5, 5}},
ImageMargins -> 0,
AspectRatio -> 1
]
];
(*-------------------- 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"],
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"], Style["disk angle 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[
Row[{"position (", Style["m", Italic], ")"}], 12],
None}, {Style["time"], 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"], Style["bob velocity 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 -> {330, 160},
ImagePadding -> {{60, 10}, {30, 25}},
ImageMargins -> {{2, 1}, {1, 1}},
Axes -> False,
FrameLabel -> plotTitle,
PlotStyle -> Red,
AspectRatio -> 0.27]
];
(*------- public ----------------------------------*)
self@makeDisplay[plotType_String,
m0_?numericPositive,
M0_?numericStrictPositive,
L0_?numericStrictPositive,
r0_?numericStrictPositive,
wasHit_?bool,
k_?numericPositive,
showPlots_?bool,
showCounters_?bool,
trace_?bool,
w_] :=
If[showPlots,
If[showCounters,
Grid[{
{makeCounters[m0, M0, L0, k]},
{makePlot[plotType, w]},
{makeDiskDiagram[k, m0, M0, L0, r0,
wasHit, {1.4 ContentSizeW, 0.515 ContentSizeH}, trace, w]}
}, Spacings -> 0]
,
Grid[{
{makePlot[plotType, w]},
{makeDiskDiagram[k, m0, M0, L0, r0,
wasHit, {1.4 ContentSizeW, 0.787 ContentSizeH}, trace, w]}
}, Spacings -> 0]
],
If[showCounters,
Grid[{
{makeCounters[m0, M0, L0, k]},
{makeDiskDiagram[k, m0, M0, L0, r0,
wasHit, {1.4 ContentSizeW, 0.93 ContentSizeH}, trace, w]}
}, Spacings -> 0]
,
Grid[{
{makeDiskDiagram[k, m0, M0, L0, r0,
wasHit, {1.4 ContentSizeW, 1.2 ContentSizeH}, trace, w]}
}, Spacings -> 0]
]
];
self
];
(*--------------------------------------------*)
(* 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];
(*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] :=
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}}
];
(*--- constant parameters size and width of display ---*)
ContentSizeW = 260;
ContentSizeH = 405;
(* objects used by the simulation.
These must be here in the initialization section *)
list = listClass[500, 0.02, 0.5];
list@add[{0, 0.05, 0.0, Pi/4, 0.09}];
solver = solverClass[];
display = displayClass[];
}
]