Guest User

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.  
  7. load  = Loader('~/Documents/fishing/SkyData')  # single instance for big files
  8. ts    = load.timescale()
  9. de421 = load('de421.bsp')
  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:
  16.     lines32 = infile.readlines()
  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:
  23.     lines33 = infile.readlines()
  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