import numpy as np
    import matplotlib.pyplot as plt

    halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
    to_degs, to_rads = 180/pi, pi/180

    e = 0.055
    a = 1.0
    mu = 1.0
    T = twopi # period for a=1 and mu=1

    # from https://space.stackexchange.com/a/43237/12102

    theta = np.linspace(-pi, pi, 501)
    E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
    t = a * np.sqrt(a/mu) * (E - e * np.sin(E))

    libration = theta - twopi * t/T

    if True:
        plt.figure()
        fs = 14
        plt.subplot(2, 2, 1)
        plt.plot(t/T, theta*to_degs)
        plt.xlabel('time (fraction of period)', fontsize=fs)
        plt.ylabel('theta (deg)', fontsize=fs)
        plt.xlim(-0.5, 0.5)
        plt.ylim(-180, 180)
        plt.subplot(2, 2, 2)
        plt.plot(theta*to_degs, t/T)
        plt.xlabel('theta (deg)', fontsize=fs)
        plt.ylabel('time (fraction of period)', fontsize=fs)
        plt.xlim(-180, 180)
        plt.ylim(-0.5, 0.5)
        plt.subplot(2, 1, 2)
        plt.plot(t/T, libration*to_degs)
        plt.xlabel('time (fraction of period)', fontsize=fs)
        plt.ylabel('apparent libration (deg)', fontsize=fs)
        plt.xlim(-0.5, 0.5)
        plt.text(0.15, 0, 'max: ' + str(round((libration*to_degs).max(), 2))  + u'\xb0', fontsize=fs)
        # https://stackoverflow.com/questions/3215168/how-to-get-character-in-a-string-in-python
        plt.show()