Guest User

Untitled

a guest
Apr 12th, 2018
132
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 skyfield.api import Loader, Topos, EarthSatellite
  4.  
  5. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]  # parentheses for Py3
  6. degs, rads = 180/pi, pi/180
  7.  
  8. load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
  9. ts   = load.timescale()
  10.  
  11. data    = load('de421.bsp')
  12. earth   = data['earth']
  13. ts      = load.timescale()
  14.  
  15. fname = "SMM tles 1989.txt"
  16.  
  17. with open(fname, 'r') as infile:
  18.     lines = infile.readlines()
  19.  
  20. L1s, L2s = lines[0::2], lines[1::2]
  21. pairs    = zip(L1s, L2s)
  22.  
  23. alts, daynums = [], []
  24. for i, (L1, L2) in enumerate(pairs):
  25.     if not i%20:
  26.         print i,
  27.     daynum = float(L1[20:32])
  28.     revs   = float(L2[52:64])
  29.     T_mins = 24.*60/revs
  30.     dayfraction = daynum % 1.0
  31.     dayinteger  = int(daynum)
  32.     minutes_epoch = 24. * 60 * dayfraction
  33.  
  34.     minutes = np.linspace(0, T_mins, 101)[:-1] + minutes_epoch
  35.     times   = ts.utc(1989, 1, dayinteger, 0, minutes)
  36.  
  37.     sat     = EarthSatellite(L1, L2)
  38.     pos     = sat.at(times).position.km
  39.     r       = np.sqrt((pos**2).sum(axis=0))
  40.     alt     = r - 6378.13
  41.     alts.append(alt)
  42.     daynums.append(daynum)
  43.  
  44. mins    = np.array([alt.min()  for alt in alts])
  45. maxs    = np.array([alt.max()  for alt in alts])
  46. means   = np.array([alt.mean() for alt in alts])
  47. daynums = np.array(daynums)
  48.  
  49. if True:
  50.     plt.figure()
  51.  
  52.     fs = 16
  53.  
  54.     plt.subplot(2, 1, 1)
  55.    
  56.     plt.plot(daynums, mins)
  57.     plt.plot(daynums, maxs)
  58.     plt.xlabel('day number in 1989', fontsize=fs)
  59.     plt.ylabel('altitude (km)', fontsize=fs)
  60.  
  61.     plt.subplot(2, 1, 2)
  62.  
  63.     plt.plot(daynums, mins)
  64.     plt.plot(daynums, maxs)
  65.     plt.plot(daynums, means, '-k', linewidth=1.5)
  66.  
  67.     plt.xlim(0, 150)
  68.     plt.ylim(425, 455)
  69.     plt.xlabel('day number in 1989', fontsize=fs)
  70.     plt.ylabel('altitude (km)', fontsize=fs)
  71.  
  72.     plt.suptitle('Solar Maximum Mission 1980-014A, 11703', fontsize=fs)
  73.  
  74.     plt.show()
RAW Paste Data