function [spikeindex, okayindex] = harv_de_spike(P911_cind, S911_cind);
% function [spikeindex, okayindex] = harv_de_spike(P911_cind, S911_cind);
%
% Finds spikes in salinity data by computing a running mean and checking for outliers,
% interpolates flagged values in density, and then backs out the salinity.
%
% Written by Miles Sundermeyer, 7/19/96

global column

reclength = size(S911_cind,1);
avgpts = 20;		% set the number of points used as to compute incremental mean
outlier = 2.;		% set the number of standard deviations to use as spike criterium
%saldiffthresh = 0.075;	% threshold value for a point to be considered a salin spike (for w)
saldiffthresh = 0.1;  % value used for u,v velocities
% Set up to iteratively weed out salinity spikes
% figure out how many columns to give our work arrays (may add one in a sec)
arryflr = avgpts * floor(reclength/avgpts);
% make work arrays so that we can get incremental averages
Salwork = reshape(S911_cind(1:arryflr),avgpts,arryflr/avgpts);
Presswork = reshape(P911_cind(1:arryflr),avgpts,arryflr/avgpts);
indexwork = reshape([1:arryflr],avgpts,arryflr/avgpts);
if ~(arryflr == reclength)	% if we have left points off, add another column
   Salwork(:,size(Salwork,2)+1) = S911_cind(reclength-avgpts+1:reclength);
   Presswork(:,size(Presswork,2)+1) = P911_cind(reclength-avgpts+1:reclength);
   indexwork(:,size(indexwork,2)+1) = (reclength-avgpts+1:reclength)';
end
flagwork=ones(size(Salwork));	% make a flag array to keep track of salinity spikes (0)

% ** Go through a loop to iteratively weed out spikes in salinity
lastpass=1;	% assign this flag to have some non-empty value
rempass=1;
while(~isempty(lastpass))
  if(1)		% SHUT OFF THIS SECTION OF CODE
   % de-trend the work array
   if rempass==1
      % first time through, just detrend and go
      Salwork = detrend(Salwork);
   else
      % successive passes, try to reduce the bias due to the already flagged points
      Salmean = sum(Salwork.*flagwork)./sum(flagwork);
      Salwork = detrend(Salwork.*flagwork + (ones(avgpts,1)*Salmean).*[~flagwork]);
   end
  end		% END SHUT OFF SECTION OF CODE
   % compute means and std deviations of the columns of this Salwork array
   Salmean = sum(Salwork.*flagwork)./sum(flagwork);
   Salstd = sqrt( sum(((Salwork-ones(avgpts,1)*Salmean).*flagwork).^2) ./ ...
						(sum(flagwork)-1) );
   Pressmean = mean(Presswork);

   % flag points where salinity exceeds some number of std deviations from mean
   % AND the "spike" is a jump of more than some threshold value
   lastpass=find( flagwork.*abs(Salwork-ones(avgpts,1)*Salmean) > ...
					abs(outlier*ones(avgpts,1)*Salstd) & ...
   				flagwork.*abs(Salwork-ones(avgpts,1)*Salmean) > ...
					saldiffthresh );

   % set the flag array to reflect the new points we found
   flagwork(lastpass)=zeros(size(lastpass));

   disp([' *** done pass ',num2str(rempass),' (pause) ***']);
   rempass=rempass+1;
   %pause
end

% Keep track of where the spikes were
spikeindex = indexwork(find(flagwork==0));
okayindex = indexwork(find(flagwork==1));


