function [lonrange, latrange] = map_range(x, y, varargin)
% function [lonrange, latrange] = map_range(x, y, ...)
%
%  Calculate a bounding box in which to map a set of points.
%     x, y are vectors of longitudes and latitudes
%     lonrange, latrange are 1x2 vectors with the minimum
%          and maximum longitude and latitude of a suitable box.
%     optional keyword/argument pairs:
%            'aspect',    2      maximum aspect ratio.
%            'stretch',   1.3    factor by which box extends beyond the data
%            'smallest',  0.5    minimum length (degrees lat) of a box side.
%            'lat_limit', 80     max north or south lat of box boundary
%            'round_to',  0      Round boundaries to multiple of this.
%                                   0 means no rounding.
%

%
%  EF 2001/08/09 (based on a simpler version for autovect.m)
%

aspect = 2;
stretch = 1.3;
lat_limit = 80;
smallest = 0.5;
round_to = 0;

for ii = 1:2:length(varargin),
   if strcmp(varargin{ii}, 'aspect'),
      aspect = varargin{ii+1};
   elseif strcmp(varargin{ii},  'stretch'),
      stretch = varargin{ii+1};
   elseif strcmp(varargin{ii},  'smallest'),
      smallest = varargin{ii+1};
   elseif strcmp(varargin{ii},  'lat_limit'),
      lat_limit = varargin{ii+1};
   elseif strcmp(varargin{ii}, 'round_to'),
      round_to = varargin{ii+1};
   else
      help map_range
      error('Unrecognized keyword')
   end
end


lonrange0 = [min(x) max(x)];
latrange0 = [min(y) max(y)];
ymid = mean(latrange0);
xmid = mean(lonrange0);
xsc = cos(ymid * pi/180);
h = diff(latrange0);
w = diff(lonrange0)*xsc;
s = max([h, w, (smallest * aspect / stretch)]);
h = max([h, s/aspect]) * stretch;
w = max([w, s/aspect]) * stretch;

lonrange = xmid + 0.5 * [-w, w] / xsc;
latrange = ymid + 0.5 * [-h, h];

latrange(1) = max([-lat_limit, latrange(1)]);
latrange(2) = min([lat_limit, latrange(2)]);

if round_to > 0,
   latrange = round(latrange / round_to) * round_to;
   lonrange = round(lonrange / round_to) * round_to;
   % Make sure the data are still inside the box.
   if lonrange(1) > lonrange0(1), lonrange(1) = lonrange(1) - round_to; end
   if latrange(1) > latrange0(1), latrange(1) = latrange(1) - round_to; end
   if lonrange(2) < lonrange0(2), lonrange(2) = lonrange(2) + round_to; end
   if latrange(2) < latrange0(2), latrange(2) = latrange(2) + round_to; end
end


