! rw1d.f90
! A simple random walk program  in 1D.
program rw1d
implicit none
integer ::    i,     N  ! number of steps
integer, dimension(8) :: seed
real :: Psx,Pdx
integer :: icount1, icount2, icount_rate, ix, irun, istep, nruns
real, dimension(:), allocatable :: rnd          ! array of random numbers
integer, dimension(:), allocatable :: x_N, x2_N ! sum of deviations and 
                                                ! squares over the runs
integer, dimension(:), allocatable :: P_N     ! final positions, sum over runs

print *, "Enter number of steps >"
read *, N
print *, "Enter number of runs >"
read *, nruns
print *, "Psx >"
read *, Psx

Pdx = 1 - Psx

open(UNIT=8, FILE='seed',ACTION='READ')
	do i=1,8
	read(8,*) seed(i)
	end do
  call random_seed(put=seed)

allocate(rnd(N))
allocate(x_N(N))
allocate(x2_N(N))
allocate(P_N(-N:N))
x_N  = 0
x2_N = 0
P_N  = 0


do irun = 1, nruns
  ix = 0 ! initial position of each run
  call random_number(rnd) ! get a sequence of random numbers
  do istep = 1, N
    if (rnd(istep) < Psx) then ! random move
      ix = ix - 1 ! left
    else
      ix = ix + 1 ! right
    end if
    x_N (istep) = x_N (istep) + ix
    x2_N(istep) = x2_N(istep) + ix**2 
  end do
  P_N(ix) = P_N(ix) + 0.5 ! accumulate (only for istep = N)
end do


print*,"# N=",N,"  nruns=",nruns

print*,' TEORICO '
print*,"# media TH = ", real(N*(Pdx-Psx)) 
print*,"# varianza TH = ", real((N*(Psx-Pdx))**2) + 4.*real(Psx)*real(Pdx)*real(N)

print*,' NUMERICO '
print*,"# media NUM = ",real(x_N(N))/nruns 
print*,"# varianza NUM = ",real(x2_N(N))/nruns - (real(x_N(N))/nruns)**2 

print*,"accuracy = ",abs( (   real(x2_N(N))/nruns - (real(x_N(N))/nruns)**2) / 	( 4.*real(Psx)*real(Pdx)*real(N) ) - 1.)



end program rw1d