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

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

%

% - compute shear profiles 

% - use only central difference

% - use 2*std editing

% 

% version 8  last change 06.08.2019

%  Martin Visbeck, LDEO, 3/7/97

% some modifications and code cleanup                          GK, 16.05.2015  2-->3

% reinstated the use of d.weight for the identification of shear pairs

%                                                              GK, 03.12.2018  3-->4 

% reorganization of the dz handling                            GK, 12.12.2018  4-->5

% added error output, catch different shear methods            GK, 22.07.2019  5-->6

% more text output                                             

% added new parameter  ps.shear_throw_out_percent              GK, 05.08.2019  6-->7

% fixed bug in shear calculation                               GK, 06.08.2019  7-->8

%======================================================================

%                    C A L C _ S H E A R 3 . M 

%                    doc: Wed Sep  4 17:55:22 2019

%                    dlm: Wed Sep  4 18:55:36 2019

%                    (c) 2019 G. Krahmann, M. Visbeck

%                    uE-Info: 400 0 NIL 0 0 72 0 2 4 NIL ofnI

%======================================================================

% CHANGES BY ANT:

%   Sep  4, 2019: - replaced messages by LDEO warn mechanism

%

% general function info

%

disp(' ')

disp(['CALC_SHEAR3:  calculate a baroclinic velocity profile based on shears only'])

%

% decide whether central differences or simple differences are to be used in

% the shear calculation

%

% 1: simple differences

% 2: central differences

%

diff_type = 2;

%

% resolution of final shear profile in meter

%

dz = ps.dz;

shear_dz = dz * diff_type;

disp(['    Averaging shear profile over ',num2str(shear_dz),' m intervals'])

%

% Discard a certain amount of data as suspected outliers that

% for relatively small numbers of values will skew the mean.

% In the old calc_shear2.m the allowed range was the std of the shears 

% times the set factor stdf. Default was stdf=2 .

% In calc_shear3.m the outlier determination is iterative and thus a bit

% safer in calculation. The allowed range is set as a fraction of the whole

% population. For gaussian distributions stdf converts to this fraction.

% As usually stdf is not varied much, we have implemented a lookup list

% here.

%

use_new_outlier = 1;

ps.shear_stdf = 4;								% arbitrary guess

stdf = ps.shear_stdf;

if use_new_outlier==1

  fracs = 1 - [31.8,13.4,4.5,1.3,0.3]/100;		

  stdfs = [1,1.5,2,2.5,3];

  [dummy,ind] = min(abs(stdf-stdfs));

  frac = fracs(ind);

end

disp(['    Maximum allowed std within calculation intervals : ',num2str(stdf)])

disp(['    Data deviating more from the median will be discarded.'])

%

% 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

%

% Apply a weight threshold to shear data.

%

% There are two variants.

% First one is used when ps.shear_throw_out_percent is not NaN.

%   This one throws out the xx% shears with the lowest correlation-derived weights.

% Second one is used when ps.shear_weightmin is not NaN.

%   This one sets a minimum weight to keep the shears.

%

% First one is now default with 10%.

%

disp(['    Correlation-derived weights range from ',num2str(min(d.weight(:))),' to ',num2str(max(d.weight(:)))])

if 0 % ~isnan(ps.shear_throw_out_percent)

  local_weight = d.weight;

  ind = find(isfinite(local_weight));

  [dummy,ind2] = sort(local_weight(ind));

  if length(ind2)>9

    local_weight(ind(ind2(1:floor(length(ind2)/10)))) = nan;

  end

elseif 0 % ~isnan(ps.shear_weightmin)

  disp(['    Minimum weight  ',num2str(ps.shear_weightmin),' for data to be used in shear calc.'])

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

else

  disp('>   No weight criterion applied to raw shear data.')

  disp('>   You should set  ps.shear_throw_out_percent  or  ps.shear_weightmin')

end

local_weight = d.weight;

disp(['    Removed ',int2str(100-sum(isfinite(local_weight))/sum(isfinite(d.weight))*100),...

  ' % of data with lowest weights from shear calculation.'])

disp(['    New weights range from ',num2str(min(local_weight(:))),' to ',num2str(max(local_weight(:)))])

%

% convert the weights to 1 and NaN

%

%local_weight = replace( local_weight, local_weight<=0, nan );

%local_weight = replace( local_weight, local_weight>0, 1 );

local_weight(find(local_weight <= 0)) = nan;

local_weight(find(local_weight > 0))  = 1;

%

% compute shear

%

% two ways are offered here

% first:   central differences for the shears

% second:  single differences 

% the first is similar to the ways of the old calc_shear2.m

