function fire_simulation ( m, n, prob_spread )

%*****************************************************************************80
%
%% fire_simulation() simulates a forest fire.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 January 2019
%
%  Author:
%
%    This version by John Burkardt
%
%  Input:
%
%    integer M, N, specifies the forest as a square of
%    M x N trees.  Default is 20x20.
%
%    real PROB_SPREAD, the probability that fire will spread from
%    a burning tree to its neighbor.  Default is 0.25.
%
  UNBURNT = 0;
  BURNING = 1;
  BURNT = 2;

  if ( nargin < 1 )
    m = 20;
    fprintf ( '\n' );
    fprintf ( 'Program will use default value of M = %d\n', m );
  elseif ( ischar ( m ) )
    m = str2num ( m );
  end

  if ( nargin < 2 )
    n = 20;
    fprintf ( 'Program will use default value of N = %d\n', n );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  if ( nargin < 3 )
    prob_spread = 0.25;
    fprintf ( 'Program will use default value of PROB_SPREAD = %g\n', prob_spread );
  elseif ( ischar ( prob_spread ) )
    prob_spread = str2num ( prob_spread );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'fire_simulation():\n' );
  fprintf ( 1, '  A probabilistic simulation of a forest fire.\n' );
  fprintf ( 1, '  The forest size is %d x %d\n', m, n );
  fprintf ( 1, '  The probability of tree-to-tree spread is %g\n', prob_spread );

  hour = 0;

  while ( true )

    if ( hour == 0 )
      forest(1:m,1:n) = UNBURNT;
      i_ignite = randi ( [ 1, m] );
      j_ignite = randi ( [ 1, n] );
      forest(i_ignite,j_ignite) = BURNING;
      fprintf ( 1, '  Lightning strike at (%d,%d)\n', i_ignite, j_ignite );
    else
      forest = forest_burns ( m, n, forest, prob_spread, hour );
    end

    forest_display ( m, n, forest, hour, i_ignite, j_ignite )

    if ( sum ( sum ( forest == BURNING ) ) == 0 )
      break
    end

    hour = hour + 1;

  end
%
%  Display the final forest state.
%
  forest_print ( m, n, forest, i_ignite, j_ignite );
%
%  Report the percentage of forest burned.
%
  total = sum ( sum ( forest == BURNT ) );
  percent = 100.0 * total / m / n;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Percentage of forest burned = %6.2f%%\n', percent );
  fprintf ( 1, '  Forest burned for %d hours\n', hour );

  return
end
function forest = forest_burns ( m, n, forest, prob_spread, hour )

%*****************************************************************************80
%
%% forest_burns() models a single time step of the burning forest.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 November 2017
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns of trees.
%
%    integer FOREST(M,N), an
%    array with an entry for each tree in the forest.
%
%    real PROB_SPREAD, the probability that the fire will 
%    spread from a burning tree to an unburnt one.
%
%  Output:
%
%    integer FOREST(M,N), an array with an entry for each tree in the forest.
%
  UNBURNT = 0;
  BURNING = 1;
  BURNT = 2;

  forest_old = forest;
%
%  Unburnted trees may begin to smolder;
%  Smoldering trees ignite;
%  Burning trees burn down;
%
  for j = 1 : n
    for i = 1 : m
      if ( forest_old(i,j) == UNBURNT )
        im1 = max ( i - 1, 1 );
        ip1 = min ( i + 1, m );
        jm1 = max ( j - 1, 1 );
        jp1 = min ( j + 1, n );
        nabe = forest_old(im1:ip1,jm1:jp1);
        nabe_num = sum ( sum ( nabe == BURNING ) );
        pnot = ( 1.0 - prob_spread ) ^ nabe_num;
        if ( pnot < rand ( ) )
          forest(i,j) = BURNING;
        end
      elseif ( forest_old(i,j) == BURNING )
        forest(i,j) = BURNT;
      end
    end
  end

  return
