Advertisement
Guest User

Untitled

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