function [top_lon, top_lat, topog] = best_topo(lonrange, latrange, varargin)
%  [top_lon, top_lat, topog] =  best_topo(lonrange, latrange, varargin);
% top_lon is [1xNlat], top_lon is [1xNlon], topog is [Nlat x Nlon]
%
% possible databases are 
%
%  data source      program                    resolution     lat range
%   --------     ------------                --------------   ----------
%  etopo5?          m_elev (from m_map package)    1deg         [-90,90]
%  smith-sandwell   ss_block                      .1deg         [-72,72]
%  terrain base     m_tbase (from m_map)           1/12deg      [-90, 90]
% 
%  
%   name             default      what
% nsub                0          subsample topog database: 
%                                   0: best guess (automatically subsample)
%                                integer: specify the amount of decimation
%                                   - decimation done AFTER loading, 
%                                     (except smithsandwell: 3 and 9)
% toponame       'smithsandwell'  choose the database to use.  options:
%                               'smithsandwell'     smith-sandwell
%                               'tbase'             terrain base
%                               'm_elev'            default m_map
%                               ''                 "give me anything"
%


%% JH 2001/12/05
%% should add warnings about how much data will be extracted

topogs = {'smithsandwell', 'm_tbase', 'm_elev'};

cfg.nsub  = 0;
cfg.toponame = 'smithsandwell';

cfg = fillstruct(cfg, varargin, 0);


guess_topo = 0;    %use a specific topo database (default was smith sandwell)
if (isempty(cfg.toponame))
   guess_topo = 1;
end


%  now try to do it:
top_lon = [];
top_lat = [];
topog = [];
got_topog = 0;


if (guess_topo | strcmp(cfg.toponame, 'smithsandwell'))
   biglonrange = sort(lonrange([1 end])) + [-.2 .2];
   biglatrange = sort(latrange([1 end])) + [-.2 .2];

   if (cfg.nsub == 0)  %try to figure it out
      max_llr = max(diff(biglonrange), diff(biglatrange));
      if max_llr  > 15, cfg.nsub = 9;
      elseif max_llr  > 5, cfg.nsub = 3;
      else cfg.nsub  = 1;
      end
   end

   [top_lon, top_lat, topog] = get_smith(biglonrange, biglatrange, cfg);
   if (~isempty(topog))
      got_topog = 1;
   end
end

%put this second in case they're asking for too much data
if (~got_topog & (guess_topo | strcmp(cfg.toponame, 'm_elev')))
   biglonrange = sort(lonrange([1 end])) +[-1 1];
   biglatrange = sort(latrange([1 end])) +[-1 1];

   [top_lon, top_lat, topog] = get_melev(biglonrange, biglatrange, cfg);
   if (~isempty(topog))
      got_topog = 1;
   end
end


if (~got_topog & (guess_topo | strcmp(cfg.toponame, 'tbase')))
   biglonrange = sort(lonrange([1 end])) + [-.2 .2];
   biglatrange = sort(latrange([1 end])) + [-.2 .2];

   [top_lon, top_lat, topog] = get_tbase(biglonrange, biglatrange, cfg);
   if (~isempty(topog))
      got_topog = 1;
   end
end


if (~isempty(topog))
   %lon, lat  : row vectors with the X and Y grids
   if (length(top_lon) > 2) & (length(top_lat) > 2)
      [topog top_lon] = m_fixlon(top_lon, top_lat', topog, lonrange);  
   else
      fprintf('WARNING: fewer than 3 points in one topographic variable\n')
      fprintf('%d points in latitude; %d points in longitude\n',...
              length(top_lat), length(top_lon))
   end
else
   fprintf('returning with no data\n')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function  [top_lon, top_lat, topog] = get_smith(biglonrange, biglatrange, cfg);

if ~exist('ss_block')
   % can't find function
   top_lon = [];
   top_lat = [];
   topog = [];
   fprintf('cannot find ss_block.\n')
   return
end

if (min(abs(biglatrange)) > 70)  
   %both latitudes are out of range
   top_lon = [];
   top_lat = [];
   topog = [];
   fprintf('latitude range is entirely out of ss_block range.\n')
   return
end
  

%trim it if the range is too big
biglatrange(1) = max(-72, biglatrange(1));
biglatrange(2) = min(72,  biglatrange(1));


if (ismember(cfg.samp, [1 3 9]))
   fprintf('getting topography (subsample by %d) from Smith-Sandwell', ...
           cfg.nsub)
   % eval is used to hide ss_block from the compiler--otherwise an
   % error will result even if this block of code is never executed:
   eval(['[top_lon, top_lat, topog]= ',...
         'ss_block(biglonrange, biglatrange, cfg.nsub);']);
else
   fprintf('getting topography (subsample by %d) from Smith-Sandwell', ...
           cfg.nsub)
   % eval is used to hide ss_block from the compiler--otherwise an
   % error will result even if this block of code is never executed:
   eval(['[top_lon, top_lat, topog]= ',...
         'ss_block(biglonrange, biglatrange);']);
   top_lon = top_lon(1:cfg.nsub:end);
   top_lat = top_lat(1:cfg.nsub:end);
   topog = topog(1:cfg.nsub:end, 1:cfg.nsub:end);
end
   
   
% ---------- end of get_smith ---------


function  [top_lon, top_lat, topog] = get_tbase(biglonrange, biglatrange, cfg);
if ~exist('m_tbase')
   top_lon = [];
   top_lat = [];
   topog = [];
   fprintf('cannot find m_tbase.\n')
   return
end

fprintf('getting topog from m_tbase\n')
[topog, top_lon, top_lat]=m_tbase([biglonrange, biglatrange]);

if (cfg.nsub == 0)  %try to figure it out
   max_llr = max(diff(biglonrange), diff(biglatrange));
   if max_llr  > 150, cfg.nsub = 9;
   elseif max_llr  > 50, cfg.nsub = 3;
   else cfg.nsub = 1;
   end
end

top_lon = top_lon(1,1:cfg.nsub:end);
top_lat = top_lat(1:cfg.nsub:end, 1);
topog = topog(1:cfg.nsub:end, 1:cfg.nsub:end);

% ---------- end of get_tbase ---------



function  [top_lon, top_lat, topog] = get_melev(biglonrange, biglatrange, cfg);
if ~exist('m_elev')
   top_lon = [];
   top_lat = [];
   topog = [];
   fprintf('cannot find m_elev.\n')
   return
end

fprintf('getting topog from m_elev\n')
[topog, top_lon, top_lat]=m_elev([biglonrange, biglatrange]);
top_lon = top_lon(1,:);
top_lat = top_lat(:,1);


if (cfg.nsub == 0)  %try to figure it out
   max_llr = max(diff(biglonrange), diff(biglatrange));
   if max_llr  > 15, cfg.nsub = 9;
   elseif max_llr  > 5, cfg.nsub = 3;
   else cfg.nsub = 1;
   end
end

top_lon = top_lon(1:cfg.nsub:end);
top_lat = top_lat(1:cfg.nsub:end);
topog = topog(1:cfg.nsub:end, 1:cfg.nsub:end);


% ---------- end of get_melev ---------

