function [d,p]=getbtrack(d,p);

% function [d,p]=getbtrack(d,p);

%

% create own bottom track in addition to the one used before

%

% force own distance to RDI bottom track

p=setdefv(p,'bottomdist',0);

% "manual" bottom track for down looker

p=setdefv(p,'btrk_mode',3);

% p.btrk_ts is in dB to detect bottom above bin1 level (for own btm track)

p=setdefv(p,'btrk_ts',10);

% p.btrk_below gives bin used below target strength maximum

p=setdefv(p,'btrk_below',0.5);

% p.btrk_range gives range of allowed bottom track ranges

p=setdefv(p,'btrk_range',[300 50]);

% maximum allowed difference between reference layer W and W bottom track

p=setdefv(p,'btrk_wlim',0.05);

p=setdefv(p,'btrk_used',0);

disp('GETBTRACK creates own bottom track in addition to RDI')

% convert echo amplitude to relative target strength

% at=0.039; % attenuation dB/m for 150 kHz

% at=0.062; % attenuation dB/m for 300 kHz

if d.down.Frequency==300,

  p=setdefv(p,'ts_att_dn',0.06);

else

  p=setdefv(p,'ts_att_dn',0.039);

end

d.tg(d.izd,:)=targ(d.ts(d.izd,:),d.zd,p.ts_att_dn);

if length(d.izu)>0

 if d.up.Frequency==300,

  p=setdefv(p,'ts_att_up',0.06);

 else

  p=setdefv(p,'ts_att_up',0.039);

 end

 d.tg(d.izu,:)=targ(d.ts(d.izu,:),d.zu,p.ts_att_up);

end

% save bin number where down looker starts

ib1=d.izd(1)-1;

nbin=length(d.izd);

disp(['  in: p.btrk_mode ',int2str(p.btrk_mode),' and p.btrk_used ',int2str(p.btrk_used)])

