function [c, h] = contour_etopo5(cont);

% Plot isobaths based on the etopo5 bathymetry.
% function [cs, h] = contour_etopo5(cont);
%   cont: a vector of bathymetric levels (use [cont cont] for a single contour)
%   c: contour matrix C as described in CONTOURC
%   h: a column vector H of handles to LINE or PATCH objects, one handle per
%      line. 
% Both C and H can be used as input to CLABEL. The UserData property of each
% object contains the height value for each contour.
%
%   Example usage:
%
% axis([120 123 -19 -16])
% cont = [100 200];
% [c, h] = contour_etopo5(cont);
% set(h, 'LineWidth', 2)
% set(h, 'LineStyle', '--')
% clabel(c)
%
%   Warning - this can be quite slow and etopo5 is often inaccurate

% $Id:$
% Copyright J. V. Mansbridge, CSIRO, Mon Sep 29 18:37:02 EST 1997

if nargin ~= 1
  help contour_etopo5
  error('must be exactly one input argument')
end

% Assume that we want bathymetric contours. Note that we multiply the etopo5
% height by -1 so that we will get positive labels.

if (min(cont) < 0)
  cont = -cont;
end

% Find 'NextPlot' values so that they can be reset later.

nextpl_gcf = get(gcf, 'NextPlot');
ax = gca;
nextpl_gca = get(ax, 'NextPlot');
xlim = get(gca, 'XLim');
ylim = get(gca, 'YLim');
set(gcf,'nextplot','add');
set(ax,'nextplot','add');

bath_name = 'etopo5';

lon_bath_index_min = floor(xlim(1)*12 + 1);
lon_bath_index_max = floor(xlim(2)*12 + 1);
lat_bath_index_min = floor(ylim(1)*12 + 1081);
lat_bath_index_max = floor(ylim(2)*12 + 1081);

lon_bath = getnc(bath_name, 'lon', lon_bath_index_min, ...
    lon_bath_index_max);
lat_bath = getnc(bath_name, 'lat', lat_bath_index_min, ...
    lat_bath_index_max);
height_bath = -getnc(bath_name, 'height', ...
    [lat_bath_index_min lon_bath_index_min], ...
    [lat_bath_index_max lon_bath_index_max]);

if length(cont) > 1
  [cs_t, h_t] = contour(lon_bath, lat_bath, height_bath, cont);
elseif length(cont) == 1
  [cs_t, h_t] = contour(lon_bath, lat_bath, height_bath, [cont cont]);
else
  error('cont has zero length')
end

set(gca, 'XLim', xlim);
set(gca, 'YLim', ylim);
set(gcf, 'NextPlot', nextpl_gcf);
set(ax, 'NextPlot', nextpl_gca);

if nargout == 1
  c = cs_t;
elseif nargout == 2  
  c = cs_t;
  h = h_t;
end
