function [ds,dr,ps,p]=getshear(d,p,ps,dr)

%

% - compute shear profiles 

% - use only central difference

% - use 2*std editing

% 

%  Martin Visbeck, LDEO, 3/7/97

% resolution of final shear profile in meter

idz=ps.dz;

disp(['GETSHEAR2: average shear profile over (ps.dz) ',num2str(idz),' [m]'])

% average over how much percent of the data

ps=setdefv(ps,'shear_stdf',2);

stdf=ps.shear_stdf;

disp([' maximum std (stdf) ',num2str(stdf),' of data '])

if nargin<4, dr.dummy=0; end

% check if only one istrument is to be used

if ps.up_dn_looker==2

 % down looker only

 d.weight(d.izu,:)=nan;

elseif ps.up_dn_looker==3

 % up looker only

 d.weight(d.izd,:)=nan;

end

% weightmin

ps=setdefv(ps,'shear_weightmin',0.1);

disp([' minimum weight  ',num2str(ps.shear_weightmin),' of data '])

is1=sum(sum(isfinite(d.weight)));

w=double(d.weight>ps.shear_weightmin);

ii=find(w==0);

w(ii)=NaN;

is2=sum(sum(isfinite(w)));

disp([' will use ',int2str(is2/is1*100),' % of data '])

% compute central shear

 % disp([' compute central diff shear '])

 iiok=1:length(d.z);

 us=[NaN*d.ru(1,iiok);diff2(d.ru(:,iiok))./diff2(d.izm);NaN*d.ru(1,iiok)].*w;

 vs=[NaN*d.rv(1,iiok);diff2(d.rv(:,iiok))./diff2(d.izm);NaN*d.rv(1,iiok)].*w;

 ws=[NaN*d.rw(1,iiok);diff2(d.rw(:,iiok))./diff2(d.izm);NaN*d.rw(1,iiok)].*w;

% loop over final depth profile

if existf(dr,'z')==1

 z=dr.z;

else

 izmax=-min(d.z);

 izmax=-maxnan(-[izmax,p.zbottom]);

 z=idz*0.5:idz:izmax;

end

ds.usm=zeros(length(z),1)*NaN;

ds.vsm=ds.usm;

ds.wsm=ds.usm;

ds.use=ds.usm;

ds.vse=ds.usm;

ds.wse=ds.usm;

ds.nn=ds.usm;

ds.z=ds.usm;

il=length(z);

iv=0;

%disp([' average shear estimates '])

for in=1:il

 i=z(in);

 iv =iv+1;

 % up and down

 i2=find(abs(-d.izm-i-idz/2) <= idz);

 i1=i2(find(isfinite(us(i2)+vs(i2))));

 ds.nn(iv)=length(i1);

 if ds.nn(iv) > 2

  usmm=median(us(i1)');

  ussd1=std(us(i1));

  ii=i1(find(abs(us(i1)-usmm)<stdf*ussd1)); 

  if length(ii)>1

   ds.usm(iv)=mean(us(ii));

   ds.use(iv)=std(us(ii));

  end

  vsmm=median(vs(i1)');

  vssd1=std(vs(i1));

  ii=i1(find(abs(vs(i1)-vsmm)<stdf*vssd1)); 

  if length(ii)>1

   ds.vsm(iv)=mean(vs(ii));

   ds.vse(iv)=std(vs(ii));

  end

  wsmm=median(ws(i1)');

  wssd1=std(ws(i1));

  ii=i1(find(abs(ws(i1)-wsmm)<stdf*wssd1)); 

  if length(ii)>1

   ds.wsm(iv)=mean(ws(ii));

   ds.wse(iv)=std(ws(ii));

  end

%  if exist('test')==1

%    hist(vs(ii),[-1:0.05:1]*.2e-1)

%    title(['z = ',num2str(z(iv))])

%    pause(0.01)

%  end

 end

% save depth vector

  ds.z(iv)=i;

% end of shear loop

end

% integrate shear profile (from bottom up)

ii=find(isnan(ds.usm));

ds.usm(ii)=0;

ii=find(isnan(ds.vsm));

ds.vsm(ii)=0;

ii=find(isnan(ds.wsm));

ds.wsm(ii)=0;

ds.ur=flipud(cumsum(flipud(ds.usm)))*idz;

ds.vr=flipud(cumsum(flipud(ds.vsm)))*idz;

ds.wr=flipud(cumsum(flipud(ds.wsm)))*idz;

ds.ur=ds.ur-mean(ds.ur);

ds.vr=ds.vr-mean(ds.vr);

ds.wr=ds.wr-mean(ds.wr);

if existf(d,'down')==1

 dz=2*abs(mean(diff(d.zd)));

 fac=1/tan(d.down.Beam_angle*pi/180)*sqrt(2)*dz;

 ds.ensemble_vel_err=ds.wse*fac;

 dr.ensemble_vel_err=ds.wse*fac;

end

if nargin>3

 dr.u_shear_method=ds.ur;

 dr.v_shear_method=ds.vr;

 dr.w_shear_method=ds.wr;

 uds=stdnan(dr.u-mean(dr.u)-ds.ur);

 vds=stdnan(dr.v-mean(dr.v)-ds.vr);

 uvds=sqrt(uds.^2+vds.^2);

 if uvds>meannan(dr.uerr)

  warn=(' increased error because of shear - inverse difference');

  disp(warn)

  if uvds>1.5*meannan(dr.uerr) & uvds>0.1

   p.warn(size(p.warn,1)+1,1:length(warn))=warn;

  end

  dr.uerr=dr.uerr/meannan(dr.uerr)*uvds/1.5;

 end

end

%--------------------------------------------------

function x = diff2(x,k,dn)

%DIFF2   Difference function.  If X is a vector [x(1) x(2) ... x(n)],

%       then DIFF(X) returns a vector of central differences between

%       elements [x(3)-x(1)  x(4)-x(2) ... x(n)-x(n-2)].  If X is a

%	matrix, the differences are calculated down each column:

%       DIFF(X) = X(3:n,:) - X(1:n-2,:).

%	DIFF(X,n) is the n'th difference function.

%	J.N. Little 8-30-85

%	Copyright (c) 1985, 1986 by the MathWorks, Inc.

if nargin < 2,	k = 1; end

if nargin < 3,	dn = 2; end

for i=1:k

	[m,n] = size(x);

	if m == 1

                x = x(1+dn:n) - x(1:n-dn);

	else

                x = x(1+dn:m,:) - x(1:m-dn,:);

	end

end