function []=bhmap(lat1,lat2,lon0);
%function []=bhmap(lat1,lat2,[lon0]);
%                 note: there are more islands on the original shore file
%                 (bh.map in the bath subdirectory) 
%                  but this one gets the jist of it

load bh1.map
if nargin==0
    latlim=input('latitude limits (42.2-42.5) ');
else
   latlim(1)=lat1; latlim(2)=lat2;
end
if nargin<3
    lon0=-70.93-(mean(latlim)-42.3)*.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(bh(1:gap(1),1),bh(1:gap(1),2),'c1')
hold on
for i=2:length(gap)
     ind=[gap(i) (gap(i-1)+1):gap(i)];
     plot(bh(ind,1),bh(ind,2),'c1')
end
clear ind i

