# Untitled

a guest
Dec 27th, 2019
118
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 Topos, Loader, EarthSatellite
4.
5. # ugly quick script to make some quick ugly plots
6.
10. earth = de421['earth']
11.
12. #1 41332U 16009A   17001.16907386  .00001861  00000-0  66700-4 0  9997
13. #2 41332  97.4831  58.2881 0024160 250.2364 109.6271 15.29806456 50281
14.
15. with open('TLEs_41332.txt', 'r') as infile:
17.     L1_32, L2_32 = lines32[0::2], lines32[1::2]
18.     pairs32 = zip(L1_32, L2_32)
19.     years32 = [2000 + float(line[18:20]) for line in L1_32]
20.     days32  = [       float(line[20:32]) for line in L1_32]
21.
22. with open('TLEs_41333.txt', 'r') as infile:
24.     L1_33, L2_33 = lines33[0::2], lines33[1::2]
25.     pairs33 = zip(L1_33, L2_33)
26.     years33 = [2000 + float(line[18:20]) for line in L1_33]
27.     days33  = [       float(line[20:32]) for line in L1_33]
28.
29. apos32, peris32, yearz32 = [], [], []
30. for i, ((L1, L2), year, day) in enumerate(zip(pairs32, years32, days32)):
31.     days  = day + np.arange(120.)/60./24.
32.     yearz32.append(year + day/365.2564)
33.     times = ts.utc(year, 1, days, 0, 0)
34.     posns = EarthSatellite(L1, L2).at(times).position.km
35.     r     = np.sqrt((posns**2).sum(axis=0))
36.     apos32.append(r.max())
37.     peris32.append(r.min())
38.     if not i%10:
39.         print(i,)
40.
41. apos32, peris32, yearz32 = [np.array(x) for x in (apos32, peris32, yearz32)]
42.
43. apos33, peris33, yearz33 = [], [], []
44. for i, ((L1, L2), year, day) in enumerate(zip(pairs33, years33, days33)):
45.     days  = day + np.arange(120.)/60./24.
46.     yearz33.append(year + day/365.2564)
47.     times = ts.utc(year, 1, days, 0, 0)
48.     posns = EarthSatellite(L1, L2).at(times).position.km
49.     r     = np.sqrt((posns**2).sum(axis=0))
50.     apos33.append(r.max())
51.     peris33.append(r.min())
52.     if not i%10:
53.         print(i,)
54.
55. apos33, peris33, yearz33 = [np.array(x) for x in (apos33, peris33, yearz33)]
56.
57. if True:
58.     Re = 6378.137
59.     plt.figure()
60.     plt.subplot(2, 1, 1)
61.     plt.plot(yearz32, apos32-Re, '.', markersize=1)
62.     plt.plot(yearz32, peris32-Re, '.', markersize=1)
63.     plt.title('41332 quick-n-dirty')
64.     plt.subplot(2, 1, 2)
65.     plt.plot(yearz33, apos33-Re, '.', markersize=1)
66.     plt.plot(yearz33, peris33-Re, '.', markersize=1)
67.     plt.title('41333 quick-n-dirty')
68.     plt.show()
RAW Paste Data