from pylab import * #graph1 = gdisplay(title='Vertical Position vs. Time', xtitle ='t (sec)', ytitle ='y (m)', background=color.white) #funct1=gcurve(color=color.white) #floor = box (pos=(0,0,0), length=.4, height=0.05, width=.4, color=color.blue) L = 1 theta0 = -pi/5 #ball2=sphere() ball = array([L*sin(theta0), -L*cos(theta0),0]) m =.1 pivot=array([0,0,0]) k=10000 dt=0.001 g = array([0,-9.8,0]) def mag(vec1): #this function takes a vector (array) and returns a magnitude return sqrt(vec1[0]**2+vec1[1]**2+vec1[2]**2) F=m*g-k*(ball-pivot)*(mag(ball-pivot)-L)/mag(ball-pivot) ballp=array([0,0,0]) t=0 #set up plotting plotx=[] ploty=[] plotx2=[] ploty2=[] #here is the stuff to plot the analytical solution-ish theta = theta0 thetadot=0 thetaddot=-mag(g)*sin(theta)/L y2=-L*cos(theta) x2=L*sin(theta) while t<3: #rate(100) F=m*g-k*(ball-pivot)*(mag(ball-pivot)-L)/mag(ball-pivot) ballp=ballp + F*dt ball = ball +ballp*dt/m #euler method for analytical pendulum thetaddot=-mag(g)*sin(theta)/L thetadot=thetadot+thetaddot*dt theta=theta+thetadot*dt y2=-L*cos(theta) x2=L*sin(theta) ploty2=ploty2+[2*y2] t = t+dt plotx=plotx+[t] ploty=ploty+[ball[1]] plot(plotx, ploty, plotx, ploty2, linewidth=3) grid(color='b', linestyle='-', linewidth=0.5) title('Springy Pendulum') xlabel('time (sec)') ylabel('y position (m)') show()