function d = edit_data(d,p)
% function d = edit_data(d,p)
%
% perform data editing (e.g. sidelobes, previous-ping interference, &c)
%
% HISTORY:
%  Jul  3, 2004: - implemented side-lobe contamination editing
%  Jul 15, 2004: - implemented PPI contamination editing
%  Jul 16, 2004: - implemented cross-contamination editing
%  Jul 18, 2004: - added parameter defaults
%  Jul 21, 2004: - moved bin masking from [loadrdi.m]
%  Jul 23, 2004: - implemented ensemble skipping
%  May 18, 2006: - incorporated simple fix for multi-ping ensembles
%		   provided by Mattew Alford
%		 - disabled edit_spike_filter by default
%		 - changed edit_spike_filter_max_curv default value

%======================================================================
%                    E D I T _ D A T A . M 
%                    doc: Sat Jul  3 17:13:05 2004
%                    dlm: Thu May 18 12:20:40 2006
%                    (c) 2004 A.M. Thurnherr
%                    uE-Info: 50 37 NIL 0 0 72 3 2 8 NIL ofnI
%======================================================================

%----------------------------------------------------------------------
% Initialize & Set Default Parameters
%----------------------------------------------------------------------

d.ts_edited = d.ts; % save for plotting/inspection

% Set to list of bins to remove from data.
p=setdefv(p,'edit_mask_dn_bins',[]);
p=setdefv(p,'edit_mask_up_bins',[]);

% Set to 1 to remove side-lobe contaminated data near seabed and
% surface.
p=setdefv(p,'edit_sidelobes',1);

% Set to 1 to implement time-domain spike filter on the data; 
% this removes interference from other acoustic instruments but,
% more importantly, can get rid of PPI when staggered pings
% are used. USING THIS FILTER WITH A CURVATURE CRITERION THAT IS
% TOO STRICT CAN SERIOUSLY DEGRADE THE VELOCITY PROFILES. THEREFORE,
% IT HAS BEEN DISABLED BY DEFAULT IN VERSION IX_3.
p=setdefv(p,'edit_spike_filter',0);

% Spike filtering is done using 2nd-difference
% peak detection in time. This parameter gives the maximum
% target-strength 2nd derivative that's allowed. Set to larger
% values to weaken the filtering. (Check figure 14 to see if
% filter is too strong or too weak.) Setting the value of this
% parameter too low will seriously degrade the velocity profiles,
% without any apparent sign of trouble. The optimal value of this
% parameter depends on the instrument type.
p=setdefv(p,'edit_spike_filter_max_curv',5);

% Set to 1 to remove data contaminated by previous-ping interference.
% NOTES: - Using the spike filter seems to work more robustly, as long
%     	   as staggered pings are used. Great care must be taken,
%	   however, to chose a good value for edit_spike_filter_max_curv.
p=setdefv(p,'edit_PPI',0);

% PPI layer thickness in meters; the value is taken directly from Eric
% Firing's default (2*clip_margin = 180m).
p=setdefv(p,'edit_PPI_layer_thickness',180);

% max distance from seabed at which PPI should be removed. This is
% an observed parameter and depends on the clarity of the water.
% Check Figure 14 to see whether this should be changed.
p=setdefv(p,'edit_PPI_max_hab',1000);

% set this vector to enable skipping of ensembles; skipping vector
% is wrapped around, i.e. [1 0] skips all odd ensembles, [0 1 0] skips
% ensembles 2 5 8 11.... This filter is useful to process the casts
% with only half the data to see whether the two halves agree, which
% probably means that the cast can be trusted. Note that if staggered
% ping setup is used to avoid PPI, the skipping vector should leave
% adjacent ensembles intact, i.e. use something like [1 1 0 0] and
% [0 0 1 1].
p=setdefv(p,'edit_skip_ensembles',[]);

%----------------------------------------------------------------------
% Bin Masking
%----------------------------------------------------------------------

