Advertisement
Guest User

Untitled

a guest
May 4th, 2019
242
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.57 KB | None | 0 0
  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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement