cc TO USE: New version, with standard choice of kmax~N^(1/(lambda-1)).
      program greedy
cc program to use (correction with ntrap2, correct # of traps).
      implicit none
      external kisset
      real*8 x(20000000),y(20000000),xprev,yprev
      real*8 p(0:30000000)
      real*8 rkiss,r,r1
      real*8 dmin,d,ddmin,xi0,yi0,theta,prod,xtemp,ytemp
      real*8 pn,norm,lambda,nmax,probm,sum
      integer seed,seed1,seed2,seed3,tmax
      integer eta(20000000),etamax,etamaxx,itest,nb,ktest
      integer iy,kisset,i,i0,ntrap,ntrap2,irun,it,im,n
      integer iangl,ipos,iruni,irunf,unit,diz,cent
      character*13 namerun
      character*40 resangl,respos
      character(len=10) nombres
      character(len=3) nbrun
      character(len=1) ccm
      character(len=1) cdm
      character(len=1) cm
      character(len=1) cc
      character(len=1) cd
      character(len=1) cu   
      nombres='0123456789'
cc ntrap traps in box of size 1*1.
cc parameters:
      ntrap=1000000
      tmax=20000
      namerun='1e6.t2e4.L3.0'
      lambda=3.0
      iruni=1
      irunf=10
cccc normalization and cut-off of  popularity distribution    
      etamaxx=int(10.0**(13./lambda))
      norm=0.
      do 6 n=1,etamaxx
       norm=norm+real(n)**(-lambda)
 6    continue
      probm=norm
      do 7 n=1,etamaxx 
       probm=probm-real(n)**(-lambda)
       nmax=real(ntrap)*probm/norm
       if(nmax.le.1.0) then
        etamax=n+1
        goto 77
       endif
 7    continue
 77   continue
      norm=0.
      do 9 n=1,etamax
       norm=norm+real(n)**(-lambda)
 9    continue
cc      write(*,*) 'etamaxx=',etamaxx, 'etamax=',etamax
      sum=0.
      p(0)=0.
      do 19 n=1,etamax
       sum=sum+real(n)**(-lambda)
       p(n)=sum/norm
 19   continue
cc      write(2,*) 'etamaxx=',etamaxx, 'etamax=',etamax
cc loop on runs and initialization:
      do 100 irun=iruni,irunf
      cent=irun/100
      diz=(irun-100*cent)/10
      unit=irun-100*cent-diz*10
      cc=nombres(cent+1:cent+1)
      cd=nombres(diz+1:diz+1)
      cu=nombres(unit+1:unit+1)
      nbrun=cc//cd//cu
      seed=irun*10
      seed1=(5*seed+2)**2
      seed2=(5*seed+1)**2
      seed3=(5*seed)**2
      iy=kisset(seed1,seed2,seed3)
      r=rkiss()
      itest=0  
      do 1 i=1,ntrap
       x(i)=rkiss()
       y(i)=rkiss()
       r=rkiss()
       if(r.le.0.5) then
        r=0.5-r
       else
        r=1.5-r
       endif
       do 17 n=1,etamax
        if(r.le.p(n).AND.r.ge.p(n-1)) then
         eta(i)=n
         goto 18
        endif
 17    continue
 18    continue    
 1    continue
      ntrap2=ntrap
cc         write(4,*) x(itest),y(itest)
cc checking distribution:
cc      write(*,*) 'ntrap2=',ntrap2
      do 192 n=1,etamax
       nb=0.
       do 193 i=1,ntrap2
        if(eta(i).eq.n) nb=nb+1
 193   continue
cc       write(3,*) n,nb
 192  continue   
ccccccccccccccccccccccc
cc open result files:
cc      resangl=
cc     &  namerun//'.'//nbrun//'.angl'
cc      iangl=1
cc      open(unit=iangl,file=resangl,
cc     &      form='formatted',status='unknown')
      respos=
     &  namerun//'.'//nbrun//'.pos'
      ipos=2
      open(unit=ipos,file=respos,
     &      form='formatted',status='unknown')
cccccccc start running:
cc starting site, i0 position of the walker
cc middle of the box:
      do 10 i=1,ntrap2
       if(x(i).ge.(0.47).AND.x(i).le.(0.53).AND
     &  .y(i).ge.(0.47).AND.y(i).le.(0.53)) then
        i0=i
	goto 11
       endif
 10   continue
 11   continue   
cc first "previous" site asigned on the same site:      
      if(i0.eq.ntrap2) then
       xprev=x(i0)
       yprev=y(i0)
      else
       xprev=x(i0)
       yprev=y(i0)
      endif
      write(ipos,*) x(i0),y(i0),eta(i0)
cccc loop on time:       
      do 2 it=1,tmax
       dmin=10.
       xi0=x(i0)
       yi0=y(i0)
cc finding site with lowest d/eta:       
       do 3 i=1,ntrap2      
        if(i.eq.i0) goto 3
        d=((xi0-x(i))**2+(yi0-y(i))**2)**.5/real(eta(i))
        if(d.le.dmin) then
	 dmin=d
	 im=i
	endif
 3     continue	
cc computing angle with previous move:
       prod=(x(im)-xi0)*(xi0-xprev)+(y(im)-yi0)*(yi0-yprev)   
       ddmin=((xi0-x(im))**2+(yi0-y(im))**2)**.5
       prod=prod/ddmin
       prod=prod/((xi0-xprev)**2+(yi0-yprev)**2)**.5
       theta=acos(prod)
cc       write(iangl,*) it, theta
cccc update:
       xtemp=x(i0)
       ytemp=y(i0)
cc remove i0 from the box
       x(i0)=100.
       y(i0)=100.
       eta(i0)=1
cc new i0:       
       i0=im      
       xprev=xtemp
       yprev=ytemp
       write(ipos,*) x(i0),y(i0),eta(i0)
 2    continue  
 100  continue
      end 
       
       
       
      function rkiss()
        real*8 rm
c       data rm / 4294967295.0d0 /
        data rm / 2147483648.0d0 /
        rkiss=dfloat(abs(kiss())+1)/rm
      return
      end

      function kiss()
c      data i,j,k/123456789,362436069,521288629/
c      data i,j,k/607504,6057433,609760129/
      common/rnd/i,j,k
      m(k,n)=ieor(k,ishft(k,n))
        i=69069*i+23606797
        j=m(m(j,17),-15)
        k=m(iand(m(k,18),2147483647),-13)
        kiss=i+j+k
      return
      end

      function kisset(ii,jj,kk)
      common/rnd/i,j,k
c       print*,' Current seeds:',i,j,k
      i=ii
      j=jj
      k=kk
      kisset=1
      return
      end
