function ising_2d_simulation ( m, n, iterations, thresh )

%*****************************************************************************80
%
%% ising_2d_simulation() carries out a 2D Ising simulation.
%
%  Discussion:
%
%    Note that, when all the cells are updated in a single cycle, there is
%    a mathematically stable checkerboard solution, in which all the reds
%    and blacks flip color repeatedly.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    21 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N: the dimensions of the region.
%
%    integer ITERATIONS: the number of iterations.
%
%    real THRESH: the threshold value.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ising_2d_simulation():\n' );
  fprintf ( 1, '  Monte Carlo simulation of a 2D Ising model.\n' );
%
%  Arguments not supplied are set to default values.
%
  if ( nargin < 1 )
    m = 10;
  elseif ( ischar ( m ) )
    m = str2num ( m );
  end

  if ( nargin < 2 )
    n = 10;
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  if ( nargin < 3 )
    iterations = 15;
  elseif ( ischar ( iterations ) )
    iterations = str2num ( iterations );
  end

  if ( nargin < 4 )
    thresh = 0.50;
  elseif ( ischar ( thresh ) )
    thresh = str2num ( thresh );
  end
%
%  Define the probability of flipping if you are in a neighborhood of
%  1, 2, 3, 4, or 5 of the same sign.
%
%  This should really be some kind of exponential dependent on "temperature".
%
  prob = [ 0.98, 0.85, 0.50, 0.15, 0.02 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The row dimension M = %d\n', m );
  fprintf ( 1, '  The column dimension N = %d\n', n );
  fprintf ( 1, '  The number of iterations taken is ITERATIONS = %d\n', iterations );
  fprintf ( 1, '  The threshhold THRESH = %f\n', thresh );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The transition probability table, based on the number of\n' );
  fprintf ( 1, '  neighbors with the same charge.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        1           2           3           4           5\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f\n', prob(1:5) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  An automatic pause of 2 seconds is imposed between \n;' );
  fprintf ( 1, '  successive displays.\n' );
%
%  Initialize C1.
%
  c1 = ising_2d_initialize ( m, n, thresh );
%
%  Do the simulation.
%
  c1 = transition ( m, n, iterations, prob, c1 );

  return
end
function c5 = ising_2d_agree ( m, n, c1 )

%*****************************************************************************80
%
%% ising_2d_agree() returns the number of neighbors agreeing with each cell.
%
%  Discussion:
%
%    The count includes the cell itself, so it is between 1 and 5.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of cells in each spatial dimension.
%
%    integer C1(M,N), an array of 1's and -1's.
%
%  Output:
%
%    integer C5(M,N), the number of neighbors that agree.
%    1, 2, 3, 4, or 5.
%
  c5 = c1 ...
     + circshift ( c1, [ -1,  0 ] ) ...
     + circshift ( c1, [ +1,  0 ] ) ...
     + circshift ( c1, [  0, -1 ] ) ...
     + circshift ( c1, [  0, +1 ] );

  i = find ( 0 < c1 );
  c5(i) = ( 5 + c5(i) ) / 2;

  i = find ( c1 < 0 );
  c5(i) = ( 5 - c5(i) ) / 2;

  return
end
function ising_2d_display ( step, m, n, c1 )

%*****************************************************************************80
%
%% ising_2d_display() displays the current Ising status.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer STEP, the step number.
%
%    integer M, N, the number of rows and columnss.
%
%    integer C1(M,N), the status of each cell:
%    -1, display as a shade of red.
%    +1, display as a shade of blue.
%
  figure ( 1 );
%
%  Clear the graphics frame.
%
  clf ( );
%
%  Determine the plot range.
%
  margin = 0.05;

  x_axes_min = 1.0 - 0.5 - margin;
  x_axes_max = m   + 0.5 + margin;
  y_axes_min = 1.0 - 0.5 - margin;
  y_axes_max = n   + 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.5, 0.5, 0.5 ];

  fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
  hold ( 'on' );
%
%  Draw a square, representing the bed,
%  with most of the length and width, centered at (I,J).
%
  for i = 1 : m
    for j = 1 : n

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

      if ( c1(i,j) == - 1 )
        rgb = [ 1.0, 0.0, 0.0 ];
      elseif ( c1(i,j) == + 1 )
        rgb = [ 0.0, 0.0, 1.0 ];
      end

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

    end
  end
%
%  Make a title.
%
  title_string = sprintf ( 'Ising charges +1/-1 on step %d', step );

  title ( title_string );
%
%  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' );
  axis ( 'off' );

  hold ( 'off' );

  pause ( 2 );

  return
end
function c1 = ising_2d_initialize ( m, n, thresh )

%*****************************************************************************80
%
%% ising_2d_initialize() initializes the Ising array.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of rows and columns.
%
%    real THRESH, the threshhold value, between 0 and 1.
%
%  Output:
%
%    integer C1(M,N), an array of 1's and -1's.
%
% rng ( 'shuffle' );

  c1 = ones ( m, n );

  r = rand ( m, n );
  
  i = find ( r <= thresh );

  c1(i) = -1;

  return
end
function ising_2d_stats ( step, m, n, c1 )

%*****************************************************************************80
%
%% ising_2d_stats() prints information about the current step.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer STEP, the step number.
%
%    integer M, N, the number of rows and columns.
%
%    integer C1(M,N), the current state of the system.
%
  if ( step == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step       Positives       Negatives\n' );
    fprintf ( 1, '             #    %%          #    %%\n' );
    fprintf ( 1, '\n' );
  end

  pos_count = sum ( sum ( 0 < c1 ) );
  neg_count = m * n - pos_count;
  pos_percent = ( 100 * pos_count ) / ( m * n );
  neg_percent = ( 100 * neg_count ) / ( m * n );

  fprintf ( 1, '  %4d  %6d  %6.2f  %6d  %6.2f\n', ...
    step, pos_count, pos_percent, neg_count, neg_percent );

  return
end
function neighbor_2d_display ( step, m, n, c1, c5 )

%*****************************************************************************80
%
%% neighbor_2d_display() displays the current Ising neighbor status.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer STEP, the step number.
%
%    integer M, N, the number of rows and columnss.
%
%    integer C1(M,N), the status of each cell:
%
%    integer C5(M,N), the number of agreeable neighbors.
%
  figure ( 2 );
%
%  Clear the graphics frame.
%
  clf ( );
%
%  Determine the plot range.
%
  margin = 0.05;

  x_axes_min = 1.0 - 0.5 - margin;
  x_axes_max = m   + 0.5 + margin;
  y_axes_min = 1.0 - 0.5 - margin;
  y_axes_max = n   + 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.5, 0.5, 0.5 ];

  fill ( [ x1, x2, x2, x1 ], [ y1, y1, y2, y2 ], rgb );
  hold ( 'on' );
%
%  Draw a square, representing the bed,
%  with most of the length and width, centered at (I,J).
%
  for i = 1 : m
    for j = 1 : n

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

      c = c1(i,j) * c5(i,j);

      if ( c == - 5 )
        rgb = [ 1.0, 0.0, 0.0 ];
      elseif ( c == - 4 )
        rgb = [ 1.0, 0.2, 0.2 ];
      elseif ( c == - 3 )
        rgb = [ 1.0, 0.4, 0.4 ];
      elseif ( c == - 2 )
        rgb = [ 1.0, 0.7, 0.7 ];
      elseif ( c == - 1 )
        rgb = [ 1.0, 0.8, 0.8 ];
      elseif ( c == 0 )
        rgb = [ 1.0, 1.0, 1.0 ];
      elseif ( c == + 1 )
        rgb = [ 0.7, 0.7, 1.0 ];
      elseif ( c == + 2 )
        rgb = [ 0.6, 0.6, 1.0 ];
      elseif ( c == + 3 )
        rgb = [ 0.4, 0.4, 1.0 ];
      elseif ( c == + 4 )
        rgb = [ 0.2, 0.2, 1.0 ];
      elseif ( c == + 5 )
        rgb = [ 0.0, 0.0, 1.0 ];
      end

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

    end
  end
%
%  Make a title.
%
  title_string = sprintf ( ...
    'Ising neighborhoods -5 to +5, step %d to step %d', step-1, step );

  title ( title_string );
%
%  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' );
  axis ( 'off' );

  hold ( 'off' );

  pause ( 2 );

  return
end
function neighbor_2d_stats ( step, m, n, c1, c5 )

%*****************************************************************************80
%
%% neighbor_2d_stats() prints neighbor statistics about the current step.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer STEP, the step number.
%
%    integer M, N, the number of rows and columns.
%
%    integer C1(M,N), the current state of the system.
%
%    integer C5(M,N), the number of agreeable neighbors.
%
  if ( step == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Step     Neighborhood Charge:\n' );
    fprintf ( 1, ...
      '           -5    -4    -3    -2    -1    +1    +2    +3    +4    +5\n' );
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '  %4d', step );

  for n = -5 : 5
    if ( n ~= 0 )
      c = sum ( sum ( c1 .* c5 == n ) );
      fprintf ( 1, '  %4d', c );
    end
  end
  fprintf ( 1, '\n' );
  
  return
end
function c1 = transition ( m, n, iterations, prob, c1 )

%*****************************************************************************80
%
%% transition() carries out a Monte Carlo simulation of a 2D Ising model.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, N, the number of cells in each spatial
%    dimension.
%
%    integer ITERATIONS, the number of iterations to carry out.
%
%    real PROB(1:5).  PROB(I) represents the probability
%    that the spin of a given cell will be reversed, given that it has I immediate
%    neighbors (including itself) with spin the same as its own.
%
%    integer C1(M,N), the current state of the system. 
%
%  Output:
%
%    integer C1(M,N).  On the current state of the system.
%
% rng ( 'shuffle' );

  step = 0;
  ising_2d_stats ( step, m, n, c1 );
  ising_2d_display ( step, m, n, c1 );

  for step = 1 : iterations
%
%  C5 contains 1 through 5, the number of cells that agree with the center cell.
%
    c5 = ising_2d_agree ( m, n, c1 );

    if ( false )
      neighbor_2d_stats ( step, m, n, c1, c5 );
    end

    if ( false )
      neighbor_2d_display ( step, m, n, c1, c5 );
    end
%
%  Determine the chances of flipping cell (I,J).
%
    threshhold = zeros ( m, n );

    for j = 1 : 5
      i = find ( c5 == j );
      threshhold(i) = prob(j);
    end

    r = rand ( m, n );

    i = find ( r < threshhold );

    c1(i) = - c1(i);

    ising_2d_stats ( step, m, n, c1 );
    ising_2d_display ( step, m, n, c1 );

  end

  return
end
 
