Guest User

Untitled

a guest
May 14th, 2018
470
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # Install with:
  2. #   conda install -c conda-forge jplephem
  3. #   conda install -c conda-forge sgp4
  4. #   conda install -c conda-forge skyfield
  5. #
  6. # Download JPL ephemeris (https://en.wikipedia.org/wiki/Jet_Propulsion_Laboratory_Development_Ephemeris) from
  7. #   wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430.bsp
  8. #
  9.  
  10. from skyfield.api import load, EarthSatellite
  11. from skyfield.timelib import Time
  12. import numpy as np
  13. import matplotlib.pyplot as plt
  14. from mpl_toolkits.mplot3d import Axes3D
  15.  
  16. def makecubelimits(axis, centers=None, hw=None):
  17.     lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
  18.     if centers == None:
  19.         centers = [0.5*sum(pair) for pair in lims]
  20.  
  21.     if hw == None:
  22.         widths  = [pair[1] - pair[0] for pair in lims]
  23.         hw      = 0.5*max(widths)
  24.         ax.set_xlim(centers[0]-hw, centers[0]+hw)
  25.         ax.set_ylim(centers[1]-hw, centers[1]+hw)
  26.         ax.set_zlim(centers[2]-hw, centers[2]+hw)
  27.         print("hw was None so set to:", hw)
  28.     else:
  29.         try:
  30.             hwx, hwy, hwz = hw
  31.             print("ok hw requested: ", hwx, hwy, hwz)
  32.  
  33.             ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
  34.             ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
  35.             ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
  36.         except:
  37.             print("nope hw requested: ", hw)
  38.             ax.set_xlim(centers[0]-hw, centers[0]+hw)
  39.             ax.set_ylim(centers[1]-hw, centers[1]+hw)
  40.             ax.set_zlim(centers[2]-hw, centers[2]+hw)
  41.  
  42.     return centers, hw
  43.  
  44. TLE = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
  45. 2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
  46. L1, L2 = TLE.splitlines()
  47.  
  48.  
  49. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
  50. degs, rads = 180/pi, pi/180
  51.  
  52. data = load('de430.bsp')
  53. ts   = load.timescale()
  54.  
  55. planets = load('de430.bsp')
  56. earth   = planets['earth']
  57.  
  58. Roadster = EarthSatellite(L1, L2)
  59.  
  60. print(Roadster.epoch.tt)
  61. hours = np.arange(0, 3, 0.01)
  62.  
  63. time = ts.utc(2018, 2, 7, hours)
  64.  
  65. Rpos    = Roadster.at(time).position.km
  66. Rposecl = Roadster.at(time).ecliptic_position().km
  67.  
  68. print(Rpos.shape)
  69.  
  70. re = 6378.
  71.  
  72. theta = np.linspace(0, twopi, 201)
  73. cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
  74. lon0 = re*np.vstack((cth, zth, sth))
  75. lons = []
  76. for phi in rads*np.arange(0, 180, 15):
  77.     cph, sph = [f(phi) for f in (np.cos, np.sin)]
  78.     lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
  79.                      lon0[1]*cph + lon0[0]*sph,
  80.                      lon0[2]) )
  81.     lons.append(lon)
  82.  
  83. lat0 = re*np.vstack((cth, sth, zth))
  84. lats = []
  85. for phi in rads*np.arange(-75, 90, 15):
  86.     cph, sph = [f(phi) for f in (np.cos, np.sin)]
  87.     lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
  88.     lats.append(lat)
  89.  
  90. if True:    
  91.     fig = plt.figure(figsize=[10, 8])  # [12, 10]
  92.  
  93.     ax  = fig.add_subplot(1, 1, 1, projection='3d')
  94.  
  95.     x, y, z = Rpos
  96.     ax.plot(x, y, z)
  97.     for x, y, z in lons:
  98.         ax.plot(x, y, z, '-k')
  99.     for x, y, z in lats:
  100.         ax.plot(x, y, z, '-k')
  101.  
  102.     centers, hw = makecubelimits(ax)
  103.  
  104.     print("centers are: ", centers)
  105.     print("hw is:       ", hw)
  106.  
  107.     plt.show()
  108.  
  109. r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
  110. alt_roadster = r_Roadster - re
  111.  
  112. if True:
  113.     plt.figure()
  114.     plt.plot(hours, r_Roadster)
  115.     plt.plot(hours, alt_roadster)
  116.     plt.xlabel('hours', fontsize=14)
  117.     plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
  118.     plt.show()
RAW Paste Data