import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as ode
import matplotlib.animation as anim

g_const=6.67408e-11
m_maa=5.9722e+24
gm=g_const*m_maa*1.e-9
r_maa=6371.
v1=np.sqrt(gm/r_maa)*1.
v1=9.
print (v1)
y0=[6578.,0.,0.,v1]


###--- POSSIBLE OPTIONS FOR INTAGRATION ---------
###--- Just uncomment the lines for the options you need 
###--- and comment out the ones you don't.

#-1-------- first option for integration ------   
#def func(y,t):
#   dy=np.zeros_like(y)
#   r3=np.sqrt(y[0]**2+y[1]**2)**3
#   dy[0]=y[2]
#   dy[1]=y[3]
#   dy[2]=-gm*y[0]/r3
#   dy[3]=-gm*y[1]/r3
#   return dy
#time=np.linspace(0,600000,num=1000)
#res=ode.odeint(func,y0,time)
#-1-------- end of the first option for integration ------   


##-2---------- second option for integration -----
def func(t,y):
   dy=np.zeros_like(y)
   r3=np.sqrt(y[0]**2+y[1]**2)**3
   dy[0]=y[2]
   dy[1]=y[3]
   dy[2]=-gm*y[0]/r3
   dy[3]=-gm*y[1]/r3
   return dy

t0=0.
dt=1.
nt=200000

r = ode.ode(func)
r.set_initial_value(y0, t0)
res = []
t = []
it=0
while r.successful() and np.sqrt(r.y[0]**2+r.y[1]**2) > r_maa and it<nt:
 it=it+1
 r.integrate(r.t+dt)
 res.append(r.y)
 t.append(r.t)
res = np.array(res)
t = np.array(t)
##-2---------- end of the second option for integration -----


###---END OF THE  POSSIBLE OPTIONS FOR INTAGRATION ---------


#*********** VISUALIZATION(ANIMATION) ************
fig,p=plt.subplots()
#p.set_xlim(-500000, 100000)
#p.set_ylim(-300000, 300000)
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)
p.plot(res[:,0],res[:,1])
ani = anim.FuncAnimation(fig, animf, np.arange(1, len(res[:,0])),interval=1,repeat=False,blit=True)
p.legend(["satelliit","maa"])
plt.show()