function [latgrat,longrat,mat] = read_topo82(latlim,lonlim)
% Purpose: read topo8.2.img global 2-minute (4 km) topography from
% satellite bathymetry. The data is available over the Internet via
% anonymous FTP from <ftp://topex.ucsd.edu/pub/global_topo_2min/>.
%
% Variables:
%   - latlim: vector of limits in latitude: [south:-90 north:90]
%   - lonlim: vector of limits in longitude: [east:0 west:360]
%   - latgart: latitude vector of elevations
%   - longart: longtitude vector of elevations
%   - map: vector of levations
%
% Global variables: DATA_DIR

% Using some global variables
global DATA_DIR

% get lonlim and latlim into the correct range
latlim = latlim(:)';
lonlim = lonlim(:)';

% the size of the matrix
nrows = 6336;
ncols = 10800;

% map lat and long limits to projected coordinates
[rlim,clim] = ll2rc(latlim,lonlim);

% trim latitudes to the limits of the map
rlim = [max([1,min(rlim)]) min([max(rlim),nrows])];

% extract the map matrix
readrows = rlim(1):rlim(2);
readcols = clim(1):clim(2);

% construct a graticule of row and column indices
[cIndGrat,rIndGrat] = meshgrid(readcols,readrows);

% wrap data that extends beyond the right edge of the map
readcols = mod(readcols,ncols); 
readcols(readcols == 0) = ncols;

% Check for requests straddling the edge of the data
indx = find(diff(readcols) < 1);

% read the data
if isempty(indx) % no straddle
	mat = readmtx(fullfile(DATA_DIR,'topo_8.2.img'),nrows,ncols,'integer*2',readrows,readcols,'ieee-be');
else
	mat1 = readmtx(fullfile(DATA_DIR,'topo_8.2.img'),nrows,ncols,'integer*2',readrows,readcols(1:indx(1)),'ieee-be');
	mat2 = readmtx(fullfile(DATA_DIR,'topo_8.2.img'),nrows,ncols,'integer*2',readrows,readcols(indx(1)+1:end),'ieee-be');
	mat = [mat1 mat2];
end

% map row and column graticule to lat and long
[latgrat,longrat] = rc2ll(rIndGrat,cIndGrat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [i,j] = ll2rc(lat,lon)
%LL2RC converts lat-long to row-column for the global satellite bathymetry

% a translated excerpt of Davis Sandwell's copyrighted code in
% ftp://topex.ucsd.edu/pub/global_topo_2min/src_img_sun/img2xyt.f


d2r =.0174533;

ncols=10800;
nrows=6336;
lat1=-72.006;
lon1=0.;
dlon=2/60.;

arg1=log(tan(d2r*(45.+lat1/2)));
arg2=log(tan(d2r*(45.+lat/2)));
i=ceil(nrows+1-(arg2-arg1)/(dlon*d2r));
j=floor((lon-lon1)/dlon+1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [lat,lon] = rc2ll(i,j)
%RC2ll converts lat-long to row-column for the global satellite bathymetry

% a translated excerpt of Davis Sandwell's copyrighted code in
% ftp://topex.ucsd.edu/pub/global_topo_2min/src_img_sun/img2xyt.f


d2r =.0174533;
ncols=10800;
nrows=6336;
lat1=-72.006;
lon1=0.;
dlon=2/60.;

arg1=d2r*dlon*(nrows-i+.5);
arg2=log(tan(d2r*(45.+lat1/2.)));
term=exp(arg1+arg2);
lat=2.*atan(term)/d2r-90.;

lon=lon1+dlon*(j-.5);
