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

G=6.674108e-11
m_maa=5.972e+24
m_kuu=7.34767309e+22
r_maa=6378.0
r_kuu=1738.1
gm=G*m_maa*1.e-9
gk=G*m_kuu*1.e-9
rmk=384399.
vk=1.022
omegak=vk/rmk
#phi0=-52.435*np.pi/180.
phi0=-60*np.pi/180.

def func(y,t):
    dy=np.zeros_like(y)
    xk=rmk*np.cos(omegak*t+phi0)
    yk=rmk*np.sin(omegak*t+phi0)
    rr=np.sqrt(y[0]**2+y[1]**2)
    r=np.sqrt((y[0]-xk)**2+(y[1]-yk)**2)
    dy[0]=y[2]
    dy[1]=y[3]
    dy[2]=-y[0]*gm/rr**3-gk*(y[0]-xk)/r**3
    dy[3]=-y[1]*gm/rr**3-gk*(y[1]-yk)/r**3
    return dy

y0=[-6471.,0.,0.,-11.01]
time=np.linspace(0.,800000.,num=100001)
res=ode.odeint(func,y0,time,rtol=1.e-9,atol=1.e-9)

xk=rmk*np.cos(omegak*time+phi0)
yk=rmk*np.sin(omegak*time+phi0)
rr=np.sqrt(res[:,0]**2+res[:,1]**2)-r_maa
vv=np.sqrt(res[:,2]**2+res[:,3]**2)
rkk=np.sqrt((res[:,0]-xk)**2+(res[:,1]-yk)**2)

#fig,[p1,p2]=plt.subplots(1,2)
#p1.plot(hh,vv*1000)
#p2.plot(hh,aa)
#p1.grid()
#p2.grid()
#plt.show()
#******** visual 2 ************************************
fig1=plt.figure(1,figsize=(10,5))
p1=plt.subplot(121)
p2=plt.subplot(122)
print(min(rkk)-r_kuu)
p1.plot(time,rkk,color="red",label="dist")
p2.plot(time,vv,color="black",label="vel")

p1.grid(color="black")
p2.grid(color="black")

p1.set_xlabel("t")
p1.set_ylabel("r")
p2.set_xlabel("t")
p2.set_ylabel("v")
p1.legend(loc='upper right')
p2.legend(loc='upper right')
plt.show()