

In [1]:

    
%matplotlib inline



In [2]:

    
from __future__ import division, print_function
from ivisual import *
from math import *
import matplotlib.pyplot as plt



In [7]:

    
#disk
M=2.94
R=0.6048/2
m1=0.52
m2=0.15*m1
I=R*R*(1/2*M+m1+m2)
d=(m1-m2)/(m1+m2)*R
print("d = ",d)

#initial angle with respect to -y axis
phi_deg=-90
phi=phi_deg*pi/180 #phi in rad

#initial angular velocity
phidot=19  #rad/s

#gravitational field in N/kg
g=9.8

t=0
dt=0.001

#position vector
r=d*vector(sin(phi),-cos(phi),0)

#period
T=0
ti=0
atpeak=False

tlist=[]
omegalist=[]
philist=[]

while t<1.25:
    rate(100)
    
    omegai=phidot
    
    #update phidotdot, phidot, phi
    phidotdot=-(M+m1+m2)*g*d/I*sin(phi)
    phidot=phidot+phidotdot*dt
    phi=phi+phidot*dt
    
    #update r
    r=R*vector(sin(phi),-cos(phi),0)
    
    #find period
    if(phidotdot>0 and (phidot-omegai)<0.001 and phidot>19 and atpeak==False):
        ti=t
        print("t_i = ",t)
        atpeak=True
    if(phidotdot>0 and (phidot-omegai)<0.001 and phidot>19 and atpeak==True):
        T=t-ti
        ti=t
        print("T = ",T)
    
    
    tlist.append(t)
    omegalist.append(phidot)
    philist.append(phi)
    
    t=t+dt









    



d =  0.223513043478
t_i =  0.077
T =  0.0
T =  0.333
T =  0.334
T =  0.333
T =  0.001



In [8]:

    
plt.title('omega vs. t')
plt.xlabel('t (s)')
plt.ylabel('omega (rad/s)')
plt.plot(tlist,omegalist,'m.')
plt.show()

Oscillating Disk