      program fract
      parameter(maxsize=10000,maxbin=1500,maxlp=2500)
      character ans,yes
      character*12 infile
      character*80 dummy
      integer v(2,maxsize)
      integer list(2,maxsize),nchain(maxsize),nplane(2*maxlp-1)
      dimension hist(maxbin)
      lp(i,j)=mod(100*i+99*j,maxlp)+maxlp
      yes='y'
c
c
c
   10 write(6,*) 'Input cluster file ?'
      read(5,11) infile
   11 format(a12)
      open(1,file=infile)
      rewind(1)
      write(6,*) '# header lines to skip ?'
      read(5,*) iskip
      do 30 i=1,iskip
         read(1,21) dummy
         write(6,22) dummy
   21    format(a80)
   22    format(1x,a80)
   30 continue
      nsize=0
      minx=100000
      miny=100000
      maxx=-1000
      maxy=-1000
      do 90 i=1,maxsize
         read(1,*,end=85) (v(j,i),j=1,2)
         if(v(1,i).gt.maxx) maxx=v(1,i)
         if(v(2,i).gt.maxy) maxy=v(2,i)
         if(v(1,i).lt.minx) minx=v(1,i)
         if(v(2,i).lt.miny) miny=v(2,i)
         nsize=nsize+1
   90 continue
   85 close(1)
      write(6,*) 'size=',nsize
      write(6,*) 'cutoff size, step (ints) ?'
      read(5,*) mbin,mdel
      mbin=min0(mbin,maxbin)
c
c
c
      do 45 k=1,mbin
         hist(k)=0.
   45 continue
      do 20 k1=1,mbin,mdel
         mx=(maxx-minx)/k1+1
         my=(maxy-miny)/k1+1
         do 25 i=1,2*maxlp-1
   25    nplane(i)=0
         do 26 j=1,maxsize
   26    nchain(j)=0
         index=0
         do 100 l=1,nsize
            do 40 i=1,mx
               minx0=minx+(i-1)*k1
               maxx0=minx0+k1
               if(v(1,l).ge.minx0.and.v(1,l).lt.maxx0) then
                  nx=i
                  go to 48
               endif
   40       continue
   48       do 60 j=1,my
               miny0=miny+(j-1)*k1
               maxy0=miny0+k1
               if(v(2,l).ge.miny0.and.v(2,l).lt.maxy0) then
                  ny=j
                  go to 68
               endif
   60       continue
   68       lp1=lp(nx,ny)
            n1=nplane(lp1)
            if(n1.eq.0) go to 125
            n2=n1
  127       if(nx.eq.list(1,n2).and.ny.eq.list(2,n2)) go to 100
            n2=nchain(n2)
            if(n2.eq.0) go to 125
            go to 127
  125       index=index+1
            list(1,index)=nx
            list(2,index)=ny
            nchain(index)=n1
            nplane(lp1)=index
            hist(k1)=hist(k1)+1.
  100    continue
   20 continue
c
c
c
      do 120 k=1,mbin,mdel
         write(2,130) k,hist(k)
  130    format(1x,i6,e15.8)
  120 continue
c
c
c
  140 write(6,*) 'another file ?'
      read(5,23) ans
   23 format(a1)
      if(ans.eq.yes) go to 10
c
c
c
      stop
      end
c
c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< RANF >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
      real function ranf(isd)
      isd=isd*1566083941+1
      ranf=float(isd)*2.328306e-10+0.5
      return
      end
