real y(4),work(27),m_maa,m_kuu
integer iwork(5)
common gm,gk,rmk,omegak,phi0
external func
open(10,file="kepler.dat")
pi=3.1415926536
G=6.67408e-11
m_maa=5.972e+24
m_kuu=7.34767309e+22
r_maa=6378.0
r_kuu=1738.1
gm=G*m_maa*1.e-9
gk=G*m_kuu*1.e-9
rmk=384399.
vk=1.022
omegak=vk/rmk

neqn=4
dt=10.
nt=80000
abserr=1.e-9
relerr=1.e-9
    
y(1)=-6471.0
y(2)=0.
y(3)=0.
y(4)=-11.01
phi0=-60.*pi/180.
!phi0=-52.7145*pi/180.
!phi0=-52.435*pi/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

xk=rmk*cos(omegak*tout+phi0)
yk=rmk*sin(omegak*tout+phi0)
rr=sqrt(y(1)**2+y(2)**2)-r_maa
vv=sqrt(y(3)**2+y(4)**2)
rkk=sqrt((y(1)-xk)**2+(y(2)-yk)**2)
write(10, '(10f15.4)')tout,y,xk,yk,rr,rkk-r_kuu,vv
if (rkk-r_kuu.lt.0.) then
print *,"crash"
stop
endif
enddo
stop
end

subroutine func(t,y,dy)
real y(4),dy(4)
common gm,gk,rmk,omegak,phi0
xk=rmk*cos(omegak*t+phi0)
yk=rmk*sin(omegak*t+phi0)
rr=sqrt(y(1)**2+y(2)**2)
r=sqrt((y(1)-xk)**2+(y(2)-yk)**2)
dy(1)=y(3)
dy(2)=y(4)
dy(3)=-y(1)*gm/rr**3-gk*(y(1)-xk)/r**3
dy(4)=-y(2)*gm/rr**3-gk*(y(2)-yk)/r**3
return
end