#! /usr/bin/env python3
#
def forest_fire ( m, n ):

#*****************************************************************************80
#
## forest_fire() simulates a fire in an mxn forest.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 February 2025
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer m, n: the row and column dimensions of the forest.
#
  import numpy as np

  step = 0

  while ( True ):

    if ( step == 0 ):
#
#  Initialize the height of the trees in the forest between 1 and 4.
#  Short trees will burn up faster.
#
      forest = np.random.randint ( low = 1, high = 5, size = [ m, n ] )
#
#  Pick a random tree, location (i,j).  
#  Reset its index to -4, meaning it was tall and is now on fire.
#
      i = np.random.randint ( low = 0, high = m )
      j = np.random.randint ( low = 0, high = n )
      forest[i,j] = -4
      print ( '  The fire has started at tree[',i,',',j,']' )

    else:

      forest = forest_update ( forest )

    forest_plot ( forest, step )

    if ( np.all ( 0 <= forest ) ):
      print ( '  Simulation completed on step ', step )
      break

    step = step + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'forest_fire():' )
  print ( '  Normal end of execution.' )

  return

def forest_plot ( forest, step ):

#*****************************************************************************80
#
## forest_plot() displays the current status of the forest.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 February 2025
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    int forest[m,n]: records the status of each cell of the forest.
#
#    int step: the current time step.
#
  import matplotlib.pyplot as plt

  m, n = forest.shape

  plt.clf ( )
  plt.axis ( 'equal' )
  plt.axis ( 'off' )
#
#  To plot our array so that (0,0) is in the top left, we need to
#  reverse the order of the I coordinates.
#
  for i in range ( 0, m ):
    for j in range ( 0, n ):

      if ( forest[i,j] == 0 ):
        rgb = [0.75,0.75,0.75]      # light gray
      elif ( forest[i,j] == 1 ):
        rgb = [0.8,1.0,0.8]         # light green
      elif ( forest[i,j] == 2 ):
        rgb = [0.7,1.0,0.7]         # medium green
      elif ( forest[i,j] == 3 ):
        rgb = [0.5,1.0,0.5]         # strong green
      elif ( forest[i,j] == 4 ):
        rgb = [0.3,1.0,0.3]         # dark green
      elif ( forest[i,j] == -1 ):
        rgb = [1.0,0.9,0.9]         # light red
      elif ( forest[i,j] == -2 ):
        rgb = [1.0,0.7,0.7]         # medium red
      elif ( forest[i,j] == -3 ):
        rgb = [1.0,0.5,0.5]         # strong red
      elif ( forest[i,j] == -4 ):
        rgb = [1.0,0.3,0.3]         # dark red

      x = [ j, j, j+1, j+1, j ]
      ir = m - 1 - i
      y = [ ir, ir+1, ir+1, ir, ir ]
      plt.fill ( x, y, color = rgb )
      plt.plot ( x, y, 'k-' )

  plt.title ( 'Forest Fire at step ' + str ( step ) )
  filename = 'forest_fire_%04d.png' % ( step ) 
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

def forest_update ( forest ):

#*****************************************************************************80
#
## forest_update() advances the forest fire one step in time.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    19 February 2025
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    int forest[m,n]: records the status of each cell of the forest.
#
#  Output:
#
#    int new_forest[m,n]: the updated status.
#
  import numpy as np

  m, n = forest.shape

  new_forest = forest.copy ( )
  
  for i in range ( 0, m ):
    for j in range ( 0, n ):
#
#  A negative value means the tree is burning.
#  Increase it by 1, so it burns towards 0.
#
      if ( forest[i,j] < 0 ):

        new_forest[i,j] = forest[i,j] + 1
#
#  A positive value means the tree is untouched.
#  Count the burning neighbors to see whether the tree should ignite.
#  This is done by making its index negative.
#
      elif ( 0 < forest[i,j] ):
#
#  Here, we locate the neighbors, but use modular arithmetic
#  so the grid wraps around at the edges.
#  
        im1 = ( i - 1 ) % m
        ip1 = ( i + 1 ) % m
        jm1 = ( j - 1 ) % n
        jp1 = ( j + 1 ) % n
#
#  Count the number of neighboring trees that are burning.
#
        s = ( forest[im1,j] < 0 ) + ( forest[ip1,j] < 0 ) \
          + ( forest[i,jm1] < 0 ) + ( forest[i,jp1] < 0 )

        ignite = np.random.random ( )

        if ( ignite < 0.25 * s ):
          new_forest[i,j] = - new_forest[i,j]
 
  return new_forest

if ( __name__ == "__main__" ):
  forest_fire ( 10, 10 )