function percolation_simulation ( m, n, p )

%*****************************************************************************80
%
%% percolation_simulation() simulates and analyzes a percolation system.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 July 2022
%
%  Author:
%
%    Ian Cooper
%    Modifications by John Burkardt
%
%  Input:
%
%    integer m, n: the number of rows and columns in the grid.
%
%    real p: the probability that a given site should be occupied.
%

%
%  Start the random number generator with a random seed.
%
% rng ( 'shuffle' );
%
%  The grid variable U is 0 or 1.
%
  r = rand ( m, n );
  u = zeros ( m, n );
  u ( r < p ) = 1.0;
%
%  Clusters
%  cls = cluster labels 1, 2, ... cluster_num
%
  [ cls, cluster_num ] = bwlabel ( u, 4 );
%
%  Number of occupied sites.
%
  nosites = sum ( sum ( u == 1.0 ) );
%
%  Probability of a site being occupied.
%
  posites = nosites / m / n;
%
%  Number of occupied sites in each cluster.
%
  cluster_size = zeros ( cluster_num, 1 );
  for c = 1 : cluster_num
    cluster_size(c) = length ( cls ( cls == c ) );
  end
%
%  Counter from 1 to max number of occupied sites in a cluster.
%
  s = 1 : max ( cluster_size );
%
%  ns = cluster size distribution.
%
  ns = zeros ( max ( cluster_size ), 1 );
  for c = s
    ns(c) = length ( cluster_size ( cluster_size == c ) );
  end
%
%  Probability that a site chosen at random is part of
%  an s-site cluster.
%
  probros = ns ./ sum ( ns );
%
%  Check for spanning clusters.
%  Does a "c" occur in column 1 and column n?
%  Does a "c" occur in row 1 and row m?
%
  isspanx = zeros ( cluster_num, 1 );
  isspany = zeros ( cluster_num, 1 );
  for c = 1 : cluster_num

    col1 = isempty ( find ( cls ( cls(:,1) == c ) ) );
    colN = isempty ( find ( cls ( cls(:,n) == c ) ) );
    if ( col1 == 0 && colN == 0 )
      isspanx(c) = 1;
    end

    row1 = isempty ( find ( cls ( cls(1,:) == c ) ) );
    rowN = isempty ( find ( cls ( cls(m,:) == c ) ) );
    if ( row1 == 0 && rowN == 0 )
      isspany(c) = 1;
    end
  end
%
%  Output
%
  fprintf ( 1, '\n' )
  fprintf ( 1, '  Total number of sites = %d\n', m * n );
  fprintf ( 1, '  Number of occupied sites = %d\n', nosites );
  fprintf ( 1, '  Observed probability of occupation = %g\n', posites );
  fprintf ( 1, '  Number of clusters = %d\n', cluster_num );
  fprintf ( 1, '  Mean cluster size = %g\n', mean ( cluster_size ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Cluster size distribution\n' );
  fprintf ( 1, '    Size  Number  Probability\n' );
  fprintf ( 1, '\n' );
  for c = s
    if ( 0 < ns(c) )
      fprintf ( 1, '       %2d  %2d  %g\n', s(c), ns(c), probros(c) );
    end
  end

  spanx_num = sum ( isspanx(1:cluster_num) );
  spany_num = sum ( isspany(1:cluster_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of horizontal spanning clusters = %d\n', spanx_num )
  fprintf ( 1, '  Number of vertical spanning clusters   = %d\n', spany_num )
%
%  Graphics.
%
  figure ( 1 );
  set ( gcf, 'units', 'normalized', 'position', [ 0.1, 0.52, 0.23, 0.32 ] );
  hold ( 'on' );
  for cx = 1 : n
    for cy = 1 : m
      h_rectangle = rectangle ( 'Position', [ cx-0.5, cy-0.5, 1, 1 ] );
      if ( u(cy,cx) == 0.0 )
        set ( h_rectangle, 'facecolor', [ 1, 1, 1 ] );
      else
        set ( h_rectangle, 'facecolor', [ 0, 0, 1 ], 'edgecolor', [ 0, 0, 1] );
      end
    end
  end

  for c = 1 : cluster_num
    [ rn, cn ] = find ( cls == c );
    text ( cn(1)-0.5, rn(1), num2str ( c, '%2.0f' ), 'color', 'w' );
  end

  title ( 'percolation occupancy' )
  axis ( 'equal' );
  axis ( 'tight' );
  axis ( 'off' );

  filename = 'occupation.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Second plot.
%
  hf2 = figure ( 2 );
  set ( gcf, 'units', 'normalized', 'position', [ 0.35,0.52,0.23,0.32 ] );
  img = label2rgb ( flipud ( cls ) );
  image ( img );
%
%  Add black grid lines to image.
%
  hold ( 'on' );
%
%  Vertical lines.
%
  for cx = 0 : n
    for cy = 0 : m

      yp = 0.5 : m + 0.5;
      xp = 0.5 + cx .* ones ( m + 1, 1 );
      plot ( xp, yp, 'k', 'linewidth', 0.5 );

    end
  end
%
%  Horizontal lines.
%
  for cx = 0 : n
    for cy = 0 : m

      xp = 0.5 : n + 0.5;
      yp = 0.5 + cy .* ones ( n + 1, 1 );
      plot ( xp, yp, 'k', 'linewidth', 0.5 );

    end
  end

  title ( 'percolation clusters' )
  axis ( 'equal' );
  axis ( 'tight' );
  axis ( 'off' );
  box ( 'on' );

  filename = 'clusters.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

