import numpy as np
    import matplotlib.pyplot as plt

    a_list = [0.723, 1.0, 1.524, 5.204]
    T_list = [a**1.5 for a in a_list]

    times = 2 * np.pi * np.linspace(0, 10, 2001)

    phases = [times/T for T in T_list]

    positions = np.array([[a * f(phase) for f in (np.cos, np.sin, np.zeros_like)]
                          for (a, phase) in zip(a_list, phases)])

    venus, earth, mars, jupiter = positions


    T_moon = 0.5
    a_moon = 1.3 / 150  # 1.3 million km divided by 1 AU

    phase_moon = times / T_moon

    moon_p = np.array([a_moon * f(phase_moon) for
                       f in (np.cos, np.sin, np.zeros_like)])
    moon_r = np.array([a_moon * f(-phase_moon) for
                       f in (np.cos, np.sin, np.zeros_like)])

    incs = [f * np.pi / 180 for f in (10, 19.5, 20.5, 10, 10)]
    for inc, thing in zip(incs, (venus, moon_p, moon_r, mars, jupiter)):
        sinc, cinc = np.sin(inc), np.cos(inc)
        thing1 = thing.copy()
        thing1[0] = cinc * thing[0] - sinc * thing[2]
        thing1[2] = sinc * thing[0] + cinc * thing[2]
        thing[[0, 2]] = thing1[[0, 2]]

    moon_p +=  earth
    moon_r +=  earth

    if True:
        xe, ye, ze = earth
        to_degs = 180 / np.pi
        for i, body in enumerate((venus, moon_p, moon_r, mars, jupiter)):
            d = body - earth
            lon = np.arctan2(d[1], d[0])
            lat = np.arctan2(d[2], np.sqrt(d[0]**2 + d[1]**2))
            bk1 = (lon[1:] > 2) * (lon[:-1] < -2)
            bk2 = (lon[:-1] > 2) * (lon[1:] < -2)
            bk = bk1 + bk2
            lon, lat = lon[:-1], lat[:-1]
            lon[bk] = np.nan
            plt.plot(to_degs * lon, to_degs * lat, linewidth=0.75)
            plt.xlim(-180, 180)
            plt.ylim(-90, 90)
            plt.gca().set_aspect('equal')
        plt.show()