Guest User

Untitled

a guest
Jun 28th, 2019
176
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.  
  4. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
  5. GMe               = 3.986E+14
  6.  
  7. fname = 'starlink 010.txt'
  8. with open (fname, 'r') as infile:
  9.     lines = infile.readlines()
  10.  
  11. L0s, L1s, L2s = lines[0::3], lines[1::3], lines[2::3]
  12. triples       = zip(L0s, L1s, L2s)
  13.  
  14. class Sat(object):
  15.     def __init__(self, triple):
  16.         L0, L1, L2 = triple
  17.         self.L0 = L0
  18.         self.L1 = L1
  19.         self.L2 = L2
  20.         self.name = L0.strip()
  21.         self.satnum   = int(L1[2:7])
  22.         self.intdesig = L1[9:17]
  23.         self.ecc      = float('0.' + L2[26:33])
  24.         self.revs     = float(L2[52:63])
  25.         self.T        = 24*3600. / self.revs
  26.         self.a        = (self.T**2 * GMe / twopi**2)**(1./3)
  27.         self.alt_km   = (self.a-6378137.) * 0.001
  28.         self.year     = 2000 + int(L1[18:20])
  29.         if self.year  >= 2057:
  30.             self.year -= 100
  31.         self.daynum   = float(L1[20:32])
  32.  
  33. sats  = [Sat(triple) for triple in triples]
  34.  
  35. objs   = [sat for sat in sats if 'object' in sat.name.lower()]
  36. debs   = [sat for sat in sats if 'deb'    in sat.name.lower()]
  37. slinks = [sat for sat in sats if 'star'   in sat.name.lower()]
  38.  
  39. obj_pairs_10   = [(sat.alt_km, sat.ecc) for sat in objs]
  40. slink_pairs_10 = [(sat.alt_km, sat.ecc) for sat in slinks]
  41. deb_pairs_10   = [(sat.alt_km, sat.ecc) for sat in debs]
  42.  
  43. # next one
  44.  
  45. fname = 'starlink 016.txt'
  46. with open (fname, 'r') as infile:
  47.     lines = infile.readlines()
  48.  
  49. L0s, L1s, L2s = lines[0::3], lines[1::3], lines[2::3]
  50. triples       = zip(L0s, L1s, L2s)
  51.  
  52. class Sat(object):
  53.     def __init__(self, triple):
  54.         L0, L1, L2 = triple
  55.         self.L0 = L0
  56.         self.L1 = L1
  57.         self.L2 = L2
  58.         self.name = L0.strip()
  59.         self.satnum   = int(L1[2:7])
  60.         self.intdesig = L1[9:17]
  61.         self.ecc      = float('0.' + L2[26:33])
  62.         self.revs     = float(L2[52:63])
  63.         self.T        = 24*3600. / self.revs
  64.         self.a        = (self.T**2 * GMe / twopi**2)**(1./3)
  65.         self.alt_km   = (self.a-6378137.) * 0.001
  66.         self.year     = 2000 + int(L1[18:20])
  67.         if self.year  >= 2057:
  68.             self.year -= 100
  69.         self.daynum   = float(L1[20:32])
  70.  
  71. sats  = [Sat(triple) for triple in triples]
  72.  
  73. objs   = [sat for sat in sats if 'object' in sat.name.lower()]
  74. debs   = [sat for sat in sats if 'deb'    in sat.name.lower()]
  75. slinks = [sat for sat in sats if 'star'   in sat.name.lower()]
  76.  
  77. obj_pairs_16   = [(sat.alt_km, sat.ecc) for sat in objs]
  78. slink_pairs_16 = [(sat.alt_km, sat.ecc) for sat in slinks]
  79. deb_pairs_16   = [(sat.alt_km, sat.ecc) for sat in debs]
  80.  
  81. if True:
  82.     fig = plt.figure()
  83.     ax1  = fig.add_subplot(2, 1, 1)
  84.     ax2  = fig.add_subplot(2, 1, 2)
  85.  
  86.     alts, eccs = np.array([zip(*obj_pairs_10)])[0]  # objects (not yet identified)
  87.     ax1.plot(alts, eccs, 'ok', markersize=3)
  88.  
  89.     alts, eccs = np.array([zip(*deb_pairs_10)])[0]
  90.     ax1.plot(alts, eccs, 'o', markeredgecolor='b', markerfacecolor='y',
  91.              markersize=8)
  92.  
  93.     alts, eccs = np.array([zip(*slink_pairs_16)])[0]  # STARLINKS (identified)
  94.     ax2.plot(alts, eccs, 'ok', markersize=3)
  95.  
  96.     alts, eccs = np.array([zip(*deb_pairs_16)])[0]
  97.     ax2.plot(alts, eccs, 'o', markeredgecolor='b', markerfacecolor='y',
  98.              markersize=8)
  99.  
  100.     ax1.set_xlabel('altitude (km)', fontsize=16)
  101.     ax1.set_ylabel('eccentricity', fontsize=16)
  102.     ax2.set_xlabel('altitude (km)', fontsize=16)
  103.     ax2.set_ylabel('eccentricity', fontsize=16)
  104.  
  105.     ax1.set_xlim(390, 570)
  106.     ax2.set_xlim(390, 570)
  107.  
  108.     ax1.set_ylim(0.0, 0.0012)
  109.     ax2.set_ylim(0.0, 0.0012)
  110.  
  111.     ax1.tick_params(axis='both', labelsize=14)
  112.     ax2.tick_params(axis='both', labelsize=14)
  113.  
  114.     ax1.set_title('"Starlink 60" plus debris around 01-06-2019 00:00 UTC', fontsize=20)
  115.     ax2.set_title('"Starlink 60" plus debris around 29-06-2019 00:00 UTC', fontsize=20)
  116.  
  117.     plt.show()
RAW Paste Data