real y(4),work(27),mu_atm,m_keha,a(4)
integer iwork(5)
common /aaa/g_rask,ro_atm_0,p1,p2,p3
common /bbb/ ar(2),at(2),aa(2)
external f
open(10,file="dim2.dat")
open(20,file="kiirendus.dat")
pi=3.1415926536
g_rask=9.814
r_keha=1.
s_keha=pi*r_keha**2
c=0.47
v_keha=4.*pi*r_keha**3/3
m_keha=15.
p_atm=100400.
t_atm=280.
mu_atm=29.
r_gaas=8314.
ro_atm_0=p_atm*mu_atm/(r_gaas*t_atm)
p1=g_rask*mu_atm/(r_gaas*t_atm)
p2=0.5*c*s_keha/m_keha
p3=g_rask*v_keha/m_keha
neqn=4
nt=4000
dt=0.001
t=0.
relerr=1.e-6
abserr=1.e-6

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

do it=1,nt
tout=t+dt

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

call f(t,y,a)
vabs=sqrt(y(3)*y(3)+y(4)*y(4))
aabs=sqrt(a(3)*a(3)+a(4)*a(4))
atau=(y(3)*a(3)+y(4)*a(4))/vabs
anorm=sqrt(aabs**2-atau**2)
radius=vabs**2/anorm
write (10,*) tout,y,ar,at,aa,a(3:4)
write (20,*) atau,anorm,aabs,radius
enddo

stop
end

subroutine f(t,y,yp)
real y(4),yp(4),m_keha
common /aaa/ g_rask,ro_atm_0,p1,p2,p3
common /bbb/ ar(2),at(2),aa(2)
yp(1)=y(3)
yp(2)=y(4)
ro=ro_atm_0*exp(-p1*y(2))
vabs=sqrt((y(3))**2+y(4)**2)
ar(1)=0.
ar(2)=-g_rask
at(1)=-p2*ro*vabs*y(3)
at(2)=-p2*ro*vabs*y(4)
aa(1)=0.
aa(2)=ro*p3
yp(3)=at(1)
yp(4)=ar(2)+at(2)+aa(2)
return
end