function bmask1 = bottommask(bmask, b, frac, extra)
% function bmask1 = bottommask(bmask, b, frac, extra)
% inputs:
% bmask and b are the output of bigtarget.m.
%   bmask:  1 at the deduced max amplitude bin, 0 elsewhere
%   b:      detrended version of amplitude with bottom exaggerated; 
%           maximum at the bottom interference bins
% frac:  the unaffected fraction of the vertical distance to
%        the bottom; it is 0.85 for 30-degree beams.  If we had good
%        pitch and roll, the correct number could be calculated for
%        each beam and ping; but no provision is made for this now.
% extra : added number of bins to lop off--typically 1 to
%        account for the blanking.
%
% output:  
% bmask1 is a mask with 1 at the  and below bottom interference
%              depth, and zero elsewhere. Includes 

% Vectorization tricks make the code a bit obscure.
%
% EF 2000/09/26
% JH 2003/10/24 added 3D


if length(size(bmask)) == 3
   [nbins, nbeams, nprofs] = size(b);
   
   bmask1 = bottommask(squeeze(bmask(:,1,:)), ...
                       squeeze(b(:,1,:)), frac, extra);
   for beami = 2:nbeams
      bmask1 = cat(3, bmask1, bottommask(squeeze(bmask(:,beami,:)), ...
                                         squeeze(b(:,beami,:)), frac, extra));
      
   end
   bmask1 = permute(bmask1,[1 3 2]);
else
   [nbins, nprofs] = size(b);
   
   indices = (1:nbins).';
   indices = indices(:, ones(nprofs, 1));
   
   bmax = max(b .* bmask);
   bmax = bmax(ones(nbins, 1), :);
   bmaxmask = cumsum((b == bmax) & bmask); % looking down each column, it
                                           % will be zero until the first 
                                           % bmax value
   bmaxmask = (bmaxmask > 0);
   b_ind = indices;
   b_ind(~bmaxmask) = NaN;
   b_ind = min(b_ind);
   
   b_ind = floor(frac * b_ind - extra);
   
   b_ind = b_ind(ones(nbins, 1), :);
   
   bmask1 = (indices >= b_ind);  % 1 below the adjusted depth of the bottom
end



