prepinv.m

%======================================================================
%                    P R E P I N V . M 
%                    doc: Wed Jan  7 16:46:29 2009
%                    dlm: Wed Sep  4 17:03:23 2019
%                    (c) 2009 A.M. Thurnherr
%                    uE-Info: 18 93 NIL 0 0 72 0 2 8 NIL ofnI
%======================================================================

% CHANGES BY ANT:
%   Jan  7, 2009: - tightened use of exist()
%   Oct 27, 2009: - color-coded down/upcast in fig. 6
%   Mar  4, 2010: - BUG: special code for shallow casts (no superensembles)
%			 was broken => removed
%   		  - BUG: down/upcast separation bombed wen no CTD data were
%			 available
%   Jun  9, 2014: - improved messages
%   Apr  1, 2016: - cosmetics
%   Sep  4, 2019: - BUG: superens std was calculated without removing means! (reported by GK)

function [di,p,d]=prepinv(d,p,dr)
% function [di,p,d]=prepinv(d,p,dr)
% LADCP-2 processing software version 3.0
% prepare for inverse slolve
%
% - average velocities within depth range
% 
%  Martin Visbeck, LDEO, 6/10/99

% average over all ensembles in the depth range avdz 
%  set default avdz to one bin length
p=setdefv(p,'avdz',medianan(abs(diff(d.izm(:,1)))));
p=setdefv(p,'oversample',1.0);
p=setdefv(p,'avpercent',100);
p=setdefv(p,'rotup2down',0);
p=setdefv(p,'offsetup2down',0);
p=setdefv(p,'tilt_weight',10);
p=setdefv(p,'avens',NaN);
p=setdefv(p,'soundcorr',1);
p=setdefv(p,'nav_error',30);
d=setdefv(d,'slon',d.time_jul*NaN);
d=setdefv(d,'slat',d.time_jul*NaN);
p=setdefv(p,'superens_std_min',d.down.Single_Ping_Err/sqrt(d.down.Pings_per_Ensemble));

disp('PREPINV: prepare data for inversion, form Super-Ensembles')
if isfinite(p.avens)
 disp([' average profiles over (p.avens) ',num2str(p.avens),' ensembles'])
else
 disp([' average profiles over (p.avdz) ',num2str(p.avdz),' meter'])
end

