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

%                    G E T D P T H I . M 

%                    doc: Wed Jan  7 16:25:26 2009

%                    dlm: Fri Mar  5 15:51:15 2010

%                    (c) 2009 A.M. Thurnherr

%                    uE-Info: 459 33 NIL 0 0 72 0 2 4 NIL ofnI

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

% CHANGES BY ANT:

%   Jan  7, 2009: - tightened use of exist()

function [d,p]=getdpthi(d,p)

% function [d,p]=getdpthi(d,p)

% LADCP-2 processing software v 7.0

%

%  -- make depth from raw data 

%  use inverse approach

% Matin Visbeck 

%  December 2002, LDEO

disp('GETDPTHI: Depth from vertical velocity inverse method')

% set default start, deepst  and end depth

p=setdefv(p,'zpar',[10 NaN 10]);

p=setdefv(p,'cut',15);

p=setdefv(p,'dzbelow',[2 -1]*medianan(abs(diff(d.zd))));

p=setdefv(p,'guessbottom',NaN);

p=setdefv(p,'zbottom',NaN);

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

p=setdefv(p,'navtime_av',2/60/24);

p.zpar=abs(p.zpar);

% remove empty profiles at beginning and end

ii=ones(size(d.rw(1,:)));

jj=find(~isnan(meannan(d.rw)));

ii(1:jj(1))=0;

ii(jj(end):end)=0;

d=cutstruct(d,ii);

% get time difference for w-values

dt=diff(d.time_jul)*24*3600;

dt=mean([dt([1,1:end]);dt([1:end,end])]);

% initialize a matrix where below bottom and above surface data

% get flagged as bad

d.izmflag = d.rw*0;

figure(4)

orient tall

% two sweeps throught the data

n2=2;

% only one sweep if ctd pressure data exist

if p.ctddepth>0, n2=1; end

for n=1:n2

  disp([' starting run ',int2str(n),' to get LADCP depth'])

  if n==1 

%  if n==1 | (~exist('ddoz') & ~exist('ddoz'))

%  first sweep used intergal of median w 

   ibad=find(p.wizr<1 | p.wizr>size(d.rw,1));

   p.wizr(ibad)=[];

   if length(p.wizr)>1

    dw=medianan(d.rw(p.wizr,:),1);

   elseif length(p.wizr)==1

    dw=d.rw(p.wizr,:);

   else

    dw=medianan(d.rw,1);

   end

%  set non finite w to zero

   ii=find(~isfinite(dw));

   dw(ii)=0;

%  integrate result

   zz=cumsum(dw.*dt);

%  make sure that start and end depth are as wanted

   zz=zz-linspace(-maxnan([0 p.zpar(1)]),-maxnan([0 p.zpar(3)])+zz(end),length(zz));

%  set maxdepth to be as requested

   if isfinite(p.zpar(2))

    zz=zz/max(zz)*p.zpar(2);

   end

  else

%  second sweep used inverse method to get z including bottom and

%  surface track data

%  "erase" subplots

  subplot(321), hold off

  subplot(322), hold off

  subplot(312), hold off

% set base matrix

   dw=d.rw(p.wizr,:)+d.izmflag(p.wizr,:);

   [d1,A1,ibot]=dinset(dw,dt);

   % set boundary conditions for inversion

   if isfinite(p.zpar(1))

     d1=[d1;p.zpar(1)];

   else 

     d1=[d1;10];

   end

   A1(length(d1),1)=1;

   if isfinite(p.zpar(2))

     d1=[d1;p.zpar(2)*10];

     A1(length(d1),ibot)=10;

   end

   if isfinite(p.zpar(3))

     d1=[d1;p.zpar(3)];

   else

     d1=[d1;10];

   end

   A1(length(d1),end)=1;

   % add surface/bottom reflections if present

   if exist('ddoz','var')

    [ld,lz]=size(A1);

    d1=[d1;ddoz*0.1];

    A2=sparse(1:length(ddoi),ddoi,0.1);

    A2(1,lz)=0;

    A1=[A1;A2];

   end

   if exist('dupz','var')

    [ld,lz]=size(A1);

    d1=[d1;dupz*0.1];

    A2=sparse(1:length(dupi),dupi,0.1);

    A2(1,lz)=0;

    A1=[A1;A2];

   end

   if exist('dbotdz','var') 

    bfac=0.1;

    [ld,lz]=size(A1);

    d1=[d1;dbotdz*bfac];

    ix=[1:length(dboti)]';

    A2=sparse([ix;ix],[dboti;ix*0+ibot],[ix*0-bfac;ix*0+bfac]);

    A2(1,lz)=0;

    A1=[A1;A2];

   end

% require z to be smooth (needed otherwise ill constrained)

   [A1,d1]=dismoo(A1,d1,0.01);