if ~isempty(p.edit_mask_dn_bins) | ~isempty(p.edit_mask_up_bins)

  nbad = 0;
  if length(d.zu) > 0
    for bi=1:length(p.edit_mask_up_bins)
      bn = length(d.zu)+1 - p.edit_mask_up_bins(bi);
      nbad = nbad + length(find(isfinite(d.weight(bn,:))));
      d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
    end
  end
  if length(d.zd) > 0
    for bi=1:length(p.edit_mask_dn_bins)
      bn = length(d.zu) + p.edit_mask_dn_bins(bi);
      nbad = nbad + length(find(isfinite(d.weight(bn,:))));
      d.weight(bn,:) = NaN; d.ts_edited(bn,:) = NaN;
    end
  end

  disp(sprintf(' bin masking               : set %d weights to NaN',nbad));

end % if bin masking enabled

%----------------------------------------------------------------------
% Side-Lobe Contamination
%----------------------------------------------------------------------

if p.edit_sidelobes

  nbad = 0;
  
  % first, the uplooker: d.z is -ve distance of ADCP from surface;
  % Cell_length is in cm, i.e. 0.015*Cell_length is 1.5 x bin size
  % in m --- the same value used by Firing's software
  
  if length(d.zu > 0)
  
    for b=1:length(d.zu)+length(d.zd)
      zlim(b,:) = (1 - cos(pi*d.up.Beam_angle/180)) * d.z ...
		- 0.015*d.up.Cell_length;
    end
    ibad = find(d.izm > zlim);
    nbad = nbad + length(find(isfinite(d.weight(ibad))));
    d.weight(ibad) = NaN; d.ts_edited(ibad) = NaN;
  
  end
  
  % now, the downlooker: p.zbottom is the +ve depth of the sea bed; therefore,
  % -d.z - p.zbottom is the -ve distance from the sea bed
  
  for b=1:length(d.zu)+length(d.zd)
    zlim(b,:) = -p.zbottom ...
	      + (1 - cos(pi*d.down.Beam_angle/180)) * (d.z+p.zbottom) ...
	      + 0.015*d.down.Cell_length;
  end
  ibad = find(d.izm < zlim);
  nbad = nbad + length(find(isfinite(d.weight(ibad))));
  d.weight(ibad) = NaN; d.ts_edited(ibad) = NaN;
  
  disp(sprintf(' side-lobe contamination   : set %d weights to NaN',nbad));

end %if p.edit_sidelobes

%----------------------------------------------------------------------
% Time-Domain Spike Filter
%----------------------------------------------------------------------

if p.edit_spike_filter

  nbad = 0;
  for b=1:length(d.zd)+length(d.zu)
    ibad = find(diff(diff(d.ts(b,:))) < -1*p.edit_spike_filter_max_curv) + 1;
    nbad = nbad + length(find(isfinite(d.weight(b,ibad))));
    d.weight(b,ibad) = NaN; d.ts_edited(b,ibad) = NaN;
  end
  disp(sprintf(' spike filter              : set %d weights to NaN',nbad));

end %if p.edit_spike_filter

%----------------------------------------------------------------------
% Previous-Ping Interference
%----------------------------------------------------------------------

