import numpy as np import matplotlib.pyplot as plt from scipy import integrate # Question 1 def f(Y, t): x,y,vx,vy = Y d31 = ((x+0.25)**2+y**2)**1.5 M = 1. return np.array([vx,vy,-M*x/d31,-M*y/d31]) # Question 2 t0 = 0 tfinal = 1.339 y0 = np.array([0.3, 0., 0.0, 0.05]) # initial position and velocity N = 5000 time = np.linspace(t0, tfinal, N) planet = integrate.odeint(f, y0, time) # Question 3 x = planet.T[0] y = planet.T[1] vx = planet.T[2] vy = planet.T[3] plt.plot(x,y) plt.savefig('planet.png') plt.close() # Question 4 perp = x*vx+y*vy P = np.mean(perp) print P # Question 5 diff = np.diag(-1.0*np.ones(N), 0) + np.diag(1.0*np.ones(N-1), 1) # Question 6 delta = time[1]-time[0] ax = np.dot(diff, vx)/delta ay = np.dot(diff, vy)/delta # Question 7 r = np.sqrt(x**2+y**2) v = np.sqrt(vx**2+vy**2) a = np.sqrt(ax**2+ay**2) # Question 8 plt.plot(r,v) plt.plot(r,a) plt.xlabel(r) plt.legend(["v", "a"]) plt.savefig('phase.png') # Question 9 fp = open("planet.txt", "w") for i in range(len(time)): print >> fp, time[i], x[i], y[i], vx[i], vy[i] fp.close() # Question 10 parallel = x*ax + y*ay print np.mean(parallel)