import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
GMe               = 3.986E+14

fname = 'starlink 010.txt'
with open (fname, 'r') as infile:
    lines = infile.readlines()

L0s, L1s, L2s = lines[0::3], lines[1::3], lines[2::3]
triples       = zip(L0s, L1s, L2s)

class Sat(object):
    def __init__(self, triple):
        L0, L1, L2 = triple
        self.L0 = L0
        self.L1 = L1
        self.L2 = L2
        self.name = L0.strip()
        self.satnum   = int(L1[2:7])
        self.intdesig = L1[9:17]
        self.ecc      = float('0.' + L2[26:33])
        self.revs     = float(L2[52:63])
        self.T        = 24*3600. / self.revs
        self.a        = (self.T**2 * GMe / twopi**2)**(1./3)
        self.alt_km   = (self.a-6378137.) * 0.001
        self.year     = 2000 + int(L1[18:20])
        if self.year  >= 2057:
            self.year -= 100
        self.daynum   = float(L1[20:32])

sats  = [Sat(triple) for triple in triples]

objs   = [sat for sat in sats if 'object' in sat.name.lower()]
debs   = [sat for sat in sats if 'deb'    in sat.name.lower()]
slinks = [sat for sat in sats if 'star'   in sat.name.lower()]

obj_pairs_10   = [(sat.alt_km, sat.ecc) for sat in objs]
slink_pairs_10 = [(sat.alt_km, sat.ecc) for sat in slinks]
deb_pairs_10   = [(sat.alt_km, sat.ecc) for sat in debs]

# next one

fname = 'starlink 016.txt'
with open (fname, 'r') as infile:
    lines = infile.readlines()

L0s, L1s, L2s = lines[0::3], lines[1::3], lines[2::3]
triples       = zip(L0s, L1s, L2s)

class Sat(object):
    def __init__(self, triple):
        L0, L1, L2 = triple
        self.L0 = L0
        self.L1 = L1
        self.L2 = L2
        self.name = L0.strip()
        self.satnum   = int(L1[2:7])
        self.intdesig = L1[9:17]
        self.ecc      = float('0.' + L2[26:33])
        self.revs     = float(L2[52:63])
        self.T        = 24*3600. / self.revs
        self.a        = (self.T**2 * GMe / twopi**2)**(1./3)
        self.alt_km   = (self.a-6378137.) * 0.001
        self.year     = 2000 + int(L1[18:20])
        if self.year  >= 2057:
            self.year -= 100
        self.daynum   = float(L1[20:32])

sats  = [Sat(triple) for triple in triples]

objs   = [sat for sat in sats if 'object' in sat.name.lower()]
debs   = [sat for sat in sats if 'deb'    in sat.name.lower()]
slinks = [sat for sat in sats if 'star'   in sat.name.lower()]

obj_pairs_16   = [(sat.alt_km, sat.ecc) for sat in objs]
slink_pairs_16 = [(sat.alt_km, sat.ecc) for sat in slinks]
deb_pairs_16   = [(sat.alt_km, sat.ecc) for sat in debs]

if True:
    fig = plt.figure()
    ax1  = fig.add_subplot(2, 1, 1)
    ax2  = fig.add_subplot(2, 1, 2)

    alts, eccs = np.array([zip(*obj_pairs_10)])[0]  # objects (not yet identified)
    ax1.plot(alts, eccs, 'ok', markersize=3)

    alts, eccs = np.array([zip(*deb_pairs_10)])[0]
    ax1.plot(alts, eccs, 'o', markeredgecolor='b', markerfacecolor='y',
             markersize=8)

    alts, eccs = np.array([zip(*slink_pairs_16)])[0]  # STARLINKS (identified)
    ax2.plot(alts, eccs, 'ok', markersize=3)

    alts, eccs = np.array([zip(*deb_pairs_16)])[0]
    ax2.plot(alts, eccs, 'o', markeredgecolor='b', markerfacecolor='y',
             markersize=8)

    ax1.set_xlabel('altitude (km)', fontsize=16)
    ax1.set_ylabel('eccentricity', fontsize=16)
    ax2.set_xlabel('altitude (km)', fontsize=16)
    ax2.set_ylabel('eccentricity', fontsize=16)

    ax1.set_xlim(390, 570)
    ax2.set_xlim(390, 570)

    ax1.set_ylim(0.0, 0.0012)
    ax2.set_ylim(0.0, 0.0012)

    ax1.tick_params(axis='both', labelsize=14)
    ax2.tick_params(axis='both', labelsize=14)

    ax1.set_title('"Starlink 60" plus debris around 01-06-2019 00:00 UTC', fontsize=20)
    ax2.set_title('"Starlink 60" plus debris around 29-06-2019 00:00 UTC', fontsize=20)

    plt.show()