function [c x y] = ca(rule, nsite, ngen, initopt)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implements 1d 2-state (k=2), range-2 (r=2) totalistic cellular 
% automata with "dead" boundary conditions.
%
% Formalism follows Wolfram, Rev. Mod Phys, v55, 601-644 (1983)
%
% Input arguments:
%
%   rule:     Update rule (0-62) as discussed in class: should be 
%             even to maintain quiescence.
%   nsite:    Number of sites (positive integer >= 5).
%   ngen:     Number of generations (positive integer >=2).
%   initopt:  Initialization option
%               == 1  ->  Set sites alive with proability 1/2.
%               ~= 1  ->  Set central site alive, others dead.
%
% Output argument:
%
%   c: (integer array, size ngen x nsite) Complete evolution history of 
%      cellular automata state for ngenerations.
%   x: (integer vector) x coordinates (site number) of all live sites.
%   y: (integer vector) y coordinates (reverse generation number) of all 
%      live sites (generation number is reversed so that when y is 
%      plotted vs x, evolution proceeds from top down).
%
% Rules that generate "interesting" spatiotemporal patterns:
%
%  2, 6, 10, 12, 14, 18, 20 (!), 22, 26, 28, 30, 34, 38, 42, 44, 46,
%  50, 52 (!)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Set output arguments to 'trivial' values for case that one or 
   % more input arguments is invalid.
   c = zeros(0,0);
   x = zeros(0,0);
   y = zeros(0,0);

   % Check input arguments for validity.
   err = 0;
   if fix(rule) ~= rule | rule < 0 | rule > 62 | (rem(rule,2) == 1)
      fprintf('ca: rule must be an even integer between 0 and 32\n');
      err = 1;
   end
   if fix(nsite) ~= nsite | nsite < 5
      fprintf('ca: nsite must be positive integer >= 5\n');
      err = 1;
   end
   if fix(ngen) ~= ngen | ngen < 2
      fprintf('ca: ngen must be positive integer >= 2\n');
      err = 1;
   end
   if fix(initopt) ~= initopt
      fprintf('ca: initopt must be an integer\n');
      err = 1;
   end
   if err 
      return
   end

   c = zeros(0,0);
   % Decode update rule
   update = zeros(1,6);
   lrule = rule;
   for i = 1:6
      update(i) = (rem(lrule,2) == 1);
      lrule = fix(lrule/2);
   end

   % Create state array and define first generation.
   c = zeros(ngen, nsite);
   if initopt == 1
      % Set sites alive with probability 1/2.
      c(1,:) = fix(1 + 2*rand(1,nsite)) - 1;
   else
      % Make central site alive and the rest dead.
      c(1,fix(nsite/2)) = 1;
   end

   % Impose 'dead' boundary conditions.
   c(1,1:2) = 0;
   c(1,nsite-1:nsite) = 0;

   % Perform the evolution ... Note the use of vector operations as 
   % well as octave/MATLAB's ability to index one array with another.
   for igen = 2 : ngen
         c(igen,3:nsite-2) = ...
            update(1 + c(igen-1,1:nsite-4) + c(igen-1,2:nsite-3) +  ...
               c(igen-1,3:nsite-2) + c(igen-1,4:nsite-1) + c(igen-1,5:nsite));
   end

   % Generate x and y coordinates of all the live sites for plotting.

   % First generate x and y arrays of the same size as the state array
   % and such that x defines the site position and y defines the generation
   % number (but reversed in sense so that when plotted, the evolution
   % proceeds from top to bottom).
   x = zeros(ngen, nsite);
   y = zeros(ngen, nsite);
   for i = 1:ngen
      x(:, i) = i;
   end
   for i = 1:ngen
      y(ngen-i+1, :) = i;
   end
   % The following statements convert x and y into 1d arrays that 
   % contain the coordinates of only the live sites.
   x = x(c == 1);
   y = y(c == 1);
end
