function handle = topomap(depth,axes,color,maxlat)
% TOPOMAP makes map of the world.
%
%        TOPOMAP(DEPTH,REGION,[COLOR]) draws a topographic profile of the chosen 
%        region at 1/6 degree resolution. REGION is [LON_MIN LON_MAX LAT_MIN
%        LAT_MAX].  COLOR is optional and must be one of the Matlab standard color
%        arguments, i.e. 'r' 'g' 'b' 'y' 'm' 'c' 'w' or 'k'. The default color 
%        is white.
%
%        TOPOMAP (called without arguments) draws the coastline of 
%        the whole world in white with the Pacific in the center.
%
%        TOPOMAP(PROJECTION), where projection is 'rec' (Mercator) 'north' (North
%        polar or 'south' (South Polar), draws the coastline in the projection 
%        specified.  TOPOMAP('rec') is identical to TOPOMAP with no arguments.
% 
%        TOPOMAP(DEPTH,PROJECTION,COLOR,[MAXLAT]), where projection is 'north' or
%        'south', draws a polar projection of the specified isobath or topographic
%        contour from the pole down to latitude MAXLAT.  MAXLAT is optional and is
%        taken to lie in the same hemisphere as the pole).  
%        
%        TOPOMAP(FILENAME,REGION,COLOR) behaves as above but loads a user-defined
%        bathymetry file.  Bathymetry files are .mat files with a single variable,
%        'dlines', which is a 2xN matrix of geographic points with longitude on
%        the first line and latitude on the second.  Contours are separated by
%        NaNs, e.g., 
%
%          [85.8320   85.8333   85.8347   85.8333       NaN   75.3333   75.2908 ...
%           28.5000   28.4982   28.5000   28.5014       NaN   36.3857   36.3333 ...]
%
% TOPOMAP in the default uses a 10-minute resolution bathymetry/topography 
%         available in 1000m intervals from -10,000m to +7000m.  This
%         bathymerty is based on etopo5.
%
%        Examples:   topomap(-1000,[-70 30 -50 0],'r')
%                    % maps the 1000m isobath in the South Atlantic region
%                    
%                    topomap(0,[-334 25 -90 90],'y')
%                    % maps the world coast in a rectlinear projection
%                    with high resolution.  Map is centered in the Pacific
%
%                    topomap(3000,'north','b',50)
%                    % maps the 3000m topography in the north polar 
%                    % region from 50N to 90N in blue.
%                    
%                    topomap(-5000)
%                    % maps the 5000m isobath in white, high resolution
%
%                    topomap('south')
%                    % maps the coast of Antarctic from 40S to 90S in white
%                    % (low resolution)
%
%                    topomap('myfilename','south','y')
%                    % maps user-generated file in south polar projection

%Notes: RP (WHOI) 6/Dec/91
%                 7/Nov/92  Changed for matlab4.0
%                 17/Oct/93 New outline file.
%                 14/Mar/94 Fixed polar projections, and changed name
%                           to 'worldmap'.
%
%  D Byrne 7 Dec 1995 Added in subsampled Etopo5 dataset.
%                     Fixed axis labelling.
%

% check inputs
if nargin == 0
  load sealevel
elseif ((nargin == 1) & isstr(depth))
  load sealevel
else
  if isstr(depth)
    arg = sprintf(['load ' depth]);
    eval(arg)
  else
    if depth < 0
      arg = sprintf('load b%4.4i',abs(depth));
      eval(arg)
    elseif depth == 0
      eval(['load sealevel'])
    else
      arg = sprintf('load t%4.4i',depth);
      eval(arg)
    end
  end
end

if nargin < 3
  color = 'w';
end

shiftx=0;
drawtwice=0;
proj='rec';

if (nargin < 4), maxlat = 40; else maxlat = abs(maxlat); end;
if (maxlat > 90); error('MAXLAT greater than 90!'); end;
if (nargin > 1),
  if (isstr(axes)),
    if (axes(1:3)=='nor'),
      proj='npl';
      axes=(90-maxlat)*[-1 1 -1 1];
    elseif (axes(1:3)=='sou'),
      proj='spl';
      axes=(90-maxlat)*[-1 1 -1 1];
    elseif (axes(1:3) == 'rec')
      axes=[-334 25 -90 90];
      shiftx=-360;
      drawtwice=1;
    else
      error('map: Unrecognized projection!');
    end;
  else         
    axes=axes(:)';
    if (max(size(axes)) ~= 4), 
      error('map: wrong number of limit args!');
    end;
    if ( (axes(4)-axes(3))>180. ),
      error('map: Lat range greater than 180 degrees');
    end;
    if ( (axes(2)-axes(1))>360. ),
      error('map: Long range greater than 360 degrees');
    else
      if (axes(1)<=-180. ),
	shiftx=-360;
	drawtwice=1;
      elseif (axes(2)>180.), 
	shiftx=360;
	drawtwice=1;
      else
	shiftx=0;
	drawtwice=0;
      end;
    end;
  end;
