Tiangong-2 orbit

class Satellite(object):
    def __init__(self, daynum):
        self.daynum = daynum

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from skyfield.api import Loader, EarthSatellite

Re = 6378.137

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

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

print earth.at(ts.now()).position.km   # for no reason at all

with open('Tiangong-2 2018 TLEs.txt', 'r') as infile:
    lines = infile.readlines()

L1s, L2s = [lines[i::2] for i in range(2)]

sats = []
for L1, L2 in zip(L1s, L2s):
    if L1[0]=='1' and L2[0]=='2':
        daynum = float(L1[20:32])
        sat = Satellite(daynum)
        sats.append(sat)
        sat.L1, sat.L2 = L1, L2
        sat.satobj = EarthSatellite(L1, L2) # add the Skyfield object

daynums     = np.array([sat.daynum for sat in sats])

minuterange = np.arange(0, 120, 0.5)
dayrange    = minuterange / (60. * 24.)

for sat in sats:
    sat.days  = sat.daynum + dayrange
    sat.tsutc = ts.utc(2018, 1, dayrange)
    sat.pos   = sat.satobj.at(sat.tsutc).position.km
    sat.r     = np.sqrt((sat.pos**2).sum(axis=0))
    sat.rmax  = sat.r.max()
    sat.rmin  = sat.r.min()

rmaxes = np.array([sat.rmax for sat in sats])
rmins  = np.array([sat.rmin for sat in sats])

if True:
    plt.figure()
    fs = 16
    plt.plot(daynums, rmaxes - Re)
    plt.plot(daynums, rmins  - Re)
    plt.plot(daynums, rmaxes - Re, '.k')
    plt.plot(daynums, rmins  - Re, '.k')
    plt.title('Tiangong-2 peri and apo (km) June 2018 so far', fontsize=fs)
    plt.xlabel('daynumber in 2018', fontsize=fs)
    plt.xlim(151, 180)
    plt.show()

if True:
    plt.figure()
    fs = 16
    plt.plot(daynums, rmaxes - Re)
    plt.plot(daynums, rmins  - Re)
    plt.plot(daynums, rmaxes - Re, '.k')
    plt.plot(daynums, rmins  - Re, '.k')
    plt.title('Tiangong-2 peri and apo (km) 2018 so far', fontsize=fs)
    plt.xlabel('daynumber in  2018', fontsize=fs)
    plt.show()