program saw
      parameter (maxedge=100,maxz=4,maxstep=100)
      implicit real*8 (a-h,o-z)
      dimension wk(maxstep),ce2(maxstep),ce4(maxstep),cr2(maxstep),
     c          cr4(maxstep),e2(maxstep),r2(maxstep)
      character lat(-maxedge:maxedge,-maxedge:maxedge)
      common /blk/ lx(0:maxstep),ly(0:maxstep),lat
      real etime,time(2)
      character*12 out
c
c
c
   10 print *,'SAW on square lattice: seed, #steps, #trials ?'
      read(5,*) isd,nstep,ntri
      isd0=isd
      nstep=min(nstep,maxstep)
      print *,'output file name ?'
      read(5,15) out
   15 format(a12)
      do 30 i=1,nstep
         wk(i)=0.d0
         ce2(i)=0.d0
         ce4(i)=0.d0
         cr2(i)=0.d0
         cr4(i)=0.d0
   30 continue
      do 60 i=-nstep,nstep
      do 60 j=-nstep,nstep
   60 lat(i,j)='0'
      lat(0,0)='1'
      lat(1,0)='1'
      lx(0)=0
      ly(0)=0
      lx(1)=1
      ly(1)=0
      mstep=1
      t0=etime(time)
      do 100 i=1,ntri
         call makesaw(nstep,mstep,e2,r2,isd)
         do 40 n=1,mstep
            wk(n)=wk(n)+1.d0
            ce2(n)=ce2(n)+e2(n)
            ce4(n)=ce4(n)+e2(n)**2
            cr2(n)=cr2(n)+r2(n)
            cr4(n)=cr4(n)+r2(n)**2
   40    continue
  100 continue
      t1=etime(time)
      ntrue=0
      do 50 i=1,nstep
         if(wk(i).eq.0) go to 70
         ntrue=ntrue+1
         v=1.d0/dfloat(ntrue)**2
         w=1.d0/wk(i)
         ce2(i)=ce2(i)*w
         sde=ce4(i)*w-ce2(i)**2
         if(sde.lt.0) sde=0.d0
         ce4(i)=dsqrt(sde)
         cr2(i)=cr2(i)*v*w
         sdr=cr4(i)*w*v**2-cr2(i)**2
         if(sdr.lt.0) sdr=0.d0
         cr4(i)=dsqrt(sdr)
   50 continue
   70 open(1,file=out)
      rewind(1)
      write(6,24)
   24 format(/1x,'steps',2x,'ms r_ee',9x,'sd',10x
     c         ,'ms r_g',10x,'sd',11x,'#wks/4')
      do 20 i=1,ntrue
         write(6,25) i,ce2(i),ce4(i),cr2(i),cr4(i),wk(i)
         write(1,25) i,ce2(i),ce4(i),cr2(i),cr4(i),wk(i)
   20 continue
   25 format(1x,i3,1x,1p,5(e12.5,2x))
      write(6,35) t1-t0,isd0,isd
   35 format(/1x,'CPU time (sec) = ',1p,e12.5,2(2x,i12))
      stop 
      end
c
c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c     subroutine "makesaw" to calculate SAW statistics on one uster
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
c
      subroutine makesaw(nstep,mstep,e2,r2,isd)
      parameter (maxedge=100,maxz=4,maxstep=100)
      implicit real*8 (a-h,o-z)
      real ranf
      dimension e2(1),r2(1)
      character lat(-maxedge:maxedge,-maxedge:maxedge)
      common /blk/ lx(0:maxstep),ly(0:maxstep),lat
      dimension ndr(2,maxz)
      data ndr /-1,0,1,0,0,-1,0,1/
c
c
c
      do 10 n=2,mstep
            lat(lx(n),ly(n))='0'
   10 continue
      e2(1)=1.d0
      r2(1)=1.d0
      ldir=1
c
c
c
      do 20 n=2,nstep
c
c
c
         nbr=0
         i1=lx(n-1)
         j1=ly(n-1)
         do 280 k=1,maxz
            if(k.eq.ldir) go to 280
            if(lat(i1+ndr(1,k),j1+ndr(2,k)).eq.'0') nbr=nbr+1
  280    continue
         if(nbr.eq.0) then
            mstep=n-1
            return
         endif
c
c NOW THE RIGHT STUFF HERE
c
         ir=int(ranf(isd)*(maxz-1))+1
         ipass=0
         do 250 idir=1,maxz
            if(idir.eq.ldir) go to 250
            ipass=ipass+1
            if(ipass.eq.ir) then
               i=lx(n-1)+ndr(1,idir)
               j=ly(n-1)+ndr(2,idir)
               if(lat(i,j).eq.'1') then
                  mstep=n-1
                  return
               endif
               lx(n)=i
               ly(n)=j
               lat(i,j)='1'
               e2(n)=dfloat(i)**2+dfloat(j)**2
               r2(n)=r2(n-1)+e2(n)
               do 260 k=1,n-1
  260          r2(n)=r2(n)+dfloat(i-lx(k))**2+dfloat(j-ly(k))**2
               ldir=2*mod(idir,2)+idir-1
               go to 20
            endif
  250    continue
c
c
c
   20 continue
      mstep=nstep
c
c
c
      return
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< RANF >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      real function ranf(isd)
      isd=isd*1566083941+1
      ranf=float(isd)*2.328306e-10+0.5
      return
      end