! Kinetically weighted 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 ksaw_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),dr(4,2),available(4)
   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
   data 1,0,-1,0,0,1,0,-1
   mat read dr
   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
         hits = 0
         for k = 1 to 4
            available(k)=1
         next k
         for k = 0 to j-1
            for l = 1 to 4
               if x(k)=x(j-1)+dr(l,1) and y(k)=y(j-1)+dr(l,2) then
                  hits=hits+1
                  available(l)=0
                  exit for
               end if
            next l
            if hits=4 then exit for
         next k
         if hits=4 then
            if plot_flag = 1 then
               plot
               set color "red"
               plot area: x(j-1)-0.4,y(j-1)-0.4;x(j-1)-0.4,y(j-1)+0.4;x(j-1)+0.4,y(j-1)+0.4;x(j-1)+0.4,y(j-1)-0.4;x(j-1)-0.4,y(j-1)-0.4
               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
            r = int(rnd*(4-hits))+1
            if r=5-hits then r=r-1
            r0 = 0
            for k = 1 to 4
               r0=r0+available(k)
               if r=r0 then
                  k0=k
                  exit for
               end if
            next k
            x(j) = x(j-1) + dr(k0,1)
            y(j) = y(j-1) + dr(k0,2)
            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("KSAW 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