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
%                               ''                  (default: try the above
%                                                       list in order
%


%% 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 = '';

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));
      aspect = diff(biglonrange)/diff(biglatrange);
      aspect = max(aspect, 1/aspect);
      if (aspect < 5) % it's squarish
         if max_llr  > 15, cfg.nsub = 9;
         elseif max_llr  > 5, cfg.nsub = 3;
         else cfg.nsub  = 1;
         end
      else
         cfg.nsub = 1;
      end
   end

   %% unnecessary??
   %if (biglonrange(2) < 0 )
   %   biglonrange = biglonrange + 360;
   %end

   [top_lon, top_lat, topog] = get_smith(biglonrange, biglatrange, cfg);
   
   if (~isempty(topog))
      got_topog = 1;
      if (length(top_lon) < 3) | (length(top_lat) < 3)
         fprintf('got some Smith Sandwell topography, but not enough\n')
         got_topog = 0;
      end
   else
      got_topog = 0;
   end
end

%use this if we have it
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


% last resort
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 (~isempty(topog))
   %lon, lat  : row vectors with the X and Y grids
   if (length(top_lon) >= 2) & (length(top_lat) >= 2)
      fprintf('running m_fixlon\n')
      [topog, top_lon] = m_fixlon(top_lon(:)', top_lat(:)',topog, lonrange);  
   else
      fprintf('WARNING: fewer than 2 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(2));


if (ismember(cfg.nsub, [1 3 9]))
   fprintf('getting topography (subsample by %d) from Smith-Sandwell\n', ...
           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\n', ...
           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
   

lonwrap = find(diff(top_lon) < 0);
if (~isempty(lonwrap))
   top_lon(lonwrap+1:end) = top_lon(lonwrap+1:end) + 360;
end

   


   
% ---------- end of get_smith ---------


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

top_lon = [];
top_lat = [];
topog = [];

if ~exist('m_tbase')
   fprintf('cannot find m_tbase.\n')
   return
end

try
   [topog, top_lon, top_lat]=m_tbase([biglonrange, biglatrange]);
   fprintf('getting topog from m_tbase\n')
catch
   fprintf('cannot get topography from tbase\n')
   return
end   

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);


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

