刚做了一个关于哈雷彗星运动的问题，给出了初位置和初速度，还有受力，求以后的运动
下面是我写的源程序
program halleycomet
implicit none
integer i
DOUBLE PRECISION,parameter :: G=6.67D-11
real*8,parameter :: M=1.99D+30
integer*8,PARAMETER :: N=20000000
integer IN
DOUBLE PRECISION T(N),X(N),Y(N),R(N),VX(N),VY(N),V(N),AX(N),AY(N),A(N)
real*8 DT
in=10000
X(1)=5.28D+12
Y(1)=0.0
R(1)=5.28D+12
VX(1)=0.0
VY(1)=-9.12D+2
V(1)=9.12D+2
AX(1)=-G*M/R(1)**2
AY(1)=0
A(1)=-G*M/R(1)**2
T(1)=0.0
DT=76*365*24*3600.0/N
DO 100 I=1,N

T(I+1)=DT*I
VX(I+1)=VX(I)+AX(I)*DT
VY(I+1)=VY(I)+AY(I)*DT
X(I+1)=X(I)+VX(I)*DT
Y(I+1)=Y(I)+VY(I)*DT
R(I+1)=(X(I+1)**2+Y(I+1)**2)**(1.0/2.0)
A(I+1)=-G*M/R(I+1)**2

AX(I+1)=A(I+1)*X(I+1)/R(I+1)
AY(I+1)=A(I+1)*Y(I+1)/R(I+1)

100 CONTINUE
open(1,file="zhaominghui.txt")
WRITE(1,"(4F18.3)")(X(I),Y(I),VX(I),VY(I),I=2,N,IN)
END PROGRAM