      PROGRAM random_walk
*     simulation of a random walk in one dimension
      IMPLICIT NONE
      REAL*8 p
      INTEGER N, itrial, ntrial, x2cum(64), xcum(64)
      CALL initial(p,N,xcum,x2cum,ntrial)
      DO itrial = 1, ntrial
         CALL walk(p,N,xcum,x2cum)
      END DO
      CALL output(N,xcum,x2cum,ntrial)
      END
C
      SUBROUTINE initial(p,N,xcum,x2cum,ntrial)
      IMPLICIT NONE
      REAL*8 p
      REAL dmy, rnd
      INTEGER istep, N, xcum(64), x2cum(64), ntrial
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'maximum number of steps N = '
      READ(*,*) N
      p = 0.5D0
      ntrial = 1000                         ! number of trials
      DO istep = 1, N
         xcum(istep) = 0                    ! not necessary in True BASIC
         x2cum(istep) = 0
      END DO
      END
C
      SUBROUTINE walk(p,N,xcum,x2cum)
      IMPLICIT NONE
      REAL*8 p
      REAL rnd
      INTEGER istep, x, N, xcum(64), x2cum(64)
      x = 0                                 ! initial position of walker for each trial
      DO istep = 1, N
         IF(rnd(0) .LE. p) THEN
            x = x + 1
         ELSE
            x = x - 1
         END IF
*        collect data after every step
         CALL mydata(x,xcum,x2cum,istep)
      END DO
      END
C
      SUBROUTINE mydata(x,xcum,x2cum,istep)
      IMPLICIT NONE
      INTEGER x, xcum(64), x2cum(64), istep
      xcum(istep) = xcum(istep) + x
      x2cum(istep) = x2cum(istep) + x*x
      END
C
      SUBROUTINE output(N,xcum,x2cum,ntrial)
      IMPLICIT NONE
      REAL*8 variance, x2bar, xbar
      INTEGER istep, N, xcum(64), x2cum(64), ntrial
      WRITE(*,'(3A13,2X,A)') '# steps', '<x>', '<x^2>', '<x^2> - <x>^2'
      WRITE(*,*)
      DO istep = 1, N
         xbar = DBLE(xcum(istep))/ntrial
         x2bar = DBLE(x2cum(istep))/ntrial
         variance = x2bar - xbar*xbar
         WRITE(*,'(I13,4F13.6)') istep, xbar, x2bar, variance
      END DO
      END
