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