      program fract
      parameter(maxsize=10000,maxbin=1500)
      character ans,yes,ans1,ans2,ans3
      character*12 infile
      character*80 dummy
      integer v(2,maxsize),numb(maxbin)
      dimension dj(maxbin),hist(maxbin),hist2(maxbin)
      dimension dm(maxbin),dm2(maxbin)
      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,'; put spheres randomly ?'
      read(5,23) ans1
      if(ans1.eq.yes) then
         write(6,*) '#spheres, seed ?'
         read(5,*) ndisk,isd
      else
         write(6,*) 'sphere at origin only ?'
         read(5,23) ans3
      endif
      write(6,*) 'cut off at edges ?'
      read(5,23) ans2
      if(ans2.ne.yes) then
         write(6,*) 'cutoff radius, step (ints) ?'
         read(5,*) maxr,mdel
      endif
c
c
c
      do 45 k=1,maxbin
      numb(k)=0
      hist(k)=0.
      hist2(k)=0.
   45 continue
      msize=nsize
      if(ans1.eq.yes) msize=ndisk
      if(ans3.eq.yes) msize=1
      do 40 i=1,msize
         i0=i
         if(ans1.eq.yes) i0=min0(int(nsize*ranf(isd))+1,nsize)
         if(ans2.eq.yes) then
            mbin=min0(maxx-v(1,i0),v(1,i0)-minx
     c            ,maxy-v(2,i0),v(2,i0)-miny,maxbin)
            mdel=1
         else
            mbin=maxr
         endif
         do 55 k=1,mbin,mdel
         numb(k)=numb(k)+1
         dj(k)=0.
   55    continue
         do 60 j=1,nsize
            nsq=(v(1,j)-v(1,i0))**2+(v(2,j)-v(2,i0))**2
            if(nsq.gt.mbin**2) go to 60
            do 80 k=1,mbin,mdel
               if(nsq.gt.k**2) go to 80
               if(nsq.eq.k**2) then
                  dj(k)=dj(k)+0.5
               else
                  dj(k)=dj(k)+1
               endif
   80       continue
   60    continue
         do 100 k=1,mbin,mdel
            hist(k)=hist(k)+dj(k)
            hist2(k)=hist2(k)+dj(k)**2
  100    continue
   40 continue
c
c
c
      do 120 k=1,maxbin
         if(numb(k).eq.0) go to 120
         dm(k)=hist(k)/float(numb(k))
         dm2(k)=hist2(k)/float(numb(k))-dm(k)**2
         if(dm2(k).gt.0) then
            dm2(k)=sqrt(dm2(k))
         else
            dm2(k)=0.
         endif
         write(2,130) k,dm(k),dm2(k),numb(k)
  130    format(1x,i6,2(e15.8,2x),i6)
  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
