function x = cellular_automaton ( cell_num, step_num )

%*****************************************************************************80
%
%% cellular_automaton() computes the evolution of a rule #30 cellular automaton.
%
%  Discussion:
%
%    Given an initial linear array of 0's and 1's, rule 30 produces a new
%    array using the rules:
%
%      111  110  101  100  011  010  001  000
%       V    V    V    V    V    V    V    V
%       0    0    0    1    1    1    1    0     
%
%    Note that there are 256 = 2^8 possible ways to fill in this output
%    chart, and that rule 30 gets its index by the fact that
%    (0,0,0,1,1,1,1,0) can be interpreted as the binary representation of 30.
%
%    For instance, if the current values of X(4), X(5) and X(6) are
%    0, 1 and 1, respectively, then the new value of X(5) will be 1.
%
%    The first and last entries of the array must be treated specially, since
%    they don't have a left or right neighbor.  One simple treatment is 
%    to assume that there are phantom neighbors whose values are both 0.
%    Another is to enforce periodic boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 December 2022
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    A New Kind of Science,
%    Wolfram Media, 2002,
%    ISBN13: 978-1579550080,
%    LC: QA267.5.C45.W67.
%
%  Input:
%
%    integer cell_num: the number of cells.
%
%    integer step_num: the number of steps.
%
%  Output:
%
%    integer x(step_num+1,cell_num): the cellular data.
%

%
%  Compute the entire X array.
%
  x = zeros ( step_num + 1, cell_num  );
%
%  Set the initial condition in row 1.
%
  mid = floor ( ( cell_num + 1 ) / 2 );
  x(1,mid) = 1;
%
%  Compute rows 2 through step_num+1 by applying rule 30.
%
  for i = 2 : step_num + 1

    for j = 2 : cell_num - 1

      if ( ( x(i-1,j-1) == 0 && x(i-1,j) == 0 && x(i-1,j+1) == 1 ) || ...
           ( x(i-1,j-1) == 0 && x(i-1,j) == 1 && x(i-1,j+1) == 0 ) || ...
           ( x(i-1,j-1) == 0 && x(i-1,j) == 1 && x(i-1,j+1) == 1 ) || ...
           ( x(i-1,j-1) == 1 && x(i-1,j) == 0 && x(i-1,j+1) == 0 ) )
        x(i,j) = 1;
      end

    end

  end

  return
end

