program cluster ! ! *********************************************************************** ! Program to generate percolation clusters only ! *********************************************************************** ! option nolet dim list(0,0),xiav(2),nn(0),latt(0,0) maxsize=250000 maxedge=500 maxnn=20000 maxs=15000 maxl=20000 input prompt "SQ Cluster: p, nedge, ntri => ": p,nedge,ntri nedge=min(nedge,maxedge) maxl=min(maxl,nedge^2) maxsize=min(maxsize,nedge^2) line input prompt "file name to record cluster sites => ": out$ line input prompt "file name to record cluster size distribution => ": dist$ mat redim latt(nedge,nedge),list(2,0 to maxsize),nn(maxnn) randomize jsuc=0 for itri=1 to ntri t0=time call make2(p,nmax,nedge,isuc,nsize,list,latt,nn,maxl,maxs,maxnn) t1=time cpu=round(t1-t0,2) if isuc=1 then jsuc=jsuc+1 if dist$ <> "" then call dist(nmax,xiav,px,nc,ng,nn,dist$,p,nedge,maxs,jsuc) end if if out$ <> "" then open #1: name out$, organization text, create newold if jsuc=1 then erase #1 else set #1: pointer end print #1: end if for i=1 to nsize print #1: list(1,i),list(2,i) next i close #1 end if call display(latt,p,ng,nsize,nc,nedge,cpu) else print "trial=",itri,", nmax=",nmax end if next itri end ! ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! subroutine "make" to make a percolation cluster ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! sub make2(p,nmax,nedge,isuc,nsize,list(,),latt(,),nn(),maxl,maxs,maxnn) ! ! ********************************************************************* ! * This subroutine generates a percolation cluster on the square * ! * lattice. * ! ********************************************************************* ! dim label(0),ndr(2,4) mat redim label(0 to maxl) data 0,1,0,-1,1,0,-1,0 mat read ndr ! ! ! for i=0 to maxl label(i)=0 next i ! !============================================================================== ! Start generating clusters in nedge x nedge grid. !============================================================================== ! ! Generate the first column from top to bottom. !------------------------------------------------------------------------------ ! In this program, column refers to (i,.) for fixed i even though this is ! actually the opposite of the order of matrix indices. This is simply ! because it is easier to interpret indices as (x,y). !------------------------------------------------------------------------------ ! ic=0 if rnd<=p then ic=ic+1 if ic>maxl then isuc=-2 exit sub end if latt(1,1)=ic label(ic)=1 else latt(1,1)=0 end if for j=2 to nedge if rnd<=p then if latt(1,j-1)>0 then latt(1,j)=latt(1,j-1) label(latt(1,j))=label(latt(1,j))+1 else ic=ic+1 if ic>maxl then isuc=-2 exit sub end if latt(1,j)=ic label(ic)=1 end if else latt(1,j)=0 end if next j ! ! Generate column 2 through column "nedge"; each one from top to bottom. ! for i=2 to nedge ! ! The top site of each column. ! if rnd<=p then if latt(i-1,1)>0 then latt(i,1)=latt(i-1,1) ip=latt(i,1) if label(ip)<0 then ip=-label(ip) label(ip)=label(ip)+1 else ic=ic+1 if ic>maxl then isuc=-2 exit sub end if latt(i,1)=ic label(ic)=1 end if else latt(i,1)=0 end if ! ! The rest of the column. ! for j=2 to nedge if rnd<=p then imx=max(latt(i-1,j),latt(i,j-1)) if latt(i,j-1)*latt(i-1,j)>0 and latt(i,j-1)<>latt(i-1,j) then latt(i,j)=min(latt(i-1,j),latt(i,j-1)) ip=imx if label(ip)<0 then ip=-label(ip) iq=latt(i,j) if label(iq)<0 then iq=-label(iq) if ip=iq then label(ip)=label(ip)+1 else jmin=min(ip,iq) jmax=max(ip,iq) label(jmin)=label(jmin)+label(jmax)+1 label(jmax)=-jmin for l=jmax+1 to ic if label(l)=-jmax then label(l)=-jmin next l end if else if latt(i,j-1)+latt(i-1,j)>0 then latt(i,j)=imx ip=imx if label(ip)<0 then ip=-label(ip) label(ip)=label(ip)+1 else ic=ic+1 if ic>maxl then isuc=-2 exit sub end if latt(i,j)=ic label(ic)=1 end if else latt(i,j)=0 end if next j ! ! ! next i ! !============================================================================== ! Finished generation. Find the cluster number of the largest cluster. !============================================================================== ! jmax=0 numb=0 id=0 for i=1 to ic if label(i)>0 then numb=numb+label(i) if label(i)>jmax then jmax=label(i) id=i end if end if next i ! !============================================================================== ! Set up list array for the maximal cluster. !============================================================================== ! ! First assign site number to each site of the cluster ! nsize=0 for i=1 to nedge for j=1 to nedge if label(latt(i,j))=-id or latt(i,j)=id then nsize=nsize+1 latt(i,j)=nsize else if latt(i,j)>0 then latt(i,j)=-latt(i,j) end if next j next i ! ! ! for i=1 to nedge for j=1 to nedge if latt(i,j)>0 then n=latt(i,j) list(1,n)=i list(2,n)=j end if next j next i ! !============================================================================== ! Finished setting up the output array. !============================================================================== !============================================================================== ! Set up "nn". !============================================================================== ! ns1=maxs for j=1 to maxnn nn(j)=0 next j for i=1 to ic ns=label(i) if ns>maxs then ns1=ns1+1 nn(ns1)=ns else if ns>0 then nn(ns)=nn(ns)+1 end if next i nmax=ns1 ! !============================================================================== ! Finished setting up the output arrays. !============================================================================== ! isuc=1 end sub ! ! ! sub dist(nmax,xiav(),px,nc,ng,nn(),dist$,p,nedge,maxs,jsuc) dim nna(8),nnb(8),nr(9) open #2: name dist$, organization text, create newold if jsuc=1 then erase #2 else set #2: pointer end end if set #2:zonewidth 6 for i=1 to 8 nna(i)=-1 nnb(i)=-1 next i ng=0 nc=0 js=0 xx=0 for i=1 to maxs if nn(i)<>0 then js=js+1 nna(js)=i nc=nc+nn(i) ng=ng+i*nn(i) xx=xx+i^2*nn(i) nnb(js)=nn(i) if js=8 then ! print #2: nna(1),nnb(1),nna(2),nnb(2),nna(3),nnb(3),nna(4),nnb(4),nna(5),nnb(5),nna(6),nnb(6),nna(7),nnb(7),nna(8),nnb(8) js=0 end if end if next i if nmax>maxs then for i=maxs+1 to nmax js=js+1 nnb(js)=1 nna(js)=nn(i) xx=xx+nn(i)^2 nc=nc+1 ng=ng+nn(i) if js=8 then ! print #2: nna(1),nnb(1),nna(2),nnb(2),nna(3),nnb(3),nna(4),nnb(4),nna(5),nnb(5),nna(6),nnb(6),nna(7),nnb(7),nna(8),nnb(8) js=0 end if next i end if if js<>0 then js1=js+1 for kk=js1 to 8 nna(kk)=-1 next kk ! print #2: nna(1),nnb(1),nna(2),nnb(2),nna(3),nnb(3),nna(4),nnb(4),nna(5),nnb(5),nna(6),nnb(6),nna(7),nnb(7),nna(8),nnb(8) end if if js=0 then js=8 px=nna(js) xiav(1)=0 xiav(2)=0 if ng<>0 then xiav(1)=xx/ng xiav(2)=(xx-px*px)/ng px=px/ng end if nna(1)=-1 ! print #2: nna(1),nnb(1),nna(2),nnb(2),nna(3),nnb(3),nna(4),nnb(4),nna(5),nnb(5),nna(6),nnb(6),nna(7),nnb(7),nna(8),nnb(8) conc=ng/nedge^2 set #2: zonewidth 12 print #2, using "#.#####^^^^": p; set #2: zonewidth 3 print #2: ", "; set #2: zonewidth 12 print #2, using "#.#####^^^^": conc; set #2: zonewidth 3 print #2: ", "; set #2: zonewidth 12 print #2, using "#.#####^^^^": px; set #2: zonewidth 3 print #2: ", "; set #2: zonewidth 12 print #2, using "#.#####^^^^": xiav(2); set #2: zonewidth 3 print #2: ", "; set #2: zonewidth 12 print #2: nc close #2 end sub ! Display results sub display(latt(,),p,n,nsize,nc,nedge,cpu) set background color "white" clear set cursor 1,1 print "p="; p; ", n="; n;", max="; nsize;", nedge=";nedge;", nc=";nc;", cpu=";cpu open #1: screen 0,1,0,0.98 set window -1,nedge+1,-1,nedge+3 set color "black" plot -0.5,-0.5; plot -0.5,nedge+0.5; plot nedge+0.5,nedge+0.5; plot nedge+0.5,-0.5; plot -0.5,-0.5 for i=1 to nedge for j=1 to nedge if key input then get key z c$=chr$(z) if c$="x" or c$="q" then exit sub end if if latt(i,j)>0 then set color "red" plot area: i-0.25,j-0.25;i-0.25,j+0.25;i+0.25,j+0.25;i+0.25,j-0.25;i-0.25,j-0.25 else if latt(i,j)<0 then set color "blue" plot area: i-0.25,j-0.25;i-0.25,j+0.25;i+0.25,j+0.25;i+0.25,j-0.25;i-0.25,j-0.25 end if next j next i get key z end sub