def deriv(X, t):
    th, thd = X
    thdd = const * np.sin(2*th)
    return np.hstack((thd, thdd))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

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

GMe   = 3.986E+14
Rc    = 1000. * (6378.137 + 400.)
const = -(3.*GMe) / (2.*Rc**3)
T     = twopi * np.sqrt(Rc**3/GMe)

theta_zeros = rads * np.arange(-80, 81, 20.0)

times   = np.linspace(0, 2*T, 1001)
minutes = times/60.
theta_orbit = twopi * times/T
theta_dot_orbit = twopi / T

answers = []
for theta in theta_zeros:
    X0 = np.hstack([theta, -theta_dot_orbit])
    answer, info = ODEint(deriv, X0, times, full_output=True)
    answers.append(answer)

print answer.shape

thetas     = [a.T[0] + theta_orbit     for a in answers]
theta_dots = [a.T[1] + theta_dot_orbit for a in answers]

answer = np.stack([a.T for a in answers])

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    for theta in thetas:
        plt.plot(minutes, degs*theta)
    plt.ylabel('theta (deg)', fontsize=16)
    plt.xlim(0, 190)

    plt.subplot(2, 1, 2)
    for theta_dot in theta_dots:
        plt.plot(minutes, degs*60.*theta_dot)
    plt.plot(minutes, np.zeros_like(minutes), '-k', linewidth=1.5)
    plt.xlabel('time (minutes)', fontsize=16)
    plt.ylabel('theta dot (deg/min)', fontsize=16)
    plt.xlim(0, 190)
    plt.show()