ccccccccccccccccccccccccc     Program 7.1     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
                      PROGRAM COMET
C
C Program for the time-dependent position of Halley's comet
C with the velocity version of the Verlet algorithm.
C
      PARAMETER (N = 20001,NP=10000,NI =200)
      DIMENSION X(N),Y(N),VX(N),VY(N),GX(N),GY(N),T(N),R(N)
      REAL KAPPA
C
C Initialization of the problem
C
      H     = 1./NP
      KAPPA = 39.478428
      T(1)  = 0.0
      X(1)  = 1.966843
      Y(1)  = 0.
      R(1)  = X(1)
      GX(1) = -KAPPA*X(1)/R(1)**3
      GY(1) = 0.
      VX(1) = 0.
      VY(1) = 0.815795
C
C Verlet algorithm for the position and velocity
C
      DO 100  I = 1, N-1
        T(I+1)  = I*H
        X(I+1)  = X(I)+H*VX(I)+H*H*GX(I)/2.
        Y(I+1)  = Y(I)+H*VY(I)+H*H*GY(I)/2.
        R(I+1)  = SQRT(X(I+1)*X(I+1)+Y(I+1)*Y(I+1))
        GX(I+1) = -KAPPA*X(I+1)/R(I+1)**3
        GY(I+1) = -KAPPA*Y(I+1)/R(I+1)**3
        VX(I+1) = VX(I)+H*(GX(I+1)+GX(I))/2.
        VY(I+1) = VY(I)+H*(GY(I+1)+GY(I))/2.
  100 CONTINUE
      WRITE (6,999) (T(I),R(I), I=1,N, NI)
  999 FORMAT (2F10.6)
      STOP
      END