function [d,p] = outlier(d,p);
% function [d,p] = outlier(d,p);
%
% check for spurious signals in the data
% a similar pinging frequency (e.g. pinger or hydrosweep)
% noise at the end of the beams etc.
%
% input  :	d		LADCP analysis data structure
%           n       factor for outlier rejection
%                   n(1)*rms(data) first sweep
%                   n(2)*rms(data) second sweep etc.
%
% output :	d		changed LADCP analysis data structure
%
% version 0.1	last change 27.6.2000

% G.Krahmann, LDEO Jun 2000
p=setdefv(p,'outlier',[4.0 3.0]);

% mark outlier for blocks of 5 minute duration
nblock=ceil(5./meannan(diff(d.time_jul)*24*60));
p=setdefv(p,'outlier_n',nblock);

n=p.outlier;
ii=find(n<0.01);
if length(ii)>1
 disp('OUTLIER:  I will not delete ALL data ')
 return
end
ii=find(n<1);
if length(ii)>1
 disp('OUTLIER: Are you sure you want to delete most data? ')
end

nblock=p.outlier_n;

rw = d.rw(d.izd,:);
rv = d.rv(d.izd,:);
ru = d.ru(d.izd,:);
if existf(d,'ts')==1
 rt = d.ts(d.izd,:);
end
dummy = rw*0;

bvel=d.bvel;
dummyb = bvel*0;
if size(dummyb,2)==4, ibvel=1; else, ibvel=0; end

si = size(dummy);
sn = ceil(si(2)/nblock);


lob=length(find(isnan(dummyb)));
lo=length(find(isnan(dummy)));

for i=1:length(n)
  % calculate anomaly fields
  rwm = medianan(rw);
  rw = rw - ones(size(rw,1),1)*rwm;
  ru = ru - ones(size(ru,1),1)*medianan(ru);
  rv = rv - ones(size(rv,1),1)*medianan(rv);
  if existf(d,'ts')
   rt = rt - ones(size(rt,1),1)*medianan(rt);
  end
  if ibvel, bvel(:,3)=bvel(:,3)-rwm'; end
  for m=1:sn
    ind = (m-1)*nblock+[1:nblock];
    ii = find( ind<=si(2) );
    ind = ind(ii);
    dummy2 = dummy(:,ind);
    rrw = rw(:,ind);
    badrw = find(abs(rrw)>n(i)*rms(rrw(:)));
    rru = ru(:,ind);
    badru = find(abs(rru)>n(i)*rms(rru(:)));
    rrv = rv(:,ind);
    badrv = find(abs(rrv)>n(i)*rms(rrv(:)));
    dummy2(badrw) = nan;
    dummy2(badru) = nan;
    dummy2(badrv) = nan;
    if existf(d,'ts')
     rrt = rt(:,ind);
     badrt = find(abs(rrt)>n(i)*rms(rrt(:)));
     dummy2(badrt) = nan;
    end
    dummy(:,ind) = dummy2;
    if ibvel
     % bottom track
     bvel(ind,1)=bvel(ind,1)-medianan(bvel(ind,1));
     bu = find(abs(bvel(ind,1))>n(i)*rms(bvel(ind,1)));
     dummyb(ind(bu),:)=nan;
     bvel(ind,2)=bvel(ind,2)-medianan(bvel(ind,2));
     bv = find(abs(bvel(ind,2))>n(i)*rms(bvel(ind,2)));
     dummyb(ind(bv),:)=nan;
     bw = find(abs(bvel(ind,3))>n(i)*rms(bvel(ind,3)));
     dummyb(ind(bw),:)=nan;
     hbot(ind)=d.hbot(ind)-medianan(d.hbot(ind));
     bh = find(abs(hbot(ind))>n(i)*rms(hbot(ind)));
     dummyb(ind(bh),:)=nan;
     if p.debug>1
      disp(['bottom track bad U ',int2str(length(bu))])
      disp(['bottom track bad V ',int2str(length(bv))])
      disp(['bottom track bad W ',int2str(length(bw))])
      disp(['bottom track bad H ',int2str(length(bh))])
     end
    end
  end
  rw = rw+dummy;
  rv = rv+dummy;
  ru = ru+dummy;
  if ibvel, d.bvel = d.bvel + dummyb; d.hbot = d.hbot + dummyb(:,1)'; end
