diff -r 000000000000 -r 0a450563f904 getbtrack.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/getbtrack.m	Tue Oct 20 16:23:49 2009 -0400
@@ -0,0 +1,316 @@
+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