function []=mbmap(lat1,lat2,lon0);
%function []=mbmap(lat1,lat2,[lon0]);
%    mass bay map
load mb1.map
if nargin==0
    latlim=input('latitude limits (41.7-42.7) ');
else
   latlim(1)=lat1; latlim(2)=lat2;
end
if nargin<3
    lon0=-70.42-(mean(latlim)-42)*.8;
end
dlon=diff(latlim)/.74314*1.74;
lonlim(1)=lon0-dlon/2;
lonlim(2)=lon0+dlon/2;
hold off
axis([lonlim(1) lonlim(2) latlim(1) latlim(2)])
plot(mb1(:,1),mb1(:,2))
hold on

