function [flagged, pfc] = aflagit(data, config, pfc)
% [flagged, pfc1] = aflagit(data,config,pfc)
% output: matrix of flags (zeros are unflagged data) with values
%         corresponding to pfc1.mask
%
%screening stuff "pfc" is a structure with "profile flag configuration"
% will output flagged=[] if size of matlab files is zero

% Jules 2000/11/27

if (isempty(data))
   flagged = [];
else
   nbins = length(data.bins);

   
   %initialize all masks
   flagnames = fieldnames(pfc.flag);
   dummy = zeros(size(data.e));
   pfc.mask.pgcutoff = dummy;  %initialize the structure if not here yet
   for flagi = 1:length(flagnames);
      pfc.mask = setfield(pfc.mask, flagnames{flagi}, dummy);
   end

   
   %%%%%%%%%%%%%%%%%%%%   amp: bottom return  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% should do amp first so the high return at the base is eliminated first
   %%% (otherwise it will hurt the ringing calculation)
   
   [bm1, b] = bigtarget(data.amp, pfc.cfg.bigtarget_amp, ...
                        pfc.cfg.bigtarget_offset, ...
                        pfc.cfg.bigtarget_avg);
   pfc.mask.maxampbin = bm1 .* pfc.flag.maxampbin;
   pfc.mask.bottom  = (bottommask(bm1, b, 0.85, 1)) .* pfc.flag.bottom;
   
   %%%%%%%%%%%%%%%%%%%%%%%% pg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % bad if pg below pgcutoff
   pfc.mask.pgcutoff( find(data.pg < pfc.cfg.pgcutoff)) = pfc.flag.pgcutoff;
   
   
   %%%%%%%%%%%%%%%%%%%% trim top bins %%%%%%%%%%%%%%%%%%%%
   %three reasons to flag them
   % (1) flag everything shallower than the shallowest location of the maximum pg 
   %           (or near maximum) - must do this calculation on pg with
   %           amp bounce from the bottom already out
   % (2) trim the top N bins
   % (3) eliminate everything with too low an amplitude in the shallowest bins
   %     (constant with speed, linearly decreasing with bin
   
   
   % start with 0 for bad data, 1 for good data. then use (~cumsum) to swap
   % (1) everything shallower than the shallowest location of max pg 
   if (pfc.cfg.maxpgbin_fudge == -1) % value of -1 disables it
      pg_is_greater = ones(size(data.pg));
   else
      % matrix of biggest pg per profile
      maskedpg = data.pg .* (~pfc.mask.pgcutoff) .* (~pfc.mask.maxampbin);
      maxpg_fM = ones(nbins,1)*(max(maskedpg)- pfc.cfg.maxpgbin_fudge);
      pg_is_greater = (maskedpg >= maxpg_fM);  %locations of highest pg
   end
   
   % (2) exclude top bins anyway - don't include them in looking for the bins
   pg_is_greater(1:pfc.cfg.trimtopbins,:) = ...
       zeros(size(pg_is_greater(1:pfc.cfg.trimtopbins,:)));
   %now get bin numbers of values greater than (max-fudge)
   binM = ([1:nbins]')*(ones(size(data.dday)));
   pg_is_greater_bins = pg_is_greater.*binM;
   pgmask1 = (~cumsum(pg_is_greater_bins));   %this is ones for flagged, 0 else
   
   % (3) 
   ampcutoffM = pfc.cfg.amp_cutoff_atbin1 - pfc.cfg.ampdecrperbin.*(binM-1);
   ampmask1 = (data.amp < ampcutoffM);
   
   pfc.mask.trimtopbins = (ampmask1 | pgmask1) * pfc.flag.trimtopbins;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bad ref layer %%%%%%%%%%%%%%%%%%%%
   reflbins = [pfc.refl_startbin : pfc.refl_endbin];

   numgoodbins = sum( ...
      (data.pg(reflbins,:) >= pfc.cfg.badpgrefl) & ...
      (pfc.mask.trimtopbins(reflbins,:) == 0) );
   
   badrefl = (numgoodbins <= pfc.cfg.badpgrefl_nbins);
   pfc.mask.badpgrefl = (ones(nbins,1)*badrefl) * pfc.flag.badpgrefl;
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% high jitter %%%%%%%%%%%%%%%%%%%%
   badjitter = (data.ref_std >= pfc.cfg.jitter_cutoff);
   pfc.mask.jitter = (ones(nbins,1)*badjitter) * pfc.flag.jitter;
   

   %%%%%%%%%%%%%%%%%%%%%%%  bad underway profiles %%%%%%%%%%%%%%%%%%%%%%
   if (isfield(pfc, 'reject_underway'))
      pfc = rmfield(pfc, 'reject_underway');
   end
   
   if (~isfield(pfc.mask, 'underway'))
      pfc.cfg.underway = 100;           %if it doesn't exist, disable it
      pfc.flag.underway    = 2^9;
   end
   
   underway = (data.smmps >= pfc.cfg.underway);
  
   pfc.mask.underway = (ones(nbins,1)*underway) * pfc.flag.underway;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%  error velocity cutoff %%%%%%%%%%%%%%%%%%%%%%%%%%%
   %plain old error velocity cutoff
   pfc.mask.ecutoff( find(abs(data.e) >= pfc.cfg.ecutoff)) = pfc.flag.ecutoff;
   
   %wire interference:
   osi = find(data.smmps < pfc.cfg.wire_onstation);
   binrange = pfc.cfg.wire_binrange(1):pfc.cfg.wire_binrange(2);
   %in these bins, and   errvel too big
   badbins =  ismember(data.bins(:), binrange(:)) * ones(size(data.dday)) & ...  
       (abs(data.e) >= pfc.cfg.wire_ecutoff);           
   badbins = slide_lrud(badbins);   %slide left, right, up, down by one
   pfc.mask.wire_ecutoff = badbins * pfc.flag.wire;
   
   %%%%%%%%%%%%%%%%%%%%%%%%% ensure enough supporting profiles %%%%%%%%%%
   % criterion will be to look for profiles with enough bins, then make
   % sure they have neighbors
   
   
   masknames = fieldnames(pfc.mask);
   initially_flagged = zeros(size(getfield(pfc.mask, masknames{1})));
   for maski = 1:length(masknames)
      initially_flagged = ...
          initially_flagged + getfield(pfc.mask, masknames{maski});
   end
   
   
   %make a 1xn vector to get the orphan mask
   rli =[ pfc.refl_startbin : pfc.refl_endbin];  %ref layer index
   orphan_profs = zeros(size(initially_flagged(1,:)));
   okprfs = (sum(~initially_flagged(rli,:)) >= pfc.cfg.orphan_numbins);
   orphan_profs(:,okprfs) = 1;     %orphan indicator: 1 is a profile
   
   oo = orphan_profs;  %store the original
   for slidei = 1:pfc.cfg.orphan_slide 
      %slide right
      padleft = zeros(1,slidei);
      ii = [1 : length(orphan_profs)-slidei];
      orphan_profs = orphan_profs + [padleft oo(ii)];
      %slide left
      padright = zeros(1,slidei);
      ii = [slidei+1 : length(orphan_profs)];
      orphan_profs = orphan_profs + [oo(ii) padright];
   end
   
   okprfs = (orphan_profs >= (1 + 2*pfc.cfg.orphan_slide));
   pfc.mask.orphan = (ones(nbins,1) * (~okprfs)) * pfc.flag.orphan;

   
   %now assign flags --------------------------------------------------
  
   
   masknames = fieldnames(pfc.mask);
   flagged = zeros(size(getfield(pfc.mask, masknames{1})));
   for maski = 1:length(masknames)
      flagged = flagged + getfield(pfc.mask, masknames{maski});
   end
   
end   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%functions

function badbins = slide_lrud(badbins);
%badbins = slide_lrud(badbins);
% 
% cavalier approach: take any bad bins and slide the bad regions
% left, right, up, and down by 1 bin each  (bitwise or the results)

%could potentially generalize this function for more sliding, and use it for 
% eliminating orphaned profiles as well.

rowpad = zeros(size(badbins(1,:)));
colpad = zeros(size(badbins(:,1)));
[ncols, nrows] = size(badbins);

badbins = badbins | ...
          [badbins(:,2:end) colpad] | ...     %slide left  [2 3 ... end X]
          [colpad badbins(:,1:end-1)] | ...   %slide right [X 1 2 ... end-1]
          [badbins(2:end,:) ; rowpad] | ...
          [rowpad ; badbins(1:end-1,:)] ;

