import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
M=14.73
m=1
g=9.81
def Atwood(y,t,M,m,g):
theta, omega, radius, Rho = y
dydt = [omega, (1/radius)*(-g*(np.sin(theta))-2*(Rho*omega)), Rho, (1/(M+m))*(m*radius*(omega**2)-M*g+m*g*(np.cos(theta)))]
return dydt
y0 = [np.pi/2 ,0 ,1 ,0]
t = np.linspace(0,100,1000000)
sul = odeint(Atwood, y0, t, args=(M,m,g))
plt.plot(sul[:,0], sul[:,1], 'r')
plt.ylabel('Angular velocity (radians per second)')
plt.xlabel('Angular displacement (radians)')
plt.title('Fig 16 - μ = 14.73')
plt.grid()
plt.show()
