!-----------------------------------------------------------------------------------------------------------------------------------
! Program to study 2D site percolation on an N = L x L square lattice.
!
! Data structure:
! ===============
!
! Sites:  coordinates: (i,j)   index:  k = (L-1)*j + i   i,j = 1, ... , L  k = 1, ... , N
!                      (i,j)   has nearest neighbors (i,mod(j+/-1+N,N)), (mod(i+/-1+N,N),j)
!
! logical lattice(L,L): A site (i,j) is occupied when lattice(i,j) = .true.
! integer parent (N)  : parent(k) contains the index of the parent of site with index k 
!         parent (k) == k iff k is the root of a cluster
! integer sz     (N)  : sz(k) contains the size of the cluster with root k. If k is not a root, sz(k) contains garbage
! real    max_cluster : The fraction of the size of the largest cluster P_N  (order parameter)
!
! function find_root(k) result (root): Returns the root of the cluster the site k belongs to. 
!                                      It recursively computes parent(parent(...parent(k))) until we reach root. 
!                                      On the way it "compresses" the path, making root to be the parent of all visited sites.
!
! subroutine unite(k1,k2): It unites the clusters where the sites with indices k1, k2 belong to.
!                          It computes the roots root1, root2 of the two clusters, and makes the parent of the root of the smaller cluster to be the root of the larger.
!                          It updates the size of the new cluster sz(root), where root is the root of the new cluster.
!                          By calling find_root(k1) and find_root(k2), there is partial root path compression.
!-----------------------------------------------------------------------------------------------------------------------------------
!Compile: gfortran -O2 perc01.f90 -o p
!Run:     echo <L> <p> <nmeas> | ./p   > perc.out
!Output:  grep ^C p.out ----> C  <P_N>  
!-----------------------------------------------------------------------------------------------------------------------------------
program percolation_2D
 use iso_fortran_env, only: int32, real64
 implicit none
 integer       , parameter :: dp = real64, i4=int32
 ! --- Arrays ---
 ! parent:  stores the index of the parent site for Union-Find
 ! sz:      stores the size  of the cluster (only valid at the root)
 ! lattice: occupancy grid
 integer(i4)               :: L               ! Lattice dimension
 integer(i4)               :: N               ! Total sites
 real   (dp)               :: p
 integer(i4), allocatable  :: parent   (:)
 integer(i4), allocatable  :: sz       (:)
 logical    , allocatable  :: lattice  (:, :)
 real(dp)   , allocatable  :: rand_grid(:, :)
 integer(i4)               :: k, i, j
 integer                   :: nmeas,imeas
 integer(i4)               :: idx, n_idx      ! linear indices
 integer(i4)               :: ni, nj          ! neighbor coordinates
 integer(i4)               :: max_cluster

 ! --- Initialization ---
 call random_init(.false.,.true.)
 
 print '(A             )',"# Enter L, p, nmeas: "
 read  *                 ,         L, p, nmeas
 print '(A,I6,F13.8,I12)',"# ",    L, p, nmeas

 N = L * L
 allocate(parent(N),sz(N),lattice(L,L),rand_grid(L,L))

 loopmeas: do imeas = 1, nmeas

 ! 1. Generate Lattice (Vectorized)
  call random_number(rand_grid)
  where (rand_grid < p)
   lattice = .true.
  elsewhere
   lattice = .false.
  end where

  ! 2. Initialize Union-Find Arrays: Initially, every occupied site is its own root with size 1.
  ! Empty sites are size 0. ??? sz is set to 1, not 0 ??
  do k       = 1, N
   parent(k) = k
   sz    (k) = 1
  end do

  ! 3. Cluster Analysis: Loop over all sites
  siteloopj: do  j = 1, L ; do i = 1, L
   if (.not. lattice(i, j)) cycle
    
   idx = (j - 1) * L + i

   ! -- Right   Neighbor (Periodic) --
   ni = mod(i, L) + 1 ; nj = j
   if (lattice(ni, nj)) then
    n_idx = (nj - 1) * L + ni
    call unite(idx, n_idx)
   end if

    ! -- Bottom Neighbor (Periodic) --
   ni = i             ; nj = mod(j, L) + 1
   if (lattice(ni, nj)) then
    n_idx = (nj - 1) * L + ni
    call unite(idx, n_idx)
   end if
  end do ; end do siteloopj

  ! 4. Find Max Cluster
  ! We only check roots of occupied sites
  max_cluster = 0
  siteloopj2:do  j = 1, L ; do i = 1, L
   if (lattice(i, j)) then
    idx = (j - 1) * L + i
    ! The array sz has correct size for sites where parent(idx) == idx . For other sites, sz(idx) has "stale" data, since we don't zero size each time we change a parent.
    if  (parent(idx) == idx        ) then
     if (sz    (idx) >  max_cluster) max_cluster = sz(idx)
    !if you wish to measure distribution of cluster sizes, then print sx(idx) only here
    end if
   end if
  end do ; end do siteloopj2
  
  print '(G28.16)', real(max_cluster, dp) / real(N, dp)
 end do loopmeas

!-----------------------------------------------------------------------------------------------------------------------------------
contains
!-----------------------------------------------------------------------------------------------------------------------------------

 !----------------------------------------------------------------------------------------------------------------------------------
 ! Recursive Find with Path Compression: This routine finds the path to the root and makes all sites to have parent(idx) = root
 !                                       Unless we call this routine, parent(idx) does not necessarily point to the root
 ! We don't call this function, unless we need to know the root, this saves time when we call unite(i_in,j_in) 
 recursive function find_root(curr_idx) result(root)
  integer(i4), intent(in) :: curr_idx
  integer(i4)             :: root

  if (parent(curr_idx) == curr_idx) then
   root                =  curr_idx
  else
   ! Path compression: set parent directly to the root
   parent   (curr_idx) = find_root(parent(curr_idx))
   root                = parent(curr_idx)
  end if
 end function find_root
 !----------------------------------------------------------------------------------------------------------------------------------
 ! Weighted Union
 subroutine unite(i_in, j_in)
  integer(i4), intent(in) :: i_in, j_in
  integer(i4)             :: root_i, root_j
  
  root_i           = find_root(i_in)
  root_j           = find_root(j_in)

  if (     root_i /=    root_j) then
   ! Attach smaller tree to larger tree to keep depth low
   if (sz (root_i) < sz(root_j)) then
    parent(root_i) =    root_j               ! root_i is not root anymore. The rest of the cluster of root_i does not know that, yet, unless we call find_root for some such site.
                                             ! still, a site in the cluster of root_i has a parent in the same cluster, and find_root will eventually lead to root_j, since parent(root_i) = root_j (restoring the path)
    sz    (root_j) = sz(root_j) + sz(root_i) ! sz(root_i) contains garbage now: array sz contains correct data only when parent(idx) == idx
   else
    parent(root_j) =    root_i               ! root_j is not root anymore
    sz    (root_i) = sz(root_i) + sz(root_j) ! sz(root_j) contains garbage now
   end if
  end if
 end subroutine unite
 !----------------------------------------------------------------------------------------------------------------------------------
end program percolation_2D
!-----------------------------------------------------------------------------------------------------------------------------------
!===================================================================================================================================
!-----------------------------------------------------------------------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos, Physics Department, National Technical University of Athens, 2025
!  konstant@mail.ntua.gr, www.physics.ntua.gr/konstant
!  
!  This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as 
!  published by the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along with this program. If not, see http://www.gnu.org/licenses
!-----------------------------------------------------------------------------------------------------------------------------------