end
function forest_display ( m, n, forest, hour, i_ignite, j_ignite )

%*****************************************************************************80
%
%% forest_display() displays the current state of the forest fire.
%
%  Discussion:
%
%    An MxN forest of trees is assumed.  FOREST(I,J) reports 
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 November 2017
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns.
%
%    integer FOREST(M,N), the status
%    of the tree at location (I,J), which is
%    0 UNBURNT;
%    1 BURNING;
%    2 BURNT;
%
%    integer HOUR, the current hour.
%
%    integer I_IGNITE, J_IGNITE, the index of the tree that
%    first caught fire.
%
  UNBURNT = 0;
  BURNING = 1;
  BURNT = 2;
%
%  Clear the graphics frame.
%
  clf ( );
%
%  Determine the plot range.
%
  margin = 0.05;

  x_axes_min = 1.0 - 0.5 - margin;
  x_axes_max = n   + 0.5 + margin;
  y_axes_min = 1.0 - 0.5 - margin;
  y_axes_max = m   + 0.5 + margin;
%
%  Fill in the background with black.
%
  x1 = x_axes_min;
  x2 = x_axes_max;
  y1 = y_axes_min;
  y2 = y_axes_max;

  rgb = [ 0.0, 0.0, 0.0 ];

  fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
  hold ( 'on' );
%
%  Draw a square, representing the tree,
%  with most of the length and width, centered at (I,J).
%  WHITE = fire start,
%  GREEN = unburnt
%  RED = burning
%  GRAY = burnt.
%
  for i = 1 : m
    for j = 1 : n

      y1 = m + 1 - ( i - 0.47 );
      y2 = m + 1 - ( i + 0.47 );
      x1 = j - 0.47;
      x2 = j + 0.47;

      if ( forest(i,j) == UNBURNT )
        rgb = [ 0.0, 1.0, 0.0 ];
      elseif ( forest(i,j) == BURNING )
        rgb = [ 1.0, 0.0, 0.0 ];
      elseif ( forest(i,j) == BURNT )
        rgb = [ 0.85, 0.85, 0.85 ];
      end

      fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );

      if ( i == i_ignite && j == j_ignite )
        xc = ( x1 + x2 ) / 2.0;
        yc = ( y1 + y2 ) / 2.0;
        plot ( [x1,x2], [y2,y1], 'r-', 'LineWidth', 3 );
        plot ( [x1,x2], [y1,y2], 'r-', 'LineWidth', 3 );
      end

    end
  end
%
%  Make a title.
%
  title_string = sprintf ( 'Forest fire at hour = %d', hour );

  title ( title_string, 'FontSize', 24 );
%
%  Choose the aspect ratio and other plot details.
%
  axis ( [ x_axes_min, x_axes_max, y_axes_min, y_axes_max] );
  axis ( 'equal' );
  axis ( 'tight' );

  hold ( 'off' );

  pause ( 1 );

  return
end
function forest_print ( m, n, forest, i_ignite, j_ignite )

%*****************************************************************************80
%
%% forest_print() prints the state of the trees in the forest.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 January 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns of trees.
%
%    integer FOREST(M,N), an array
%    with an entry for each tree in the forest.
%
%    integer I_IGNITE, J_IGNITE, the location of the start 
%    of the fire.
%
  BURNT = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map of fire damage.\n' );
  fprintf ( 1, '  Fire started at "*".\n' );
  fprintf ( 1, '  Burned trees are indicated by ".".\n' );
  fprintf ( 1, '  Unburned trees are indicated by "X".\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m
    fprintf ( 1, '  ' );
    for j = 1 : n
      if ( i == i_ignite && j == j_ignite )
        fprintf ( 1, '*' );
      elseif ( forest(i,j) == BURNT )
        fprintf ( 1, '.' );
      else
        fprintf ( 1, 'X' );
      end
    end
    fprintf ( 1, '\n' );
  end

  return
end

