Untitled

a guest
Jan 22nd, 2019
123
0
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.
29. for theta in theta_zeros:
30.     X0 = np.hstack([theta, -theta_dot_orbit])
31.     answer, info = ODEint(deriv, X0, times, full_output=True)
33.
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.
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()