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

    two = """1 48413U 21038BN  21132.52488535  .00253791  39528-4  39172-3 0  9996
    2 48413  53.0519   5.2321 0014496 204.7636 155.2701 16.03554980  1249"""

    L1, L2 = two.splitlines()

    earth_radius_km = 6378.137

    load  = Loader('~/Documents/fishing/SkyData')  # single instance for big files

    sat = EarthSatellite(L1, L2, '48413')
    ts = load.timescale()
    epoch = sat.epoch.tt
    t = ts.tt_jd(np.linspace(epoch, epoch + 0.08, 1000))

    g = sat.at(t)
    valid = [m is None for m in g.message]

    timez = t.utc_datetime()
    altitude = np.where(valid, g.distance().km - earth_radius_km, np.nan)
    i_apo, i_peri = [np.argmax(a) for a in (altitude, -altitude)]
    apo, peri = [altitude[i] for i in (i_apo, i_peri)]
    t_apo, t_peri = [timez[i] for i in (i_apo, i_peri)]
    apo_label, peri_label = [str(round(alt, 1)) for alt in (apo, peri)]
    a_approx = 0.5 * (apo + peri)
    print('a_approx: ', a_approx)

    fig, ax = plt.subplots(1, 1)

    ax.set_title(sat.name + ' altitude @ epoch' + sat.epoch.ut1_strftime())

    ax.plot(timez, altitude)

    ax.text(timez[20], altitude[0], 'epoch', verticalalignment='center', fontsize=12)
    ax.plot(timez[:1], altitude[:1], 'or')

    ax.text(t_apo, apo, apo_label, verticalalignment='bottom', fontsize=12)
    ax.plot(t_apo, apo, 'ok')

    ax.text(t_peri, peri, peri_label, verticalalignment='top', fontsize=12)
    ax.plot(t_peri, peri, 'ok')

    ax.plot(timez, np.full_like(altitude, a_approx), '-k')
    a_approx_text = str(round(a_approx, 1))
    ax.text(timez[0], a_approx, a_approx_text,
            verticalalignment='bottom', fontsize=12)

    plt.show()