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