!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: comet.f90                                           !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: July 7, 2006                                            !
!                                                                         !
!  (1) This F90 program is created for the book, "An Introduction to      !
!      Computational Physics, 2nd Edition," written by Tao Pang and       !
!      published by Cambridge University Press on January 19, 2006.       !
!                                                                         !
!  (2) No warranties, express or implied, are made for this program.      !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM COMET
!
! An example to study the time-dependent position and
! velocity of Halley's comet via the Verlet algorithm.
!
  IMPLICIT NONE
  INTEGER, PARAMETER :: N = 20000, M = 200
  INTEGER :: I
  REAL :: H = 2.0/(N-1), K = 39.478328, R2, R3
  REAL, DIMENSION (N) :: T, X, Y, R, VX, VY, GX, GY
!
! Initialization of the problem
!
  T(1)  = 0.0
  X(1)  = 1.966843
  Y(1)  = 0.0
  R(1)  = X(1)
  VX(1) = 0.0
  VY(1) = 0.815795
  GX(1) = -K/(R(1)*R(1))
  GY(1) = 0.0
!
! 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
    R2 = X(I+1)*X(I+1)+Y(I+1)*Y(I+1)
    R(I+1)  = SQRT(R2)
    R3 = R2*R(I+1)
    GX(I+1) = -K*X(I+1)/R3
    GY(I+1) = -K*Y(I+1)/R3
    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, M)
END PROGRAM COMET