import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos, EarthSatellite

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]  # parentheses for Py3
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts   = load.timescale()

data    = load('de421.bsp')
earth   = data['earth']
ts      = load.timescale()

fname = "SMM tles 1989.txt"

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

L1s, L2s = lines[0::2], lines[1::2]
pairs    = zip(L1s, L2s)

alts, daynums = [], []
for i, (L1, L2) in enumerate(pairs):
    if not i%20:
        print i,
    daynum = float(L1[20:32])
    revs   = float(L2[52:64])
    T_mins = 24.*60/revs
    dayfraction = daynum % 1.0
    dayinteger  = int(daynum)
    minutes_epoch = 24. * 60 * dayfraction

    minutes = np.linspace(0, T_mins, 101)[:-1] + minutes_epoch
    times   = ts.utc(1989, 1, dayinteger, 0, minutes)

    sat     = EarthSatellite(L1, L2)
    pos     = sat.at(times).position.km
    r       = np.sqrt((pos**2).sum(axis=0))
    alt     = r - 6378.13
    alts.append(alt)
    daynums.append(daynum)

mins    = np.array([alt.min()  for alt in alts])
maxs    = np.array([alt.max()  for alt in alts])
means   = np.array([alt.mean() for alt in alts])
daynums = np.array(daynums)

if True:
    plt.figure()

    fs = 16

    plt.subplot(2, 1, 1)

    plt.plot(daynums, mins)
    plt.plot(daynums, maxs)
    plt.xlabel('day number in 1989', fontsize=fs)
    plt.ylabel('altitude (km)', fontsize=fs)

    plt.subplot(2, 1, 2)

    plt.plot(daynums, mins)
    plt.plot(daynums, maxs)
    plt.plot(daynums, means, '-k', linewidth=1.5)

    plt.xlim(0, 150)
    plt.ylim(425, 455)
    plt.xlabel('day number in 1989', fontsize=fs)
    plt.ylabel('altitude (km)', fontsize=fs)

    plt.suptitle('Solar Maximum Mission 1980-014A, 11703', fontsize=fs)

    plt.show()