function [maxamp, maxampbin, b, threshM] = ampscales(amp, varargin);
% [maxamp, maxampbin, b, threshM]  = amp_scales(amp, varargin);
% input amp (AGC) or sound (dB)
% return maxamp and maxampbin%   both are (nbeams X nprofs)
% (also return b and threshM from bigstarget)
%
% NOTE: Acoustic interference from another instrument may cause failure.
% Recommend using clean_amp first, and then using amp_scales.
%
%
% varargin: name, value
%
%  NAME                 DEFAULT            DESCRIPTION
% ---------           ----------         ---------------------
% thresh                [15 10]             size of amplitude bump for a 
%                                         shallow bottom, (like old waterfall
%                                         plots).  
%                      scalar             - constant for all bins, all beams
%                      vector             - same for all beams, linearly 
%                                              interpolate the vector in 
%                                              depth for a vertical profile 
%                                              (first number is shallowest)
%                     structure           - field names are "beam1", "beam2",
%                                             (etc), one per beam.  Each field
%                                             is treated as scalar or vector
%                                             as above.  Use this option
%                                             if one beam is noisy and has 
%                                             low threshold of detection
%
% nside                  1                run a box filter in vertical with 
%                                              this halfwidth 
% offset                 3                number of bins away (in the 
%                                              vertical) to look to define
%                                              the bump (roughly)
% nanmask               (none)            mask of bins to ignore when finding 
%                                              max amp bin

%JH, 2003/10/13

cfg.thresh = [15 10];
cfg.nside = 1;
cfg.offset = 3;
cfg.nanmask = [];

cfg = fillstruct(cfg, varargin);
   

if length(size(amp)) == 3
   [nbins, nbeams, nprofs]=size(amp);
else
   [nbins, nprofs]=size(amp);
   nbeams = 1;
end


%initialize for use or bailout
maxampbin = NaN*zeros(nbeams,nprofs);
maxamp    = NaN*zeros(nbeams,nprofs);
bb  = NaN*amp;


if (0)
%%%%% in case we're dealing with sdb and nans, fill them in
if (length(find(isnan(amp))) > 0)
   if nbeams == 1
      amp = fillends(amp, 'dir', 'col', 'fill',1);
   else
      filledamp = [];
      for beami = 1:nbeams
         eval(sprintf('a = squeeze(amp(:,%d,:));', beami))
         a = fillends(a, 'dir', 'col', 'fill',1);
         filledamp = cat(3, filledamp, a);
      end
      amp = permute(filledamp, [1,3,2]);
   end
end
end


[bmask, b, threshM]=bigstarget(amp,...
                      'nanmask', cfg.nanmask,...
                      'thresh', cfg.thresh,...
                      'nside', cfg.nside',...
                      'offset', cfg.offset);

bmask1 = bottommask(bmask, b, 1, 0);


%% beam1
if nbeams == 1
   goodi = find(sum(bmask1) ~= zeros(1,nprofs));
   maxampbin(goodi) = nbins - sum(bmask1(:,goodi)) + 1;
   
   for jj=1:length(goodi)
      aprof = amp(:,  goodi(jj));
      maxamp(goodi(jj)) = aprof(maxampbin(goodi(jj)));
   end   


else
   for beami = 1:nbeams
      bb1 = squeeze(bmask1(:,beami,:));
      goodi = find(sum(bb1) ~= zeros(1,nprofs));
      maxampbin(beami,goodi) =  nbins - sum(bb1(:,goodi)) + 1;
      %should vectorize this
      for jj=1:length(goodi)
         aprof = squeeze(amp(:, beami, goodi(jj)));
         maxamp(beami,goodi(jj)) = aprof(maxampbin(beami,goodi(jj)));
      end   
   end
end

