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

gamma=6.67408e-11
m_maa=5.9722e+24
gm=gamma*m_maa*1.e-9
r_maa=6378.
g_rask=9.814

m_keha=2800.
r_keha=2.
c_keha=0.47
s_keha=np.pi*r_keha**2

t_ohu=300.
p_ohu_0=101000.
mu_ohu=29.
r_gaas=8314.
ro_ohu_0=p_ohu_0*mu_ohu/(r_gaas*t_ohu)

p1=mu_ohu*g_rask*1000./(r_gaas*t_ohu)
p2=500.*c_keha*s_keha/m_keha

def func(y,t):
   yp=np.zeros_like(y)
   r=np.sqrt(y[0]**2+y[1]**2)
   h=r-r_maa
   ro_ohu=ro_ohu_0*np.exp(-p1*h)
   v=np.sqrt(y[2]**2+y[3]**2)
   r3=r**3
   yp[0]=y[2]
   yp[1]=y[3]
   yp[2]=-gm*y[0]/r3-p2*ro_ohu*v*y[2]
   yp[3]=-gm*y[1]/r3-p2*ro_ohu*v*y[3]
   return yp

v0=7.7884079
alpha=2.*np.pi/180.
vx=-v0*np.sin(alpha)
vy=v0*np.cos(alpha)

y0=[6578.,0.,vx,vy]
time=np.linspace(0,760,num=10000)
res=ode.odeint(func,y0,time)

r=np.sqrt(res[:,0]**2+res[:,1]**2)
hh=r-r_maa
vv=np.sqrt(res[:,2]**2+res[:,3]**2)
r3=r**3
ro=ro_ohu_0*np.exp(-p1*hh)
ax=-gm*res[:,0]/r3-p2*ro*vv*res[:,2]
ay=-gm*res[:,1]/r3-p2*ro*vv*res[:,3]
aa=np.sqrt(ax**2+ay**2)*1000
#np.savetxt ("out.dat",hh,vv)
print (max(aa),min(vv*1000*3.6))
#******* visual 1*************************************

#fig,[p1,p2]=plt.subplots(1,2)
#p1.plot(hh,vv*1000)
#p2.plot(hh,aa)
#p1.grid()
#p2.grid()
#plt.show()
#******** visual 2 ************************************
fig=plt.subplots()

p1=plt.subplot(121)
p2=plt.subplot(122)

p1.set_xlabel("aeg")
p1.set_ylabel("kiirus")
p1.plot(hh,vv,color='blue')
p1.tick_params(axis='y',labelcolor='blue')
p1.grid()

p2.set_ylabel("kiirendus")
p2.plot(hh,aa,color='red')
p2.tick_params(axis='y',labelcolor='red')
p2.grid()
plt.show()