      parameter(md=4,mxs=2000,mxs1=2001,mlen=31,max=25)
      implicit real*8 (v)
      dimension v3(md),p0(mlen)
      integer coord(mxs1,md),mstep(mlen),dim,seed,dir(mxs),maxc(md)
      character*35 title1,title2
      character*36 titlea,titleb
      character ans,yes,null,more,amark,agrid,ascale
      character*3 mark
      data null/x'00'/,more/x'02'/
      mark='$1'//null
      yes='y'
c
c
c
   80 write(6,11)
   11 format(1x,'type of walk?'//
     x' -1 = RW, -2 = BRW, -3 = RLF, -4 = BRLF'/
     y' 1 = SAW, 2 = BSAW, 3 = NALF, 4 = BNALF, 5 = PALF, 6 = BPALF'/)
      read(5,*) iwalk
      if(iwalk.eq.1) then
          write(6,51)
   51     format(1x,'seed, dim, #steps, #trials ?')
          read(5,*) seed,dim,nstep,ntri
      elseif(iwalk.eq.2) then
          write(6,52)
   52     format(1x,'seed, dim, pg, #steps, #trials ?')
          read(5,*) seed,dim,pg,nstep,ntri
      elseif(iwalk.eq.3.or.iwalk.eq.5) then
          write(6,53)
   53     format(1x,'seed, dim, mu, #steps, #trials ?')
          read(5,*) seed,dim,dmu,nstep,ntri
      elseif(iwalk.eq.4.or.iwalk.eq.6) then
          write(6,54)
   54     format(1x,'seed, dim, pg, mu, #steps, #trials ?')
          read(5,*) seed,dim,pg,dmu,nstep,ntri
      elseif(iwalk.eq.-1) then
          write(6,41)
   41     format(1x,'seed, dim, #steps ?')
          read(5,*) seed,dim,nstep
      elseif(iwalk.eq.-2) then
          write(6,42)
   42     format(1x,'seed, dim, pg, #steps ?')
          read(5,*) seed,dim,pg,nstep
      elseif(iwalk.eq.-3) then
          write(6,43)
   43     format(1x,'seed, dim, mu, #steps ?')
          read(5,*) seed,dim,dmu,nstep
      elseif(iwalk.eq.-4) then
          write(6,44)
   44     format(1x,'seed, dim, pg, mu, #steps ?')
          read(5,*) seed,dim,pg,dmu,nstep
      else
          write(6,55)
   55     format(1x,'invalid type, try again')
          go to 80
      endif
      if(dim.gt.md.or.nstep.gt.mxs) go to 80
      if(iwalk.lt.0) ntri=1
c
c
c
      if(iabs(iwalk).ge.3) then
          lb=2
          a=float(lb)**dmu
          p1=(a-1.)/a
          p0(1)=p1
          mstep(1)=1
          do 72 i=2,max
          p1=p1/a
          p0(i)=p0(i-1)+p1
          mstep(i)=lb*mstep(i-1)
   72     continue
          mstep(max+1)=lb*mstep(max)
      endif
c
c
c
      if(iabs(iwalk).eq.3.or.iwalk.eq.-1) then
          np=2*dim
      else
          np=2*dim-2
          if(iwalk.eq.1.or.iwalk.eq.5) pg=float(np)/float(np+1)
      endif
c
c
c
      do 27 i=1,nstep+1
      do 27 j=1,dim
   27 coord(i,j)=0
c
c
c
      do 10 itri=1,ntri
c
c
c
      iseed=seed
      do 20 i=2,dim
      coord(1,i)=0
   20 coord(2,i)=0
      coord(1,1)=0
      if(iabs(iwalk).le.2) then
          coord(2,1)=1
      else
          p=ranf(seed)
          do 71 i=1,max
          if(p.le.p0(i)) go to 73
   71     continue
          i=max+1
   73     coord(2,1)=mstep(i)
      endif
c
      imod=1
      jd=1
      dir(1)=1
      v1=0.d0
      do 21 i=1,dim
      maxc(i)=iabs(coord(2,i))
   21 v3(i)=dfloat(coord(2,i))
c
c
c
      do 30 i=3,nstep+1
      p=ranf(seed)
      if(iwalk.eq.-1.or.iabs(iwalk).eq.3) then
              ip=int(float(np)*p)
              jd=ip/2+1
              imod=2*mod(ip,2)-1
      elseif(p.lt.pg) then
              ip=int(float(np)*p/pg)
              kd=ip/2+1
              if(kd.ge.jd) kd=kd+1
              jd=kd
              imod=2*mod(ip,2)-1
      endif
c
c
c
      do 31 j=1,dim
      if (j.eq.jd) then
          istep=1
          if(iabs(iwalk).ge.3) then
              p=ranf(seed)
              do 81 k=1,max
              if(p.le.p0(k)) go to 83
   81         continue
              k=max+1
   83         istep=mstep(k)
          endif
	  coord(i,j)=coord(i-1,j)+imod*istep
      else
          coord(i,j)=coord(i-1,j)
      endif
   31 continue
      dir(i-1)=jd