% remove bottom "bad" track values 
% This can be because of range 
if isfinite(p.zbottom)
 p=setdefv(p,'btrk_range',[300 50]);
 % backward compatible
 if length(p.btrk_range)<2, p.btrk_range(2)=0; end

 ii=find(abs(d.z+p.zbottom)>max(p.btrk_range) & isfinite(d.bvel(:,1)') );
 if length(ii)>0
  disp([' discarded ',int2str(length(ii)),...
    ' bottom tracks velocities because of height above bottom > ',int2str(max(p.btrk_range))])
  d.bvel(ii,:)=NaN;
  d.bvels(ii,:)=NaN;
 end

 ii=find(abs(d.z+p.zbottom)<min(p.btrk_range) & isfinite(d.bvel(:,1)') );
 if length(ii)>0
  disp([' discarded ',int2str(length(ii)),...
    ' bottom tracks velocities because of height above bottom < ',int2str(min(p.btrk_range))])
  d.bvel(ii,:)=NaN;
  d.bvels(ii,:)=NaN;
 end

 ii=find(abs(d.z+p.zbottom-d.hbot)> 100 );
 if length(ii)>0
  disp([' discarded ',int2str(length(ii)),...
    ' bottom distances because of depth difference > 100'])
  d.hbot(ii)=NaN;
 end
else
 d.bvel(:,:)=NaN;
 d.bvels(:,:)=NaN;
end


% reduce weight for large tilts
if p.tilt_weight>0 & existf(d,'tilt_weight')~=1
 fac=1.-tanh(d.tilt/p.tilt_weight)/2;
 d.weight=d.weight.*meshgrid(fac,d.weight(:,1));
 disp([' reduce weight for larger tilts 0.5 at ',...
        num2str(p.tilt_weight),' degree'])
 d.tilt_weight=p.tilt_weight;
end

% prepare for heading averaging
i=sqrt(-1);
u1d=exp(-i*(d.hdg(1,:))*pi/180); 
if length(d.izu)>1
 u1u=exp(-i*(d.hdg(2,:))*pi/180); 
 % get mean heading offset from COMPASS comparison
 hoff=compoff(u1d,u1u);
 % hoff=angle(u1d/u1u)*180/pi;
 p.up_dn_comp_off=hoff;
 disp([' mean heading offset from compasses = ',num2str(hoff),' deg'])
 
 % check tilt sensors
 % rotate by compass offset
  diary off
  hoff2=fminsearch('checktilt',0,[],[d.rol(2,:);d.pit(2,:);d.rol(1,:);d.pit(1,:)]);
  diary on
  disp([' mean heading offset from pitch/roll = ',num2str(hoff2),' deg'])
  p.up_dn_pit_rol_comp_off=hoff2;
  [d.rol(3,:),d.pit(3,:)]=uvrot(d.rol(2,:),d.pit(2,:),-hoff2);
  p.up_dn_rol_off=mean(d.rol(1,:)-d.rol(3,:));
  p.up_dn_pit_off=mean(d.pit(1,:)-d.pit(3,:));
 
% plot compass and tilt meter differences
  figure(6)
  clf
  orient tall

% find downcast/upcast separation
  [btm,btmi] = min(d.z);

% plot compass comparison between up and down instrument
  da=-angle(u1d)+angle(u1u*exp(i*hoff*pi/180));
  ii=find(da>pi);
  da(ii)=da(ii)-2*pi;
  ii=find(da<-pi);
  da(ii)=da(ii)+2*pi;
  subplot(311)
  plot(d.hdg(1,1:btmi),da(1:btmi)*180/pi,'r.')
  hold on
  plot(d.hdg(1,btmi:end),da(btmi:end)*180/pi,'b.')
  ylabel('Heading difference up-down')
  title([' Heading offset : ',num2str(hoff),' (r/b: down-/upcast)'])
  xlabel('Heading down')
  grid
  dhmax=sort(abs(da*180/pi));
  dhmax=dhmax(fix(end*0.95));
  if dhmax > 15 & ~existf(d,'hrot')
   warn=[' Large compass deviation: ',num2str(dhmax)];
   p.warn(size(p.warn,1)+1,1:length(warn))=warn;
  end
  dhmax=max([10 dhmax*1.3]); 
  axis([0 360 [-1 1]*dhmax])

% plot tilt meter difference
  tiltmax=p.tiltmax(1);
  subplot(312)
  pa=d.pit(1,:)-d.pit(3,:)-p.up_dn_pit_off;
  plot(d.pit(1,1:btmi),pa(1:btmi),'r.')
  hold on
  plot(d.pit(1,btmi:end),pa(btmi:end),'b.')
  ylabel('Pitch difference up-down')
  xlabel('Pitch of down instrument')
  title([' Pitch offset : ',num2str(p.up_dn_pit_off)])
  grid;
  axis([[-1 1]*tiltmax [-1 1]*5])
 
  subplot(313)
  ra=d.rol(1,:)-d.rol(3,:)-p.up_dn_rol_off;
  plot(d.rol(1,1:btmi),ra(1:btmi),'r.')
  hold on
  plot(d.rol(1,btmi:end),ra(btmi:end),'b.')
  ylabel('Roll difference up-down')
  xlabel('Roll of down instrument')
  title(['Roll offset : ',num2str(p.up_dn_rol_off)])
  grid;
  axis([[-1 1]*tiltmax [-1 1]*5])
 
  streamer([p.name,' Figure 6']);
  pause(0.01)
end


% offset upward looking ADCP to downward looking ADCP
if (p.offsetup2down~=0 & length(d.zd)~=length(d.ru(:,1)) & exist('dr','var'))
  if p.rotup2down==2
   % will not rotate to match velocities and correct offset
   p.rotup2down=1;
  end
  i=sqrt(-1);
  disp(' remove first guess ocean velocity from raw data')
  % OCEAN velocity
  [ib,it]=size(d.ru);
  z=-d.izm+d.ru*0;
  ii=find(z>=min(dr.z) & z<=max(dr.z));
  l.uoce(:,1)=interp1q(dr.z,dr.u,z(ii));
  l.uoce(:,2)=interp1q(dr.z,dr.v,z(ii));
  [prof,bin]=meshgrid(1:it,1:ib);
  l.ru=full(sparse(bin(ii),prof(ii),l.uoce(:,1)));
  l.rv=full(sparse(bin(ii),prof(ii),l.uoce(:,2)));
  l.ru(ib,it)=NaN;
  l.rv(ib,it)=NaN;
  ii=find(~(z>=min(dr.z) & z<=max(dr.z)));
  l.ru(ii)=NaN;
  l.rv(ii)=NaN;

  uu=medianan(d.ru(d.izu,:)+i*d.rv(d.izu,:)...
             -l.ru(d.izu,:)-i*l.rv(d.izu,:)+d.weight(d.izu,:)*0,2);
  ud=medianan(d.ru(d.izd,:)+i*d.rv(d.izd,:)...
             -l.ru(d.izd,:)-i*l.rv(d.izd,:)+d.weight(d.izd,:)*0,2);
  clear l

 ii=find(~isfinite(uu+ud));
 uu(ii)=0;
 ud(ii)=0;
 uoff=(ud-uu)*p.offsetup2down; 
 uoffm(d.izu,:)=meshgrid(uoff/2,d.izu);
 uoffm(d.izd,:)=meshgrid(-uoff/2,d.izd);
 d.ru=d.ru+real(uoffm); 
 d.rv=d.rv+imag(uoffm); 
 if existf(d,'bvel')
  d.bvel(:,1)=d.bvel(:,1)+real(-uoff/2)';
  d.bvel(:,2)=d.bvel(:,2)+imag(-uoff/2)';
 end
 if existf(d,'uoff')
  d.uoff=d.uoff+uoff;
 else
  d.uoff=uoff;
 end
 % estimate tilt error that could explain difference
 % beam velocity ~ W
 wbeam=(abs(d.wm)*cos(d.down.Beam_angle*pi/180));
 ii=find(wbeam<0.4);
 wbeam(ii)=NaN;
 % uerror for 1 degree tilt error
 du=wbeam/sin((d.down.Beam_angle+0.5)*pi/180)-...
    wbeam/sin((d.down.Beam_angle-0.5)*pi/180);
 % linear projection to estimate tilt error
 d.tilterr=abs(d.uoff)./abs(du);

 figure(10), clf
 orient tall
 subplot(311)
 plot(real(d.uoff))
 uoffav=boxav(d.uoff',20);
 hold on
 plot([0.5:length(uoffav)]*20,real(uoffav),'-r')
 grid
 axis tight
 ax=axis;
 ax(3:4)=[-1 1]*max(medianan(abs(real(uoffav)))*6,0.1);
 axis(ax)

 title('U offset [m/s]')
 subplot(312)
 plot(imag(d.uoff))
 hold on
 plot([0.5:length(uoffav)]*20,imag(uoffav),'-r')
 grid
 axis tight
 ax=axis;
 ax(3:4)=[-1 1]*max(medianan(abs(imag(uoffav)))*6,0.1);
 axis(ax)
 title('V offset [m/s]')

 subplot(313)
 plot(d.tilterr)
 hold on
 tilterrav=boxav(d.tilterr',20);
 plot([0.5:length(tilterrav)]*20,tilterrav,'-r')
 axis tight
 ax=axis;
 ax(3:4)=[0 1]*medianan(d.tilterr)*6;
 if isfinite(sum(ax))
  ax(4)=max(3,ax(4));
  axis(ax)
 end
 grid
 title('Tilt error [degree] consistent with offset')
  
 streamer([p.name,' Figure 10']);
 pause(0.01)

 disp(' adjusted for velocity offset in up and down looking ADCP')

end

if existf(p,'drot')
  if ~isfinite(p.drot)
   warn=[' magnetic deviation given is NAN '];
   p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
  end
else
   warn=[' NO magnetic deviation given '];
   p.warnp(size(p.warnp,1)+1,1:length(warn))=warn;
end
  

 % check if rotated already

% rotate upward looking ADCP to downward looking ADCP
if (p.rotup2down~=0 & length(d.zd)~=length(d.ru(:,1)))

  da=angle(u1d)-angle(u1u*exp(-i*hoff*pi/180));
  ii=find(da>pi); da(ii)=da(ii)-2*pi; ii=find(da<-pi); da(ii)=da(ii)+2*pi;
  d.diff_hdg=da;

  u1uc=exp(-i*(d.hdg(2,:)-hoff)*pi/180);
  hrotcomp=angle(u1uc./u1d)*180/pi;
  d.rot_comp=hrotcomp;

  % get mean heading offset from velocities
  if exist('dr','var')
   disp(' remove first guess ocean velocity from raw data')
   % OCEAN velocity
   [ib,it]=size(d.ru);
   z=-d.izm+d.ru*0;
   ii=find(z>=min(dr.z) & z<=max(dr.z));
   l.uoce(:,1)=interp1q(dr.z,dr.u,z(ii));
   l.uoce(:,2)=interp1q(dr.z,dr.v,z(ii));
   [prof,bin]=meshgrid(1:it,1:ib);
   l.ru=full(sparse(bin(ii),prof(ii),l.uoce(:,1)));
   l.rv=full(sparse(bin(ii),prof(ii),l.uoce(:,2)));
   l.ru(ib,it)=NaN;
   l.rv(ib,it)=NaN;
   ii=find(~(z>=min(dr.z) & z<=max(dr.z)));
   l.ru(ii)=NaN;
   l.rv(ii)=NaN;

   uu=medianan(d.ru(d.izu,:)+i*d.rv(d.izu,:)...
              -l.ru(d.izu,:)-i*l.rv(d.izu,:)+d.weight(d.izu,:)*0,2);
   ud=medianan(d.ru(d.izd,:)+i*d.rv(d.izd,:)...
              -l.ru(d.izd,:)-i*l.rv(d.izd,:)+d.weight(d.izd,:)*0,2);
   clear l

  else
   iz=1:4;
   uu=medianan(d.ru(d.izu(iz),:)+i*d.rv(d.izu(iz),:)+d.weight(d.izu(iz),:),1);
   ud=medianan(d.ru(d.izd(iz),:)+i*d.rv(d.izd(iz),:)+d.weight(d.izd(iz),:),1);
  end
%  try to take speed into account
  hrotvel=angle(uu./ud)*180/pi;
  d.rot_vel=hrotvel;

%  decide which to use 
  if p.rotup2down==1
   hrot=hrotcomp;
   disp(' rot up2down use mean up/down compass')
  elseif p.rotup2down==2
   hrot=hrotvel;
   disp(' rot up2down use velocities')
  elseif p.rotup2down==3
   hrot=hrotcomp;
   disp(' rot up2down use down compass')
  elseif p.rotup2down==4
   hrot=hrotcomp;
   disp(' rot up2down use up compass')
  else
   hrot=d.ru(1,:)*0;
   disp(' not sure what you want rot=0')
  end  

% plot what I am about to do
  figure(5)
  clf
  orient tall
  plot([hrotcomp-90;hrotvel-180;hrot]')
  hold on
  ax=axis;
  plot(ax(1:2),[0 0],'-k',ax(1:2),[0 0]-90,'-k',ax(1:2),[0 0]-180,'-k') 
  text(1,20,'rotation used ')
  text(1,-60,'rotation compass')
  text(1,-150,'rotation velocity')
  ax(3)=-250;
  axis(ax)
  ylabel(' heading correction for uplooking')
  streamer([p.name,' Figure 5']);
  pause(.01)

 if existf(d,'hrot')
  disp(' rotated earlier, use difference ')
  oldrot=d.hrot;
 else
  oldrot=0;
 end

 % apply rotation
  hrotm=meshgrid(hrot-oldrot,d.izd);
   % dont rotate NaNs
  ii=find(~isfinite(hrotm));
  hrotm(ii)=meannan(hrot-oldrot);
 % rotate both by half difference
  if p.rotup2down==4
   [ru,rv]=uvrot(d.ru(d.izd,:),d.rv(d.izd,:),-hrotm);
   d.ru(d.izd,:)=ru;
   d.rv(d.izd,:)=rv;
   [bu,bv]=uvrot(d.bvel(:,1),d.bvel(:,2),-(hrot-oldrot)');
   d.bvel(:,1)=bu;
   d.bvel(:,2)=bv;
  elseif p.rotup2down~=3
   [ru,rv]=uvrot(d.ru(d.izd,:),d.rv(d.izd,:),-hrotm/2);
   d.ru(d.izd,:)=ru;
   d.rv(d.izd,:)=rv;
   [bu,bv]=uvrot(d.bvel(:,1),d.bvel(:,2),-(hrot-oldrot)'/2);
   d.bvel(:,1)=bu;
   d.bvel(:,2)=bv;
  end
 
  hrotm=meshgrid(hrot-oldrot,d.izu);
   % dont rotate NaNs
  ii=find(~isfinite(hrotm));
  hrotm(ii)=meannan(hrot);
 % rotate both by half difference
  if p.rotup2down==4
   [ru,rv]=uvrot(d.ru(d.izu,:),d.rv(d.izu,:),hrotm);
   d.ru(d.izu,:)=ru;
   d.rv(d.izu,:)=rv;
  elseif p.rotup2down~=3
   [ru,rv]=uvrot(d.ru(d.izu,:),d.rv(d.izu,:),hrotm/2);
   d.ru(d.izu,:)=ru;
   d.rv(d.izu,:)=rv;
  end

  d.hrot=hrot;

end


% sound speed correction
if p.soundcorr==1
   if existf(d,'p')~=1
     d.p=press(abs(d.z));
   end
   if existf(d,'ss')~=1
    disp(' make soundspeed based on pressure and ADCP temp')
    d.ss=sounds(d.p,d.temp(1,:),34.5);
   end
   if d.soundc==0
    disp(' correct velocities for sound speed ')
    sc=meshgrid(d.ss./d.sv(1,:),d.izd);
    d.ru(d.izd,:)=d.ru(d.izd,:).*sc;
    d.rv(d.izd,:)=d.rv(d.izd,:).*sc;
    d.rw(d.izd,:)=d.rw(d.izd,:).*sc;
    if length(d.zd)~=length(d.ru(:,1))
     sc=meshgrid(d.ss./d.sv(1,:),d.izu);
     d.ru(d.izu,:)=d.ru(d.izu,:).*sc;
     d.rv(d.izu,:)=d.rv(d.izu,:).*sc;
     d.rw(d.izu,:)=d.rw(d.izu,:).*sc;
    end
    d.soundc=1;
   else
    disp(' will not correct for sound speed twice')
   end
end
 

if (isnan(p.avdz) | p.avdz<=0 ) & (isnan(p.avens) | p.avens < 2)
 disp(' avdz=NAN  => No pre-averaging done !!!')
 di.ru=d.ru;
 di.rv=d.rv;
 di.ruvs=d.ru*0+p.superens_std_min;
 di.rw=d.rw;
 di.re=d.re;
 di.ts=d.ts;
 di.tg=d.tg;
 di.weight=d.weight;
 di.bvel=d.bvel';
 di.hbot=d.hbot;
% make up std
 di.bvels=d.bvel'*0+p.single_ping_accuracy;
 di.hdg=d.hdg;
 di.pit=d.pit;
 di.rol=d.rol;
 di.temp=d.temp;
 di.tsd=d.ts(d.izd(2),:);
 di.tsd_out=d.ts(d.izd(end),:);
 di.dtiv=d.z*0+1;
 di.time_jul=d.time_jul;
 di.z=d.z;
 di.izm=d.izm;
 di.slat=d.slat;
 di.slon=d.slon;
 di.izd=d.izd;
 di.izu=d.izu;

else

% default reference velocity bin range bin 2 + 3
 if length(d.izd)>2
   izr=[min(d.izd)+[1:2]];
 end
 if length(d.izu)>2
  izr=[izr,max(d.izu)-[1:2]];
 end
 d=setdefv(d,'izr',izr);
 di.izr=d.izr;
 izr=di.izr;

% remove reference velocity and then average ensembles
  disp(' remove reference velocity and average ensembles ')
  ilast=1;
  il=length(d.izm);
  im=0;
  ibin=1:p.nbins;

% big loop
  while ilast<il 
   im=im+1;
   i=ilast;
   if p.avens>0

%  fixed number of ensembles
    i1=ilast+[1:p.avens];
   else

%  set up index for ensembles to average based on depth
    ii=find(abs(d.izm(1,(ilast+1):il)-d.izm(1,ilast))>p.avdz);
    if length(ii)<1, ii=il-ilast; end
    i1=ilast+[1:ii(1)];
   end

   i1l=length(i1)/2*p.oversample;
   i1=round(mean(i1)+[-i1l:i1l]);
   ii=find(i1<1); i1(ii)=[];
   ii=find(i1>il); i1(ii)=[];
   if length(i1)==1, i1=[i1 i1]; end
   ilast=max(i1);
  

   w=d.weight(izr,i1)*0+1;
   w2=d.weight(:,i1)*0+1;

% U
   ur=medianan(d.ru(izr,i1).*w);
   ruav=meannan(ur);
   i3=find(isnan(ur));
   ur(i3)=i3*0;
   iav=round(length(ur)/200*p.avpercent);
   ur=meshgrid(ur,ibin);
   di.ru(:,im)=medianan([d.ru(:,i1).*w2-ur]',iav)'+ruav;
   rus=stdnan([d.ru(:,i1).*w2-ur]')';
% V
   vr=medianan(d.rv(izr,i1).*w);
   rvav=meannan(vr);
   i3=find(isnan(vr));
   vr(i3)=i3*0;
   iav=round(length(vr)/200*p.avpercent);
   vr=meshgrid(vr,ibin);
   di.rv(:,im)=medianan([d.rv(:,i1).*w2-vr]',iav)'+rvav;
% estimate mean STD of U and V
   di.ruvs(:,im)=sqrt(rus.^2+stdnan([d.rv(:,i1).*w2-vr]')'.^2);
% W
   wr=medianan(d.rw(izr,i1).*w);
   rwav=meannan(wr);
   i3=find(isnan(wr));
   wr(i3)=i3*0;
   iav=round(length(wr)/200*p.avpercent);
   wr=meshgrid(wr,ibin);
   di.rw(:,im)=medianan([d.rw(:,i1).*w2-wr]',iav)'+rwav;

%EA
   di.re(:,im)=meannan(d.re(:,i1)')';

%TS
   di.ts(:,im)=meannan(d.ts(:,i1)')';
   di.tg(:,im)=meannan(d.tg(:,i1)')';

% weight
   di.weight(:,im)=meannan(d.weight(:,i1)')';

% bottom track
   di.bvel(:,im)=meannan(d.bvel(i1,:))';
   bvel=d.bvel(i1,:);
%  remove mean vertical velocity from bottom track w prior to STD 
   bvel(:,3)=bvel(:,3)-wr(1,:)';
   di.bvels(:,im)=stdnan(bvel)';
%  distance of bottom
   di.hbot(im)=meannan(d.hbot(i1));

% bin depth
   di.izm(:,im)=mean(d.izm(:,i1)')';

% heading
   di.hdg(1,im)=-angle(mean(u1d(i1)))*180/pi;
   if length(d.izu)>1
     di.hdg(2,im)=-angle(mean(u1u(i1)))*180/pi;
   end

% pitch and roll
   di.pit(:,im)=mean(d.pit(:,i1),2);
   di.rol(:,im)=mean(d.rol(:,i1),2);

% target strength
   di.tsd(im)=mean(d.ts(d.izd(2),i1),2);
   di.tsd_out(im)=mean(d.ts(d.izd(end),i1),2);

% target strength
   di.temp(im)=mean(d.temp(i1));

% ships position
   di.slon(im)=medianan(d.slon(i1));
   di.slat(im)=medianan(d.slat(i1));

% number of ensembles
   di.dtiv(im)=length(i1);

% time
   di.time_jul(im)=meannan(d.time_jul(i1));

% depth
   di.z(im)=meannan(d.z(i1));
% end of big loop
end

% adjust compass to RDI definition [0-360]
di.hdg=di.hdg+(di.hdg<0)*360;

% remove outlier
di.izd=d.izd;
di.izu=d.izu;
[di,p]=outlier(di,p);

% remove bottom track data with single ping or large wstd
if isfinite(p.zbottom)
 % remove single ping bottom track ensembles
 ii=find(prod(di.bvels(1:3,:))==0);
 disp([' found ',int2str(length(ii)),' bottom track std==0 set to 0.1 m/s'])
 di.bvels(:,ii)=0.1;

 ii=find(di.bvels(3,:)>0);
 if length(ii)>0
  disp([' found ',int2str(length(ii)),' finite bottom track ensembles'])
  p=setdefv(p,'btrk_wstd',median(di.bvels(3,ii))*2);
  ii=find(di.bvels(3,:)>p.btrk_wstd | di.bvels(3,:)==0);
  disp([' discarded ',int2str(length(ii)),...
    ' bottom tracks velocities because of wstd  > ',num2str(p.btrk_wstd)])
  di.bvel(:,ii)=NaN;
  di.bvels(:,ii)=NaN;
 else
  ii=find(isfinite(d.bvel(3,:)));
  disp([' found no valid bottom track ensemble from ',int2str(length(ii)),...
        ' finite raw bottom tracks '])
 end 
end


% remove ensembles without velocity data
ii=find(~isfinite(maxnan(di.ru)));

if length(ii)>0
 disp([' removed ',int2str(length(ii)),...
        ' non finite super ensembles'])
 di.rw(:,ii)=[];
 di.ru(:,ii)=[];
 di.rv(:,ii)=[];
 di.ruvs(:,ii)=[];
 di.bvel(:,ii)=[];
 di.bvels(:,ii)=[];
 di.weight(:,ii)=[];
 di.izm(:,ii)=[];
 di.z(ii)=[];
 di.hdg(:,ii)=[];
 di.pit(:,ii)=[];
 di.rol(:,ii)=[];
 di.time_jul(ii)=[];
 di.dtiv(ii)=[];
 di.hbot(ii)=[];
 di.tsd(ii)=[];
 di.temp(ii)=[];
 di.slon(ii)=[];
 di.slat(ii)=[];
end

% check for small std's
% blank out data using weight 
di.ruvs=di.ruvs+di.weight*0;
ii=find(di.ruvs==0);
di.weight(ii)=nan;
disp([' set ',int2str(length(ii)),' weight values to nan  because super ensemble std =0 '])
di.ruvs=di.ruvs+di.weight*0;
ii=find(di.ruvs<p.superens_std_min);
di.ruvs(ii)=p.superens_std_min;
disp([' set ',int2str(length(ii)),' values to minimum super ensemble std ',...
     num2str(p.superens_std_min)]);

disp([' reduced profile length = ',int2str(length(di.z)),' super-ensemble bins'])
if length(di.z)<5
	error('not enough data to process station ')
end
end

% time difference in seconds
dtt=diff(di.time_jul)*24*3600;
di.dt=mean([dtt([1,1:end]);dtt([1:end,end])]);

% set positons
p=setdefv(p,'poss',[0 0 0 0]);
p=setdefv(p,'pose',[0 0 0 0]);

slat = p.poss(1)+p.poss(2)/60.0;
elat = p.pose(1)+p.pose(2)/60.0;
slon = p.poss(3)+p.poss(4)/60.0;
elon = p.pose(3)+p.pose(4)/60.0;

% compute ships drift
sjul =julian(p.time_start);
ejul =julian(p.time_end);

p.dt_profile = (ejul-sjul)*24*3600;

dlat = elat - slat;
dlon = elon - slon;

p.lat= ( slat + elat ) /2.0;
p.lon= ( slon + elon ) /2.0;

p.xdisp = dlon * cos(p.lat*pi/180) * 60.0 * 1852.0;
p.ydisp = dlat * 60.0 * 1852.0;

p.uship= p.xdisp / p.dt_profile;
p.vship= p.ydisp / p.dt_profile;


%--------------------------------------------------
function  hoff=compoff(u1,u2)
% compute mean compass offset 

h1=-angle(u1)*180/pi;
h1=h1+(h1<0)*360;
nhead=36;
dhead=360/2/nhead;
h0=linspace(5,355,nhead);

for i=1:nhead
 ii=find(abs(h1-h0(i))<=dhead);
 if length(ii)>1
  u1a(i)=mean(u1(ii));
  u2a(i)=mean(u2(ii));
 else
  u1a(i)=NaN;
  u2a(i)=NaN;
 end
end

ii=find(isfinite(u1a+u2a));

if length(ii)>0
 hoff=angle(u1a(ii)/u2a(ii))*180/pi;
else
 hoff=0;
end

return

%-----------------------------------------------------
function y=boxav(x,n)
% boxaverage and subsample
[ld,lv]=size(x);
in=fix(ld/n);

for i=1:lv
 xm=reshape(x(1:(in*n),i),n,in);
 y(:,i)=meannan(xm)';
end