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,ro_atm_0,p1,p2,p3
external func

open(20,file="maandumine.dat")
pi=3.1415926536
c_keha=0.47
r_keha=2.
m_keha=2800.
p_atm=101000.
mu_atm=29.
r_gaas=8314.
t_atm=300.

g_const=6.67408e-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.
ro_atm_0=p_atm*mu_atm/(r_gaas*t_atm)
call rask(g_const,r_maa,m_maa,0.,g_rask)
p1=mu_atm*g_rask*1.e3/(r_gaas*t_atm)
p2=s_keha*c_keha*500./m_keha
p3=v_keha*1.e-3/m_keha

!print *,g_rask

neqn=4
dt=0.1
nt=2000000
abserr=1.e-9
relerr=1.e-9
v0=7.7884079
alpha=2.
h0=200.

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
ntlv=300
nlv=0.
    
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.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

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

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

enddo

stop
end

subroutine func(t,y,dy)
implicit none
real y(4),dy(4),gm,r_maa,ro_atm_0,p1,p2,p3
real h,rr,vv,g_rask,m_maa,t,g_const,ro
common gm,g_const,r_maa,m_maa,ro_atm_0,p1,p2,p3
rr=sqrt(y(1)**2+y(2)**2)
vv=sqrt(y(3)**2+y(4)**2)
h=rr-r_maa
ro=ro_atm_0*exp(-p1*h/(1+h/r_maa))
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