function theHandle = subplotsafe(m, n, ind, jnd)

% subplotsafe -- Subplot via Matlab indexing.
%  subplotsafe(m, n, ind, jnd) returns the handle
%   associated with the subplot that spans
%   the (ind, jnd) indices in a template
%   of size [m n] overall.  If such a subplot
%   already exists, having been created with
%   this routine, it is made current.
%   Note: we use regular Matlab array
%   indexing here, rather than the counter-
%   intuitive scheme of the "subplot" function.  
%  subplotsafe (no arguments) demonstrates itself.
 
% Copyright (C) 2001 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 18-Oct-2001 11:11:04.
% Updated    18-Oct-2001 16:22:05.

if nargin < 1
	help(mfilename)
	m = 5;
	n = 2;
	ind = 4:5;
	jnd = 1:2;
	feval(mfilename, m, n, ind, jnd)
	s = [' ## ' mfilename '(' int2str(m) ', ' int2str(n) ', ' ...
			mat2str(ind) ', ' mat2str(jnd) ')'];
	disp(s)
	return
end

if nargin < 4, jnd = 1; end

x = zeros(n, m);
x(:) = 1:prod(size(x));
x = x';

% Create invisible figure of subplots.

fig = figure('Visible', 'off', 'Units', 'normalized');

h = [];   % Handles.
k = 0;    % Counter.
kmax = prod([size(ind) size(jnd)]);
pos = zeros(kmax, 4);

for j = jnd
	for i = ind
		k = k+1;
		h(k) = subplot(m, n, x(i, j));
		pos(k, :) = get(h(k), 'Position');
	end
end
	
delete(fig)

% Measure the extent of the request.

pos(:, 3:4) = pos(:, 3:4) + pos(:, 1:2);
pos(:, 1:2) = -pos(:, 1:2);
pos = max(pos);
pos(:, 1:2) = -pos(:, 1:2);
pos(:, 3:4) = pos(:, 3:4) - pos(:, 1:2);

% Activate the existing subplot or create a new one.
%  Careful: the Matlab "plot" command erases the
%  "Tag" property of the "gca", so we cannot search
%  for it here.

oldUnits = get(gcf, 'Units');   % Save units.
set(gcf, 'Units', 'normalized')   % Set to normalized units.

ax = findobj(gcf, 'Type', 'axes', 'Position', pos);

if any(ax)
	h = ax(1);
else
	h = axes('Position', pos, 'Tag', mfilename);   % Create.
end

axes(h)

set(gcf, 'Units', oldUnits)   % Restore units.

if nargout > 0
	theHandle = h;
end
