!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 7.1   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                       !
! Please Note:                                                          !
!                                                                       !
! (1) This computer program is written by Tao Pang in conjunction with  !
!     his book, "An Introduction to Computational Physics," published   !
!     by Cambridge University Press in 1997.                            !
!                                                                       !
! (2) No warranties, express or implied, are made for this program.     !
!                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM COMET
!
! Program for the time-dependent position of Halley's comet
! with the velocity version of the Verlet algorithm.
! Copyright (c) Tao Pang 1997.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=20001,NP=10000,NI =200
  INTEGER :: I
  REAL :: H,KAPPA
  REAL, DIMENSION (N) :: X,Y,VX,VY,GX,GY,T,R
!
! Initialization of the problem
!
  H     = 1.0/NP
  KAPPA = 39.478428
  T(1)  = 0.0
  X(1)  = 1.966843
  Y(1)  = 0.0
  R(1)  = X(1)
  GX(1) = -KAPPA*X(1)/R(1)**3
  GY(1) = 0.0
  VX(1) = 0.0
  VY(1) = 0.815795
!
! Verlet algorithm for the position and velocity
!
  DO I = 1, N-1
    T(I+1)  = I*H
    X(I+1)  = X(I)+H*VX(I)+H*H*GX(I)/2.0
    Y(I+1)  = Y(I)+H*VY(I)+H*H*GY(I)/2.0
    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.0
    VY(I+1) = VY(I)+H*(GY(I+1)+GY(I))/2.0
  END DO
  WRITE (6,"(2F16.8)") (T(I),R(I), I=1,N,NI)
END PROGRAM COMET