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']);