Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib.pyplot as plt
- from mpl_toolkits.mplot3d import Axes3D
- linez = []
- JDs, positions, velocities = [], [], []
- names = ["Sun", "Mercury", "Venus", "Earth", "Moon", "Mars",
- "Jupiter", "Saturn", "Uranus", "Neptune", "PlutoBary",
- "Voyager 1", "Voyager 2"]
- for name in names:
- fname = "horizons_results_1dayLONG " + name + ".txt"
- print "trying: ", fname
- with open(fname, "r") as infile:
- lines = infile.read().splitlines() # [-631:-34]
- iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
- iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]
- print iSOE, iEOE, lines[iSOE], lines[iEOE]
- linez.append(lines)
- lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
- JD = np.array([float(line[0]) for line in lines])
- pos = np.array([[float(item) for item in line[2:5]] for line in lines])
- vel = np.array([[float(item) for item in line[5:8]] for line in lines])
- JDs.append(JD)
- positions.append(pos)
- velocities.append(vel)
- JDzero = min([JD[0] for JD in JDs])
- JDoffs = [JD[0] - JDzero for JD in JDs]
- JDs[-2] = np.hstack((np.zeros(3170),JDs[-2]))
- JDs[-1] = np.hstack((np.zeros(3154),JDs[-1]))
- JDs[-2][:3170] = np.nan
- JDs[-1][:3154] = np.nan
- positions[-2] = np.vstack((np.zeros(3*3170).reshape(-1, 3), positions[-2]))
- positions[-1] = np.vstack((np.zeros(3*3154).reshape(-1, 3), positions[-1]))
- positions[-2][:3170] = np.nan
- positions[-1][:3154] = np.nan
- velocities[-2] = np.vstack((np.zeros(3*3170).reshape(-1, 3), velocities[-2]))
- velocities[-1] = np.vstack((np.zeros(3*3154).reshape(-1, 3), velocities[-1]))
- velocities[-2][:3170] = np.nan
- velocities[-1][:3154] = np.nan
- if True:
- fig = plt.figure(figsize=[12, 10]) # [12, 10]
- ax = fig.add_subplot(1, 1, 1, projection='3d')
- for pos in positions:
- x, y, z = pos.T
- ax.plot(x, y, z)
- plt.show()
- if True:
- ep, sp, V1p, V2p = [positions[i] for i in [3, 0, -2, -1]]
- def angsep(a, b, c):
- ba, ca = b-a, c-a
- rba, rca = np.sqrt( (ba**2).sum(axis=1) ), np.sqrt( (ca**2).sum(axis=1) )
- nba, nca = ba/rba[:, None], ca/rca[:, None]
- dot = (nba * nca).sum(axis=1)
- return np.arccos(dot)
- an1 = angsep(V1p, sp, ep)
- an2 = angsep(V2p, sp, ep)
- JD = JDs[0]
- JD0 = JD[0]
- span = JD[-1] - JD[0]
- years = 1969.0 + (2031.0 - 1969.0) * (JD-JD0) / span
- if False:
- plt.figure()
- plt.subplot(2,1,1)
- plt.plot(years, (180./np.pi)*an1)
- plt.plot(years, (180./np.pi)*an2)
- plt.xlim(1975, 2020)
- plt.ylim(0, 20.)
- plt.subplot(2,1,2)
- plt.plot(years, (180./np.pi)*an1)
- plt.plot(years, (180./np.pi)*an2)
- plt.xlim(1975, 2020)
- plt.ylim(0, 0.5)
- plt.show()
- if True:
- JDnow = 2457591.0 # UTC noon July 21, 2016
- iJD = np.argmax(JD>JDnow)
- pstoday = np.array(positions)[:, :iJD]
- maxes = np.nanmax(pstoday.reshape(-1, 3), axis=0)
- mins = np.nanmin(pstoday.reshape(-1, 3), axis=0)
- centerz = np.array([0.0, 0.5 * (maxes + mins)[1], 0.0])
- hw = np.nanmax( np.abs(pstoday-centerz).reshape(-1, 3), axis=0).max()
- diffs = maxes - mins
- maxdiff = diffs.max()
- hw = 0.5 * maxdiff
- xlim, ylim, zlim = [(cen-hw, cen+hw) for cen in centerz]
- if True:
- nframes = 60
- theta = 0.5 * np.pi * np.linspace(-1, 1, nframes)
- things = 0.5 * (1. + np.sin(theta))
- more_things = np.hstack((15*[0], things, 15*[1]))
- for i, thing in enumerate(more_things[::10]):
- fig = plt.figure(figsize=[12, 10]) # [12, 10]
- ax = fig.add_subplot(1, 1, 1, projection='3d')
- az = 30
- el = 50 * thing
- ax.view_init(el, az)
- for pos in pstoday:
- x, y, z = pos.T
- ax.plot(x, y, z)
- ax.set_xlim(xlim)
- ax.set_ylim(ylim)
- ax.set_zlim(zlim)
- for pos in pstoday:
- x, y, z = pos.T
- ax.plot(x, y, zlim[0]*np.ones_like(z), '-', color='0.7')
- ax.set_xlim(xlim)
- ax.set_ylim(ylim)
- ax.set_zlim(zlim)
- for pos in pstoday:
- x, y, z = pos[-1].T
- ax.scatter(x, y, z, 'ok')
- ax.set_xlim(xlim)
- ax.set_ylim(ylim)
- ax.set_zlim(zlim)
- ax.set_axis_off()
- plt.savefig('VoyagerQQQ_1_2_' + str(1000+i)[1:] )
- plt.close()
- # plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement