real y(4),work(27),mu_atm,m_keha,m_maa,a(4)
integer iwork(5)
common gm,r_maa,ro_atm_0,p1,p2,p3
external func
open(20,file="fall.dat")
pi=3.1415926536
c_keha=0.47
r_keha=2.8
m_keha=2800.
p_atm=101000.
mu_atm=29.
r_gaas=8314.
t_atm=300.
g_rask=9.814

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

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)
p1=mu_atm*g_rask*1.e3/(r_gaas*t_atm)
p2=s_keha*c_keha*500./m_keha
p3=g_rask*v_keha*1.e-3/m_keha

neqn=4
dt=1.
nt=200000
abserr=1.e-8
relerr=1.e-8
v0=7.7884079
alpha=0.5
      
y(1)=6578.0
y(2)=0.
y(3)=-v0*sin(pi*alpha/180.)
y(4)=v0*cos(pi*alpha/180.)
t=0
iflag=1
    
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.5.)then
r_keha=20.
s_keha=pi*r_keha**2
p2=c_keha*s_keha*500./m_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)
real y(4),dy(4)
common gm,r_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)
dy(1)=y(3)
dy(2)=y(4)
dy(3)=-y(1)*gm/rr**3-p2*ro*vv*y(3)+p3*ro*y(1)/rr
dy(4)=-y(2)*gm/rr**3-p2*ro*vv*y(4)+p3*ro*y(2)/rr
return
end