import numpy as np
import matplotlib.pyplot as plt

# Constants
m1 = 5.972E24   # kg mass Earth
m2 = 100.       # kg mass of astronaut
mu = m1*m2/(m1+m2)
G = 6.674E-11   # m^3 kg^-1 s^-2  gravitational constant
l = 5.2E12     # kg m^2 s^-1 angular momentum of astronaut (assuming ISS conditions)
rE = 6.371E6    # m radius of Earth
d = 4.00E5      # m typical ISS orbit height

# Plot Earth radius
r = rE*np.ones(100)
phi = 2 * np.pi * np.linspace(0, 1, 100)
plt.polar(phi, r, 'g--')

# Plot ISS circular orbit
r = (rE+d)*np.ones(100)
plt.polar(phi, r, 'r')

# Plot astronaut orbit
# Calculate eccentricity
dp = 2000.      #  kg m s^-1
s0 = G*m1*m2*mu/l**2
eps = dp/l/s0
print 'Eccentricity: ', eps

# Calculate r(\phi) assuming phase \phi_0=0
r = (rE+d)/(1-eps*np.cos(phi-np.pi/2))
plt.polar(phi, r, 'b')

# Decorations
plt.xlabel('$\phi$')
plt.ylabel('$r$ [m]', labelpad=25)
plt.ylim([0.97*rE, 1.07*rE])
plt.savefig('orbit.png')