! 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
