function jumping_bean_simulation ( bean_num, step_num )

%*****************************************************************************80
%
%% jumping_bean_simulation() simulates the movement of jumping beans.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 July 2023
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of beans = %d\n', bean_num );
  fprintf ( 1, '  Number of steps = %d\n', step_num );
%
%  Initialize positions and temperatures.
%
  bean_x = rand ( bean_num, 1 );
  bean_y = rand ( bean_num, 1 );
  bean_t = 15.0 + 20.0 * rand ( bean_num, 1 );

  clf ( );
  hold ( 'on' );
  for step = 1 : step_num

    bean_x_old = bean_x;
    bean_y_old = bean_y;

    for bean = 1 : bean_num

      x = bean_x(bean);
      y = bean_y(bean);
%
%  Adjust temperature of bean.
%
      g = ground_temperature ( x, y );
      t = bean_t(bean);
      if ( t < g )
        bean_t(bean) = bean_t(bean) + 1;
      elseif ( g < t )
        bean_t(bean) = bean_t(bean) - 1;
      end
      t = bean_t(bean);
%
%  Decide whether to jump.
%
      if ( rand ( ) <= jump_prob ( t ) )
        x = mod ( x + 0.10 * randn ( ), 1.0 );
        y = mod ( y + 0.10 * randn ( ), 1.0 );
        bean_x(bean) = x;
        bean_y(bean) = y;
      end

    end
%
%  Display current position.
%
    plot ( bean_x_old, bean_y_old, '.', 'color', [0.5, 0.5, 0.5 ], 'markersize', 30 );
    plot ( bean_x, bean_y, 'r.', 'markersize', 30 );
    for bean = 1 : bean_num
      plot ( [ bean_x_old(bean), bean_x(bean) ], [ bean_y_old(bean), bean_y(bean) ], ...
        'color', [0.5, 0.5, 0.5 ] );
    end
    plot ( bean_x, bean_y, 'r.', 'markersize', 30 );
    grid ( 'on' );
    axis ( 'square' );
    axis ( [ 0.0, 1.0, 0.0, 1.0 ] );

    pause ( 2 );

  end

  hold ( 'off' );
  
  return
end
function t = ground_temperature ( x, y )

%*****************************************************************************80
%
%% ground_temperature() returns t(x,y), the ground temperature at a point.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 August 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real x, y: the current coordinates of the bean.
%
%  Output:
%
%    real t: the ground temperature at this point.
%
  if ( x < 0.5 || y < 0.5 )
    t = 25.0;
  else
    t = 15.0;
  end

  return
end
function p = jump_prob ( t )

%*****************************************************************************80
%
%% jump_prob() returns the jumping probability as a function of temperature.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 August 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real t: the local ground temperature.
%
%  Output:
%
%    real p: the probability that the bean will jump in the next time interval.
%
  p = ( t - 15.0 ) / ( 30.0 - 15.0 );
  p = min ( p, 1.0 );
  p = max ( p, 0.0 );

  return
end


