implicit none
real y(4),work(27),mu_atm,m_keha,m_maa,a(4)
real t,tout,relerr,abserr,v0,alpha,c_keha,r_keha,p_atm
real r_gaas,s_keha,ro_atm_0,gm,t_atm,p1,p2,p3,pi,h,h0
real v_keha,g_rask,dt,g_const,aa,ro,r_maa,vv,r,nlv,ntlv
integer iwork(5),i,nt,neqn,iflag
common gm,g_const,r_maa,m_maa,p2,p3
external func

open(20,file="maandumine.dat")
pi=3.1415926536
c_keha=0.47
!                C*S/(2*m) 
r_keha=1.55554 !0.000638  time=400s
!r_keha=4.91904 !0.00638   time=630s
!r_keha=10.9993  !0.0319    time=1130s

m_keha=2800.
g_const=6.67430e-11
m_maa=5.9722e+24
gm=g_const*m_maa*1.e-9
r_maa=6371.

s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
p2=s_keha*c_keha*500./m_keha
p3=v_keha*1.e-3/m_keha

neqn=4
dt=0.1
nt=2000000
abserr=1.e-9
relerr=1.e-9
v0=7.620
!v0=8.0
alpha=3.
h0=91.5
      
y(1)=r_maa+h0
y(2)=0.
y(3)=-v0*sin(pi*alpha/180.)
y(4)=v0*cos(pi*alpha/180.)
t=0
tout=0
iflag=1

nlv=0.
ntlv=100

do i=1,nt
tout=t+dt
call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)

if(iflag.ne.2) then
iflag=2
endif

vv=sqrt(y(3)**2+y(4)**2)
r=sqrt(y(1)**2+y(2)**2)
h=r-r_maa
call func(t,y,a)
aa=sqrt(a(3)**2+a(4)**2)

if(h.le.0.)then
stop
endif

if (h.le.10..and.nlv<ntlv)then
nlv=nlv+1.
r_keha=2.+18.*nlv/ntlv
s_keha=pi*r_keha**2
p2=c_keha*s_keha*500./m_keha
!print *,nlv,r_keha
endif

write(20, '(10f15.5)') tout,y,h,vv*1000,aa*1000

enddo

stop
end

subroutine func(t,y,dy)
use COESA_module
implicit none
real*8 hc,ro
real y(4),dy(4),gm,r_maa,p2,p3
real h,rr,vv,g_rask,m_maa,t,g_const
common gm,g_const,r_maa,m_maa,p2,p3
rr=sqrt(y(1)**2+y(2)**2)
vv=sqrt(y(3)**2+y(4)**2)
h=rr-r_maa

!****** COESA subroutine for air density calculation ******
hc=h*1000.
ro=COESA_density(hc)
!**********************************************************

call rask(g_const,r_maa,m_maa,h,g_rask)
dy(1)=y(3)
dy(2)=y(4)
dy(3)=-y(1)*gm/rr**3-p2*ro*vv*y(3)+p3*g_rask*ro*y(1)/rr
dy(4)=-y(2)*gm/rr**3-p2*ro*vv*y(4)+p3*g_rask*ro*y(2)/rr
return
end

subroutine rask(g_const,r_maa,m_maa,h,g_rask)
implicit none
real g_const,r_maa,m_maa,h,g_rask
g_rask=g_const*m_maa*1.e-6/(r_maa+h)**2
return
end