elseif nargin == 1
  if isstr(depth),
    if (depth(1:3)=='nor'),
      proj='npl';
      axes=(90-maxlat)*[-1 1 -1 1];
    elseif (depth(1:3)=='sou'),
      proj='spl';
      axes=(90-maxlat)*[-1 1 -1 1];
    elseif (axes(1:3) == 'rec')
      axes=[-334 25 -90 90];
      shiftx=-360;
      drawtwice=1;
    else
      error('map: Unrecognized projection!');
    end;
  else
    axes=[-334 25 -90 90];
    shiftx=-360;
    drawtwice=1;
  end
else % nargin is 0: default plots Mercator coast
  axes=[-334 25 -90 90];
  shiftx=-360;
  drawtwice=1;
end;

if (proj=='rec'),
  h=plot(dlines(1,:)+shiftx,dlines(2,:),'Color',color,'Erasemode','none');
elseif (proj=='npl'),
  kk=find(dlines(2,:)<maxlat);
  xx=(90-dlines(2,:)).*cos(dlines(1,:)*pi/180-pi/2);
  yy=(90-dlines(2,:)).*sin(dlines(1,:)*pi/180-pi/2);
  xx(kk)=NaN*kk;
  h=plot(xx,yy,'Color',color,'Erasemode','none');
  set(gca,'aspect',[1 1]);
elseif (proj=='spl'),
  kk=find(dlines(2,:)>-maxlat);
  xx=(90+dlines(2,:)).*cos(-dlines(1,:)*pi/180+pi/2);
  yy=(90+dlines(2,:)).*sin(-dlines(1,:)*pi/180+pi/2);
  xx(kk)=NaN*kk;
  h = plot(xx,yy,'Color',color,'Erasemode','none');
  set(gca,'aspect',[1 1]);
end;

if (drawtwice),
  h2 = line(dlines(1,:),dlines(2,:),'Linestyle','-',...
      'Color',color,'Erasemode','none');
  h = [h(:); h2(:)];
end;

if (proj=='rec'),
  axis(axes); % new method for xlabels which is robust
  xl = get(gca,'xtick'); [Nn,Mm] = size(xl);
  zz=rem(xl+540,360)-180; zz = zz(:);
  nxl = []; 
  for kk=1:Mm,
    ixl = num2str(abs(zz(kk)));
    s = size(ixl); s = s(2);
    for ii=1:3-s
      ixl = [' ' ixl];
    end
    if zz(kk) < 0, 
      nxl = str2mat(nxl,[ixl 'W']);
    elseif zz(kk) > 0  
      nxl = str2mat(nxl,[ixl 'E']);
    else
      nxl = str2mat(nxl,[ixl ' ']);
    end;
  end;
  nxl = nxl(2:Mm,:);
%  set(gca,'Xticklabels',nxl)
  yl=get(gca,'Ytick');
 [Mm,Nn]=size(yl);
 nyl = [];
 for kk=1:Nn,
    iyl = num2str(abs(yl(kk)));
    s = size(iyl); s = s(2);
    for ii=1:3-s
      iyl = [' ' iyl];
    end
   if (yl(kk) < 0)
     nyl = str2mat(nyl,[iyl 'S']);
   elseif (yl(kk) > 0)
     nyl = str2mat(nyl,[iyl 'N']);
   else
     nyl = str2mat(nyl,[iyl ' ']);
   end;
 end;
 nyl = nyl(2:Nn,:);
% set(gca,'Yticklabels',nyl)
 set(gca,'Ytick',yl);  % Need this to reset position limits

else % polar projection
  axis(axes);
  set(gca,'visible','off');
  for kk=10:10:(90-maxlat),
    xx=kk*cos([0:10:360]*pi/180);
    yy=kk*sin([0:10:360]*pi/180);
    line(xx,yy,'color','w','linestyle',':');
  end;
  xx=(90-maxlat)*cos([0:5:360]*pi/180);
  yy=(90-maxlat)*sin([0:5:360]*pi/180);
  line(xx,yy,'color','w','linestyle','-');
  for kk=10:10:(90-maxlat-2),
    if (proj=='npl'), 
      text(0,kk,sprintf('%2.0fN',90-kk),'horizontal','center');
    else
      text(0,kk,sprintf('%2.0fS',90-kk),'horizontal','center');
    end;
  end;
  for kk=0:30:359,
    xx=[10 90-maxlat]*cos(kk*pi/180);
    yy=[10 90-maxlat]*sin(kk*pi/180);
    line(xx,yy,'color','w','linestyle',':');
    if (kk<=180),nyl=sprintf('%3.0fE',kk);
    else nyl=sprintf('%3.0fW',abs(360-kk)); 
    end;
    if (proj=='npl'), 
      text(yy(2),-xx(2),nyl,'rotation',kk,'horizontal','center','vertical','top');
    else
      text(yy(2),xx(2),nyl,'rotation',-kk,'horizontal','center','vertical','bottom');
    end;
  end;
end;

if nargout == 1
  handle = h;
end
clear ans arg axes color depth dlines drawtwice 
clear kk lh maxlat nyl proj shiftx str xx yy
return
