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