import numpy as np
    import matplotlib.pyplot as plt
    GMe = 3.986E+14 # m^3/s^2
    Re  = 6378137.  # meters

    halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
    degs, rads        = 180/pi, pi/180
    Nhalf             = 200
    theta             = np.linspace(0, twopi, 2*Nhalf+1)
    peri, apo         = 1.1, 6
    a, ep             = 0.5*(apo + peri), (apo-peri)/(apo+peri)
    rcirc             = 1.1
    rellipse          = a * (1-ep**2)/ (1-ep*np.cos(theta))
    x0, y0            = [1.0   * f(theta) for f in (np.cos, np.sin)]
    xc, yc            = [rcirc   * f(theta) for f in (np.cos, np.sin)]
    xe, ye            = [rellipse* f(theta) for f in (np.cos, np.sin)]

    if True:
        plt.plot(x0, y0, '-b', linewidth=2)
        plt.plot(xc, yc, '-k')
        plt.plot(xe, ye, '-k')
        plt.plot(xe[::Nhalf], ye[::Nhalf], 'ok', markersize=8)
        plt.xlim(-2, 7)
        plt.show()

    rrel = np.arange(60) + 1.0
    r    = Re * rrel

    vorb    = np.sqrt(GMe/r)
    vesc    = np.sqrt(2) * vorb
    aa      = 0.5*(r + Re*peri)
    vapo    = np.sqrt(GMe*(2./r - 1./aa))
    vexcrit = vapo + vesc

    if True:
        plt.figure()
        things = (vesc, vapo, vexcrit)
        names  = ('vesc', 'vapo', 'vexcrit')
        rplot = r * 0.001 # km
        ilab  = 25
        for thing, name in zip(things, names):
            plt.plot(rplot, thing)
            plt.text(rplot[ilab], thing[ilab]+100, name, fontsize=14)
        plt.ylim(None, 6000)
        plt.xlabel('r_apo (km)', fontsize=18)
        plt.ylabel('velocity (m/s)', fontsize=18)
        for v in (2000, 3000, 4000):
            plt.plot(rplot, v*np.ones_like(rplot), '-k', linewidth=0.5)

        plt.show()