c     
c     
c
      if(iwalk.lt.0) go to 90
      if(iwalk.le.4) then
         do 32 j=1,i-1
         do 33 k=1,dim
   33    if(coord(j,k).ne.coord(i,k)) go to 32
         go to 10
   32    continue
         go to 90
      else
         ir=dir(i-1)
         do 35 j=1,i-4
         jr=dir(j)
         do 36 k=1,dim
         if(k.eq.jr.or.k.eq.ir) go to 36
         if(coord(j,k).ne.coord(i-1,k)) go to 35
   36    continue
         if(jr.ne.ir) then
            j1=coord(i,ir)-coord(j,ir)
            j2=coord(i-1,ir)-coord(j,ir)
            i1=coord(j+1,jr)-coord(i,jr)
            i2=coord(j,jr)-coord(i,jr)
            if((j1.eq.0.or.j2.eq.0).or.(isign(j1,j2).ne.j1)) then
               if((i1.eq.0.or.i2.eq.0).or.(isign(i1,i2).ne.i1)) go to 10
            endif
         else
            len=max0(coord(i,ir),coord(i-1,ir),coord(j+1,jr),
     c               coord(j,jr))-min0(coord(i,ir),coord(i-1,ir),
     c               coord(j+1,jr),coord(j,jr))
            if(len.le.(iabs(coord(i,ir)-coord(i-1,ir))
     c                +iabs(coord(j+1,jr)-coord(j,jr)))) go to 10
         endif
   35    continue
      endif
c
c
c
   90 v1old=v1
      do 34 id=1,dim
      v3(id)=v3(id)+dfloat(coord(i,id))
      if(iabs(coord(i,id)).gt.maxc(id)) maxc(id)=iabs(coord(i,id))
   34 v1=v1+dfloat(coord(i,id))**2
c
c
c
   30 continue
      go to 40
c
c
c
   10 continue
      write(6,45) ntri
   45 format(1x,'Failed after ',i4,' trials.')
      go to 200
c
c
c
   40 v2=0.d0
      do 50 id=1,dim
   50 v2=v2+v3(id)**2
      rg=(dfloat(nstep)*v1-v2)**0.5/dfloat(nstep)
      r=dsqrt(v1-v1old)
      write(6,91) r,rg,iseed,(maxc(id),id=1,dim)
   91 format(1x,'r = ',f12.4,', rg= ',f12.4,', seed=',i20/
     c1x,'abs(maxcoord) = ',7i10)
c     write(6,92) (coord(i,1),coord(i,2),i=1,nstep+1)
c  92 format(1x,2i10)
c
c
c
   70 if(dim.ge.3) then
         write(6,26) dim
   26    format(1x,'projection onto which 2 axes from 1 - ',i1,' ?')
         read(5,*) m1,m2
      else
         m1=1
         m2=2
      endif
      write(6,14)
   14 format(1x,'title ?')
      read(5,15) title1
   15 format(a35)
      titlea=title1//null
      write(6,19)
   19 format(1x,'subtitle ?')
      read(5,15) title2
      titleb=title2//null
      write(6,17) 
   17 format(1x,'display 6x6 grid ?')
      read(5,5) agrid
      write(6,24)
   24 format(1x,'set scale ?')
      read(5,5) ascale
      if(ascale.eq.yes) then
         write(6,18)
   18    format(1x,'max +window size ?')
         read(5,*) win
      else
         win=amax0(maxc(m1),maxc(m2))
      endif
      grid=0.33333333*win
      write(6,23)
   23 format(1x,'mark end points ?')
      read(5,5) amark
c
c
c
      call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
      call pltlun(19,.true.,.false.)
      call pltlfn('graph.out')
      call plots
      call plot(0.,3.5,-3)
      if(agrid.eq.yes) then
         call plot(3.,1.,3)
         call plot(3.,7.,2)
         call plot(4.,7.,3)
         call plot(4.,1.,2)
         call plot(5.,1.,3)
         call plot(5.,7.,2)
         call plot(6.,7.,3)
         call plot(6.,1.,2)
         call plot(7.,1.,3)
         call plot(7.,7.,2)
         call plot(8.,7.,3)
         call plot(8.,1.,2)
         call plot(9.,1.,3)
         call plot(9.,7.,2)
         call plot(3.,7.,2)
         call plot(3.,6.,3)
         call plot(9.,6.,2)
         call plot(9.,5.,3)
         call plot(3.,5.,2)
         call plot(3.,4.,3)
         call plot(9.,4.,2)
         call plot(9.,3.,3)
         call plot(3.,3.,2)
         call plot(3.,2.,3)
         call plot(9.,2.,2)
         call plot(9.,1.,3)
         call plot(3.,1.,2)
      endif
c
c     Draw a cross at the center.
c
      call plot(6.,4.,-3)
      call plot(-0.1,-0.1,3)
      call plot(0.1,0.1,2)
      call plot(0.1,-0.1,3)
      call plot(-0.1,0.1,2)
      call plot(0.,0.,3)
      x1=0.
      y1=0.
c
c
c
      do 60 i=2,nstep+1
      x=float(coord(i,m1))/grid
      y=float(coord(i,m2))/grid
      if((ascale.ne.yes).or.
     c   ( (abs(x).le.3..and.abs(x1).le.3.)
     c     .and.(abs(y).le.3..and.abs(y1).le.3.) )) then
         call plot(x1,y1,3)
         call plot(x,y,2)
         if(amark.eq.yes) then
            call symbol(x,y,0.1,mark,0.)
         endif
      endif
      x1=x
      y1=y
   60 continue
c
c
c
      call plot(-6.,-4.,-3)
      call symbol(2.5,9.,0.3,titlea,0.)
      call symbol(2.5,8.5,0.2,titleb,0.)
      call plot(0.,0.,-3)
      call plot(0.,0.,999)
      read(5,5) ans
c
c
c
      write(6,13) more
   13 format(1x,a1)
  200 write(6,16)
   16 format(1x,'another view ?')
      read(5,5) ans
      if(ans.eq.yes) go to 70
      write(6,62)
   62 format(1x,'another walk ?')
      read(5,5) ans
    5 format(a1)
      if(ans.eq.yes) go to 80
c
c
c
      stop
      end
c
c
c
      function ranf(i)
      i=i*1566083941+1
      ranf=real(i)*2.3283064e-10+0.5
      return
      end
