real work(27),y(4),m_maa,m_keha,a(4)
integer iwork(5)
common gm
external func
open(10,file="kepler.dat")

relerr=1.e-7
abserr=1.e-7
m_maa=5.9722e+24
m_keha=2000.
r_maa=6378.
G=6.67408e-11
gm=G*m_maa*1.e-9
t=0.
y(1)=6578.
y(2)=0.
y(3)=0.
y(4)=9.8337
nt=200000
dt=1.

t=0.
iflag=1
neqn=4
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)
rr=sqrt(y(1)**2+y(2)**2)
ekin=0.5*m_keha*vv*vv*1.e6
epot=-G*m_keha*m_maa*1.e-3/rr
etot=ekin+epot
call func(t,y,a)
write(10,'(100e20.10)') tout,y,a(3:4),ekin,epot,etot
enddo

stop
end

subroutine func(t,y,dy)
common gm
real y(4),dy(4)
dy(1)=y(3)
dy(2)=y(4)
dy(3)=-gm*y(1)/sqrt(y(1)**2+y(2)**2)**3
dy(4)=-gm*y(2)/sqrt(y(1)**2+y(2)**2)**3
return
end