end  


d.weight(d.izd,:)=d.weight(d.izd,:)+dummy;
d.rw(d.izd,:) = d.rw(d.izd,:)+dummy;
d.ru(d.izd,:) = d.ru(d.izd,:)+dummy;
d.rv(d.izd,:) = d.rv(d.izd,:)+dummy;

disp([' Outlier discarded ',int2str(length(find(isnan(dummy)))-lo),' bins down looking'])
if ibvel
disp([' Outlier discarded ',int2str(length(find(isnan(dummyb)))-lob),' bottom track'])
end

if size(d.izd)~=size(d.rw,1)

rw = d.rw(d.izu,:);
rv = d.rv(d.izu,:);
ru = d.ru(d.izu,:);
dummy = rw*0;
si = size(dummy);
sn = ceil(si(2)/nblock);

lo=length(find(isnan(dummy)));


for i=1:length(n)
  % calculate anomaly fields
  rw = rw - ones(size(rw,1),1)*medianan(rw);
  ru = ru - ones(size(ru,1),1)*medianan(ru);
  rv = rv - ones(size(rv,1),1)*medianan(rv);
  for m=1:sn
    ind = (m-1)*nblock+[1:nblock];
    ii = find( ind<=si(2) );
    ind = ind(ii);
    dummy2 = dummy(:,ind);
    rrw = rw(:,ind);
    badrw = find(abs(rrw)>n(i)*rms(rrw(:)));
    rru = ru(:,ind);
    badru = find(abs(rru)>n(i)*rms(rru(:)));
    rrv = rv(:,ind);
    badrv = find(abs(rrv)>n(i)*rms(rrv(:)));
    dummy2(badrw) = nan;
    dummy2(badru) = nan;
    dummy2(badrv) = nan;
    dummy(:,ind) = dummy2;
  end
  rw = rw+dummy;
  rv = rv+dummy;
  ru = ru+dummy;
end  

d.weight(d.izu,:)=d.weight(d.izu,:)+dummy;
d.rw(d.izu,:) = d.rw(d.izu,:)+dummy;
d.ru(d.izu,:) = d.ru(d.izu,:)+dummy;
d.rv(d.izu,:) = d.rv(d.izu,:)+dummy;

disp([' Outlier discarded ',int2str(length(find(isnan(dummy)))-lo),' bins up looking'])


end

return

%
function y = rms(x,dim)
%RMS    root-mean square. For vectors, RMS(x) returns the standard
%       deviation.  For matrices, RMS(X) is a row vector containing
%       the root-mean-square of each column. The difference to STD is
%       that here the mean is NOT removed.  
%       RMS(X,DIM) returns the root-mean-square of dimension DIM
%
%	See also STD,COV.
%

%       Uwe Send, IfM Kiel, Apr 1992
% added NaN handling   Gerd Krahmann, IfM Kiel, Oct 1993, Jun 1994
% removed bug in NaN handling   G.Krahmann, Aug 1994
% added compatibility to MATLAB 5	G.Krahmann, LODYC Paris, Jul 1997


  if nargin<2
    dim=min(find(size(x)>1));
  end

  if all(isnan(x))
    y=nan;
    return
  end    

  x = shiftdim(x,dim-1);
  s = size(x);
  so = s(1);
  s(1) = 1;

  for n = 1:prod(s)
    good = find(~isnan(x(:,n)));
    if ~isempty(good)
      y(1,n) = norm( x(good,n) ) / sqrt(length(good));
    else
      y(1,n) = NaN;
    end
  end
  y = reshape(y,s);
  y = shiftdim( y, ndims(x)-dim+1 );