(*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[]; } ]