function [bmask, b] = bigtarget(a, thresh, offset, nside, mfilt)
% [bmask, b] = bigtarget(a, thresh, offset, nside)
% inputs :
%  a      : amplitude matrix (numdepths x numprofiles)
%  thresh : average of diffs between amp peak and value at offset bins shallower
%         : and deeper than max amp bin (like old offset in waterfall plots)
%  offset : how many bins (shallower and deeper than max amp bin) to determine
%         : amp bump
%  nside  : half-width of the boxcar (nside = 1, 2, or 3 is for smoothing
%         : over 3, 5 or 7 points respectively); nside = 0 is no smoothing
%  mfilt  : run a horizontal median filter of this length, (no halfwidth)
%         :  mfilt is optional; 
%           
%  suggest offset = 5, nside = 2,  thresh about 25, mfilt = 3
%
% outputs:
% bmask   : 1 at max amp bin, 0 otherwise
% b       : detrended amplitude, smoothed away from highest amplitudes 
%           (values mostly near 0, but high at the 'bottom')
%
           
           
[nd, np] = size(a);   % maybe extend later to handle all 4 beams at once

if (nargin == 5)
   if ~isodd(mfilt)
      fprintf('ERROR in bigtarget: median filter must be odd')
      bmask = [];
      b = [];
      return
   end
   run_mfilt = 1;
   mfilt_tol = 10;  %AGC units -- toss out outliers  (should be a variable?)
else
   run_mfilt = 0;
end

ii = 1:nd;

%run a boxcar to smooth the amplitude
switch nside
   case 1
      aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :)) / 3;
   case 2
      aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :) + ...
                a([1, 1, 1:(nd-2)], :) + a([3:nd, nd, nd], :)) / 5;
   case 3
      aa = (a + a([1, 1:(nd-1)], :) + a([2:nd, nd], :) + ...
                a([1, 1, 1:(nd-2)], :) + a([3:nd, nd, nd], :) + ...
                a([1, 1, 1, 1:(nd-3)], :) + a([4:nd, nd, nd, nd], :)) / 7;
   case 0
      aa = a;
   otherwise
      error('nside must be in range 0-3');
end

%taper for endpoints
wu = ones(nd, 1);
wu(1:offset) = (0:(offset-1)).' / offset;
wu(nd:-1:(nd - offset + 1)) = 2 - wu(1:offset);
wd = flipud(wu);
%propagate to the right for nprofiles
wu = wu(:, ones(np, 1));
wd = wd(:, ones(np, 1));


iu = [ones(1, offset), 1:(nd - offset)];
id = [(1 + offset):nd, nd * ones(1, offset)];

%  unsmoothed version at the apex and smoothed versions
% at the shallower and deeper ends
b = 0.5 * a - .25 * (wu .* aa(iu, :) + wd .* aa(id, :));

% Ignore the top (offset + nside) bins
topmask = ones(nd, np);
topmask(1:(offset+ nside), :) = 0;


if (run_mfilt)
   if (length(b(1,:)) >= mfilt)
      for ii=1:length(b(:,1))
         b(ii,:) = medfilt(b(ii,:), mfilt, mfilt_tol);
      end
   end
end




bmask = (b > thresh) & topmask;




