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
% 2003/02/21 - adding rzapbins and badtimes; moved flag definition to here
%  from aflagit_setup.m

%set the flags
pfc.flag.pgcutoff           = 2^0;
pfc.flag.trimtopbins        = 2^1;
pfc.flag.badpgrefl          = 2^2;
pfc.flag.maxampbin          = 2^3;     %flag for maxamp bin
pfc.flag.bottom             = 2^4;     %flag for simulating bottom-editing
pfc.flag.ecutoff            = 2^5;
pfc.flag.orphan             = 2^6;
pfc.flag.wire               = 2^7;
pfc.flag.jitter             = 2^8;     %flag for high jitter
pfc.flag.underway           = 2^9;     %flag to flag underway data
pfc.flag.rzapbins           = 2^10;    %flag for manually-edited bins
                                       %used in manflags function (gautoedit)
                                       % and alistit.m
pfc.flag.badtimes           = 2^11;    %flag for manually-editedprofiles
                                       %used in manflags function (gautoedit)
                                       % and alistit.m
pfc.flag.glitch             = 2^12;    %flag glitches in time, but only if 
                                       % accelerating.  bins are independent.
                                       
                                       
%special case:  need to have a way to return flags
if (nargin == 0)
   flagged = pfc.flag;
   fprintf('special case: no arguments returns with flags\n')
   return
end


%% now continue 

if (isemptys(data))
   flagged = [];
   return
else
   nbins = length(data.bins);
end
   

%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;

% not written to disk: used to show effect of editing
% also used to mask bottom when doing other editing
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.bottom);
   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

if (~isemptys(data.smmps))
   underway = (data.smmps >= pfc.cfg.underway);
   pfc.mask.underway = (ones(nbins,1)*underway) * pfc.flag.underway;
else
   pfc.mask.underway = zeros(nbins,length(data.dday));
end



%%%%%%%%%%%%%%%%%%%%%%%%%  error velocity cutoff %%%%%%%%%%%%%%%%%%%%%%%%%%%
%plain old error velocity cutoff
pfc.mask.ecutoff( find(abs(data.e) >= pfc.cfg.ecutoff)) = pfc.flag.ecutoff;

%wire interference:
binnums = pfc.cfg.wire_binrange(1):pfc.cfg.wire_binrange(2);
%on station and in these bins, errvel too big
testbins =  ismember(data.bins(:), binnums(:)) * ones(size(data.dday));
if ~isemptys(data.smmps)
   onstation = repmat(data.smmps,nbins,1) < pfc.cfg.wire_onstation;
else
   onstation = zeros(nbins,length(data.dday));
end

bad_e =  abs(data.e) >= pfc.cfg.wire_ecutoff;

%slide left, right, up, down by one

badbins = slide_lrud(testbins & onstation & bad_e);
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;

%%%%%%%%%%%%%%%%%%%%%%%% glitch %%%%%%%%%%%%%%%%%%%%%%
%% glitch is 0 if bad
%% do this INCLUDING application of orphans first

pfc.cfg.topglitchbin = max(1,pfc.cfg.topglitchbin);
pfc.cfg.topglitchbin = min(length(data.bins),pfc.cfg.topglitchbin);

glitch_is_0 = ones(size(data.uabs_rawnav));
for ibin=pfc.cfg.topglitchbin:length(data.bins)
   goodi = find(~isnan( ...
      dotmmask(pfc.mask.orphan(ibin,:).*data.uabs_rawnav(ibin,:))));
   
      
   if (length(goodi) >= 3)
      pm =  glitch(data.uabs_rawnav(ibin,goodi),pfc.cfg.glitchfactor) .* ...
               glitch(data.vabs_rawnav(ibin,goodi),pfc.cfg.glitchfactor);
   else
      pm=ones(size(goodi));
   end
   if (~isemptys(pm))
      pmi = find(pm==0);
      if ~isemptys(pmi)
         glitch_is_0(ibin,goodi(pmi))=0;
      end
   end
end

glitchmask = ~glitch_is_0;
degpermin0 = [0 diff(data.heading)./diff(data.dday*24*60)];
degpermin1 = [diff(data.heading)./diff(data.dday*24*60) 0];
dpm = max(abs(degpermin0), abs(degpermin1));


turnmask = (abs(dpm) >= pfc.cfg.agility);                   %turning
if ~isemptys(data.rawmps)
   fastmask = (data.rawmps >= pfc.cfg.wire_onstation);         %fast
else
   fastmask = zeros(1, length(data.dday));
end
   
spinmask = (abs(dpm) >= pfc.cfg.degperminute);              %very strong turning 
turning = (turnmask & fastmask) | spinmask;
accmask = repmat(turning, length(data.bins),1);

pfc.mask.glitch =  pfc.flag.glitch*(accmask&glitchmask);


%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

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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,:)] ;

