1000 program cluster
1010 !
1020 ! ***********************************************************************
1030 !     Program to generate percolation clusters by breadth first search
1040 ! ***********************************************************************
1050 !
1060 option nolet
1070 dim list(0,0),iflag(2)
1080 maxsize=50000
1090 maxz=4
1100 maxperim=100000
1101 maxlp=10000
1110 mat redim list(2,0 to maxsize)
1120 input prompt "SQ Cluster:  p, n, ntri => ": p,n,ntri
1130 line input prompt "file name to record cluster sites => ": out$
1140 randomize
1150 for itri=1 to ntri
1160    t0=time
1170    call bmake2(p,n,m2,iflag(),isuc,list(,),maxsize,maxz,maxperim,maxlp)
1180    t1=time
1190    cpu=round(t1-t0,2)
1200    if isuc=1 then
1202       set background color "white"
1204       clear
1206       set cursor 1,1
1208       print "p="; p; ", n="; m2,", CPU=",cpu;
1220       print ", flag=",iflag(1),iflag(2)
1230       if out$ <> "" then
1240          open #1: name out$, organization text, create newold
1250 	      erase #1
1260          for i=1 to m2
1270             print #1: list(1,i),list(2,i)
1280          next i
1290          close #1
1300       end if
1310       xmin=list(1,1)
1320       xmax=list(1,1)
1330       ymin=list(2,1)
1340       ymax=list(2,1)
1350       for i=2 to m2      
1360           if list(1,i)<xmin then xmin=list(1,i)
1370           if list(1,i)>xmax then xmax=list(1,i)
1380           if list(2,i)<ymin then ymin=list(2,i)
1390           if list(2,i)>ymax then ymax=list(2,i)
1400       next i
1410       call display(list(,),m2,xmin,xmax,ymin,ymax)
1420       stop
1430    else
1440   	   print "trial#",itri,", size=",m2,", flag=",iflag(1),iflag(2)
1450    end if
1460 next itri
1470 end
1471
1472 def lp(i,j,maxlp)
1474    lp=mod(100*i+99*j,maxlp)+maxlp
1476 end def
1480 
1490 ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1500 !     "bmake2" to make a percolation cluster by breadth first
1510 ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1520 
1530 sub bmake2(p,n,m2,iflag(),isuc,list(,),maxsize,maxz,maxperim,maxlp)
1540 dim r1(2),r(2),dr(0,0),perim(0,0)
1542 dim sitepl(0),perimpl(0),sitech(0),perimch(0)
1544 mat redim sitepl(2*maxlp-1),perimpl(2*maxlp-1),sitech(maxsize),perimch(maxperim)
1550 mat redim dr(2,maxz),perim(2,maxperim)
1552 declare def lp
1560 !
1570 !-----------------------------------------------------------------------
1580 !    The following block specifies lattice geometry for SQ.
1590 !-----------------------------------------------------------------------
1600 !
1610 data 0,0,1,-1,1,-1,0,0 
1620 mat read dr
1630 !
1640 !-----------------------------------------------------------------------
1650 !                   Initialize exit conditions
1660 !-----------------------------------------------------------------------
1670 !
1680 lim=maxsize
1690 if n<0 then lim=min(-n,maxsize)
1700 nok=min(abs(n),lim)
1710 !
1720 !-----------------------------------------------------------------------
1730 !             Put origin on stack
1740 !-----------------------------------------------------------------------
1750 !
1760 list(1,1)=0
1770 list(2,1)=0
1778 !
1780 !------------------------------------------------------------------------
1782 !             Initialize planes and chains for hash coding
1784 !------------------------------------------------------------------------
1786 !
1791 for i=1 to 2*maxlp-1
1792    sitepl(i)=0
1793    perimpl(i)=0
1794 next i
1795 sitepl(lp(0,0,maxlp))=1
1796 for i=1 to maxsize
1797    sitech(i)=0
1798 next i
1799 for i=1 to maxperim
1800    perimch(i)=0
1801 next i
1805 !
1810 !-----------------------------------------------------------------------
1820 !     m1 is the lower pointer where link check is being done and
1830 !     m2 is the upper pointer where the cluster has grown up to.
1840 !     m3 is the upper index for perim and
1850 !-----------------------------------------------------------------------
1860 !
1870 m1=1
1880 m2=1
1890 m3=0
1900 !
1910 !-----------------------------------------------------------------------
1920 !          No neighbor yet for the origin.  Set flags to OK for now.
1930 !-----------------------------------------------------------------------
1940 !
1950 iflag(1)=0
1960 iflag(2)=0
1970 !
1980 !======================================================================
1990 !               START GENERATION OF CLUSTER
2000 !======================================================================
2010 !
2020 do
2030    r(1)=list(1,m1)
2040    r(2)=list(2,m1)
2070 !
2080 !
2090 !
2100    for icord=1 to maxz
2110       r1(1)=r(1)+dr(1,icord)
2120       r1(2)=r(2)+dr(2,icord)
2130 !
2140 !-----------------------------------------------------------------------
2150 !     Check if this neighbor is already on the stack
2160 !-----------------------------------------------------------------------
2170 !
2172       lp1=lp(r1(1),r1(2),maxlp)
2174       n1=sitepl(lp1)
2176       if n1=0 then 2320
2178       n2=n1
2180       if r1(1)=list(1,n2) and r1(2)=list(2,n2) then 2780
2182       n2=sitech(n2)
2184       if n2=0 then 2320
2186       goto 2180
2260 !
2270 !----------------------------------------------------------------------
2280 !     If any of the limits has been reached, then no more addition of
2290 !     sites is done.
2300 !----------------------------------------------------------------------
2310 !
2320       if (iflag(1)+iflag(2))>0 then exit do
2330 !
2340 !------------------------------------------------------------------------
2350 !     Check if this site has already been decided to be a perimeter site.
2360 !------------------------------------------------------------------------
2370 !
2372       n3=perimpl(lp1)
2374       if n3=0 then 2520
2376       n4=n3
2378       if r1(1)=perim(1,n4) and r1(2)=perim(2,n4) then 2780
2380       n4=perimch(n4)
2382       if n4=0 then 2520
2384       goto 2378
2460 !
2470 !-----------------------------------------------------------------------
2480 !     Generate new site with probability p.
2490 !-----------------------------------------------------------------------
2500 !
2520       if rnd>p then
2540          m3=m3+1
2550          perim(1,m3)=r1(1)
2560          perim(2,m3)=r1(2)
2562          perimpl(lp1)=m3
2564          perimch(m3)=n3
2570          if m3=maxperim then iflag(2)=1
2580       else
2600              
2610 !
2620 !------------------------------------------------------------------------
2630 !     Add the new site on 'list'
2640 !------------------------------------------------------------------------
2650 !
2660          m2=m2+1
2670          list(1,m2)=r1(1)
2680          list(2,m2)=r1(2)
2682          sitepl(lp1)=m2
2684          sitech(m2)=n1
2690          if m2=lim then iflag(1)=1
2720       end if
2730 !
2740 !-----------------------------------------------------------------------
2750 !     Advance the lower pointer m1 to check the next site on 'list'.
2760 !-----------------------------------------------------------------------
2770 !
2780    next icord
2790    m1 = m1+1
2800 loop while m1 <= m2
2810 !
2820 !=======================================================================
2830 !                    GENERATION ATTEMPT IS FINISHED
2840 !=======================================================================
2850 !
2860 if m2>=nok then
2870    isuc=1
2880 else
2890    isuc=0
2900 end if
2910 !
2920 !
2930 !
2940 end sub
2950 
2960 ! Display results
2970 sub display(list(,),n,xmin,xmax,ymin,ymax)
2996    open #1: screen 0,1,0,0.98
3000    set window xmin-1,xmax+1,ymin-1,ymax+1
3010    set color "black"
3020    plot xmin-0.5,ymin-0.5;
3030    plot xmin-0.5,ymax+0.5;
3040    plot xmax+0.5,ymax+0.5;
3050    plot xmax+0.5,ymin-0.5;
3060    plot xmin-0.5,ymin-0.5
3090    set color "red"
3100    for i=1 to n
3110       plot area: list(1,i)-0.25,list(2,i)-0.25;list(1,i)-0.25,list(2,i)+0.25;list(1,i)+0.25,list(2,i)+0.25;list(1,i)+0.25,list(2,i)-0.25;list(1,i)-0.25,list(2,i)-0.25
3120    next i
3122    close #1
3130    get key z
3140 end sub