import numpy as np
import matplotlib.pyplot as plt

GMe      = 3.986004418E+14
oneday   = 23*3600 + 56*60. + 4.09
twopi    = 2. * np.pi
onethird = 1./3.

fname191 = 'All TLEs 2018 day 191.txt'
fname192 = 'All TLEs 2018 day 192.txt'

with open (fname191, 'r') as infile:
    lines191 = infile.readlines()

with open (fname192, 'r') as infile:
    lines192 = infile.readlines()

L1s, L2s = lines191[0::2], lines191[1::2]
pairs191 = zip(L1s, L2s)

L1s, L2s = lines192[0::2], lines192[1::2]
pairs192 = zip(L1s, L2s)

data191  = []
for L1, L2 in pairs191:
    year = float(L1[18:20])
    day  = float(L1[20:32])
    nperday = float(L2[52:63])
    T       = oneday/nperday
    a       = (T**2 * GMe / twopi**2)**onethird
    data191.append((year, day, nperday, T, a))

data192 = []
for L1, L2 in pairs192:
    year = float(L1[18:20])
    day  = float(L1[20:32])
    nperday = float(L2[52:63])
    T       = oneday/nperday
    a       = (T**2 * GMe / twopi**2)**onethird
    data192.append((year, day, nperday, T, a))

data191 = np.array(zip(*data191))
data192 = np.array(zip(*data192))

a, b = np.histogram(data191[1], bins=np.arange(0, 201, 0.2))
c, d = np.histogram(data192[1], bins=np.arange(0, 201, 0.2))

if True:
    plt.figure()
    # plt.subplot(2, 1, 1)
    plt.plot(d[1:], c, '-r')
    plt.plot(b[1:], a, '-b')
    plt.yscale('log')
    plt.ylim(0.1, None)
    plt.xlim(175, 198)
    plt.title('epochs (0.2 day bins)')
    plt.suptitle('TLEs issued on days 191 and 192 of 2018', fontsize=18)
    plt.show()

if True:
    # plt.subplot(2, 1, 2)
    plt.plot(data192[1], data192[4]*0.001, '.r', markersize=2)
    plt.plot(data191[1], data191[4]*0.001, '.b', markersize=3)
    plt.xlim(175, 198)
    plt.ylim(None, 120000)
    plt.title('semimajor axis (km) versus epoch')
    plt.suptitle('TLEs issued on days 191 and 192 of 2018', fontsize=18)
    plt.show()