import numpy as np import matplotlib.pyplot as plt from scipy import integrate g = 9.8 # kg m s^-2 beta = 0.05 # kg m^-1 gamma = 0.0001 m = 100. # kg # Question 1 def Stokes(eqs, t, beta, g, m): z = eqs[0] v = eqs[1] dz = -v dv = g - (beta/m)*v return (dz, dv) # Question 2 initial = (39000., 0.) time = np.arange(0, 300, 1) FelixStokes = integrate.odeint(Stokes, initial, time, args = (beta, g, m)) # Question 3 def Rayleigh(eqs, t, beta, g, m, gamma): z = eqs[0] v = eqs[1] dz = -v dv = g - (beta/m)*v - gamma*v**2 return (dz, dv) FelixRayleigh = integrate.odeint(Rayleigh, initial, time, args = (beta, g, m, gamma)) # Question 4 plt.plot(time, FelixStokes.T[1]) plt.plot(time, FelixRayleigh.T[1]) plt.legend(('Stokes','Rayleigh')) plt.xlabel('time') plt.ylabel('V') plt.savefig('FelixV.png') plt.close() # Question 5 plt.plot(time, FelixStokes.T[0]) plt.plot(time, FelixRayleigh.T[0]) plt.legend(('Stokes','Rayleigh')) plt.xlabel('time') plt.ylabel('Z') plt.savefig('FelixZ.png') plt.close() # Question 6 def ground(time, Z): for i in range(len(Z)): if Z[i] < 0: t = (time[i]*Z[i-1]-time[i-1]*Z[i])/(Z[i-1]-Z[i]) return t print ground(time, FelixStokes.T[0]), ground(time, FelixRayleigh.T[0]) # Question 7 print np.max(FelixStokes.T[1]), np.max(FelixRayleigh.T[1])