! from:
!http://www.physics.nau.edu/~hart/classes/550_spring_2004/hw_schedule/solutions/perc.f90
module perc
implicit none
private
public perccheck, cluster_update
contains
function perccheck(grid)
integer, intent(out) :: perccheck
integer, intent(in) :: grid(:,:)
integer i, N
perccheck = 0
N = size(grid,1)
do i = 1, maxval(grid)
if(any(grid(2,2:N-1)==i) .and. &
any(grid(N-1,2:N-1)==i) .and. &
any(grid(2:N-1,2)==i) .and. &
any(grid(2:N-1,:N-1)==i)) then ! ith cluster touches all sides--it percolates
perccheck = i
exit
endif
enddo
end function perccheck
subroutine cluster_update(x,y,grid)
integer, intent(in) :: x, y
integer, intent(inout) :: grid(:,:)
integer :: label, i, j
i = x + 1 ! offset the x and y values to compensate for the fact
j = y + 1 ! that there is an extra site on each edge
! Check each of the four neighbors of the given point to see if
! clusters need to be combined and relabeled
if(grid(i+1,j)/=0)then ! neighboring site is occupied
label = min(grid(i,j),grid(i+1,j)) ! select the smaller cluster number
where(grid == max(grid(i,j),grid(i+1,j))) ! relabel all the sites with the
grid = label ! higher cluster number
endwhere
endif
if(grid(i-1,j)/=0)then ! neighboring site is occupied
label = min(grid(i,j),grid(i-1,j)) ! select the smaller cluster number
where(grid == max(grid(i,j),grid(i-1,j))) ! relabel all the sites with the
grid = label ! higher cluster number
endwhere
endif
if(grid(i,j+1)/=0)then ! neighboring site is occupied
label = min(grid(i,j),grid(i,j+1)) ! select the smaller cluster number
where(grid == max(grid(i,j),grid(i,j+1))) ! relabel all the sites with the
grid = label ! higher cluster number
endwhere
endif
if(grid(i,j-1)/=0)then ! neighboring site is occupied
label = min(grid(i,j),grid(i,j-1)) ! select the smaller cluster number
where(grid == max(grid(i,j),grid(i,j-1))) ! relabel all the sites with the
grid = label ! higher cluster number
endwhere
endif
end subroutine cluster_update
end module perc
program percolation
use perc
use random_stuff
implicit none
integer, parameter :: L=22 ! size of grid
integer, dimension(0:L+1,0:L+1) :: site = 0 ! grid of sites
real :: p ! occupation ratio
real :: ratio ! ratio of sites in percolating cluster vs. number of total occupied sites
real :: rnd(2) ! 2 random numbers for generating random x, y positions
real :: invL ! inverse box edge length
integer :: i, j, loop ! general loop counters
integer :: x, y ! x, y positions (integer)
integer :: nsites = 0 ! number of occupied sites
integer :: iclust = 1 ! current cluster number
integer :: PercClusLab ! Cluster number (label) of the percolated cluster
call set_random_seed(0) ! "Randomly" set the random seed
invL = 1/real(L)
open(11, file="percavg.data", status = "old", position = "append")
do loop = 1, 20
site = 0
nsites = 0
iclust = 1
do ! Loop until percolation happens
call random_number(rnd) ! Select two random numbers
x = int(rnd(1)*L)+1 ! Change to random numbers between 1 and L (x, y positions)
y = int(rnd(2)*L)+1
if(site(x,y)==0) then ! it's unoccupied
site(x,y) = iclust ! mark the site as occupied
nsites = nsites + 1 ! keep track of the number of occupied sites
! check to see if this new site is already part of a cluster
call cluster_update(x,y,site) ! relabel any clusters as necessary
iclust = iclust + 1
if(perccheck(site)/=0) exit ! Stop adding sites if a cluster has percolated
endif
enddo
p = real(nsites)/L**2
ratio = count(site==perccheck(site))/real(nsites)
print *, "occupancy: ", p
print *, "ratio in percolating cluster: ", ratio
write(11,*) invL, p, ratio
enddo
close(11)
! Write cluster data to file
PercClusLab = site(x,y)
open(10,file="perc.data"); open(11,file="spanningcluster.data")
do i=1,L; do j=1,L
if(site(i,j)/=0) write(10,'2i10') j,i
if(site(i,j)==site(x,y)) write(11,'2i10') j,i ! write sites for spanning cluster only
enddo;enddo
close(10); close(11)
end program percolation