function nbodyplot(t, r, plotoptions, markersize)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function nbodyplot(t, r, plotoptions)
%
% Plots particle trajectories from N-body calculations.
% 
% Assumes motion is restricted to x-y plane.
% 
% INPUT ARGUMENTS
%
%   t (real 1D array, length nt): Contains output times.  At the current
%     time, nt must be >= 2 (see note below).
%
%   r (real 3D array, dimensions np x 3 x nt): (x,y,z) coordinates of
%     all particles at all time steps. np is the number of particles. 
%     nt is the number of time steps stored. 
%     
%     (r(i,1,n), r(i,2,n), r(i,3,n)) are the x, y and z coordinates, 
%     respectively, of the i-th particle at the n-th time step.
% 
%     NOTE: In MATLAB/octave it is apparently impossible to define 
%     a 3D array with a third dimension of length 1: for example,
%     although zeros will happily appear to create one if you type 
%     'zeros(2,2,1)', the result will actually be a 2D array with 
%     size 2 x 2.  Although nt = 1 could be implemented as a special 
%     case here, I have not yet done so.  Thus the restriction to 
%     nt >= 2.
%
%   plotoptions (OPTIONAL, string column vector, arbitrary length): 
%     Contains sequence of plot options for successive curves generated 
%     by function (see 'help plot' for more details). 
% 
%     Example: (reproduces default)
% 
%        plotoptions = ...
%           [ 'ok';  'ob';  'or';  'og';  'om';  'oc' ];
%
%     Important: plotoptions MUST be defined as a column vector
%     since attempted definition as row vector, e.g.
% 
%        plotoptions = ...
%           [ 'ok',  'ob',  'or',  'og',  'om',  'oc' ];
% 
%     results in concantenation, so the above is equivalent to
%
%        plotoptions = [ 'okoborogomoc' ];
% 
%     If plotoptions is the empty vector, then the default
%     options are used
% 
%   markersize (OPTIONAL, integer scalar): marker size for plots
%
% OUTPUT ARGUMENTS
%
%   None
% 
% EXAMPLES
% 
%   % Plot using default colors and style
%   nbodyplot(t, r);
%
%   % Plot using red, green, blue lines cyclically for particles 
%   % 1, 2, 3, ...
%   nbodyplot(t, r, ['r', 'g', b']);
%
%   % Plot using default colors / style with markersize 10
%   nbodyplot(t, r, [], 10);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % Function name for diagnostic messages ...
   S = 'nbodyplot';

   % Default plot option array
   plotoptions_default = ...
      [ 'ok';  'ob';  'or';  'og';  'om';  'oc' ];

   % Default markersize
   markersize_default = 2;

   % BEGIN argument processing / checking ...
   err = 0;
   np = 0;
   nt = 0;

   % Need at least two arguments ...  
   if nargin < 2
      fprintf('%s: Must be called with at least 2 arguments.\n', S);
      fprintf('%s: Type ''help %s'' for usage information.\n', S, S);
      err = 1;
      return;
   end

   % Check that t is a vector; if it is, set nt to its length,
   % and ensure that nt >= 2 ...
   if ~ isvector(t) 
      fprintf('%s: First argument [t] must be a vector.\n', S);
      err = 1;
   else
      nt = length(t);
      if nt < 2 
         fprintf('%s: Number of time steps [length(t)=%d] must be >= 2.\n', S, nt);
         err = 1;
      end
   end

   % Check dimensions of r ...
   szr = size(r);
   szszr = size(size(r));
   rankr = szszr(2);
   if rankr ~= 3
      fprintf('%s: Second argument [r] must be a 3D array.\n', S);
      err = 1;
   else
      if szr(2) ~= 3
         fprintf('%s: Second dimension of second argument [r] must be 3.\n', S);
         err = 1;
      else
         % Determine number of particles from dimensions of r and 
         % ensure that length of 3rd dimension is nt ...
         np = szr(1);
         if szr(3) ~= nt 
            fprintf('%s: Third dimension of second argument [r] ', S);
            fprintf('should be %d but is %d.\n', nt, szr(3));
            err = 1;
         end
      end
   end

   if nargin < 3
      % Use default plotoptions not specified in arg list ...
      plotoptions = plotoptions_default;
   elseif length(plotoptions) == 0
      % ... or if supplied arg is empty
      plotoptions = plotoptions_default;
   else
      % Check that plotoptions is an array of strings
      for i = 1 : length(plotoptions)
         if ~ ischar(plotoptions(i))
            fprintf('%s: Third argument [plotoptions] must be an array of strings.\n', S);
            break;
            err = 1;
         end
      end
   end
   npo = length(plotoptions);

   % Use default markersize if not specifed in arg list
   if nargin < 4
      markersize = markersize_default;
   else
      markersize = round(markersize);
   end

   if err 
      fprintf('\n%s: One or more errors detected in input arguments.\n', S);
      fprintf('%s: Aborting without generating plots.\n', S);
      fprintf('%s: Type ''help %s'' for usage information.\n', S, S);
      return
   end

   clf;
   for ip = 1  : np
      % Clear figure, and enable multiplot mode ...
      hold on;
      % Extract x and y coordinates, ensuring that slices are converted to vectors
      x = r(ip,1,:);
      x = x(:);
      y = r(ip,2,:);
      y = y(:);
      % ... and generate plot ... 
      ipo = mod(ip,npo);
      if ipo == 0
         ipo = npo;
      end
      plotoptions(ipo,:);
      plot(x, y, plotoptions(ipo,:), 'markersize', markersize);
   end
end
