function a = avg(instclass, r, a, p)
% r is a raw data structure; a is initialized for one profile
% p is a structure of parameters that control the averaging

%JH cut this out of bb/load_raw.m so I could use it to write nb/load_raw
%     (it may be too specific)

enu = beam_enu(r.vel, r.heading + p.h_align, 0, 0, 'error', 1, ...
               'geometry', p.geometry);

no_vel_mask = isnan(squeeze(enu(:, 1, :)));

ev_mask = abs(enu(:, 4, :)) > p.max_ev;
% Calculate other masks and combine with logical or.

bmask = zeros(size(no_vel_mask));

%%% bad  implimention of  bottom tracking test
%%% (os) if r.config.bt.BP > 0,
if ~(isnan(sum(r.bt.pg(:,1)))) %no bottom tracking: mask
   btd = min(r.bt.range);
   dd = r.depth.' - r.config.tr_depth;
   bmask = mask_below(dd, btd);
else
   for ii = 1:4
       [bm1, b] = bigtarget(squeeze(r.amp(:, 1, :)), ...
            p.bigtarget_amp, p.bigtarget_offset, p.bigtarget_avg);
       bmask = bmask | bottommask(bm1, b, 0.85, 1);
   end
end

bad_mask = squeeze(ev_mask) | bmask;
% The following could probably be made more efficient by
% simply passing the bad mask into refavg, and optimizing
% refavg to use a single set of masks and counts for all
% iterations.

% maybe make refavg also operate
% on all 4 components at once.

% Eliminate from the reference layer levels with insufficient
% samples.
pg = 100 *  sum(~(bad_mask | no_vel_mask), 2)  / a.num_pings;
ipg = find(pg(p.ilayer) > p.min_pg);

if length(ipg) < 2,  % No decent ref layer--just average.
   for ii = 1:4
      uu = squeeze(enu(:, ii, :));
      uu(bad_mask) = NaN;
      a.vel(:,ii) = mstdgap(uu, 2);
   end
else
   for ii = 1:4
      uu = squeeze(enu(:, ii, :));
      uu(bad_mask) = NaN;
      [uu, ulayer] = refavg(uu, p.ilayer(ipg), 5);
      a.vel(:,ii) = uu + mstdgap(ulayer); % simple mean for now
   end
end

a.pg(:,:) = 0;
a.pg(:, 4) = round(pg);

a.sw = round(mean(r.sw, 3)); % Or, could use the masks to
                               % get the average of the values actually used.

a.amp = round(mean(r.amp, 3));


%------------------------------
function mask = mask_below(depth, bottom)
np = length(bottom);
nd = length(depth);
depth = depth(:, ones(np, 1));
bottom = bottom(ones(nd, 1), :);

mask = depth > bottom;


