! Self-avoiding walk in 2 dimensions
! Modified from program to accompany "Computational Physics" by N. Giordano and H. Nakanishi
! Copyright Prentice Hall 1997, 2006
! H. Nakanishi

program saw_2d
   option nolet
   library "sgfunc.trc"
   library "sglib.trc"
   randomize
   dim x2ave(1),n_walks(1)
   call initialize(x2ave,n_walks,n_steps,cutoff)
   ymax = n_steps^0.75  ! set up window for plotting
   open #1: screen 0,1,0,0.9
   set window -2*ymax,2*ymax,-2*ymax,2*ymax
   set background color "white"
   call calculate(x2ave,n_walks,n_steps,maxstep,cutoff)
   call display(x2ave,maxstep,n_walks(1))
   close #1
end

! initialize variables
! n_walks = number of walkers   n_steps = number of steps taken by each walker
sub initialize(x2ave(),n_walks(),n_steps,cutoff)
   input prompt "maximum number of steps per walk => ": n_steps
   input prompt "number of walks at step 1 => ": n_walks(1)
   input prompt "fraction of surviving walks for statistics cutoff  => ": cutoff
   print "Hit p to pause, anything else to stop plotting/watching for key input."
   mat redim x2ave(n_steps),n_walks(n_steps)
end sub

! do the calculation here
! x2ave(n) contains the average of r^2 at step n
! n_walks(n) = total number of walkers still alive at n-th step
! n_steps = max number of steps taken by each walker
sub calculate(x2ave(),n_walks(),n_steps,maxstep,cutoff)
   dim x(0),y(0)
   mat redim x(0 to n_steps),y(0 to n_steps)
   randomize
   plot_flag = 1
   cross_flag = 1
   for j=2 to n_steps
      n_walks(j)=0
      x2ave(j)=0
   next j
   x(0) = 0
   y(0) = 0
   for i = 1 to n_walks(1)
      if plot_graph = 1 then
         set color "green"
         plot area: -0.4,-0.4;-0.4,0.4;0.4,0.4;0.4,-0.4;-0.4,-0.4
         set color "black"
         plot 0,0;
      end if
      for j = 1 to n_steps
         r = rnd
         if r <= 0.25 then
            x(j) = x(j-1) + 1
            y(j) = y(j-1)
         else if r <= 0.5 then
            x(j) = x(j-1) - 1
            y(j) = y(j-1)
         else if r <= 0.75 then
            x(j) = x(j-1)
            y(j) = y(j-1) + 1
         else
            x(j) = x(j-1)
            y(j) = y(j-1) - 1
         end if
         found=0
         for k = 0 to j-2
            if x(j)=x(k) and y(j)=y(k) then 
               found=1
               exit for
            end if
         next k
         if found=1 then
            if plot_flag = 1 then
               set color "red"
               plot x(j),y(j)
               set color "black"
               if cross_flag = 1 and key input then
                  get key z
                  c$ = chr$(z)
                  if c$="s" or c$="p" then
                     get key z
                  else
                     cross_flag=0
                     plot_flag=0
                  end if
               end if
            end if
            exit for
         else
            if j>1 then n_walks(j)=n_walks(j)+1
            x2ave(j) = x2ave(j) + x(j)^2 + y(j)^2
            if plot_flag = 1 and j < n_steps then
               plot x(j),y(j);
            else if plot_flag = 1 then
               plot x(j),y(j)
            end if
         end if
         if key input then ! after this just accumulate <x^2> for later display
            get key z
            c$ = chr$(z)
            if c$="s" or c$="p" then
               get key z
            else
               cross_flag = 0
               plot_flag = 0
            end if 
         end if
      next j
      clear
   next i
   for i = 1 to n_steps               ! normalize x2ave when finished
      if n_walks(i) > cutoff*n_walks(1) then
         x2ave(i) = log(x2ave(i) / n_walks(i))
      else
         exit for
      end if
      maxstep=i
   next i
end sub

! display results now
sub display(x2ave(),n,n_starts)
   dim t(0)  ! dummy array for plotting
   call setcanvas("white")
   set color "black"
   mat redim t(n),x2ave(n)
   for i = 1 to n
      t(i) = log(i)
   next i
   call settitle("SAW in d=2")
   call sethlabel("ln N")
   call setvlabel("ln <R^2>")
   call datagraph(t,x2ave,3,0,"black")
   call addlsgraph(t,x2ave,1,"red")   ! compute and plot least squares fit
   call fitline(t,x2ave,m,b)
   set cursor 4,20
   print "slope = ";m;", starts =";n_starts;", max =";n
   get key z
end sub