%

% central differences

if diff_type==2

  local_weight = [repmat(nan,1,size(local_weight,2));diff2(local_weight)+1;repmat(nan,1,size(local_weight,2))];

  ushear = [NaN*d.ru(1,:);diff2(d.ru(:,:))./diff2(d.izm);NaN*d.ru(1,:)].*local_weight;

  vshear = [NaN*d.rv(1,:);diff2(d.rv(:,:))./diff2(d.izm);NaN*d.rv(1,:)].*local_weight;

  wshear = [NaN*d.rw(1,:);diff2(d.rw(:,:))./diff2(d.izm);NaN*d.rw(1,:)].*local_weight;

  zshear = -d.izm;

% single differences

else

  ushear = diff( d.ru.*local_weight )./diff(d.izm);

  vshear = diff( d.rv.*local_weight )./diff(d.izm);

  wshear = diff( d.rw.*local_weight )./diff(d.izm);

  zshear = -(d.izm(1:end-1,:)+d.izm(2:end,:))/2;

end

ds.ushear = ushear;

ds.vshear = vshear;

ds.wshear = wshear;

ds.zshear = zshear;

%

% set depth levels

%

z = dr.z;

%

% prepare shear solution result vectors

%

ds.usm = repmat(nan,length(z),1);

ds.vsm = ds.usm;

ds.wsm = ds.usm;

ds.usmd = ds.usm;

ds.vsmd = ds.usm;

ds.use = ds.usm;

ds.vse = ds.usm;

ds.wse = ds.usm;

ds.nn = ds.usm;

ds.z = z;

%

% loop over depth levels and calculate the average shear at that level

%

% in the case of central differences this is oversampled here

% but by sticking with the same resolution it makes the results easier

% to work with

%

for n=[1:length(z)]

  i1 = find( ( abs( zshear - z(n) ) <= shear_dz/2 ) & isfinite( ushear + vshear ) );

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

  if ds.nn(n) > 2

    % two ways to select outliers

    % first:   select all that are beyond a fixed range around the median

    % second:  iteratively reject the worst (largest distance from mean)

    %          until a fixed fraction is rejected

    % the second is usually the safer calculation but is a bit slower

    if 1

      usmm = median( ushear(i1) );

      ussd1 = std( ushear(i1) );

      vsmm = median( vshear(i1) );

      vssd1 = std( vshear(i1) );

      wsmm = median( wshear(i1) );

      wssd1 = std( wshear(i1) );

      ii1 = i1( find(abs(ushear(i1)-usmm)<stdf*ussd1) );

      ii2 = i1( find(abs(vshear(i1)-vsmm)<stdf*vssd1) );

      ii3 = i1( find(abs(wshear(i1)-wsmm)<stdf*wssd1) );

    else

      [dummy,ii1] = meanoutlier(ushear(i1),frac);

      [dummy,ii2] = meanoutlier(vshear(i1),frac);

      [dummy,ii3] = meanoutlier(wshear(i1),frac);

      ii1 = i1(ii1);

      ii2 = i1(ii2);

      ii3 = i1(ii3);

    end

    % two ways of calculating the mean and std of the selected shears

    % first:  if there is a rejected one in any of u,v,w shears then use it

    %         for non of the calculations

    % second: if there is a rejected one in any of u,v,w shears then use it

    %         only in u,v or w calculations

    % the second one is the one used by the old calc_shear2.m

    % but to me this does not make sense, GK May 2015

    if 1

      dummy = zeros(size(ushear));

      dummy(ii1) = 1;

      dummy(ii2) = dummy(ii2)+1;

      dummy(ii3) = dummy(ii3)+1;

      ii = find(dummy==3);

      if length(ii)>1

        ds.usm(n) = mean(ushear(ii));

        ds.usmd(n) = median(ushear(ii));

        ds.use(n) = std(ushear(ii));

        ds.ii(n) = length(ii);

        ds.vsm(n) = mean(vshear(ii));

        ds.vsmd(n) = median(vshear(ii));

        ds.vse(n) = std(vshear(ii));

        ds.wsm(n) = mean(wshear(ii));

        ds.wsmd(n) = median(wshear(ii));

        ds.wse(n) = std(wshear(ii));

        % debugging plot 

        if 0

          figure(3)

          clf

          subplot(3,1,1)

          hist(ushear(ii),30)

          hold on

          ax = axis;

          plot([1,1]*mean(ushear(ii)),ax(3:4),'r')

          plot([1,1]*median(ushear(ii)),ax(3:4),'m')

          ind = find(dr.z==z(n));

          if ind<length(dr.u)

            plot([1,1]*(dr.u(ind-1)-dr.u(ind+1))/20,ax(3:4),'g')

          end

          title(int2str(z(n)))

          subplot(3,1,2)

          hist(vshear(ii),30)

          hold on

          ax = axis;

          plot([1,1]*mean(vshear(ii)),ax(3:4),'r')

          plot([1,1]*median(vshear(ii)),ax(3:4),'m')

          ind = find(dr.z==z(n));

          if ind<length(dr.u)

            plot([1,1]*(dr.v(ind-1)-dr.v(ind+1))/20,ax(3:4),'g')

          end

          subplot(3,1,3)

          hist( zshear(ii)-z(n), 30 )

          pause

        end 

      end

    else

      if length(ii1)>1

        ds.usm(n) = mean(ushear(ii1));

        ds.usmd(n) = median(ushear(ii1));

        ds.use(n) = std(ushear(ii1));

      end

      if length(ii2)>1

        ds.vsm(n) = mean(vshear(ii2));

        ds.vsmd(n) = median(vshear(ii2));

        ds.vse(n) = std(vshear(ii2));

      end

      if length(ii3)>1

        ds.wsm(n) = mean(wshear(ii3));

        ds.wsmd(n) = median(wshear(ii3));

        ds.wse(n) = nstd(wshear(ii3));

      end

    end

  end

