function [map,maplegend] = tbase(scalefactor,latlim,lonlim)
%TBASE  Read TerrainBase Global 5-Min digital terrain data.
%
%  [Z, REFVEC] = TBASE(SCALEFACTOR) reads the data for the entire world,
%  downsampling the data by the scale factor. The result is returned as a
%  regular data grid and referencing vector. Elevations and depths are
%  given in meters above or below mean sea level.
%
%  [Z, REFVEC] = TBASE(SCALEFACTOR, LATLIM, LONLIM) reads the data for the
%  part of the world within the latitude and longitude limits. The limits
%  must be two-element vectors in units of degrees. Longitude limits can be
%  defined in the range [-180 180] or [0 360]. For example, 
%  lonlim = [170 190] returns data centered on the dateline, while 
%  lonlim = [-10  10] returns data centered on the prime meridian.
%
%  For details on locating TBASE data for download over the Internet, see
%  the following documentation at the MathWorks web site: 
%
%  <a href="matlab: 
%  web('http://www.mathworks.com/support/tech-notes/2100/2101.html#terrainbase') 
%  ">http://www.mathworks.com/support/tech-notes/2100/2101.html</a>
%
%  See also GTOPO30, ETOPO, USGSDEM.

%  Copyright 1996-2004 The MathWorks, Inc.
%  $Revision: 1.1.6.3 $ $Date: 2004/12/18 07:46:36 $
%  Written by:  A. Kim, W. Stumpf

%  Binary data file (byteorder - little endian)
%  Data arranged in W-E columns (-180 to 180) by N-S rows (90 to -90).
%  Elevation in meters

if nargin==1
	subset = 0;
elseif nargin==3
	subset = 1;
else
	error('Incorrect number of arguments')
end

sf = scalefactor;
dcell = 5/60;			% 5 minute grid
shift = 0;

if ~subset

%  Check to see if scalefactor fits matrix dimensions
	if mod(1080,sf)~=0 | mod(4320,sf)~=0
		error('Scalefactor must divide evenly into 1080 rows and 4320 columns.')
	end

else

%  Check to see if data needs to be shifted (0 to 2pi)
	if lonlim(2)>180
		shift = 1;
	end

%  Check lat and lon limits
	errnote = 0;
	if latlim(1)>latlim(2)
		warning('1st element of latlim must be greater than 2nd')
		errnote = 1;
	end
	if lonlim(1)>lonlim(2)
		warning('1st element of lonlim must be greater than 2nd')
		errnote = 1;
	end
	if latlim(1)<-90 | latlim(2)>90
		warning('latlim must be between -90 and 90')
		errnote = 1;
	end
	if ~shift & (lonlim(1)<-180 | lonlim(2)>180)
		warning('lonlim must be between -180 and 180')
		errnote = 1;
	end
	if shift & (lonlim(1)<0 | lonlim(2)>360)
		warning('lonlim must be between 0 and 360')
		errnote = 1;
	end
	if errnote
		error('Check limits')
	end

%  Convert lat and lon limits to row and col limits
	if latlim(2)==90
		rowlim(1) = 1;
	else
		rowlim(1) = floor(-12*(latlim(2)-90)) + 1;
	end
	if latlim(1)==-90
		rowlim(2) = 2160;
	else
		rowlim(2) = ceil(-12*(latlim(1)-90));
	end
	if ~shift
		lon0 = -180;
	else
		lon0 = 0;
	end
	if (~shift & lonlim(1)==-180) | (shift & lonlim(1)==0)
		collim(1) = 1;
	else
		collim(1) = floor(12*(lonlim(1)-lon0)) + 1;
	end
	if (~shift & lonlim(2)==180) | (shift & lonlim(2)==360)
		collim(2) = 4320;
	else
		collim(2) = ceil(12*(lonlim(2)-lon0));
	end

end

%  Read TBASE binary image file
fid = fopen('tbase.bin','rb','ieee-le');
if fid==-1
	error('tbase.bin file not found')
end
if ~subset
	firstrow = 0;
	lastrow = 2*sf*floor(2159/sf)*4320;
	colindx = 1:sf:sf*floor(4319/sf)+1;
	maptop = 90;
	mapleft = -180;
else
	firstrow = 2*(rowlim(1)-1)*4320;
	lastrow = firstrow + 2*sf*floor((rowlim(2)-rowlim(1))/sf)*4320;
	colindx = collim(1):sf:collim(1)+sf*floor((collim(2)-collim(1))/sf);
	maptop = 90 - dcell*(rowlim(1)-1);
	mapleft = -180 + dcell*(collim(1)-1);
	if shift
		mapleft = dcell*(collim(1)-1);
	end
end
srow = firstrow:2*sf*4320:lastrow;			% start row position indicators
rows = length(srow);
%  Read from bottom to top of map (first row of matrix is bottom of map)
for m=rows:-1:1
%	if mod(m,10)==0
%		disp(m)
%	end
	fseek(fid,srow(m),'bof');
	temp = fread(fid,[1 4320],'int16');
	if shift
		temp = [temp(2161:4320) temp(1:2160)];
	end
	map(rows-m+1,:) = temp(colindx);
end
fclose(fid);
cellsize = 1/(sf*dcell);
maplegend = [cellsize maptop mapleft];