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



#****************************
#****** parameetrid  **************
#****************************************

xm0=-6471.
ym0=0.
vmx0=0.
vmy0=-11.01

phi=-60.
#phi=-52.435
#phi=-54

animkiirus=100
#*****************************************
#***********************************
#****************************
asa=100000
aeg=800000.

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

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=[xm0,ym0,vmx0,vmy0]
time=np.linspace(0.,aeg,num=asa)
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)

c_text=""
for i in np.arange(0,len(time)):
    if rkk[i] <= r_kuu:
      tanim=i
      c_text="Kukkus Kuule!!!"
      break
    tanim=i
if c_text == "":
    for i in np.arange(0,len(time)):
        if rr[i] <= 0.:
          tanim=i
          c_text="Kukkus Maale!!!"
          break
        tanim=i
#print rkk[tanim]
    

#*******************************************

fig=plt.figure(2,figsize=(15,10))
p1=plt.subplot(131)

p1.plot(time[0:tanim]/3600,rkk[0:tanim]-r_kuu,color="red",label="Kuu ja mooduli vaheline kaugus")
p1.grid(color="black")
p1.legend(loc='upper right')
p1.set_xlabel("Aeg [tund]")
p1.set_ylabel("kaugus [km]")

p2=plt.subplot(132)
p2.plot(time[0:tanim]/3600.,vv[0:tanim],color="black",label="mooduli kiirus")
p2.set_xlabel("Aeg [tund]")
p2.set_ylabel("Kiirus [km/s]")
p2.grid(color="black")
p2.legend(loc='upper right')

p4=plt.subplot(133)
p4.plot(res[0:tanim,0],res[0:tanim,1],color="red",label="moodul")
p4.plot(xk[0:tanim],yk[0:tanim],color="blue",label="Kuu")
plt.xlim(-600000, 600000)
plt.ylim(-600000, 600000)
p4.grid()

if c_text == "":
    kaugus_text="Minimaalne kaugus="+str('%15.5f' % min(rkk[0:tanim]-r_kuu) )+' km'
else:
    kaugus_text="Minimaalne kaugus= 0.0 km"

p4.text(0.02,0.93,c_text,transform=p4.transAxes,color='red')
p4.text(0.02,0.2,kaugus_text,transform=p4.transAxes,color='black')
p4.set_aspect("equal")
p4.grid(color="black")
p4.legend(loc='upper right')

#******** visual ************************************
fig1=plt.figure(1,figsize=(10,10))
p3=plt.subplot(111)

plt.xlim(-600000, 600000)
plt.ylim(-600000, 600000)
p3.grid()

p3.set_aspect("equal")
p3.set_xlabel("x-koordinaat [km]")
p3.set_ylabel("y-koordinaat [km]")

maa=plt.Circle((0,0),6378,color="green")
ohu=plt.Circle((0,0),6478,color="blue",fill=False)
p3.add_artist(maa)
p3.add_artist(ohu)

sat, = p3.plot([], [], '.',color="red",label="Moodul")
kuu, = p3.plot([], [], 'o',color="blue",label="Kuu")
vel_txt = p3.text(0.02,0.95,'',transform=p3.transAxes)
dist_txt = p3.text(0.02,0.91,'',transform=p3.transAxes)

def animf(i):
    kx = (xk[i])
    ky = (yk[i])
    sx = (res[i,0])
    sy = (res[i,1])
    sat.set_data(sx, sy)
    kuu.set_data(kx, ky)
    vel_txt.set_text('Kiirus = %15.5f km/s' % (vv[i]))
    dist_txt.set_text('Kuu-moodul kaugus = %20.10f km' % (rkk[i]-r_kuu))
    return sat,kuu,vel_txt,dist_txt

ani = anim.FuncAnimation(fig1, animf, np.arange(1, tanim,animkiirus),interval=1,repeat=False,blit=True)
p3.legend(loc='upper right')

plt.show()