Guest User

Untitled

a guest
Jan 22nd, 2019
35
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2.  
  3.    
  4. def deriv(X, t):
  5.     th, thd = X
  6.     thdd = const * np.sin(2*th)
  7.     return np.hstack((thd, thdd))
  8.  
  9. import numpy as np
  10. import matplotlib.pyplot as plt
  11. from scipy.integrate import odeint as ODEint
  12.  
  13. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
  14. degs, rads = 180/pi, pi/180
  15.  
  16. GMe   = 3.986E+14
  17. Rc    = 1000. * (6378.137 + 400.)
  18. const = -(3.*GMe) / (2.*Rc**3)
  19. T     = twopi * np.sqrt(Rc**3/GMe)
  20.  
  21. theta_zeros = rads * np.arange(-80, 81, 20.0)
  22.  
  23. times   = np.linspace(0, 2*T, 1001)
  24. minutes = times/60.
  25. theta_orbit = twopi * times/T
  26. theta_dot_orbit = twopi / T
  27.  
  28. answers = []
  29. for theta in theta_zeros:
  30.     X0 = np.hstack([theta, -theta_dot_orbit])
  31.     answer, info = ODEint(deriv, X0, times, full_output=True)
  32.     answers.append(answer)
  33.  
  34. print answer.shape
  35.  
  36. thetas     = [a.T[0] + theta_orbit     for a in answers]
  37. theta_dots = [a.T[1] + theta_dot_orbit for a in answers]
  38.  
  39. answer = np.stack([a.T for a in answers])
  40.  
  41. if True:
  42.     plt.figure()
  43.     plt.subplot(2, 1, 1)
  44.     for theta in thetas:
  45.         plt.plot(minutes, degs*theta)
  46.     plt.ylabel('theta (deg)', fontsize=16)
  47.     plt.xlim(0, 190)
  48.  
  49.     plt.subplot(2, 1, 2)
  50.     for theta_dot in theta_dots:
  51.         plt.plot(minutes, degs*60.*theta_dot)
  52.     plt.plot(minutes, np.zeros_like(minutes), '-k', linewidth=1.5)
  53.     plt.xlabel('time (minutes)', fontsize=16)
  54.     plt.ylabel('theta dot (deg/min)', fontsize=16)
  55.     plt.xlim(0, 190)
  56.     plt.show()
RAW Paste Data