BeginPackage["HoppingHoop`"]; HoppingHoop::usage= "HoppingHoop[initialAngle_, initialAngularSpeed_] calculates and plots the motion of the asymmetrically weighted hoop described in the accompanying article."; Options[HoppingHoop]= { initialHoopTranslationSpeed -> noSlipping, inertiaParameter -> 0.01, massRatio -> 0.95, coefficientOfStaticFriction -> 1., coefficientOfKineticFriction -> 0.8, isAirborneTolerance -> 0.01, slipTimeStep -> 0.002, visualizationType -> strobedHoops, numberOfStrobeFlashes -> 40 }; initialHoopTranslationSpeed::usage= "The option initialHoopTranslationSpeed specifies initial speed of horizontal translation of the center of the hoop. The default value is chosen so that the hoop rolls without slipping."; inertiaParameter::usage= "The option inertiaParameter is a dimensionless number, generally quite small, which is related to the moment of rotational inertia of the attached object about its own center of mass. The limit inertiaParameter -> 0 corresponds to an attached point mass and is non-physical."; massRatio::usage= "The option massRatio specifies the ratio of the mass of the object attached to the rim of the hoop to the combined mass of the object plus the hoop itself. The limit massRatio -> 1 corresponds to a completely massless hoop and is non-physical."; coefficientOfStaticFriction ::usage= "The option coefficientOfStaticFriction specifies the coefficient of static friction between the hoop and the supporting surface. This quantity is always greater than the coefficient of kinetic friction."; coefficientOfKineticFriction::usage= "The option coefficientOfKineticFriction specifies the coefficient of kinetic friction between the hoop and the supporting surface. This quantity is always less than the coefficient of static friction."; isAirborneTolerance::usage= "The option isAirborneTolerance specifies the tolerance to be used in determining whether the hoop is in contact with the supporting surface: hoop and surface are considered to be in contact as long as the distance between the bottom of the hoop and the surface does not exceed isAirborneTolerance."; slipTimeStep::usage= "The option slipTimeStep specifies the time step to be used in computing the motion of the hoop during the phase in which it slips as it rolls."; visualizationType::usage= "The option visualizationType specifies how the motion of the hoop is plotted. Allowed values are 'strobedHoops', 'schematic', 'animation', and 'none'."; strobedHoops::usage= "Setting visualizationType -> strobedHoops generates a representation of the motion of the hoop similar to a stroboscopic photograph, showing images of the hoop at numberOfStrobeFlashes equally spaced times."; schematic::usage= "Setting visualizationType -> schematic generates a stroboscopic representation of the motion of the hoop similar to that produced by visualizationType -> strobedHoops, but with less detail. Shown is the trajectory of the center of mass of the system, with superimposed dots indicating the respective positions of the attached object and the center of the hoop at numberOfStrobeFlashes equally spaced times."; animation::usage= "Setting visualizationType -> animation generates a movie of the hoop with numberOfStrobeFlashes frames."; none::usage= "Setting visualizationType -> none suppresses the graphic of the hoop."; plotSoln::usage= "The function plotSoln[] uses the results of function HoppingHoop to generate a visual representation of the motion of the asymmetrically weighted hoop."; Options[plotSoln]= { visualizationType -> strobedHoops, numberOfStrobeFlashes -> 40 }; numberOfStrobeFlashes::usage= "The option numberOfStrobeFlashes specifies the number of individual images of the hoop appearing in the strobedHoops, schematic, and animation visualizations of the motion."; angleOnsetOfSlipping::usage= "The function angleOnsetOfSlipping[initialAngle_, initialAngularSpeed_, massRatio_, inertiaParameter_, coefficientOfStaticFriction_] returns the value in radians of the critical angle at which the hoop begins to slip." xSlip::usage= "The horizontal position of the center of mass at the onset of slippage."; ySlip::usage= "The vertical height of the center of mass at the onset of slippage."; qSlip::usage= "The critical angle of onset of slippage."; tSlip::usage= "The time at which the hoop begins to slip as it rolls."; xHop::usage= "The horizontal position of the center of mass at the instant the hoop loses contact with the supporting surface."; yHop::usage= "The vertical height of the center of mass at the instant the hoop loses contact with the supporting surface."; qHop::usage= "The angle with the vertical of the line segment connecting the center of the hoop with the attached weight at the instant the hoop loses contact with the supporting surface."; tHop::usage= "The time at which the hoop loses contact with the supporting surface."; tImpact::usage= "The time at which the hoop regains contact with the supporting surface."; Begin["`Private`"]; HoppingHoop[(q0_)?NumericQ, (qDot0_)?NumericQ, opts___] := ( {initialVxHoop, e, h, staticFrictionCoeff, kineticFrictionCoeff, tolerance, dt} = {initialHoopTranslationSpeed, inertiaParameter, massRatio, coefficientOfStaticFriction, coefficientOfKineticFriction, isAirborneTolerance, slipTimeStep} /. {opts} /. Options[HoppingHoop]; (* Definitions and housekeeping *) g = 9.8; a = 1; c = a/g*qDot0^2; Print["c = ", c]; cos0 = Cos[q0]; nCyc[x_] = 1/(2 + e*h + 2*h*x)^2* (4*g + 4*g*e*h - 2*g*h^2 - 2*c*g*h^2 + g*e^2*h^2 - g*e*h^3 - c*g*e*h^3 + 8*g*h*x - 4*c*g*h*x + 4*g*e*h^2*x - 4*c*g*e*h^2*x - c*g*e^2*h^3*x + 10*g*h^2*x^2 - 2*c*g*h^2*x^2 + 3*g*e*h^3*x^2 - c*g*e*h^3*x^2 + 4*g*h^3*x^3 - 2*g*h^3*cos0 - 2*c*g*h^3*cos0 - 4*g*h^2*x*cos0 - 4*c*g*h^2*x*cos0 - 2*g*e*h^3*x*cos0 - 2*c*g*e*h^3*x*cos0 - 2*g*h^3*x^2*cos0 - 2*c*g*h^3*x^2*cos0); (* Main logic *) Catch[ If[TrueQ[initialVxHoop == noSlipping], isNotSlipping = True; initialVxHoop = a*qDot0, If[NumericQ[initialVxHoop], isNotSlipping = initialVxHoop == a*qDot0, Return[Print["The option initialHoopTranslationSpeed must be either ", "a numeric value, or the value 'noSlipping'. ", "Execution terminated."]]]]; (* end IF TrueQ *) Clear[CMSoln, qSoln]; If[isNotSlipping, (* THEN *) If[qDot0 == 0, If[q0 == 0 || q0 == Pi, Return[Print["The hoop is in a position of equilibrium.", "It is, and will remain, motionless."]]]; rollingFromRest[q0], rollingWithoutSlipping[q0, qDot0]]; If[NumericQ[qSlip], slipping[xSlip, ySlip, qSlip, vxSlip, vySlip, qDotSlip, nCyc[Cos[qSoln[tSlip - dt]]] - g, nCyc[Cos[qSlip]], tolerance, dt]; freeFall[]], (* ELSE it's initally slipping, so ... *) tSlip = 0; qSlip = q0; qDotSlip = qDot0; xSlip = h*a*Sin[q0]; ySlip = a*(1 + h*Cos[q0]); vxSlip = initialVxHoop + h*a*Cos[q0]*qDot0; vySlip = -h*a*Sin[q0]*qDot0; nSlip = -(-g + a*qDot0^2*h*Cos[q0])/ (1 + 1/(-1 - e*h + h^2)* h*Sin[q0]* (-h*Sin[q0] + Sign[a*qDot0 - initialVxHoop]* (1 + h*Cos[q0])*kineticFrictionCoeff)); qDblDotSlip = nSlip*(h*Sin[q0] - kineticFrictionCoeff*Sign[a*qDot0 - initialVxHoop]* (1 + h*Cos[q0]))/ (a*(1 + e*h - h^2)); qDotPrior = qDotSlip - dt*qDblDotSlip; qPrior = q0 - dt/2*(qDot0 + qDotPrior); priorAyCM = -g + (g - a*qDotPrior^2*h*Cos[qPrior])/ (1 + (h*Sin[qPrior]* (-h*Sin[qPrior] + (1 + h*Cos[qPrior])* Sign[a*qDot0 - initialVxHoop]* kineticFrictionCoeff))/ (-1 - e*h + h^2)); slipping[xSlip, ySlip, qSlip, vxSlip, vySlip, qDotSlip, priorAyCM, nSlip, tolerance, dt]; freeFall[] (* end of the ELSE clause *) ]; (* end IF isNotSlipping *) plotSoln[opts]; ] (* extent of CATCH *) ); (* end function HoppingHoop *) (****************************************************************) rollingWithoutSlipping[(initialAngle_)?NumericQ, (initialAnglrSpeed_)?NumericQ] := Module[{c, thetaDot, theta}, cos0 = Cos[initialAngle]; c = a/g*initialAnglrSpeed^2; thetaDot = (Sqrt[g/a]*Sqrt[2*c + c*e*h + 2*h*cos0 + 2*c*h*cos0 - 2*h*Cos[q[t]]])/ Sqrt[2 + e*h + 2*h*Cos[q[t]]]; cycloid[q_] := a*{q - initialAngle + h*Sin[q], 1 + h*Cos[q]}; qSlip = angleOnsetOfSlipping[initialAngle, initialAnglrSpeed, h, e, staticFrictionCoeff]; Which[(* Case 1 of WHICH *) NumericQ[qSlip], theta = q /. NDSolve[{Derivative[1][q][t] == thetaDot, q[0] == initialAngle}, q, {t, 0, 500}, StoppingTest -> q[t] > qSlip][[1]]; tSlip = t /. FindRoot[theta[t] == qSlip, {t, theta[[1,1,2]]}][[1]]; Print[StringForm["\!\(t\_Slip\) = `` seconds; \!\(\[Theta]\_Slip\) = \ \ `` degrees", NumberForm[tSlip, 4], NumberForm[qSlip/Degree, 4]]]; qSoln[tt_ /; 0 <= tt <= tSlip] := theta[tt]; CMSoln[tt_ /; 0 <= tt <= tSlip] := cycloid[theta[tt]]; {xSlip, ySlip} = cycloid[qSlip]; {vxSlip, vySlip} = D[cycloid[theta[t]], t] /. t -> tSlip; qDotSlip = thetaDot /. q[t] -> qSlip, (* Case 2 of WHICH *) qSlip == neverSlips, theta = q /. NDSolve[{Derivative[1][q][t] == thetaDot, q[0] == initialAngle}, q, {t, 0, 500}, StoppingTest -> q[t] > initialAngle + 2*Pi][[1]]; tImpact = theta[[1,1,2]]; Print["Hoop never slips. Integration terminated at t = ", tImpact, " seconds."]; qSoln[tt_ /; 0 <= tt <= tImpact] := theta[tt]; CMSoln[tt_ /; 0 <= tt <= tImpact] := cycloid[theta[tt]], (* Case 3 of WHICH *) True, Throw[Print["No numerical solution."]] ] (* end WHICH *) ] (* end module RollingWithoutSlipping *) (****************************************************************) rollingFromRest[(initialAngle_)?NumericQ] := Module[{thetaDot, thetaDblDot, theta}, cos0 = Cos[initialAngle]; thetaDot = (Sqrt[g/a]*Sqrt[2*h*cos0 - 2*h*Cos[q[t]]])/ Sqrt[2 + e*h + 2*h*Cos[q[t]]]; thetaDblDot = (g*h/a)Sin[q[t]](2+e*h+2*h*cos0)/(2+e*h+2*h*Cos[q[t]])^2; cycloid[q_] := a*{q - initialAngle + h*Sin[q], 1 + h*Cos[q]}; qSlip = angleOnsetOfSlipping[initialAngle, 0, h, e, staticFrictionCoeff]; Which[(* Case 1 of WHICH *) NumericQ[qSlip], theta = q /. NDSolve[{Derivative[2][q][t] == thetaDblDot, q[0] == initialAngle, Derivative[1][q][0] == 0}, q, {t, 0, 500}, StoppingTest -> q[t] > qSlip][[1]]; tSlip = t /. FindRoot[theta[t] == qSlip, {t, theta[[1,1,2]]}][[1]]; Print[StringForm["\!\(t\_Slip\) = `` seconds; \!\(\[Theta]\_Slip\) = \ \ `` degrees", NumberForm[tSlip, 4], NumberForm[qSlip/Degree, 4]]]; qSoln[tt_ /; 0 <= tt <= tSlip] := theta[tt]; CMSoln[tt_ /; 0 <= tt <= tSlip] := cycloid[theta[tt]]; {xSlip, ySlip} = cycloid[qSlip]; {vxSlip, vySlip} = D[cycloid[theta[t]], t] /. t -> tSlip; qDotSlip = thetaDot /. q[t] -> qSlip, (* Case 2 of WHICH *) qSlip == neverSlips, theta = q /. NDSolve[{Derivative[2][q][t] == thetaDblDot, q[0] == initialAngle, Derivative[1][q][0] == 0}, q, {t, 0, 500}, StoppingTest -> q[t] > initialAngle + 2*Pi][[1]]; tImpact = theta[[1,1,2]]; Print["Hoop never slips. Integration terminated at t = ", tImpact, " seconds."]; qSoln[tt_ /; 0 <= tt <= tImpact] := theta[tt]; CMSoln[tt_ /; 0 <= tt <= tImpact] := cycloid[theta[tt]], (* Case 3 of WHICH *) True, Throw[Print["No numerical solution."]] ] (* end WHICH *) ] (* end module RollingFromRest *) (****************************************************************) angleOnsetOfSlipping[(initialAngle_)?NumericQ, (initialAngularSpeed_)?NumericQ, (h_)?NumericQ, (e_)?NumericQ, (coefficientOfStaticFriction_)?NumericQ] := Module[{c, cos0, slipEqn, angles0ToPi}, If[ !(NumericQ[g] && NumericQ[a]), Return[Print["Assign numerical values to the constants ", "HoppingHoop`Private`g and HoppingHoop`Private`a!"]]]; c = a/g*initialAngularSpeed^2; cos0 = Cos[initialAngle]; slipEqn = (1 - x^2)*h^2* (2 + 6*x*h + e*h + 4*x^2*h^2 + 3*x*e*h^2 - 2*(1 + c)*cos0*h*(1 + x*h + e*h) - c*(2 + e*h)*(1 + x*h + e*h))^2 == coefficientOfStaticFriction^2* (-e^2*h^2 - 4*x^3*h^3 + e*h*(-4 + (1 + c)*h^2) + 2*(-2 + (1 + c)*h^2 + (1 + c)*cos0*h^3) + x*h*(2 + e*h)* (-4 + 2*cos0*h + c*(2 + 2*cos0*h + e*h)) + x^2*h^2*(-10 + 2*cos0*h - 3*e*h + c*(2 + 2*cos0*h + e*h)))^2; angles0ToPi = ArcCos[Cases[x /. NSolve[slipEqn, x], _Real]]; If[angles0ToPi == {}, neverSlips, Min[Select[Union[angles0ToPi, 2*Pi - angles0ToPi, 2*Pi + angles0ToPi, 4*Pi - angles0ToPi], #1 >= initialAngle & ]]] ] (* end module angleOnsetOfSlipping *) (****************************************************************) slipping[X0_, Y0_, theta0_, Vx0_, Vy0_, thetaDot0_, AyPrior_, n0_, tolerance_, dt_] := (x[0] = X0; y[0] = Y0; th[0] = theta0; vx[0] = Vx0; vy[0] = Vy0; thDot[0] = thetaDot0; n = n0; ay[-1] = AyPrior; xctr[i_] = -a*h*Sin[th[i]] + x[i]; yctr[i_] = -a*h*Cos[th[i]] + y[i]; vxctr[i_] = -a*h*Cos[th[i]]*thDot[i] + vx[i]; vyctr[i_] = a*h*Sin[th[i]]*thDot[i] + vy[i]; ayctr[i_] = -g + n + (n*h*Sin[th[i]]* (-h*Sin[th[i]] + (1 + h*Cos[th[i]])* fSign[-a*h*Sin[th[i]] + x[i], thDot[i]]* kineticFrictionCoeff))/ (-1 - e*h + h^2) + a*h*Cos[th[i]]*thDot[i]^2; fSign[hoopVx_, w_] := If[Positive[a*w - hoopVx], Plus[1], -1]; For[i = 0, yctr[i] - a <= tolerance, i++, nRevised = (2*(a + a*h*Cos[th[i]] + 1/2*dt^2* (g - a*h*Cos[th[i]]*thDot[i]^2) - dt*(a*h*Sin[th[i]]*thDot[i] + vy[i]) - y[i]))/ (dt^2*(1 + (h*Sin[th[i]]* (-h*Sin[th[i]] + (1 + h*Cos[th[i]])* fSign[-a*h*Sin[th[i]] + x[i], thDot[i]]* kineticFrictionCoeff))/ (-1 - e*h + h^2))); If[yctr[i] + vyctr[i]*dt + ayctr[i]*dt^2/2 < a, n = Max[0, nRevised]]; ax[i] = kineticFrictionCoeff*n*fSign[vxctr[i], thDot[i]]; ay[i] = n - g; thDblDot[i] = n/(a*(1 - h^2 + h*e))* (h*Sin[th[i]] + kineticFrictionCoeff*(1 + h*Cos[th[i]])); vx[i + 1] = vx[i] + ax[i]*dt; vy[i + 1] = vy[i] + ay[i]*dt; thDot[i + 1] = thDot[i] + thDblDot[i]*dt; x[i + 1] = x[i] + (vx[i] + vx[i + 1])*dt/2; y[i + 1] = y[i] + (vy[i] + vy[i + 1])*dt/2; th[i + 1] = th[i] + (thDot[i] + thDot[i + 1])*dt/2; n = Max[0, n + (ay[i] - ay[i - 1])]]; nsteps = i; tHop = tSlip + nsteps*dt; {xHop, yHop} = {x[nsteps], y[nsteps]}; {vxHop, vyHop} = {vx[nsteps], vy[nsteps]}; {qHop, qDotHop} = {th[nsteps], thDot[nsteps]}; i =. ; CMSoln[tt_ /; Inequality[tSlip, LessEqual, tt, Less, tHop]] = {Interpolation[Table[{j*dt, {x[j], vx[j]}}, {j, 0, nsteps}]][tt - tSlip], Interpolation[Table[{j*dt, {y[j], vy[j]}}, {j, 0, nsteps}]][tt - tSlip]}; qSoln[tt_ /; Inequality[tSlip, LessEqual, tt, Less, tHop]] = Interpolation[Table[{j*dt, {th[j], thDot[j]}}, {j, 0, nsteps}]][tt - tSlip]; Print[nsteps, " time steps in the slipping phase."]; ); (* end function Slipping *) (****************************************************************) freeFall := ( CMHop = {xHop, yHop}; vCMHop = {vxHop, vyHop}; CMSoln[tt_ /; tt >= tHop] := CMHop + vCMHop*(tt - tHop) + {0, -g}/2*(tt - tHop)^2; qSoln[tt_ /; tt >= tHop] := qHop + qDotHop*(tt - tHop); ctrHoop[tt_] := CMSoln[tt] - h*a*{Sin[qSoln[tt]], Cos[qSoln[tt]]}; rootInterval = dt; accuracy = 1./10^6; bottomHoop[tt_] := ctrHoop[tt][[2]] - a; For[t1 = tHop, Positive[bottomHoop[t2 = t1 + rootInterval]], t1 += rootInterval]; While[t2 - t1 > accuracy, tmidpt = (t1 + t2)/2; If[Positive[bottomHoop[tmidpt]], t1 = tmidpt, t2 = tmidpt]]; tImpact = tmidpt; ); (* end function freeFall *) (****************************************************************) plotSoln[opts___] := ( {plotType, ntSteps} = {visualizationType, numberOfStrobeFlashes} /. {opts} /. Options[plotSoln]; ctrHoop[tt_] := CMSoln[tt] - h*a*{Sin[qSoln[tt]], Cos[qSoln[tt]]}; rMass[t_] := ctrHoop[t] + a*{Sin[qSoln[t]], Cos[qSoln[t]]}; deltaT = N[tImpact/ntSteps]; ptsCM = CMSoln /@ Range[0, tImpact, deltaT]; ptsCtrHoop = ctrHoop /@ Range[0, tImpact, deltaT]; ptsMass = rMass /@ Range[0, tImpact, deltaT]; hoops = (Circle[#1, a] & ) /@ ptsCtrHoop; Switch[plotType, strobedHoops, lines = {{GrayLevel[0.8], RGBColor[1, 0.4, 0.4], Line /@ Transpose[{ptsCtrHoop, ptsMass}]}, {AbsolutePointSize[4], RGBColor[1, 0, 0], Point /@ ptsMass}, {AbsolutePointSize[2], RGBColor[1, 0, 0], (If[#1[[2]] > a + tolerance, {RGBColor[1, 1, 0], PointSize[0.01], Point[#1]}, Point[#1]] & ) /@ ptsCtrHoop}}; xRange = {Min[#]-a, Max[#]+a}& [First /@ ptsCtrHoop]; yRange = {Min[#]-a, Max[#]+a}& [Last /@ ptsCtrHoop]; ParametricPlot[CMSoln[tt], {tt, 0, tImpact}, Compiled -> False, Prolog -> {hoops, {AbsoluteThickness[2], RGBColor[1, 1, 0], (Circle[#1, a] & ) /@ Select[ptsCtrHoop, #1[[2]] > 1 + tolerance & ]}}, Epilog -> {lines, AbsoluteThickness[1.5], RGBColor[1, 1, 0], Line[{{xRange[[1]], aa = yRange[[2]] + 0.018}, {xRange[[2]], aa}}], AbsoluteThickness[0.5], Line[{{xRange[[1]], 0}, {xRange[[2]], 0}}]}, AspectRatio -> Automatic, Background -> RGBColor[0.3, 0.4, 1], DefaultColor -> GrayLevel[0], PlotRange -> {{xRange[[1]]-0.2, xRange[[2]]+0.3}, {-0.1, yRange[[2]]*1.1}}, Axes -> None, Frame -> True], (* end code for plotType strobedHoops *) schematic, lines = Graphics[{{GrayLevel[0.8], Line /@ Transpose[{ptsCtrHoop, ptsMass}]}, {AbsolutePointSize[4], Point /@ ptsMass}, {AbsolutePointSize[2], Point /@ ptsCtrHoop}}]; Show[lines, ParametricPlot[CMSoln[tt], {tt, 0, tImpact}, Compiled -> False, Axes -> False, AspectRatio -> Automatic, PlotRange -> All, DisplayFunction -> Identity]], (* end code for plotType schematic *) animation, nColors = 16; arcAngle = 2*Pi/nColors; hoopThickness = 0.05; bkgrndColor = Hue[0.6, 0.2, 0.999]; ground = Graphics[{Line[{{-a, 0}, {5*a, 0}}], Line[{{0, 2.5*a}, {0, 0}}]}]; hoop := Graphics[{Table[{Hue[iColor/nColors], Disk[ctrHoop[ti], a, arcAngle*{iColor - 1, iColor} - qSoln[ti]]}, {iColor, 1, nColors}], {AbsolutePointSize[10], Point[rMass[ti]]}, Circle[ctrHoop[ti], a], {bkgrndColor, Disk[ctrHoop[ti], a - hoopThickness]}, Circle[ctrHoop[ti], a - hoopThickness], {Dashing[{0.03, 0.01}], Line[{ctrHoop[ti] + {0, a}, ctrHoop[ti]}]}, {Dashing[{0.03, 0.01}], Line[{ctrHoop[ti], rMass[ti]}]}}]; points := Graphics[{AbsolutePointSize[5], {RGBColor[1, 0, 0.5], Point /@ Take[ptsCM, i]}, {RGBColor[0, 0, 1], Point /@ Take[ptsCtrHoop, i]}}]; Do[ti = i*deltaT; Show[hoop, ground, If[i >= 1, ParametricPlot[{CMSoln[tt], ctrHoop[tt]}, {tt, 0, ti}, PlotStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}}, Compiled -> False, DisplayFunction -> Identity, AspectRatio -> Automatic], {}], points, Axes -> False, AspectRatio -> Automatic, PlotRange -> All, Background -> bkgrndColor], {i, 0, ntSteps}], (* end DO *) (* end code for plotType animation *) none, Null, (* Do not generate a plot. *) _, Print["Allowed values of visualizationType are strobedHoops, schematic, ", "animation, and none. No plot generated."] ] (* end SWITCH *) ); (* end function plotSoln *) End[]; EndPackage[];