      PROGRAM planet                        ! planetary motion
      IMPLICIT NONE
      REAL*8 GM, dt, pos(2), t, vel(2)
      INTEGER IR_, counter, nshow
      LOGICAL Ltmp1_, TBkeyinput
      CALL GWopen(IR_, 0)
      CALL initial(pos,vel,t,GM,dt,nshow)
      CALL output(pos,t)
      counter = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL Euler(pos,vel,t,GM,dt)
         counter = counter + 1
         IF(mod(counter,nshow) .EQ. 0) THEN
            CALL output(pos,t)              ! draw position of "earth"
         END IF
         Ltmp1_ = .NOT. TBkeyinput()
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(pos,vel,t,GM,dt,nshow)
      IMPLICIT NONE
      REAL*8 pi, pos(2), vel(2), t, GM, dt
      REAL r, xwin, ywin
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, nshow
      GM = 4.0D0*pi**2                      ! GM constant in astronomical units
      t = 0
      WRITE(*,'(A,$)') 'time step = '
      READ(*,*) dt
      WRITE(*,'(A,$)') 'number of time steps between plots = '
      READ(*,*) nshow
      WRITE(*,'(A,$)') 'initial x position = '
      READ(*,*) pos(1)
      pos(2) = 0                            ! initial y-position
      vel(1) = 0                            ! initial x-velocity
      WRITE(*,'(A,$)') 'initial y-velocity = '
      READ(*,*) vel(2)
*     assumed maximum value of semi-major axis
      r = 2*pos(1)
      CALL compute_aspect_ratio(r,xwin,ywin) ! in library
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      END
C
      SUBROUTINE Euler(pos,vel,t,GM,dt)
      IMPLICIT NONE
      REAL*8 accel(2), r2, r3, pos(2), vel(2), t, GM, dt
      INTEGER i
      r2 = pos(1)*pos(1) + pos(2)*pos(2)
      r3 = r2*sqrt(r2)
      DO i = 1, 2
         accel(i) = -GM*pos(i)/r3
         vel(i) = vel(i) + accel(i)*dt
         pos(i) = pos(i) + vel(i)*dt
      END DO
      t = t + dt
      END
C
      SUBROUTINE output(pos,t)
      IMPLICIT NONE
      REAL*8 pos(2), t
      REAL radius
      INTEGER IR_
      IF(t .EQ. 0) THEN
         CALL GWclear(IR_, -1)
         radius = 0.05                      ! radius of "sun" (not to scale)
         CALL GWsetpen(IR_, 13, -1, -1, -1)
         CALL GWsetmrk(IR_, 6, 2*radius, 13, -1, 4)
         CALL GWputmrk(IR_, 0.0, 0.0)
      END IF
      CALL GWsetpxl(IR_, REAL(pos(1)), REAL(pos(2)), 16)
      END