if p.edit_PPI

  % NB: at present, PPI filtering is only implemented for the downlooker
  
  nbad = 0;
  
  % calc ping-intervals; dt(1) contains the difference (in seconds)
  % between the first two ensembles (t(2) - t(1)).
  
  dt = diff(d.time_jul)*86400;

  % adjust for multi-ping ensembles; THIS FIX, PROVIDED BY MATTHEW ALFORD,
  % DOES NOT WORK FOR IRREGULAR PINGING RATES (e.g. 10s ensembles of 3
  % pings each with 1s between pings)

  dt = dt/d.down.Pings_per_Ensemble;
  
  % use the mean sound speed below the approximate expected PPI depth;
  % this is anal but not very expensive; using 1500m/s would be nearly
  % as good.
  
  if isfield(d,'ctdprof_z') & isfield(d,'ctdprof_ss')
    guess_z = -p.zbottom + 1500 * meannan(dt) / 2;
    SS = meannan(d.ctdprof_ss(find(d.ctdprof_z > -guess_z)));
    if ~isfinite(SS), SS = 1500; end
  else
    SS = 1500;
  end
  
  % calculate the depth limits to remove the PPI for all (but the first)
  % ensembles; the beam-angle limits were found to be too conservative
  % and were replaced by a nominal layer_tickness.
  
  %PPI_min_beam_angle = 0.0 * d.down.Beam_angle;
  %PPI_max_beam_angle = 1.2 * d.down.Beam_angle;
  %PPI_max_z = -p.zbottom + SS * dt/2 * cos(pi/180 * PPI_min_beam_angle);
  %PPI_min_z = -p.zbottom + SS * dt/2 * cos(pi/180 * PPI_max_beam_angle);
  
  PPI_hab = SS * dt/2 * cos(pi/180 * d.down.Beam_angle);
  PPI_hab(find(PPI_hab > p.edit_PPI_max_hab)) = inf;
  
  PPI_min_z = -p.zbottom + PPI_hab - p.edit_PPI_layer_thickness / 2;
  PPI_max_z = PPI_min_z + p.edit_PPI_layer_thickness;
  
  % remove the contaminated data from the downlooker bins
  
  for b=length(d.zu)+1:length(d.zu)+length(d.zd)
    ibad = find(d.izm(b,2:end) > PPI_min_z & d.izm(b,2:end) < PPI_max_z) + 1;
    nbad = nbad + length(find(isfinite(d.weight(b,ibad))));
    d.weight(b,ibad) = NaN; d.ts_edited(b,ibad) = NaN;
  end
  
  disp(sprintf(' previous-ping interference: set %d weights to NaN',nbad));

end %if p.edit_PPI

%----------------------------------------------------------------------
% Ensemble Skipping
%----------------------------------------------------------------------

if ~isempty(p.edit_skip_ensembles)

  nskipped = 0;
  iskip = [];
  for i=1:length(p.edit_skip_ensembles)
    if p.edit_skip_ensembles(i)
      iskip = [iskip i:length(p.edit_skip_ensembles):length(d.time_jul)];
    end
  end

  for b=1:length(d.zd)+length(d.zu)
    nskipped = nskipped + length(find(isfinite(d.weight(b,iskip))));
    d.weight(b,iskip) = NaN; d.ts_edited(b,iskip) = NaN;
  end

  disp(sprintf(' ensemble skipping         : set %d weights to NaN',nskipped));

end % if p.edit_skip_ensembles enabled

%----------------------------------------------------------------------
% Plot Results of Editing
%----------------------------------------------------------------------

bin_no = [0];
if length(d.zu) > 0, bin_no = [-length(d.zu):1 bin_no]; end
if length(d.zd) > 0, bin_no = [bin_no 1:length(d.zd)]; end

figure(14);
clf;
orient landscape;
colormap([[1 1 1]; jet(128)]);

subplot(2,1,1);
imagesc(1:size(d.ts,2),bin_no,...
	[d.ts(1:length(d.zu),:); ...
	 ones(1,size(d.ts,2))*NaN; ...
	 d.ts(size(d.ts,1)-length(d.zd)+1:end,:)...
        ]);
csc = caxis;
xlabel('Ensemble #');
ylabel('Bin #');
title('Before Data Editing');

subplot(2,1,2);
imagesc(1:size(d.ts,2),bin_no,...
	[d.ts_edited(1:length(d.zu),:); ...
	 ones(1,size(d.ts,2))*NaN; ...
	 d.ts_edited(size(d.ts,1)-length(d.zd)+1:end,:)...
        ]);
csc = caxis;
xlabel('Ensemble #');
ylabel('Bin #');
title('After Data Editing');

streamer([p.name,'  Figure 14']);