! Random walk in 1 dimension
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program walk_1d
   option nolet
   library "sgfunc*","sglib*" 
   randomize
   dim x2ave(1000)    ! keep x^2 averages here
   call initialize(x2ave,n_walks,n_steps)
   call calculate(x2ave,n_walks,n_steps)
   call display(x2ave)
end
! initialize variables
! n_walks = number of walkers   n_steps = number of steps taken by each walker
sub initialize(x2ave(),n_walks,n_steps)
   n_steps = 100
   n_walks = 200
   mat redim x2ave(n_steps)
   ymax = sqr(n_steps)  ! set up window for plotting
   clear
   set window 0,n_steps,-2*ymax,2*ymax
   plot 0,0;n_steps,0
end sub
! do the calculation here
! x2ave(n) contains the average of x^2 at step n
! n_walks = total number of walkers   
! n_steps = number of steps taken by each walker
sub calculate(x2ave(),n_walks,n_steps)
   plot_flag = 1
   for i = 1 to n_walks
      x = 0                 ! current location of the walker
      for j = 1 to n_steps
         if rnd < 0.5 then
            x = x + 1
         else               ! just use an else statement 
            x = x - 1       ! DO NOT generate a new value using rnd
         end if
         if plot_flag = 1 then plot j,x      ! plot each walk until a keystroke is hit
         if key input then        ! after this just accumulate <x^2> for later display
            plot_flag = 0
            get key z
         end if
         x2ave(j) = x2ave(j) + x^2
      next j
   next i
   for i = 1 to n_steps               ! normalize x2ave when finished
      x2ave(i) = x2ave(i) / n_walks
   next i
end sub
! display results now
sub display(x2ave())
   dim t(0)  ! dummy array for plotting
   n = size(x2ave)
   mat redim t(n)
   for i = 1 to n
      t(i) = i
   next i
   call datagraph(t,x2ave,1,0,"black")
   call addlsgraph(t,x2ave,1,"black")   ! compute and plot least squares fit
   call fitline(t,x2ave,m,b)
   set cursor 4,12
   print "slope = ";m
end sub