function shade_nan(x, y, z, shade_color)
 
% SHADE_NAN: replaces NaNs in a 2-d diagram with coloured patches.
% function shade_nan(x, y, z, shade_color)
%  INPUT:
% x, y & z describe position in a way suitable for a call to
% contour(x, y, z), i.e., x and y may be vectors or a matrix and z must be a
% matrix of compatible size.
% shade_color is a 3 element vector giving the rgb values of the colour to be
% used to patch the positions of the NaNs.

% $Id: shade_nan.m,v 1.2 1998/12/11 07:15:30 mansbrid Exp $
% Copyright J. V. Mansbridge, CSIRO, Fri Dec 11 16:07:43 EST 1998
  
% Check the input arguments. At this stage determine whether x and y are
% vectors or not.

if nargin ~= 4
  error('There should be 4 input arguments')
end

if (min(size(x)) == 1)
  x_is_vector = 1;
else
  x_is_vector = 0;
end

num_rows = size(z, 1);
num_cols = size(z, 2);
if x_is_vector
  x = x(:);
  y = y(:);
  if (num_rows ~= length(y))
    error('the number of rows in z should be the length of y')
  end
  if (num_cols ~= length(x))
    error('the number of columns in z should be the length of x')
  end
end

ff_nan = find(isnan(z));
if (~isempty(ff_nan))
  
  % Identify the regions to be covered by patches. When there are NaNs at
  % interior points the patches will go all the way to the neighbouring
  % points. When there are NaNs on the boundaries it is necessary to stop the
  % patches from going outside the boundaries.

  row = mod(ff_nan, num_rows);
  col = floor(ff_nan/num_rows) + 1;
  col_m = col - 1;
  ff = find(col_m < 1);
  if (~isempty(ff))
    col_m(ff) = repmat(1, size(ff));
  end;
  col_p = col + 1;
  ff = find(col_p > num_cols);
  if (~isempty(ff))
    col_p(ff) = repmat(num_cols, size(ff));
  end;
  row_m = row - 1;
  ff = find(row_m < 1);
  if (~isempty(ff))
    row_m(ff) = repmat(1, size(ff));
  end;
  row_p = row + 1;
  ff = find(row_p > num_rows);
  if (~isempty(ff))
    row_p(ff) = repmat(num_rows, size(ff));
  end;

  % Specify the areas to be patched. This is done differently for vectors and
  % matrices.
  
  if x_is_vector
    x_bound = [x(col_m) x(col_m) x(col_p) x(col_p)]';
    y_bound = [y(row_m) y(row_p) y(row_p) y(row_m)]';
  else
    ff_bl = row_m + num_rows*(col_m - 1);
    ff_tl = row_p + num_rows*(col_m - 1);
    ff_br = row_m + num_rows*(col_p - 1);
    ff_tr = row_p + num_rows*(col_p - 1);
    x_bound = [x(ff_bl) x(ff_tl) x(ff_tr) x(ff_br)]';
    y_bound = [y(ff_bl) y(ff_tl) y(ff_tr) y(ff_br)]';
  end
  
  % Draw the patches.
  
  patch(x_bound, y_bound, shade_color, 'EdgeColor', shade_color);
end