end

%

% a debugging figure

%

if 0

sfigure(3);

clf

orient tall

subplot(1,2,1)

plot(ushear,zshear,'b.','markersize',3)

hold on

plot(ds.usm,ds.z,'r')

plot(ds.usmd,ds.z,'k')

inv_shear_u = -diff(dr.u-mean(dr.u))/dz;

plot(inv_shear_u,(z(1:end-1)+z(2:end))/2,'g')

set(gca,'ydir','reverse')

subplot(1,2,2)

plot(vshear,zshear,'b.','markersize',3)

hold on

plot(ds.vsm,ds.z,'r')

plot(ds.vsmd,ds.z,'k')

inv_shear_v = -diff(dr.v-mean(dr.v))/dz;

plot(inv_shear_v,(z(1:end-1)+z(2:end))/2,'g')

set(gca,'ydir','reverse')

sfigure(2)

end

%

% integrate shear profile (from bottom up)

%

%%ds.usm = replace( ds.usm, isnan(ds.usm), 0 );

%%ds.vsm = replace( ds.vsm, isnan(ds.vsm), 0 );

%%ds.wsm = replace( ds.wsm, isnan(ds.wsm), 0 );

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

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

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

if length(ind)/length(ds.usm)*100>5

  disp(['>   Found ',num2str(length(ind)/length(ds.usm)*100),...

    '% Nan in shear data. Integration result might be problematic.'])

end

if 1

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

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

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

else

  ds.ur = flipud(cumsum(flipud(ds.usmd)))*dz;

  ds.vr = flipud(cumsum(flipud(ds.vsmd)))*dz;

  ds.wr = flipud(cumsum(flipud(ds.wsmd)))*dz;

end

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

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

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

%

% This is a peculiar place for the single ping error estimate. But 

% as it is based on the variability in the data itself, it makes sense.

% The assumption is that there should be basically zero shear in the

% vertical velocities. At least it is so small as to be not detectable

% here. Thus any variability in the vertical shear is caused by the

% errors/noise of the measurement. Together with an angular conversion

% factor this gives an error/noise value for the horizontal velocities.

%

if ~isfield(d,'zd')

  dz2 = abs(mean(diff(d.z)));

else

  dz2 = diff_type*abs(mean(diff(d.zd)));

end

if isfield(d,'down')

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

else

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

end

ds.ensemble_vel_err = ds.wse*fac;

dr.ensemble_vel_err = ds.wse*fac;

%

% store result and give text output

%

dr.u_shear_method = ds.ur;

dr.v_shear_method = ds.vr;

dr.w_shear_method = ds.wr;

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

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

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

disp(['    Inversion average error : ',num2str(mean( dr.uerr ) ),' m/s'])

if uvds>mean(dr.uerr)*1.5

  error_increase_factor = 1/mean(dr.uerr)*uvds/1.5;

  warn = ('>   Increasing error estimate because of elevated shear - inverse difference');

  disp(warn)

  disp(['>   by a factor of ',num2str(error_increase_factor)])

  disp(['>   std of difference between regular and shear profile : ',num2str(uvds),' m/s'])

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

  dr.uerr = dr.uerr * error_increase_factor;

end

disp(['    Final average error : ',num2str(mean( dr.uerr ) ),' m/s'])

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

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