function flat_pcolor(x, y, c)
% function flat_pcolor(x, y, c)
% where x and y are vectors, makes a pseudocolor plot of the matrix c
% on the grid defined by x and y.  This uses pcolor to produce a
% simple color image in the way that one
% would intuitively expect pcolor to work.  The problem with pcolor is
% that if shading is 'flat' or 'faceted' then the value of c(i,j) is
% used to draw a rectangle where the bottom left hand corner is at
% (x(j), y(i)) and the top right hand corner is at (x(j+1), y(i+1)).
% Accordingly, the last row and column of c are not used.
% flat_pcolor draws uses (x(j), y(i)) as the center of the rectangle
% and all rows and columns of c are used.
% Note that if shading is 'interp' then pcolor does a slow, bilinear
% interpolation between data points as would be expected but this
% means that a single NaN fills creates 4 blank spaces.
% Note also that if x or y is a scalar then the axis that it is used
% for will go from 0 to 1.  This is useful for drawing a legend.

%     Copyright J. V. Mansbridge, CSIRO, june 28 1993

[m, n] = size(c);

% Create interpolated and extrapolated vectors or matrices in the x
% and y directions.

if nargin == 1
  
  xs = 1:(n + 1);
  ys = 1:(m + 1);
  
elseif nargin == 3
  
  [mx, nx] = size(x);
  if nx == 1
    if mx == 1    % scalar case
      xs = [ 0 1 ];
    else          % column vector case
      x0 = 2*x(1) - x(2);
      xf = 2*x(mx) - x(mx-1);
      xs1 = [ x; xf ];
      xs = [ x0; x ];
      xs = (xs + xs1)/2;
    end
  else
    if mx == 1    % row vector case
      x0 = 2*x(1) - x(2);
      xf = 2*x(nx) - x(nx-1);
      xs1 = [ x xf ];
      xs = [ x0 x ];
      xs = (xs + xs1)/2;
    else          % matrix case
      xs = [ x   , x(:, nx);
	     x(mx, :), 
    end
  end

% Create an interpolated and extrapolated vector in the y direction.

[m, n] = size(y);
if n == 1
  if m == 1
    ys = [ 0 1 ];
  else
    y0 = 2*y(1) - y(2);
    yf = 2*y(m) - y(m-1);
    ys1 = [ y; yf ];
    ys = [ y0; y ];
    ys = (ys + ys1)/2;
  end
else
  y0 = 2*y(1) - y(2);
  yf = 2*y(n) - y(n-1);
  ys1 = [ y yf ];
  ys = [ y0 y ];
  ys = (ys + ys1)/2;
end

end

% Create an extended version of the c matrix with an extra row and
% column at the ends.

[m, n] = size(c);
cs = zeros(m+1, n+1);
cs(1:m, 1:n) = c;
cs(m+1,:) = cs(m,:);
cs(:, n+1) = cs(:, n);

% Call pcolor to produce the picture.

pcolor(xs, ys, cs);
shading('flat')
