2D and 3D simulation of a pendulum that swings on a hoop which rotates with varying ω
by Nasser Abbasi 7/14/07
see http://12000.org/my_courses/FULLERTON_COURSES/summer_2007/Math_503/HW/HW6/index.htm
for analytical solution.

getBobCoordinates2D[θ_,r_]:=Module[{q,x,y,θNew},{
If[θ>=0,θNew=Mod[θ,2 Pi],{θNew=-θ;θNew=Mod[θNew,2 Pi];θNew=-θNew;θNew=θNew+2 Pi}];
q=First[getQuadrent[θNew]];

If[q==1,{x=r Sin[θNew],y=-r Cos[θNew]}];
If[q==2,{θNew=θNew-Pi/2;x=r Cos[θNew],y=r Sin[θNew]}];
If[q==3,{θNew=θNew-Pi;x=-r Sin[θNew],y=r Cos[θNew]}];
If[q==4,{θNew=θNew-(3/2)Pi;x=-r Cos[θNew],y=-r Sin[θNew]}];
{x,y}}
]

getBob[ϕ_,r_,θ_,is3D_]:=Module[{normalizedPhi,q,coords,l},
{
If[is3D,{coords={r Cos[θ] Sin[ϕ],r Sin[θ]Sin[ϕ],r Cos[ϕ]};l=Line[{{0,0,0},coords}];{Sphere[coords,.1 r],l}
},Line[{{0,0},First@getBobCoordinates2D[ϕ,r]}]]
}
]

getQuadrent[a_]:=Module[{q},
{If[a≥0&&a≤Pi/2,q=1];
If[a>Pi/2&&a≤Pi,q=2];
If[a>Pi&&a≤(3/2)Pi,q=3];
If[a>(3/2)Pi&&a<2Pi,q=4];
q}
]

getCurrentCircle[θ_,r_,is3D_]:=Module[{ϕ},
If[is3D,Table[{r Cos[θ]Sin[ϕ],r Sin[θ] Sin[ϕ],r Cos[ϕ]},{ϕ,0.,2π,0.01 π}],
Circle[{0,0},r]]
]

getXYLines[or_,r_]:=Module[{},Union[{Dashed,Line[{or,{-2 r,0}}]},{Dashed,Line[{or,{2 r,0}}]},{Dashed,Line[{or,{0,2 r}}]},{Dashed,Line[{or,{0,-2 r}}]}]]

getXYZLines[or_,r_]:=Module[{},
{ {Dashed,Line[{or,{-2 r,0,0}}]},{Dashed,Line[{or,{2 r,0,0}}]},
{Dashed,Line[{or,{0,2 r,0}}]},{Dashed,Line[{or,{0,-2 r,0}}]},
{Red,Line[{or,{0,0,2 r}}]},{Red,Line[{or,{0,0,-2 r}}]}}
]

process2[s_,ϕ_,ts_,θ_,r_,is3D_,ω_]:=Module[{p,ϕNew,axesLines=0,or3={0,0,0},or2={0,0},t=0,bobLine=0,eq},
If[is3D,{axesLines=getXYZLines[or3,r];ϕNew=π-ϕ;Graphics3D[{axesLines,{Thick,Line[getCurrentCircle[Mod[θ,2Pi],r,is3D]]},getBob[ϕNew,r,θ,is3D]},Boxed->False, ImageSize->{250,250},
ViewPoint->{1.391,-3.001,.713},ViewAngle->20Degree,PlotRange->{{-2r,2r},{-2r,2r},{-2r,2r}}]
},
{axesLines=getXYLines[or2,r];ϕNew=ϕ;p=Plot[(180./Pi) s[t],{t,0,If[ω==0 Or ω<0.5,10,3((2. π)/ω)]},AxesOrigin->{0,0},AxesLabel->{"time (sec)","ϕ(t) in degrees"}];
GraphicsRow[{p,Graphics[{axesLines,getCurrentCircle[0,r,is3D],getBob[ϕNew,r,θ,is3D],PointSize[.03],Red,Point[ getBobCoordinates2D[ϕNew,r] ]},PlotRange->{{-2r,2r},{-2r,2r}},ImageSize->{200,200}]}]
}
]
]

process[ω_,r_,ϕ0_,v0_,is3D_]:=Module[{s,g=9.8,eq,period,t,ϕ,maxSteps,currentTime,maxTime},{
If[ω≠0,period=((2. π)/ω),period=100.];
maxTime=100 period; (*seconds*)
eq= ϕ''[t]+Sin[ϕ[t]]((g/r)-(ω^2)Cos[ϕ[t]]);

s=ϕ/.First[NDSolve[{eq==0,ϕ[0]==ϕ0 Degree,ϕ'[0]==v0 Degree},ϕ,{t,0,maxTime},
MaxSteps->Infinity
]];

Dynamic[Animate[
First@ process2[   s, s[currentTime],currentTime,If[is3D,ω currentTime,0],r,is3D,ω],{{currentTime,0,"current time in sec"},0,maxTime,period/1000.,AnimationRate->0.05 period  If[ω==0,.05,ω],RefreshRate->1,AnimationRunning->True}]]
}
]

Manipulate[
First@ process[ω ,r,ϕ0,v0,is3D],
{{ω,1,"ω in rad/sec. Angular speed of hoop"},0,50,.5,ContinuousAction->True,Appearance->"Labeled"},
{{r,1,"R in meters. Radius of hoop"},1,10,Appearance->"Labeled"},{{ϕ0,45,"initial position in degrees ϕ(0)"},0,180,1,Appearance->"Labeled"},
{{v0,0,"initial bob angular speed in degree/sec ϕ'(0)"},0,100,1,Appearance->"Labeled"},
{{is3D,False,"display type"},{False->"2D",True->"3D"}},
FrameLabel->"Pendulum on rotating hoop: ϕ''[t]+Sin[ϕ[t]]((g/r)-(ω^2)Cos[ϕ[t]])",AutorunSequencing->{{1,300},SaveDefinitions->True}
]

u=Export["hoop.swf",u]