!-----------------------------------------------------------------------------------------------------------------------------------
! Program to study 2D site percolation on an N = L x L square lattice. It also  computes the distribution of cluster sizes.
!
! 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.
!-----------------------------------------------------------------------------------------------------------------------------------
! This code also computes the distribution of cluster sizes.
!
! S_avg: The average size of finite clusters (i.e. excluding the largest one)
!
! cluster     (1:ncluster ): contains the sites of the lattice that are the roots of the clusters. The ordering is done in descending cluster size. cluster(1) has the root of the largest cluster.
! cluster     (1:ncluster2): contains the  roots of the clusters whose size is larger than 1
! cluster     (ncluster+1:): contains garbage
! 
! clustersizes(1:ncluster ): clustersizes(k) contains the size of the cluster with root cluster(k).  k1 > k2 => clustersizes(k1) >= clustersizes(k2) 
! 
! histogram(s)             : The fraction of clusters of size s.
!
! subroutine print_clusters: Prints occupied sites' coordinates (i,j), together with the "color" of the cluster they belong to.
!                            The color of the cluster is an integer that is increasing as the size of the lattice is decreasing. (color(cluster(k)) = k)
!                            The largest cluster has color = 1.
!
! Sorting of cluster sizes is done using shell sort.
!-----------------------------------------------------------------------------------------------------------------------------------
!Compile: gfortran -O2 perc03.f90 -o p
!Run:     echo <L> <p> <nmeas> | ./p   > perc.out
!Output:  grep ^C p.out ----> C  <P_N>  <S_avg> <ncluster> <ncluster2>
!         grep ^H p.out ----> H  <size> <histogram(size)>
!         grep ^P p.out ----> P  imeas  i  j  color             (i.e. the site (i,j) and the color of the cluster that it belogs to)
!-----------------------------------------------------------------------------------------------------------------------------------
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
 !----------------------------------------------------------------------------------------------------------------------------------
 integer(i4), allocatable  :: cluster(:)       ! Stores indices of root sites
 integer(i4), allocatable  :: clustersizes(:)  ! Stores sizes of these roots
 real   (dp), allocatable  :: histogram(:)
 integer(i4)               :: ncluster         ! Pointer to the last root (Total count)
 integer(i4)               :: ncluster2        ! Pointer to the last cluster with size > 1
 integer(i4)               :: tmp_idx, tmp_sz, gap, k_sort, m_sort
 logical                   :: swapped
 real(dp)                  :: S_avg, S_num, S_denom

 ! --- 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),cluster(N), clustersizes(N),histogram(N))
 histogram = 0.0_dp

 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 Roots, Collect, and Sort
  ncluster    = 0
  ncluster2   = 0
  max_cluster = 0
  
  !---------------------------------------------------------------------------------------------------------------------------------
  ! Step A: Collect all valid roots
  siteloopj2: do j = 1, L    ; do i = 1, L
   if (lattice(i, j)) then
    idx                     = (j - 1) * L + i
           
    ! Check if this site is a root:
    if (parent(idx)        == idx) then
     ncluster               = ncluster + 1
     cluster     (ncluster) = idx
     clustersizes(ncluster) = sz(idx)
              
     ! Track max_cluster:
     if (sz(idx) > max_cluster) max_cluster = sz(idx)
    end if ! if (parent(idx)        == idx)
   end if  ! if (lattice(i, j))
  end do                    ; end do siteloopj2

  !---------------------------------------------------------------------------------------------------------------------------------
  ! Step B: Sort arrays in DESCENDING order of size (Shell Sort)
  ! We sort 'clustersizes' and swap 'cluster' accordingly to keep them synced.
  gap                       = ncluster / 2
  sortloop: do while (gap   > 0)
   ksortloop: do k_sort     = gap + 1, ncluster
    tmp_sz                  = clustersizes       (k_sort)
    tmp_idx                 = cluster            (k_sort)
    m_sort                  =                     k_sort
        
    do while (    m_sort > gap .and. clustersizes(m_sort - gap) < tmp_sz)
     clustersizes(m_sort)   =        clustersizes(m_sort - gap)
     cluster     (m_sort)   =        cluster     (m_sort - gap)
     m_sort                 =                     m_sort - gap
    end do
        
    clustersizes(m_sort)    = tmp_sz
    cluster     (m_sort)    = tmp_idx
   end do ksortloop
   gap                      = gap / 2
  end  do  sortloop

  !---------------------------------------------------------------------------------------------------------------------------------
  ! --- Step C: Set Pointers ---
  ! Find ncluster2: The index of the last cluster with size > 1.
  ! Since list is DESCENDING (e.g., 100, 50, 2, 1, 1), we just count from start.
  ncluster2    = 0
  if (ncluster > 0) then
   max_cluster = clustersizes(1) ! First element is now guaranteed max
      
   do k = 1, ncluster
    if (clustersizes(k) > 1) then
     ncluster2 = k
    else
     exit ! As soon as we hit a '1', we know the rest are also '1'
    end if
   end do

  end if ! if (ncluster > 0)
  
  !---------------------------------------------------------------------------------------------------------------------------------
  ! --- Step D: Measurements of distribution of sizes and average size
  do k = 1, ncluster
   histogram(clustersizes(k)) = histogram(clustersizes(k)) + 1
  end do

  !S_avg has the average size of FINITE clusters, i.e. omiting the largest one. Trivial clusters must be included.
  !S(p) ~ (p-pc)^{-γ}
  !Probability that a  randomly chosen occupied site (belonging to a finite cluster) belongs to a cluster of size s is:
  !  w(s) = s . n(s)/(sum_s' s' . n(s'))                            (Sums exclude the largest cluster)
  !    S  = sum_s s . w(s) = sum_s s^2 . n(s) / sum_s s . n(s)
  S_avg   = 0.0_dp
  S_num   = sum( real(clustersizes(2:ncluster),dp)**2) !if ncluster = 1, then sum returns 0 for empty slices
  S_denom = sum( real(clustersizes(2:ncluster),dp)   ) !avoid integer overflow: convert to real before taking the square
  if(S_denom > 0.0_dp) S_avg = S_num/S_denom

  !---------------------------------------------------------------------------------------------------------------------------------
  ! --- Print results:
  print                      '(A,2G28.16,2I12)','C ',  real(max_cluster ,dp)/real(N,dp),S_avg,ncluster,ncluster2
  ! --- Print configuration:
  if(imeas <= 50) call print_clusters
 end do loopmeas

 do k = 1, N
  if(histogram(k) > 0) print '(A,I12,G28.16  )','H ',k,real(histogram(k),dp)/real(nmeas,dp)
 end do

!-----------------------------------------------------------------------------------------------------------------------------------
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
 !----------------------------------------------------------------------------------------------------------------------------------
 ! Print configuration:
 subroutine     print_clusters
  integer(i4), allocatable :: color(:)
  integer(i4)              :: k, i, j, idx, root

  if(ncluster < 1 ) return

  allocate(color(N))

  color  = 0
  do k   = 1, ncluster
   color(cluster(k)) = k    ! give a color to each cluster in descending size order
  end do

  do j   = 1, L ; do i = 1, L
   if(lattice(i,j))then
    idx  = (j-1)*L + i
    root = find_root(idx)
    print '(A,4I12)','P ',imeas,i,j,color(root)
   end if
  end do      ; end do

 end subroutine print_clusters
 !----------------------------------------------------------------------------------------------------------------------------------
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
!-----------------------------------------------------------------------------------------------------------------------------------