% solve for best depth time series

   zz=lesqchol(d1,A1)';

  end

  % set depth below which to delete all data 

  dzbelow=p.dzbelow(n);

  if n==1

  % get sound speed time series

   if ~existf(d,'ss')

    if existf(d,'ctd_ss')

      disp(' take soundspeed from CTD time series')

      d.ss=d.ctd_ss;

    elseif existf(d,'ctdprof_ss')

     disp(' take soundspeed from CTD profile ')

     zctd=d.ctdprof_z;

     zctd(1)=-1e5;

     zctd(end)=1e5;

     d.ss=interp1q(zctd,d.ctdprof_ss,zz')';

    else

      if existf(d,'ctd_temp')

        disp(' make soundspeed based on CTD pressure and temp')

        pp=press(abs(d.z));

        d.ss=sounds(pp,d.ctd_temp,34.5);

      else

        disp(' make soundspeed based on pressure and ADCP temp')

        pp=press(abs(zz));

        d.ss=sounds(pp,d.temp(1,:),34.5);

      end

    end

   end

   % sound speed correction

   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 existf(d,'hbot')

     d.hbot=d.hbot.*sc(1,:);

     d.bvel(:,1:3)=d.bvel(:,1:3).*sc(1:3,:)';

     if existf(d,'bvel_rdi')

      d.bvel_rdi(:,1:3)=d.bvel_rdi(:,1:3).*sc(1:3,:)';

     end

     if existf(d,'bvel_own')

      d.bvel_own(:,1:3)=d.bvel_own(:,1:3).*sc(1:3,:)';

     end

    end

    if length(d.zd)~=length(d.ru(:,1))

     sc=meshgrid(d.ss./d.sv(2,:),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

    if existf(d,'hsurf')

     d.hsurf=d.hsurf.*sc(1,:);

    end

    d.soundc=1;

   else

    disp(' will not correct for sound speed twice')

   end

  end

  if p.ctddepth==0

   % save results only if CTD-depth was not available

   d.z=-zz;

   p.ladcpdepth=2;

   disp(' use LADCP depth from constrained integrated W ')

  else

   d.z_ladcp=-zz;

   dz=d.z_ladcp-d.z;

   ii=find(isfinite(dz));

   p.ladcpr_CTD_depth_std=[mean(dz), std(dz)];

   disp(' use CTD time series depth, will not do depth inversion ')

   disp([' LADCP minus CTD depth mean: ',num2str(p.ladcpr_CTD_depth_std(1)),...

          '  std: ',num2str(p.ladcpr_CTD_depth_std(2))]);

   p.ladcpdepth=0;

  end

  [p.maxdepth,ibottom]=max(-d.z);

  % plot near surface LADCP data 

  [dum,ibot]=max(-d.z);

  % for very shallow stations turn of surface detection

  if dum<100, p.surfdist=0; disp(' shallow station no surface detection '),end

  % first down cats

  ii=1:ibot;

  iok=ii(find(d.z(ii)>-200 & d.z(ii)<-30 ));

  if length(iok)>2

   iok2=1:iok(end);

   subplot(321)

   plot(iok2,d.z(iok2),'.r',iok2,iok2*0,'-k'),

   hold on

   title('start depth (.red) (.blue) Surface distance')

   xlabel('time in ensembles')

   ylabel('depth in meter')

   axis tight

   ax=axis;

   ax(4)=max(ax(4),20);

   axis(ax);

  else

   disp('no surface reflections down-cast')

  end

  % then up cast

  ii=ibot:length(d.z);

  iok=ii(find(d.z(ii)>-200 & d.z(ii)<-30 ));

  if length(iok)>2

   iok2=iok(1):length(d.z);

   subplot(322)

   plot(iok2,d.z(iok2),'.r',iok2,iok2*0,'-k'),

   hold on

   title('end depth')

   xlabel('time in ensembles')

   ylabel('depth in meter')

   axis tight

   ax=axis;

   ax(4)=max(ax(4),20);

   axis(ax);

  else

   disp('no surface reflections up-cast')

  end

  % extract surface distance to find start depth

  if existf(d,'hsurf') & p.surfdist

   if sum(isfinite(d.hsurf))>10

    [zmax,ibot]=max(-d.z);

   % start depth

    ii=find(isfinite(d.hsurf(1:ibot)));

   % limit to surface detections within 30m of answer

    iok=ii(find(d.z(ii)>-200 & d.z(ii)<-40 & abs(d.hsurf(ii)+d.z(ii))<30));

    if ~isempty(iok)

     ddoz=d.hsurf(iok)';

     ddoi=iok';

    % temporary plot of surface detection

      iok2=1:iok(end);

      subplot(321)

      plot(iok,d.hsurf(iok)+d.z(iok),'.',iok2,d.z(iok2),'.r'),

      ax=axis;

      ax(4)=max(ax(4),20);

      axis(ax);

      axis tight

    else

     disp(' no surface reflections down-cast')

    end

   % surface distance to find end depth

    ii=find(isfinite(d.hsurf(ibot:end)))+ibot-1;

   % limit to surface detections within 30m of answer

    iok=ii(find(d.z(ii)>-200 & d.z(ii)<-40 & abs(d.hsurf(ii)+d.z(ii))<30));

    if ~isempty(iok)

     dupz=d.hsurf(iok)';

     dupi=iok';

    % temporary plot of surface detection

      iok2=iok(1):length(d.z);

      subplot(322)

      plot(iok,d.hsurf(iok)+d.z(iok),'.',iok2,d.z(iok2),'.r'),

      ax=axis;

      ax(4)=max(ax(4),20);

      axis(ax);

      axis tight

    else

     disp(' no surface reflections up-cast')

    end

   end 

  end 

  disp([' maximum depth from int W is :',int2str(max(zz))])

  disp([' should be                   :',int2str(p.zpar(2))])

  % bottom depth

  subplot(312)

  iok=find((max(-d.z)+d.z)<200);

  iok2=iok(1):iok(end);

  plot(iok2,d.z(iok2),'.r'),

  if ~isfinite(p.zbottom)

   if sum(isfinite(d.hbot))>10

    % look for bottom only close to deepest CTD depth

    iok=find((max(-d.z)+d.z)<200 & d.hbot>0 );

    if length(iok)>2

      % fit polynom to bottom depth time series

      if n==n2

      % use deepest point to set bottom depth

       if isfinite(p.guessbottom)

	  zbottom=p.guessbottom;

       else

          zbottom= median(d.hbot(iok)-d.z(iok));

       end

       zbottomerr= zbottom-(d.hbot(iok)-d.z(iok)) ;

       [dum,is]=sort(abs(zbottomerr));

       is=is(1:fix(length(is)/2));

       c=polyfit(iok(is),d.hbot(iok(is))-d.z(iok(is)),1);

       zbottomerr= polyval(c,iok)-(d.hbot(iok)-d.z(iok)) ;

       [dum,is]=sort(abs(zbottomerr));

       is=is(1:fix(length(is)/2));

       iok=iok(find(abs(zbottomerr)<2*std(zbottomerr(is)) | abs(zbottomerr)<30 ));

       c=polyfit(iok,d.hbot(iok)-d.z(iok),2);

       zbottomerr= polyval(c,iok)-(d.hbot(iok)-d.z(iok)) ;

       zbottom=polyval(c,ibottom);

       p.zbottom=zbottom;

      else

       if isfinite(p.guessbottom)

	  zbottom=p.guessbottom;

       else

          zbottom= median(-d.hbot(iok)-d.z(iok));

       end

       zbottomerr= zbottom-(d.hbot(iok)-d.z(iok)) ;

       [dum,is]=sort(abs(zbottomerr));

       is=is(1:fix(length(is)/2));

       zbottom=medianan(d.hbot(iok(is))-d.z(iok(is)));

       zbottomerr= zbottom-(d.hbot(iok)-d.z(iok)) ;

       [dum,is]=sort(abs(zbottomerr));

       is=is(1:fix(length(is)/2));

       iok=iok(find(abs(zbottomerr)<2*std(zbottomerr(is)) | abs(zbottomerr)<50 ));

       c=polyfit(iok,d.hbot(iok)-d.z(iok),1);

       % save bottom distances for inversion

       dbotdz=(d.hbot(iok)-polyval(c,iok)-d.z(ibottom))';

       dboti=iok';

      end

      p.zbottomerror = medianan(abs(zbottomerr)); 

    % temporary plot of bottom detection

      iok2=iok(1):iok(end);

      plot(iok,-d.hbot(iok)+d.z(iok),'.',iok2,d.z(iok2),'.r'),

      hold on

      plot(iok2,-zz(iok2),'-k')

      hold on, plot(iok2,iok2*0-p.zbottom,'--k') 

      hold on, plot(iok2,-polyval(c,iok2),'-b')

      title('bottom (--k) LADCP depth (-k) bottom distance (.b)')

      xlabel('time in ensembles')

      ylabel('depth in meter')

      if n==n2

 % remove suspicious bottom track data

       axis tight

       ax=axis;

       text(ax(1)+abs(diff(ax(1:2)))*.15,ax(3)+abs(diff(ax(3:4)))*.8,...

       ['bottom at: ',int2str(p.zbottom),' [m]   ADCP was ',...

         int2str(p.zbottom-max(-d.z)),' m above bottom'])

       ibad=1:length(d.hbot);

       % good data are

       ibad(iok)=[];

       d.hbot(ibad)=NaN;

       d.bvel(ibad,:)=NaN;

      else

      end

    else

      zbottom=NaN;

      p.zbottomerror = nan;

    end

   else

    zbottom=NaN;

    p.zbottomerror=NaN;

   end

  else

    disp(['  bottom preset at ',int2str(p.zbottom)])

    zbottom=p.zbottom;

    p.zbottomerror=0;

  end

   % check if bottom is shallower that maxctd-depth an

    if ((zbottom-p.maxdepth<-(p.maxdepth*0.01+10) & isfinite(zbottom)) |...

       	p.zbottomerror > 20 )

      disp('  no bottom found')

      disp(['   given maximum profile depth : ',int2str(p.maxdepth)])

      disp(['   extracted bottom depth      : ',int2str(zbottom)])

      disp(['        bottom depth error     : ',int2str(p.zbottomerror)])

      p.zbottom=NaN;

    elseif n==n2

      p.zbottom=zbottom;

      disp(['  bottom found at ',int2str(p.zbottom),' +/- ',...

                                 int2str(p.zbottomerror),' m'])

      if (p.zbottom<p.maxdepth)

        disp('  extracted bottom within 20m above given maximum profile depth')

      end

    end

  pause(0.1)

  % assign a depth to each bin

  [izm1,izm]=meshgrid([fliplr(d.zu),-d.zd],d.z);

  izm1=izm1';

  if d.soundc==1

  % make sound speed correction for depth vector

   sc=meshgrid(d.ss./d.sv(1,:),d.izd);

   izm1(d.izd,:)=izm1(d.izd,:).*sc;

   if length(d.zu)>0

    sc=meshgrid(d.ss./d.sv(2,:),d.izu);

    izm1(d.izu,:)=izm1(d.izu,:).*sc;

   end

   disp(' correct bin length for sound speed')

  end

  % add the two parts

  d.izm=izm'+izm1;

  % flag all data below bottom as bad

  if ~isnan(p.zbottom)

    ii = find(d.izm<-p.zbottom-dzbelow);

    d.izmflag(ii)=NaN; 

  end

  % flag all data close to the surface as bad

  if length(d.zu)>0

   ii = find(d.izm>-(d.zu(2)-d.zu(1))/2);

   d.izmflag(ii)=NaN; 

  end

end 

% set velocities deeper than bottom to NaN

bad = find( isnan(d.izmflag) & isfinite(d.ru) );

if ~isempty(bad) & isfinite(p.zbottom)

  disp([' removing ',int2str(length(bad)),...

	' values below bottom'])

 %d.ru = d.ru+d.izmflag;

 %d.rv = d.rv+d.izmflag;

 %d.rw = d.rw+d.izmflag;

 d.weight = d.weight + d.izmflag;

end

% compute pressure from depth

d.p=press(abs(d.z));

if existf(d,'wctd')==1

 d.wm=d.wctd;

else

 d.wm=-gradient(d.z)./gradient(d.time_jul*24*3600);

end

if d.z(1)<-50, 

 warn=[' first LADCP depth is ',int2str(d.z(1))];

 disp(warn)

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

end

if d.z(end)<-50, 

 warn=[' last LADCP depth is ',int2str(d.z(end))];

 disp(warn)

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

end

% cut raw data to only include profile

i1=find(d.z(1:ibottom)>-p.cut);

i2=find(d.z(ibottom:end)>-p.cut)+ibottom-1;

if length(i1)==0

 i1=1;

end

if length(i2)==0

 i2=length(d.z);

end

ii=d.z*0;

ic=i1(end):i2(1);

ii(ic)=1;

if (sum(ii)~=length(ii)) & p.cut>0 

 disp(' remove data at begining and end of cast')

 disp(' adjust start and end time ')

 d=cutstruct(d,ii);

 p.zpar([1 3])=p.cut;

 subplot(321), ax=axis; plot([1 1]*ic(1),ax(3:4),'--k')

 subplot(322), ax=axis; plot([1 1]*ic(end),ax(3:4),'--k')

 p.time_start=gregoria(d.time_jul(1));

 p.time_end=gregoria(d.time_jul(end));

 if existf(d,'slon')

  % average over first p.navtime_av days

  ii=find(d.time_jul<(d.time_jul(1)+p.navtime_av));

  slon=median(d.slon(ii));

  slat=median(d.slat(ii));

  if isfinite(slon+slat)

   p.poss=[fix(slat), (slat-fix(slat))*60, fix(slon), (slon-fix(slon))*60];

  end

  ii=find(d.time_jul>(d.time_jul(end)-p.navtime_av));

  elon=median(d.slon(ii));

  elat=median(d.slat(ii));

  disp(' find new start end position')

 elseif existf(p,'poss')

  % don't have time variable navigation data, adjust end position

  % to account for shorter cast time

  slat=p.poss(1)+p.poss(2)/60; 

  slon=p.poss(3)+p.poss(4)/60;