if p.btrk_mode>=1

    if p.btrk_ts>0

      disp(' using increased bottom echo amplitudes to create bottom track')

      % don't use first bin

      fitb1=1;

      zd=d.zd(fitb1:end);

      ead=d.tg(d.izd(fitb1:end),:);

      % fit parabola to locate bottom

      [zmead,mead,imead]=localmax2(zd',ead);

      imead=imead+fitb1-1;

      dts=mead-ead(1,:);

      % decide which bin to use for bottom velocities

      dz=abs(diff(d.zd(1:2)));

      % imeadbv=round(imead+p.btrk_below);

      imeadbv=round((zmead-d.zd(1))/dz+1+p.btrk_below);

      if p.btrk_used==1

       if p.bottomdist==0

        % check RDI bottom track only if non zero

        ii=find(d.hbot<min(p.btrk_range));

        if length(ii)>0

         disp([' found ',int2str(length(ii)),' bottom depth below btrk_range ',...

            int2str(min(p.btrk_range))])

         d.bvel(ii,:)=nan;

         d.hbot(ii)=nan;

        end

        ii=find(d.hbot>max(p.btrk_range));

        if length(ii)>0

         disp([' found ',int2str(length(ii)),' bottom depth above btrk_range ',...

            int2str(max(p.btrk_range))])

         d.bvel(ii,:)=nan;

         d.hbot(ii)=nan;

        end

       end

       % save RDI bottom track

       d.bvel_rdi=d.bvel;

       d.hbot_rdi=d.hbot;

       ii=find(isfinite(d.bvel(:,1)+d.bvel(:,2)));

       if length(ii)<10, 

        disp(' found less than 10 RDI bottom track values, try own')

        p.btrk_used=0; 

       end

      end

      disp([' use ',num2str(p.btrk_below),...

         ' bins below maximum target strength for own bottom track velocity'])

      % locate acceptable bottom tracks (don't accept first/last bin)

      ii=find(dts>p.btrk_ts & ...

              zmead>min(p.btrk_range) & zmead<max(p.btrk_range) & ...

               imeadbv<(nbin-1) & imeadbv>fitb1);

      if length(ii)>0

      % save bottom distance data

       d.hbot_own=d.hbot+NaN;

       d.hbot_own(ii)=zmead(ii);

      % force bottom distance if RDI mode fails to report distance

       if p.bottomdist | p.btrk_mode==2 | p.btrk_used~=1

        d.hbot=d.hbot_own;

        disp([' created ',int2str(length(ii)),...

		' bottom distances'])

        if p.bottomdist, p.btrk_used = 12; end

       else

        disp([' created ',int2str(length(ii)),...

		' bottom distances keeping original'])

       end

       % make bottom velocity data

       bv=d.bvel+NaN;

       for j=1:length(ii)

          ji=ii(j);

          bv(ji,1)=medianan(d.ru(ib1+imeadbv(ji)+[-1,0,0,1],ji));

          bv(ji,2)=medianan(d.rv(ib1+imeadbv(ji)+[-1,0,0,1],ji));

          bv(ji,3)=medianan(d.rw(ib1+imeadbv(ji)+[-1,0,0,1],ji));

          bv(ji,4)=medianan(d.re(ib1+imeadbv(ji)+[-1,0,0,1],ji));

       end

       % check for W-bot 

       wref=medianan(d.rw(d.izd,:),2);

       ii=find(abs(wref'-bv(:,3))>p.btrk_wlim);

       disp([' removed ',int2str(length(ii)),...

           ' bottom track profiles W_btrk - W_ref difference > ',num2str(p.btrk_wlim)])

       bv(ii,:)=nan;

       % check for outlier

       bv=boutlier(bv,d.hbot_own,p);

       d.bvel_own=bv;

       ii=find(isfinite(bv(:,1)+bv(:,2)));

       if (p.btrk_used~=1 & p.btrk_used~=12) | p.btrk_mode==2

        p.btrk_used = 2;

        d.bvel=d.bvel_own;

        disp([' created ',int2str(length(ii)),...

		' bottom track data from normal velocities'])

       else

        disp([' created ',int2str(length(ii)),...

		' bottom track velocities keeping original'])

       end

      else

        if (p.btrk_used~=1 & p.btrk_used~=12)

         p.btrk_used = -1;

        end

        disp(' did not find any bottom echos to create own bottom track ')

      end

    else

      disp(' no valid own bottom track. Increase target strength difference ? ')

    end

else

 disp('force no bottom track data ')

 p.btrk_used = -1;

 d.bvel=d.bvel+nan;

 d.hbot=d.hbot+nan;

end

% summary output

disp([' out: p.btrk_mode ',int2str(p.btrk_mode),' and p.btrk_used ',int2str(p.btrk_used)])

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

function [bvel,p] = boutlier(bvel,hbot,p);

% function [bvel,p] = boutlier(bvel,hbot,p);

%

%

% input  : bvel bottom track velocity

%           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

n=p.outlier;

nblock=p.outlier_n;

dummyb = bvel*0;

if size(dummyb,2)~=4, disp(' no data '), return, end

si = size(dummyb);

sn = ceil(si(1)/nblock);

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

for i=1:length(n)

  for m=1:sn

    ind = (m-1)*nblock+[1:nblock];

    ii = find( ind<=si(1) );

    ind = ind(ii);

    % bottom track

     bvelt(ind,1)=bvel(ind,1)-medianan(bvel(ind,1));

     bu = find(abs(bvelt(ind,1))>n(i)*rms(bvelt(ind,1)));

     dummyb(ind(bu),:)=nan;

     bvelt(ind,2)=bvel(ind,2)-medianan(bvel(ind,2));

     bv = find(abs(bvelt(ind,2))>n(i)*rms(bvelt(ind,2)));

     dummyb(ind(bv),:)=nan;

     bw = find(abs(bvel(ind,3))>n(i)*rms(bvel(ind,3)));

     dummyb(ind(bw),:)=nan;

     hbot(ind)=hbot(ind)-medianan(hbot(ind));

     bh = find(abs(hbot(ind))>n(i)*rms(hbot(ind)));

     dummyb(ind(bh),:)=nan;

  end

  bvel = bvel + dummyb; 

  hbot = hbot + dummyb(:,1)'; 

end  

disp([' boutlier removed ',int2str(length(find(isnan(dummyb)))-lob),...

      ' bottom track velocities '])

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

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

function [ts,bcs]=targ(ea,dis,at,bl,eas,ap)

% function [ts,bcs]=targ(ea,dis,at,bl,eas,ap)

% Target strength of EA for volume scatterer

% ea = echoamp in  dB

% dis = distance in  m

% at = attenuation dB/m

% bl = pulse/bin legth in  m

% eas = source level

% ap = aperature in degree

% M. Visbeck 2004

% make distance matrix if needed

[lr,lc]=size(ea);

if size(dis,2)==1 | size(dis,1)==1

 dis=dis(:);

 dis=repmat(dis,[1,lc]);

end

if nargin<3

 at=0.039; % attenuation dB/m for 150 kHz

 at=0.062; % attenuation dB/m for 300 kHz

end

%binlength

if nargin<4 , bl=median(abs(diff(dis(:,1)))); end

% source level in dB

if nargin <5, eas=100; end

% beam aperature in DEGREE convert to radian

if nargin <6, ap=2; end

al=ap*pi/180; 

% radius of top and bottom of each bin

r1=tan(al)*(dis-bl/2);

r2=tan(al)*(dis+bl/2);

% ensonified volume 

v=pi*bl/3*(r1.^2+r2.^2+r1.*r2);

% transmission loss

tl=20*log10(dis)+at*dis;

% target strength

ts=ea-eas+2*tl-10*log10(v);

if nargout>1

 %backscatter cross section

 bcs=exp(ts./10);

end