program cluster c c *********************************************************************** c Program to generate percolation clusters by breadth first search c *********************************************************************** c parameter (maxsize=50000,maxz=4,maxchem=1501,maxperim=100000) parameter (maxlp=10000) common list(2,0:maxsize) real etime,time(2) integer iflag(2) character*12 out 10 write(6,120) 120 format(' SQ Cluster: p, n, ntri, isd ?') read(5,*) p,n,ntri,isd write(6,121) 121 format(' file name to record cluster sites ?') read(5,122) out 122 format(a12) do 300 itri=1,ntri isd0=isd t0=etime(time) call bmake2(p,n,m2,iflag,isuc,isd) t1=etime(time) if(isuc.eq.1) then write(6,200) m2,t1-t0,iflag,isd0,isd 200 format(1x,'size=',i6,', CPU=',1p,e10.3,', flags=',2i1,2x,2i12) open(1,file=out) rewind(1) do 250 i=1,m2 250 write(1,251) list(1,i),list(2,i) 251 format(1x,2(i6,1x)) close(1) stop else write(6,520) itri,m2,iflag 520 format(' trial#',i3,', size=',i5,', flag=',2i1) endif 300 continue stop end c c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< c "bmake2" to make a percolation cluster by breadth first c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c subroutine bmake2(p,n,m2,iflag,isuc,isd) parameter (maxsize=50000,maxz=4,maxperim=100000) parameter (maxlp=10000) integer r1(2),r(2),dr(2,maxz),iflag(2) integer perim(2,maxperim) common list(2,0:maxsize) integer sitepl(2*maxlp-1),perimpl(2*maxlp-1) c,sitech(maxsize),perimch(maxperim) c c----------------------------------------------------------------------- c Function "lp" gives the plane number c----------------------------------------------------------------------- c lp(i,j)=mod(100*i+99*j,maxlp)+maxlp c c----------------------------------------------------------------------- c The following block specifies lattice geometry for SQ. c----------------------------------------------------------------------- c data dr/0,1,0,-1,1,0,-1,0/ c c----------------------------------------------------------------------- c Initialize exit conditions c----------------------------------------------------------------------- c ncord=4 lim=maxsize if(n.lt.0) lim=min0(-n,maxsize) nok=min0(iabs(n),lim) c c----------------------------------------------------------------------- c Put origin on stack c----------------------------------------------------------------------- c list(1,1)=0 list(2,1)=0 list(1,0)=-1000000 list(2,0)=-1000000 c c----------------------------------------------------------------------- c Initialize planes and chains for hash coding c----------------------------------------------------------------------- c do 130 i=1,2*maxlp-1 sitepl(i)=0 perimpl(i)=0 130 continue sitepl(lp(0,0))=1 do 140 i=1,maxsize sitech(i)=0 140 continue do 150 i=1,maxperim perimch(i)=0 150 continue c c----------------------------------------------------------------------- c m1 is the lower pointer where link check is being done and c m2 is the upper pointer where the cluster has grown up to. c m3 is the upper index for perim and c----------------------------------------------------------------------- c m1=1 m2=1 m3=0 c c----------------------------------------------------------------------- c No neighbor yet for the origin. Set flags to OK for now. c----------------------------------------------------------------------- c iflag(1)=0 iflag(2)=0 c c====================================================================== c START GENERATION OF CLUSTER c====================================================================== c 300 continue r(1)=list(1,m1) r(2)=list(2,m1) mk2=m2 mk3=m3 c c c do 400 icord=1,ncord r1(1)=r(1)+dr(1,icord) r1(2)=r(2)+dr(2,icord) c c----------------------------------------------------------------------- c Check if this neighbor is already on the stack c----------------------------------------------------------------------- c lp1=lp(r1(1),r1(2)) n1=sitepl(lp1) if(n1.eq.0) go to 425 n2=n1 427 if(r1(1).eq.list(1,n2).and.r1(2).eq.list(2,n2)) go to 400 n2=sitech(n2) if(n2.eq.0) go to 425 go to 427 c c---------------------------------------------------------------------- c If any of the limits has been reached, then no more addition of c sites is done. c---------------------------------------------------------------------- c 425 if((iflag(1)+iflag(2)).gt.0) go to 500 c c------------------------------------------------------------------------ c Check if this site has already been decided to be a perimeter site. c------------------------------------------------------------------------ c n3=perimpl(lp1) if(n3.eq.0) go to 430 n4=n3 437 if(r1(1).eq.perim(1,n4).and.r1(2).eq.perim(2,n4)) go to 400 n4=perimch(n4) if(n4.eq.0) go to 430 go to 437 c c----------------------------------------------------------------------- c Generate new site with probability p. c----------------------------------------------------------------------- c 430 if(ranf(isd).le.p) go to 440 m3=m3+1 perim(1,m3)=r1(1) perim(2,m3)=r1(2) perimpl(lp1)=m3 perimch(m3)=n3 if(m3.eq.maxperim) iflag(2)=1 go to 400 c c------------------------------------------------------------------------ c Add the new site on 'list' c------------------------------------------------------------------------ c 440 m2=m2+1 list(1,m2)=r1(1) list(2,m2)=r1(2) sitepl(lp1)=m2 sitech(m2)=n1 if(m2.eq.lim) iflag(1)=1 c c----------------------------------------------------------------------- c Advance the lower pointer m1 to check the next site on 'list'. c----------------------------------------------------------------------- c 400 continue if(m1.eq.m2) go to 500 m1=m1+1 go to 300 c c======================================================================= c GENERATION ATTEMPT IS FINISHED c======================================================================= c 500 if(m2.ge.nok) then isuc=1 else isuc=0 endif c c c return end c c c real function ranf(iseed) iseed=iseed*1566083941+1 ranf=float(iseed)*2.328306e-10+0.5 return end