#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as ode
import matplotlib.animation as anim

g=6.67408e-11
m=5.9722e+24
gm=g*m*1.e-9
r0=6578.
v1=np.sqrt(gm/r0)*1.3
print (v1)

def func(y,t):
   yp=np.zeros_like(y)
   r3=np.sqrt(y[0]**2+y[1]**2)**3
   yp[0]=y[2]
   yp[1]=y[3]
   yp[2]=-gm*y[0]/r3
   yp[3]=-gm*y[1]/r3
   return yp
   
#y0=[6578.,0.,0.,v1]
y0=[6578.,0.,0.,v1]
time=np.linspace(0,60000,num=1000)
res=ode.odeint(func,y0,time)

#*********** visualization ************
fig=plt.figure()
p = plt.subplot()
p=plt.axes(xlim=(-50000, 10000), ylim=(-30000, 30000))
p.grid()
p.set_aspect("equal")
sat, = p.plot([], [], '.',color="red")

def animf(i):
    sx = [res[i,0]]
    sy = [res[i,1]]
    sat.set_data(sx, sy)
    return sat,

maa=plt.Circle((0,0),6378,color="green")
p.add_artist(maa)
ani = anim.FuncAnimation(fig, animf, np.arange(1, len(res[:,0])),interval=10,repeat=False,blit=True)
plt.show()