implicit none
real y(4),work(27),mu,m_keha,dy(4)
real g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0,vtuul,t,tout,dt
real pi,r_keha,c,s_keha,v_keha,p_atm,relerr,abserr
real vabs,amod,atau,anorm,r,r_gaas,ekin,epot,etot
real g_rask,t_atm,ro
integer it,nt,iflag, iwork(5),neqn
common p1,p2,p3,ro_atm_0,g_const,r_maa,m_maa
external func
open(10,file="dim2.dat")
open(20,file="kiirendus.dat")
pi=3.1415926536
r_keha=0.5
c=0.47
m_keha=5.
m_maa=5.97219e+24
g_const=6.6743015e-11
r_maa=6371000.
p_atm=100400.
t_atm=280.
mu=29.
r_gaas=8314.
s_keha=pi*r_keha**2
v_keha=4.*pi*r_keha**3/3
ro_atm_0=p_atm*mu/(r_gaas*t_atm)

call rask(g_const,r_maa,m_maa,0.,g_rask)
p1=g_rask*mu/(r_gaas*t_atm)
p2=0.5*c*s_keha/m_keha
p3=v_keha/m_keha

neqn=4
nt=10000
dt=0.001
t=0.
relerr=1.e-8
abserr=1.e-8

y(1)=0.
y(2)=0.
y(3)=10.
y(4)=10.
iflag=1

do it=1,nt
tout=t+dt

call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
if (iflag /= 2)then
iflag=2
endif
if (y(2) <0.) then
stop
endif

call func(t,y,dy)
vabs=sqrt(y(3)*y(3)+y(4)*y(4))
amod=sqrt(dy(3)*dy(3)+dy(4)*dy(4))
atau=(y(3)*dy(3)+y(4)*dy(4))/vabs
anorm=sqrt(amod**2-atau**2)
r=vabs**2/anorm

ekin=0.5*m_keha*vabs**2
call rask(g_const,r_maa,m_maa,y(2),g_rask)
epot=m_keha*g_rask*y(2)
etot=ekin+epot

write (10,'(15f15.6)')tout,y,dy(3:4),ekin,epot,etot
write (20,'(15f15.6)') amod,atau,anorm,r
enddo

stop
end

subroutine func(t,y,dy)
implicit none
real y(4),dy(4),g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0
real ro,vabs,t,g_rask
common p1,p2,p3,ro_atm_0,g_const,r_maa,m_maa
dy(1)=y(3)
dy(2)=y(4)
call rask(g_const,r_maa,m_maa,y(2),g_rask)
ro=ro_atm_0*exp(-p1*y(2))
vabs=sqrt(y(3)**2+y(4)**2)
dy(3)=-p2*ro*vabs*y(3)
dy(4)=-g_rask-p2*ro*vabs*y(4)+p3*ro*g_rask
return
end

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