# Untitled

a guest
May 4th, 2019
107
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. import matplotlib.pyplot as plt
3. from mpl_toolkits.mplot3d import Axes3D
4.
5. linez = []
6.
7. JDs, positions, velocities = [], [], []
8.
9. names = ["Sun", "Mercury", "Venus", "Earth", "Moon", "Mars",
10.          "Jupiter", "Saturn", "Uranus", "Neptune", "PlutoBary",
11.          "Voyager 1", "Voyager 2"]
12.
13. for name in names:
14.
15.     fname = "horizons_results_1dayLONG " + name + ".txt"
16.
17.     print "trying: ", fname
18.
19.     with open(fname, "r") as infile:
20.
21.         lines = infile.read().splitlines()  # [-631:-34]
22.
23.         iSOE = [i for i, line in enumerate(lines) if "\$\$SOE" in line][0]
24.         iEOE = [i for i, line in enumerate(lines) if "\$\$EOE" in line][0]
25.
26.         print iSOE, iEOE, lines[iSOE], lines[iEOE]
27.
28.         linez.append(lines)
29.         lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
30.         JD  = np.array([float(line[0]) for line in lines])
31.         pos = np.array([[float(item) for item in line[2:5]] for line in lines])
32.         vel = np.array([[float(item) for item in line[5:8]] for line in lines])
33.
34.         JDs.append(JD)
35.         positions.append(pos)
36.         velocities.append(vel)
37.
38. JDzero = min([JD[0] for JD in JDs])
39.
40. JDoffs = [JD[0] - JDzero for JD in JDs]
41.
42. JDs[-2] = np.hstack((np.zeros(3170),JDs[-2]))
43. JDs[-1] = np.hstack((np.zeros(3154),JDs[-1]))
44. JDs[-2][:3170] = np.nan
45. JDs[-1][:3154] = np.nan
46.
47. positions[-2] = np.vstack((np.zeros(3*3170).reshape(-1, 3), positions[-2]))
48. positions[-1] = np.vstack((np.zeros(3*3154).reshape(-1, 3), positions[-1]))
49. positions[-2][:3170]  = np.nan
50. positions[-1][:3154]  = np.nan
51.
52. velocities[-2] = np.vstack((np.zeros(3*3170).reshape(-1, 3), velocities[-2]))
53. velocities[-1] = np.vstack((np.zeros(3*3154).reshape(-1, 3), velocities[-1]))
54. velocities[-2][:3170]  = np.nan
55. velocities[-1][:3154]  = np.nan
56.
57.
58. if True:
59.     fig = plt.figure(figsize=[12, 10])  # [12, 10]
60.     ax  = fig.add_subplot(1, 1, 1, projection='3d')
61.
62.     for pos in positions:
63.         x, y, z = pos.T
64.         ax.plot(x, y, z)
65.     plt.show()
66.
67.
68. if True:
69.
70.     ep, sp, V1p, V2p = [positions[i] for i in [3, 0, -2, -1]]
71.
72.     def angsep(a, b, c):
73.
74.         ba, ca = b-a, c-a
75.         rba, rca = np.sqrt( (ba**2).sum(axis=1) ), np.sqrt( (ca**2).sum(axis=1) )
76.         nba, nca = ba/rba[:, None], ca/rca[:, None]
77.         dot = (nba * nca).sum(axis=1)
78.
79.         return np.arccos(dot)
80.
81.     an1 = angsep(V1p, sp, ep)
82.     an2 = angsep(V2p, sp, ep)
83.
84. JD  = JDs[0]
85. JD0 = JD[0]
86.
87. span = JD[-1] - JD[0]
88. years = 1969.0 + (2031.0 - 1969.0) * (JD-JD0) / span
89.
90.
91. if False:
92.     plt.figure()
93.     plt.subplot(2,1,1)
94.     plt.plot(years, (180./np.pi)*an1)
95.     plt.plot(years, (180./np.pi)*an2)
96.     plt.xlim(1975, 2020)
97.     plt.ylim(0, 20.)
98.     plt.subplot(2,1,2)
99.     plt.plot(years, (180./np.pi)*an1)
100.     plt.plot(years, (180./np.pi)*an2)
101.     plt.xlim(1975, 2020)
102.     plt.ylim(0, 0.5)
103.     plt.show()
104.
105.
106. if True:
107.
108.     JDnow = 2457591.0  # UTC noon July 21, 2016
109.
110.     iJD = np.argmax(JD>JDnow)
111.
112.     pstoday = np.array(positions)[:, :iJD]
113.
114.     maxes = np.nanmax(pstoday.reshape(-1, 3), axis=0)
115.     mins  = np.nanmin(pstoday.reshape(-1, 3), axis=0)
116.
117.     centerz = np.array([0.0, 0.5 * (maxes + mins)[1], 0.0])
118.
119.     hw = np.nanmax(  np.abs(pstoday-centerz).reshape(-1, 3), axis=0).max()
120.
121.     diffs = maxes - mins
122.     maxdiff = diffs.max()
123.     hw = 0.5 * maxdiff
124.     xlim, ylim, zlim = [(cen-hw, cen+hw) for cen in centerz]
125.
126.
127. if True:
128.
129.     nframes = 60
130.
131.     theta  = 0.5 * np.pi * np.linspace(-1, 1, nframes)
132.     things = 0.5 * (1. + np.sin(theta))
133.
134.     more_things = np.hstack((15*[0], things, 15*[1]))
135.
136.     for i, thing in enumerate(more_things):
137.
138.         fig = plt.figure(figsize=[12, 10])  # [12, 10]
139.         ax  = fig.add_subplot(1, 1, 1, projection='3d')
140.
141.         az = 30
142.         el = 50 * thing
143.         ax.view_init(el, az)
144.
145.         for pos in pstoday:
146.             x, y, z = pos.T
147.             ax.plot(x, y, z)
148.         ax.set_xlim(xlim)
149.         ax.set_ylim(ylim)
150.         ax.set_zlim(zlim)
151.
152.         for pos in pstoday:
153.             x, y, z = pos.T
154.             ax.plot(x, y, zlim[0]*np.ones_like(z), '-', color='0.7')
155.         ax.set_xlim(xlim)
156.         ax.set_ylim(ylim)
157.         ax.set_zlim(zlim)
158.
159.         for pos in pstoday:
160.             x, y, z = pos[-1].T
161.             ax.scatter(x, y, z, 'ok')
162.         ax.set_xlim(xlim)
163.         ax.set_ylim(ylim)
164.         ax.set_zlim(zlim)
165.
166.         ax.set_axis_off()
167.
168.         plt.savefig('VoyagerQQQ_1_2_' + str(1000+i)[1:] )
169.
170.         plt.close()
171.         # plt.show()
RAW Paste Data