function [aa, r] =  avg(instclass, r, a, p)
% [a, r] = avg(instclass, r, a, p)
% inputs:
% r is a raw data structure; read by load_raw or by get_rndata 
%         (with 'return3d') with ancillary information
% a is initialized for one profile
% p is a structure of parameters that control the averaging and editing
%      p.editcfg is a structure that overrides pingedit defaults 
%        (see get_pedefaults)
%      
%      p.pingedit    level of pingediting to apply
%                     'none'  (or 0)        no pingediting at all
%                     'ampbt'               pre-ping edit = amp and bottom only
%                     'amp'                 pre-ping edit = amp only
%                     'bt'                  pre-ping edit =  bottom only
%                     'pings'               water only (no amp or bottom)
%                     'all'   (or 1)        amp and bottom, and water
%                     'debug'               return both answers
%                                             'none' and 'all' as
%                                             a.vel_noedit, a.vel_pedit
%                                         (a.vel is a.vel_pingedit)
%     p.use_refavg    1  (or 0)              use refavg? or (mstdgap??) 
%                                                  subsets of editing
%                                       
%
% (these are used to determine which shipu and shipv are used in the average:
%      (only take profiles with sufficient PG in the specfied layer)
%     p.min_reflpg      minimum PG cutoff integer[0-100]  (default = 20)
%                       for reference layer (NOITE: only used if 
%                       pingedit == 'none' or pingedit == 'debug'
%      p.ilayer         layer to use min_reflpg as a test for shipuv and
%                       for reference  layer averaging (default = 1:10)
%
%   p.beamangle         if present, overrides value in r.config
%   p.showverbose       plot diagnostics  [0 or 1]        (default  = 0)
%   p.instrument        used only to select specific irregular transformations
%                              (two exist: "kmson", "kmos")
%outputs:
% a: averaged parameters with pingediting (if invoked).
% fields for load_raw:
%     lastdday, depth, config, bt,
%     vel, amp, pg, cor, sw     
% additional for load_raw matdata
%     lastdday, lastlon, lastlat, num_pings, 
%     meanheading, peanpitch, meanroll, meantemperature, meansoundspeed
%     ens_mean.lon, ens_mean.lat, ens_mean.uship, ens_mean.vship
% additional fields for debugging (if requested)
%     vel_noedit (vel_pedit) % same as vel, but edited as stated
% 

%JH history: - cut this out of bb/load_raw.m to write nb/load_raw
%          - versions propagated to rdi (for os and bb), rdi/@nb, sontek/@son
% 2002/08/30: consolidating in rawadcp
% 2003/12/16: taking out pingedit, prepedit, streamlining get_pedefaults


%%% DEAL WITH THIS!!!
%fprintf('--------- replacing struct a in avg.m\n')
%a = struct;

cfg.min_reflpg = 20;
cfg.showverbose = 0;
cfg.use_refavg = 1;
cfg.pingedit = [];
cfg.instrument = [];

% fill the rest of the defaults
fnames = fieldnames(cfg);
for fnamei = 1:length(fnames)
   if isfield(p, fnames{fnamei})
      cfg = setfield(cfg, fnames{fnamei}, getfield(p, fnames{fnamei}));
   end
end


% get defaults for ping editing 
editcfg = get_pedefaults(instclass, r.config.pingtype);
if isfield(p, 'editcfg')
   editcfg = fillstruct(editcfg, p.editcfg);
end
cfg.editcfg = editcfg;

% fix other things with defaults that depend on size of data.cuv
%% this has to be here and in avg in case pingedit is called NOT from avg
if (isstr(cfg.editcfg.rl_bins) )
   if strcmp(cfg.editcfg.rl_bins, 'all')
      cfg.editcfg.rl_bins = [1:length(r.vel(:,1))]';
   else
      fprintf('string %s not allowed for editcfg "rl_bins" field\n', ...
              editcfg.rl_bins)
      return
   end
end


if ~isfield(p, 'beamangle');
   p.beamangle = r.config.beamangle;
end


if ~isstr(cfg.pingedit)
   if cfg.pingedit == 0
      cfg.pingedit = 'none';
      fprintf('setting pingedit flag to ''none''.\n')
   elseif cfg.pingedit == 1
      cfg.pingedit = 'all';
      fprintf('setting pingedit flag to ''all''.\n')
   else
      help(mfilename)
      a=[];r=[];
      fprintf('incorrect pingedit flag to avg.m -- use ''all'' or ''none''.\n')
      return
   end
end

%%%% rotate the data into earth coordinates if necessary
needs_beamenu =  ~isfield(r, 'enu') &  ~strcmp(r.config.coordsystem, 'earth');

if (needs_beamenu)   
   %get umeas, vmeas to do the editing in geographic coordinates
   enu =  beam_enu(instclass, r.vel,...
                   r.heading + p.h_align, 0, 0,...
                   'instrument', cfg.instrument,...
                   'convex', r.config.convex,...                   
                   'error', 1, 'beamangle', p.beamangle);
   r.emeas = squeeze(enu(:,4,:));
   r.bt.enu = ...
       beam_enu(instclass, r.bt.vel,...
                r.heading + p.h_align, 0, 0,...
                'instrument', cfg.instrument,...
                'convex', r.config.convex,...
                'error', 1, 'beamangle', p.beamangle);
else
   if (length(r.dday) > 1) & (sum(abs(r.heading)) == 0)
      fprintf('WARNING: %d pings averaged, all have heading=0!!!\n', ...
              length(r.dday))
   end
   if isfield(r, 'enu') %data were converted after read (vel is still beam vels)
      enu = r.enu;
   else
      enu = r.vel;
   end
   r.emeas = squeeze(enu(:,4,:));
end

[nbins, nbeams, nprofs] = size(r.vel);
r.enu = enu;
r.umeas = squeeze(enu(:,1,:));
r.vmeas = squeeze(enu(:,2,:));
r.cuv = r.umeas + i*r.vmeas;


if ismember(cfg.pingedit, {'ampbt', 'all', 'debug', 'bt', 'amp'})
   [r, bt_bmask, bt_wmask, hiamp] = ...
       prepedit(instclass, r, cfg.editcfg, cfg.showverbose, cfg.pingedit);
   %showverbose plots to figure(1)
   % this is now in prepedit:   r.bmask = bt_bmask | bt_wmask | hiamp;
   if strcmp(cfg.pingedit, 'amp') % didn't do bt
      bt_bmask = zeros(size(r.cuv));
      bt_wmask = zeros(size(r.cuv));
      r.bmask  = zeros(size(r.cuv));
   end
   if strcmp(cfg.pingedit, 'bt') % didn't do amp
      hiamp  = zeros(size(r.cuv));
      r.cleaned_amp = r.amp;
   end
else
   bt_bmask = zeros(size(r.cuv));
   bt_wmask = zeros(size(r.cuv));
   hiamp  = zeros(size(r.cuv));
   r.bmask  = zeros(size(r.cuv));
   r.cleaned_amp = r.amp;
end
r.bt_bmask = bt_bmask;
r.bt_wmask = bt_wmask;
r.hiamp = hiamp;
r.editcfg = cfg.editcfg;
%% set this so codas will always look for BT data

a.config.bt_was_on = 1;


%---------------------------------

%% apply pingediting to measured velocity if requested


if ismember(cfg.pingedit, {'all', 'debug', 'pings'})
   [instpflag, vprofile, trefl,  resid, binmasks] =  ...
       pingedit(instclass, r, 'editcfg', cfg.editcfg, ...
                'showverbose', cfg.showverbose,...
                'pflag', r.bmask);
   r.instpflag = instpflag;
   r.binmasks = binmasks;
else
   instpflag = zeros(size(r.cuv));
end


% on all 4 components at once.

if ~isfield(a,'num_pings')
   a.num_pings = nprofs;
end


%% return pingedited or unedited (whichever was requested) as a.vel
%% also return as named fields: a.vel_pedit a.vel_noedit 
%%         if requested (for debugging)

%% these are earth coordinates velocities
%% -- using "vel" so we can loop through indices. rename below
 
if ismember(cfg.pingedit, {'all','debug', 'pings'})
   bad_mask = r.bmask | instpflag | isnan(squeeze(r.vel(:,1,:)));
   for ii = 1:4
      uu = squeeze(enu(:, ii, :));
      uu(bad_mask) = NaN; 
      if (cfg.use_refavg)
         [uu, ulayer] = refavg(uu, cfg.editcfg.rl_bins, 5);
         a.enu(:,ii) = uu + mstdgap(ulayer);
      else
         a.enu(:,ii) = mstdgap(uu, 2);
      end
   end
   pingedit_pg = 100 *  sum(~bad_mask, 2)  / a.num_pings;
   a.enu_pedit = a.enu;
end

if ismember(cfg.pingedit, {'none', 'ampbt', 'debug'})
   % Eliminate from the reference layer levels with insufficient samples.
   % NOTE: this does not use a short run of good deep data
   % (i.e. distribution of pings not considered)
   bad_mask = isnan(squeeze(r.vel(:,1,:)));
   noedit_pg = 100 *  sum(~bad_mask')  / a.num_pings;
   ipg = find(noedit_pg(cfg.editcfg.rl_bins) > cfg.min_reflpg);
   if length(ipg) < 2,  % No decent ref layer--just average.
      for ii = 1:4
         uu = squeeze(enu(:, ii, :));
         uu(bad_mask) = NaN;
         a.enu(:,ii) = mstdgap(uu, 2);
      end
   else
      for ii = 1:4
         uu = squeeze(enu(:, ii, :));
         uu(bad_mask) = NaN;
         if (cfg.use_refavg)
            [uu, ulayer] = refavg(uu, cfg.editcfg.rl_bins(ipg), 5);
            a.enu(:,ii) = uu + mstdgap(ulayer); % simple mean for now
         else
            a.enu(:,ii) = mstdgap(uu, 2);
         end
      end
   end
   a.enu_noedit = a.enu;
end
 
a.pg(:,:) = 0;
if ismember(cfg.pingedit, {'none', 'ampbt', 'debug'})
   a.enu = a.enu_noedit;
   a.pg(:, 4) = round(noedit_pg);
   a.amp = round(mean(r.amp, 3));
end
if ismember(cfg.pingedit, {'all', 'debug', 'pings'})
   a.enu = a.enu_pedit;
   a.pg(:, 4) = round(pingedit_pg);
   a.amp = round(mean(r.cleaned_amp, 3));
end


% output average lon, lat; average shipu, shipv; 
% average oceanu, oceanv (was get_ensmean)
% use PG in reflayer to screen shipu


% NOTE: there's alot of code for debugging in this loop.
% should probably make it optional

goodpos_mask = ~isnan(sum(r.nav.txy2));  %ok positions
posi = find(goodpos_mask);
if (length(posi) > 1)
   a.lastdday = r.nav.txy2(1,max(posi));  
   a.lastlon     = r.nav.txy2(2,max(posi));
   a.lastlat     = r.nav.txy2(3,max(posi));
   
   [junk, x, y] =  ll2mps(r.nav.txy2(1,:), r.nav.txy2(2,:), r.nav.txy2(3,:));

   % calculate ship speed from the good position values
   %% (now stored in "get_goodshipuv.m") replaced with:
   uship = diff(x)./diff(86400*r.nav.txy2(1,:));
   uship = [uship(1) uship];
   vship = diff(y)./diff(86400*r.nav.txy2(1,:));
   vship = [vship(1) vship];
   
   ens_mean.lon = mstdgap(r.nav.txy2(2,:));   
   ens_mean.lat = mstdgap(r.nav.txy2(3,:));
   ens_mean.uship = mstdgap(uship);   
   ens_mean.vship = mstdgap(vship);   
   a.ens_mean = ens_mean;
   % end of get_goodshipuv replacement
   
else
   fprintf('---> no good positions in this ensemble\n')
   a.lastdday = NaN;
   a.lastlon     = NaN;
   a.lastlat     = NaN;
   uship = NaN*ones(size(r.dday));
   vship = NaN*ones(size(r.dday));
   ens_mean.lon = NaN;
   ens_mean.lat = NaN;
   if ismember(cfg.pingedit, {'none', 'ampbt'})
      ens_mean.uship = NaN;
      ens_mean.vship = NaN;
   else
      ens_mean.uship =  NaN; 
      ens_mean.vship =  NaN;
   end   
   a.ens_mean = ens_mean;
end


%make an entry with uabs, vabs, (and w, and e) so we can used getfstats
% the entries are all the same, but the flags are used differently in
% getfstats, so we need different names for the structures

if (0)
   %sontek doesn't have this  (is this necessary??)
   if (isa(instclass, 'son'))
      if (~isfield(r, 'sw'))
         a.sw = 0*a.pg;
      else
         r.sw(r.sw==255) = NaN;
         for beamnum = 1:4
            a.sw(:,beamnum) = mstdgap(squeeze(r.sw(:,beamnum,:))')';
         end
      end
   end
end


%extra goodies for codas database
% temperature, heading, pitch, roll
a.meantemperature         = mstdgap(r.temperature);
a.meansoundspeed          = mstdgap(r.soundspeed);
a.meanpitch               = mstdgap(r.pitch);
a.meanroll                = mstdgap(r.roll);
unwrapped_head = 180*(unwrap(r.heading*pi/180))/pi;
a.meanheading             = mstdgap(unwrapped_head);


%% extra goodies for the mat files
% sw and cor

a.lastdday = r.corr_dday(end);
%bottom track: expects u,v,w,e
a.bt = mk_bt(r,1);

if (cfg.showverbose)
   f1;
   titlestr = [' BT bottom mask'];
   plotbot(5,r,bt_bmask,titlestr);

   titlestr = ['waterping bottom mask'];
   plotbot(6,r,bt_wmask,titlestr);   

   r.vprofile = vprofile;
   r.trefl = trefl;
   r.resid = resid;
   %%% keyboard is contained within "avg_showverbose" f2 and f3
   avg_showverbose(instclass, a, r);
end


aa.enu                 = a.enu;
aa.amp                 = a.amp;
aa.pg                  = a.pg;
if isfield(a,'cor'), aa.cor                 = a.cor; end
if isfield(a,'sw'),  aa.sw                  = a.sw;  end
aa.lastdday            = a.lastdday;
aa.depth               = a.depth;
aa.config              = a.config;
aa.bt                  = a.bt;
aa.umeas               = a.enu(:,1);
aa.vmeas               = a.enu(:,2);
aa.num_pings           = a.num_pings;
aa.lastdday            = a.lastdday;
aa.lastlon             = a.lastlon;
aa.lastlat             = a.lastlat;
aa.meanheading         = a.meanheading;
aa.meanpitch           = a.meanpitch;
aa.meanroll            = a.meanroll;
aa.meantemperature     = a.meantemperature;
aa.meansoundspeed      = a.meansoundspeed;
aa.ens_mean            = a.ens_mean;
