# Untitled

a guest
Apr 12th, 2018
168
0
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.
10.
12. earth   = data['earth']
14.
15. fname = "SMM tles 1989.txt"
16.
17. with open(fname, 'r') as infile:
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()